An Investigation of Wind Direction and Speed in a Featured Wind Farm Using Joint Probability Distribution Methods

Abstract

Wind direction and speed are both crucial factors for wind farm layout; however, the relationship between the two factors has not been well addressed. To optimize wind farm layout, this study aims to statistically explore wind speed characteristics under different wind directions and wind direction characteristics. For this purpose, the angular–linear model for approximating wind direction and speed characteristics were adopted and constructed with specified marginal distributions. Specifically, Weibull–Weibull distribution, lognormal–lognormal distribution and Weibull–lognormal distribution were applied to represent the marginal distribution of wind speed. Moreover, the finite mixture of von Mises function (FVMF) model was used to investigate the marginal distribution of wind direction. The parameters of those models were estimated by the expectation–maximum method. The optimal model was obtained by comparing the coefficient of determination value (R²) and Akaike’s information criteria (AIC). In the numerical study, wind data measured at a featured wind farm in north China was adopted. Results showed that the proposed joint distribution function could accurately represent the actual wind data at different heights, with the coefficient of determination value (R²) of 0.99.

1. Introduction

With the gradually increasing consumption of nonrenewable energy, such as oil, coal, and natural gas, renewable energy has become the hope of humankind to solve the current energy shortage and future energy crisis [1,2]. Wind energy is a well-known renewable energy source and has been used worldwide due to its advantages of clean, renewable, and large reserves [3,4,5]. People capture wind energy through wind turbines, and wind speed and wind direction characteristics are therefore crucial for calculating wind power and the design and arrangement of wind turbines [6]. The wind speed probability density distribution can be used to estimate wind power density [7,8] and can also be applied to determine the most dominant wind direction [9,10,11]. In order to optimize wind farm layout in complex terrain, the joint probability distribution of wind direction and wind speed have gained great attention [12,13,14] as it can be used to obtain wind speed characteristics under different wind directions in complex terrain. Wind speed characteristics under different wind directions play a significant role in arrangement of wind turbines. It should be noted that wind speed with the highest frequency and wind speed that can capture the maximum wind energy are different. To obtain maximum energy output, the two wind speeds should be as close as possible [9]. In addition, the joint probability distribution is important in the construction of the measure–correlate–predict method and the calculation of structural wind load [15,16,17].

Up to now, a large number of probability models have been used to fit wind velocity distribution [18,19,20,21]. However, there is no fixed optimal model for different regions. Masseran et al. [22] investigated eight probability models to describe the wind speed probability distribution: Weibull, Rayleigh, lognormal, exponential, inverse Gaussian, gamma, Burr, and inverse gamma. Jaramillo et al. [23] found that the Weibull–Weibull distribution was more appropriate than the two-parameter Weibull distribution for regions where wind speed presented a bimodal probability distribution. Kollu et al. [24] used a mixture of three probability distributions—Weibull–lognormal distribution, Weibull–extreme value, and lognormal–extreme value distribution—to model wind speed distribution. The result showed that the mixed model had good fitting on both multipeak and single-peak wind speed distribution. Moreover, it is known that wind direction distribution plays a crucial role in determining the most dominant wind direction [11]. However, compared to wind speed, there are few models that are suitable for fitting wind direction distribution [25]. This can be attributed to two main difficulties. Firstly, the wind direction is generally divided into 16 intercardinal directions: north (N), north-northeast (NNE), northeast (NE), east-northeast (ENE), east (E), east-southeast (ESE), southeast (SE), south-southeast (SSE), south (S), south-southwest (SSW), southwest (SW), west-southwest (WSW), west (W), west-northwest (WNW), northwest (NW), and north-northwest (NNW). Second, the wind direction is a circular variable. Masseran et al. [11] manifested that a finite mixture of von Mises function (FVMF) could be applied to describe the actual data of wind direction in Malaysia. As a result, it was found that more than 90% of wind directions could be explained by FVMF models. So far, FVMF models have been widely applied to model the angular variable, which plays an important role in many fields, such as wind energy [26,27,28,29], environment [30], economics [31], and video analysis [32].

In modeling the joint probability function of wind direction and speed, Johnson and Wehrly [12] constructed four kinds of models, and the fourth model was a combination of wind direction and speed distribution. Carta et al. [6] derived joint distribution with isotropic Gaussian model, angular–linear model, and anisotropic Gaussian model. The conclusion was that the angular–linear model had a better degree of fit in the case of the Canary Islands. Erdem and Shi [13] constructed bivariate joint distribution based on angular–linear model, Farlie–Gumbel–Morgenstern model, and anisotropic lognormal approaches model. The results indicated that angular–linear model had a good fit to the measured wind data in North Dakota, USA. Han et al. [17] proposed two nonparametric models to fit joint speed and direction distributions and claimed that the novel methods performed better than parametric distributions. Up to now, very few publications have utilized joint probability distributions to simultaneously approximate wind speed and direction.

For a more accurate description of complex wind characteristics, the wind speed distribution is described by a mixture of three probability distributions: Weibull–Weibull (W–W), lognormal–lognormal (L–L), and Weibull–lognormal (W–L); the FVMF models are used to investigate the wind direction characteristics, and the joint probability distribution of wind direction and speed is obtained with angular–linear model. The suitability of the distribution is judged from the coefficient of determination R².

2. Wind Data Measurement

The wind data were collected for a wind farm in Taonan, China (45°24′18.16″ N, 122°20′22.24″ E) from April 2013 to March 2017. The actual wind direction and speed data were taken at a height of 70 m and 50 m above ground level and stored with a time interval of 10 min.

3. Models

3.1. Joint Distribution

In this paper, we used the angular–linear model proposed by Johnson and Wehrly [12] to construct the joint probability distribution of wind direction and speed. The joint probability distribution of wind direction and speed is as follows (Equations (1) and (2)):

$$f_{V,\Theta }=2\pi g\left(\zeta \right)f_V\left(v\right)f_{\Theta }\left(\theta \right); 0\le \theta <2\pi ; -\infty \le v<\infty$$

(1)

where $f_V\left(v\right)$ is the probability distribution of the wind speed; $f_{\Theta }\left(\theta \right)$ is the probability distribution of the wind direction; $g(\cdot )$ is the probability distribution of the angular variable ζ. The angular variable ζ is as follows:

$$\zeta =\{\begin{array}{c}2\pi \left[F_V\left(v\right)-F_{\Theta }\left(\theta \right)\right],F_V\left(v\right)\ge F_{\Theta }\left(\theta \right)\\ 2\pi \left[F_V\left(v\right)-F_{\Theta }\left(\theta \right)\right]+2\pi ,F_V\left(v\right)<F_{\Theta }\left(\theta \right)\end{array}$$

(2)

where $F_V\left(v\right)$ is the cumulative probability distribution of the wind speed; $F_{\Theta }\left(\theta \right)$ is the cumulative probability distribution of the wind direction.

3.2. Wind Speed Distribution

For the probability distribution of wind speed, we used W–W, given by Equation (3), L–L, given by Equation (4), and W–L, given by Equation (5):

$$f_V\left(v\right)={{\displaystyle \sum }}_{m=1}^{M=2}{w}_m\frac{{\beta }_m}{{\alpha }_m}{(\frac{v}{{\alpha }_m})}^{{\beta }_m-1}exp\left[-{\left(\frac{v}{{\alpha }_m}\right)}^{{\beta }_m}\right],v>0;{\alpha }_m,{\beta }_m>0$$

(3)

where $w_m$ are weight coefficients that add to one; the parameter ${\alpha }_m$ is a shape parameter; ${\beta }_m$ is a scale parameter with the same units as the wind speed.

$$f_V\left(v\right)={{\displaystyle \sum }}_{m=1}^{M=2}\frac{w_m}{v{\sigma }_m\sqrt{2\pi }}exp[\frac{-{\left(ln\left(v\right)-{\mu }_m\right)}^2}{2{\sigma }_m^2}], v>0; {\sigma }_m>0$$

(4)

where $w_m$ are weight coefficients that add to one; the parameter ${\mu }_m$ is the mean of the associated normal distribution; ${\sigma }_m$ is standard deviation of the associated normal distribution.

$$f_V\left(v\right)=w_1\frac{\beta }{\alpha }{(\frac{v}{\alpha })}^{\beta -1}exp\left[-{\left(\frac{v}{\alpha }\right)}^{\beta }\right]+w_2\frac{1}{v\mathsf{\sigma }\sqrt{2\pi }}exp[\frac{-{\left(ln\left(v\right)-\mu \right)}^2}{2{\sigma }^2}], v>0$$

(5)

where $w_1$ and $w_2$ are weight coefficients that add to one.

3.3. Circular Variable Distribution

For the probability distribution of wind direction, we used a finite mixture of von Mises function. The FVMF is a mixture distribution that consists of Mth single von Mises distribution. The finite mixture von Mises function is as follows (Equations (6) and (7)):

$$f_{\Theta }\left(\theta \right)={{\displaystyle \sum }}_{m=1}^M\frac{w_m}{2\pi I_0\left(k_m\right)}exp\left[k_m\cos \left(\theta -{\mu }_m\right)\right], 0\le \theta <2\pi ; k_m\ge 0;0\le {\mu }_m<2\pi$$

(6)

where $w_m$ are weight coefficients that add to one; the parameter ${\mu }_m$ is the mean wind direction; and the parameter $k_m$ is the concentration parameter. Here, $I_0\left(k_m\right)$ is the modified Bessel function of the first kind and order zero [33]. $I_0\left(k_m\right)$ is as follows:

$$I_0\left(k_m\right)=\frac{1}{\sqrt{2\pi }}{{\displaystyle \int }}_0^{2\pi }\exp \left[k_m\cos x\right]dx={{\displaystyle \sum }}_{k=0}^{\infty }\frac{1}{{\left(k!\right)}^2}{(\frac{k_m}{2})}^{2k}$$

(7)

For the probability distribution of the angular variable ζ, the two component mixtures von Mises function was used [6].

3.4. Cumulative Distribution Function

The cumulative probability distribution of the wind speed can be derived through Equation (8). The cumulative probability distribution of the wind direction can be derived through Equation (9).

$$F_V\left(v\right)={{\displaystyle \int }}_0^v{f}_V\left(v\right)dv$$

(8)

$$F_{\Theta }\left(\theta \right)={{\displaystyle \int }}_0^{\theta }{f}_{\Theta }\left(\theta \right)d\theta$$

(9)

4. Parameter Estimation

In this paper, we chose the expectation–maximization method to estimate the parameters of W–W, L–L, W–L, and FVMF [34].

4.1. Parameter Estimation of W–W Distribution

If we let $b_m={\beta }_m$ and $c_m={{\alpha }_m}^{{\beta }_m}$, the W–W in Equation (3) can also be written as follows (Equation (10)):

$$f_V\left(v\right)={{\displaystyle \sum }}_{m=1}^{M=2}{w}_m\frac{b_m}{c_m}{v}^{b_m-1}exp\left[-\frac{v^{b_m}}{c_m}\right]$$

(10)

The W–W parameters ($w_m$, $b_m$, and $c_m$) can be estimated through the actual wind speed data. The W–W parameters can be derived by Equations (11)–(14):

$$p_n=p\left(m|v_n;w_m,b_m,c_m\right)=\frac{w_m\frac{b_m}{c_m}{v_n}^{b_m-1}exp[-\frac{{v_n}^{b_m}}{c_m}]}{{{\displaystyle \sum }}_{m=1}^{M=2}{w}_m\frac{b_m}{c_m}{v_n}^{b_m-1}exp[-\frac{{v_n}^{b_m}}{c_m}]}$$

(11)

$${\widehat{w}}_m=\frac{1}{N}{{\displaystyle \sum }}_{n=1}^N{p}_n$$

(12)

$$\frac{{{\displaystyle \sum }}_{n=1}^N{p}_n}{{\widehat{b}}_m}+{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)p_n\right]-\frac{{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)v_n^{{\widehat{b}}_m}{p}_n\right]}{{\widehat{c}}_m}=0$$

(13)

$${{\displaystyle \sum }}_{n=1}^N\left(v_n^{{\widehat{b}}_m}{p}_n\right)-{\widehat{c}}_m{{\displaystyle \sum }}_{n=1}^N{p}_n=0$$

(14)

where N is the number of measured wind speed data; $p_n=p\left(m|v_n;w_m,{\mu }_m,{\sigma }_m\right)$ and $p_n$ is the distribution of hidden variables in the expectation–maximization method [11].

4.2. Parameter Estimation of L–L Distribution

The L–L parameters ($w_m$, ${\mu }_m$, and ${\sigma }_m$) can be determined through the actual wind speed data. The L–L parameters can be derived by Equations (15)–(18):

$$p_n=p\left(m|v_n;w_m,{\mu }_m,{\sigma }_m\right)=\frac{\frac{w_m}{v_n{\sigma }_m\sqrt{2\pi }}exp[\frac{-{\left(ln\left(v_n\right)-{\mu }_m\right)}^2}{2{\sigma }_m^2}]}{{{\displaystyle \sum }}_{m=1}^M\frac{w_m}{v_n{\sigma }_m\sqrt{2\pi }}exp[\frac{-{\left(ln\left(v_n\right)-{\mu }_m\right)}^2}{2{\sigma }_m^2}]}$$

(15)

$${\widehat{w}}_m=\frac{1}{N}{{\displaystyle \sum }}_{n=1}^N{p}_n$$

(16)

$${\widehat{\mu }}_m=\frac{{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)p_n\right]}{{{\displaystyle \sum }}_{n=1}^N{p}_n}$$

(17)

$${\widehat{\sigma }}_m={\{\frac{{{\displaystyle \sum }}_{n=1}^N[{\left(ln{v}_n-{\widehat{\mu }}_m\right)}^2{p}_n]}{{{\displaystyle \sum }}_{n=1}^N{p}_n}\}}^{\frac{1}{2}}$$

(18)

where N is the number of measured wind speed data; $p_n=p\left(m|v_n;w_m,{\mu }_m,{\sigma }_m\right)$ and $p_n$ is the distribution of hidden variables in the expectation–maximization method [11].

4.3. Parameter Estimation of W–L Distribution

If we let $b={\beta }_1$ and $c={{\alpha }_1}^{{\beta }_1}$, the W–W in Equation (5) can also be written as follows (Equation (19)):

$$f_V\left(v\right)=w_1\frac{b}{c}{v}^{b-1}exp[-\frac{v^b}{c}]+w_2\frac{1}{v\mathsf{\sigma }\sqrt{2\pi }}exp[\frac{-{\left(ln\left(v\right)-\mu \right)}^2}{2{\sigma }^2}]$$

(19)

The W–L parameters ($w_1$, $w_2$, $b$, $c$, $\mu$, and $\sigma$) can be determined through the actual wind speed data. The W–L parameters can be derived by Equations (20)–(26):

$$p\left(m|v_n;\Theta \right)=p\left(m|v_n;w_1,w_2,b,c,\mu ,\sigma \right)=\{\begin{array}{c}\frac{w_1\frac{b}{c}{v}_n^{b-1}exp[-\frac{v_n^b}{c}]}{f_V\left(v_n\right)} ,m=1\\ \frac{w_2\frac{1}{v_n\sigma \sqrt{2\pi }}exp[\frac{-{\left(ln\left(v_n\right)-\mu \right)}^2}{2{\sigma }^2}]}{f_V\left(v_n\right)},m=2\end{array}$$

(20)

$${\widehat{w}}_1=\frac{1}{N}{{\displaystyle \sum }}_{n=1}^{N}p\left(m=1|v_n;\Theta \right)$$

(21)

$${\widehat{w}}_2=\frac{1}{N}{{\displaystyle \sum }}_{n=1}^{N}p\left(m=2|v_n;\Theta \right)$$

(22)

$$\frac{{{\displaystyle \sum }}_{n=1}^{N}p\left(m=1|v_n;\Theta \right)}{\widehat{b}}+{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)p\left(m=1|v_n;\Theta \right)\right]-\frac{{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)v_n^{\widehat{b}}p\left(m=1|v_n;\Theta \right)\right]}{\widehat{c}}=0$$

(23)

$${{\displaystyle \sum }}_{n=1}^N\left(v_n^{\widehat{b}}p\left(m=1|v_n;\Theta \right)\right)-\widehat{c}{{\displaystyle \sum }}_{n=1}^{N}p\left(m=1|v_n;\Theta \right)=0$$

(24)

$$\widehat{\mu }=\frac{{{\displaystyle \sum }}_{n=1}^N\left[ln\left(v_n\right)p\left(m=2|v_n;\Theta \right)\right]}{{{\displaystyle \sum }}_{n=1}^{N}p\left(m=2|v_n;\Theta \right)}$$

(25)

$$\widehat{\sigma }={\{\frac{{{\displaystyle \sum }}_{n=1}^N[{\left(ln{v}_n-\widehat{\mu }\right)}^{2}p\left(m=2|v_n;\Theta \right)]}{{{\displaystyle \sum }}_{n=1}^{N}p\left(m=2|v_n;\Theta \right)}\}}^{\frac{1}{2}}$$

(26)

where N is the number of measured wind speed data; $\Theta =\left\{w_1,w_2,b,c,\mu ,\sigma \right\}$ and $p\left(m|v_n;\Theta \right)$ is the distribution of hidden variables in the expectation–maximization method [11].

4.4. Parameter Estimation of FVMF

The FVMF in Equation (6) can also be written as follows (Equation (27)):

$$f_{\Theta }\left(\theta \right)={{\displaystyle \sum }}_{m=1}^M\frac{w_m}{2\pi I_0\left(k_m\right)}exp\left\{k_m[\cos \left(\theta \right)\cos \left({\mu }_m\right)+\sin \left(\theta \right)\sin \left({\mu }_m\right)\right\}$$

(27)

where $\sin ({\mu }_m)$ and $\cos ({\mu }_m)$ are considered two parameters, and ${\sin }^2({\mu }_m)+{\cos }^2({\mu }_m)=1$.

The FVMF parameters {$w_m$, $\sin \left({\mu }_m\right)$, $\cos \left({\mu }_m\right)$, $k_m$} can be determined through the actual wind direction data [25]. The FVMF parameters can be derived by Equations (28)–(32):

$$p_n=p\left(m|{\theta }_n;w_m,\sin \left({\mu }_m\right),\cos \left({\mu }_m\right),k_m\right)=\frac{\frac{w_m}{2\pi I_0\left(k_m\right)}exp\left\{k_m[\cos \left({\theta }_n\right)\cos \left({\mu }_m\right)+\sin \left({\theta }_n\right)\sin \left({\mu }_m\right)\right\}}{{{\displaystyle \sum }}_{m=1}^M\frac{w_m}{2\pi I_0\left(k_m\right)}exp\left\{k_m[\cos \left({\theta }_n\right)\cos \left({\mu }_m\right)+\sin \left({\theta }_n\right)\sin \left({\mu }_m\right)\right\}}$$

(28)

$${\widehat{w}}_m=\frac{1}{N}{{\displaystyle \sum }}_{n=1}^N{p}_n$$

(29)

$${\sin }^2\left({\widehat{\mu }}_m\right)=\frac{{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\sin \left({\theta }_n\right)\right]\}}^2}{{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\sin \left({\theta }_n\right)\right]\}}^2+{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\cos \left({\theta }_n\right)\right]\}}^2}$$

(30)

$${\cos }^2\left({\widehat{\mu }}_m\right)=\frac{{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\cos \left({\theta }_n\right)\right]\}}^2}{{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\cos \left({\theta }_n\right)\right]\}}^2+{\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\sin \left({\theta }_n\right)\right]\}}^2}$$

(31)

$$\frac{I_1({\widehat{k}}_m)}{I_0({\widehat{k}}_m)}=\frac{\cos \left({\widehat{\mu }}_m\right)\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\cos \left(x_n\right)\right]\}+\sin \left({\widehat{\mu }}_m\right)\{{{\displaystyle \sum }}_{n=1}^N\left[p_n\sin \left(x_n\right)\right]\}}{{{\displaystyle \sum }}_{n=1}^N{p}_n}$$

(32)

where N is the number of measured wind direction data; $p_n$ is the distribution of hidden variables in the expectation–maximization method [11]; $I_1({\widehat{k}}_m)$ is the modified Bessel function of the first kind and order one. Here, the derivative of $I_0\left(k_m\right)$ with respect to $k_m$ is $I_1\left(k_m\right)$. Besides, to evaluate the goodness of fit for the theoretical distribution [33], the coefficient of determination (R²) and Akaike’s information criteria (AIC) are given by Equations (33) and (34):

$$R^2=\frac{{{\displaystyle \sum }}_{n=1}^N{({\widehat{F}}_n-\overline{F})}^2}{{{\displaystyle \sum }}_{n=1}^N{({\widehat{F}}_n-\overline{F})}^2+{{\displaystyle \sum }}_{n=1}^N{(F_n-{\widehat{F}}_n)}^2}$$

(33)

where ${\widehat{F}}_n$ are estimated cumulative probabilities; $\overline{F}={{\displaystyle \sum }}_{n=1}^N{\widehat{F}}_n/N$; $F_n$ are empirical cumulative probabilities.

$$AIC=-2log\left(L\right)+2k$$

(34)

where $L$ is the likelihood; $k$ is the number of parameter in the fitted model.

5. Wind Power Density

Wind power density can be used to evaluate the potential of wind energy in a region. The mean wind power density [35,36,37] can be expressed as follows (Equation (35)):

$$\overline{E}=\frac{1}{2N}{{\displaystyle \sum }}_{n=1}^N\rho v_n^3$$

(35)

where $\overline{E}$ represent the mean wind power density (W/m²); N is the number of records in the set period; $v_i$ is the wind speed of the ith record (m/s); $\rho$ is the density of air (kg/m³). In this paper, the value of ρ was 1.225 kg/m³. The probability distribution of wind power density at different wind speed can be determined using Equation (36):

$$P\left(v\right)=\frac{1}{2}\rho f_V\left(v\right)v^3$$

(36)

In addition, the wind power density can be estimated by Equation (37):

$$E_e={{\displaystyle \int }}_0^{\infty }\frac{1}{2}\rho f_V\left(v\right)v^{3}dv$$

(37)

The joint probability distribution of wind direction and speed can be used to evaluate the wind speed distribution at different wind directions, given by Equation (38):

$$f_{\left(\theta 1-\theta 2\right)}\left(v\right)={{\displaystyle \int }}_{\theta 1}^{\theta 2}{f}_{V,\Theta }d\theta$$

(38)

Moreover, two useful wind speeds parameters, namely the most probable wind speed and the wind speed carrying maximum energy, can be derived by the wind speed distribution and wind power density distribution [38,39]. The most probable wind speed is the wind speed of maximum point in the wind speed distribution. The wind speed carrying maximum energy is the wind speed of maximum point in the wind power density distribution [40].

6. Results and Discussion

Figure 1 was drawn by analyzing the data at a height of 70 m, and it demonstrates the wind speed distributions (Figure 1a) and the probability distribution of wind power density (Figure 1b). The parameters in the wind speed distributions, the corresponding R² coefficient, and AIC are shown in

Table 1. Fitted parameters in wind speed distributions.

Distribution	m	${\mathit{w}}_{\mathit{m}}$	${\mathit{\alpha }}_{\mathit{m}}/{\mathit{\mu }}_{\mathit{m}}$	${\mathit{\beta }}_{\mathit{m}}\left(\mathit{m}/\mathit{s}\right)/{\mathit{\sigma }}_{\mathit{m}}$	R2	AIC
Weibull–Weibull	1	0.5382	7.4895	2.4639	0.9942	1073705
	2	0.4617	6.0398	1.5558
Lognormal–Lognormal	1	0.7794	1.8437	0.4342	0.9945	1079551
	2	0.2206	0.8139	0.8520
Weibull–Lognormal	1	0.7594	6.1016	1.7790	0.9943	1073344
	2	0.2405	2.0413	0.3324

. It can be seen from Figure 1a that W–L distribution was better than W–W and L–L distribution for the histogram of wind speed data, which is also clear from in Table 1. Therefore, we chose W–L distribution as the marginal distribution of wind speed to construct the joint distribution. Figure 1b illustrates the probability density distribution of wind power density. The estimated wind power density was 270.79 W/m², 307.34 W/m², and 271.97 W/m² corresponding to W–W, L–L, and W–L, respectively. The mean wind power density was 270.53 W/m².

Figure 2 shows that the wind direction rose at a height of 70 m (Figure 2a) and the wind energy rose at a height of 70 m (Figure 2b). As shown in Figure 2, wind direction 1, 8, and 13 (N, SSE, and W) had a large wind direction frequency of 9.28%, 7.66%, and 10.41%, respectively, while wind direction 12 (WSW) had a large wind power density frequency (Figure 2b) of 14.92%. Therefore, in order to investigate their internal correspondence, we needed to study the wind speed characteristics under different wind directions and investigate their characteristics.

Figure 3 was obtained by numerical calculation of the circular variable, which was measured at a height 70 m. Figure 3a shows the wind direction distributions, and Figure 3b depicts the probability distribution of the angular variable ζ. It can be seen from the histogram of the measured wind direction data in Figure 3a that there were up to four main wind directions; therefore, the FVMF with M = 2, M = 3, and M = 4 was chosen to fit the wind direction distribution. Figure 3b illustrates that there were two distinct peaks in the histogram of the angular variable ζ; Therefore, the FVMF with M = 2 was selected to construct the probability distribution of the angular variable ζ. The parameters in the wind direction distributions, the corresponding R² coefficient, and AIC are shown in

Table 2. Fitted parameters in wind direction models.

Distribution	m	${\mathit{w}}_{\mathit{m}}$	${\mathit{k}}_{\mathit{m}}$	$\sin \left({\mathit{\mu }}_{\mathit{m}}\right)$	$\cos \left({\mathit{\mu }}_{\mathit{m}}\right)$	R2	AIC
FVMF with M = 2	1	0.8604	0.7638	−0.9315	0.3639	0.9890	740279
	2	0.1396	6.7758	0.4512	−0.8924
FVMF with M = 3	1	0.2414	5.2855	−0.9904	−0.1379	0.9883	737192
	2	0.3587	2.7048	−0.3228	0.9465
	3	0.3999	1.6529	0.3105	−0.9506
FVMF with M = 4	1	0.1999	5.2058	0.4319	−0.9019	0.9864	735538
	2	0.1567	8.3011	−0.1414	0.9900
	3	0.5610	1.5859	−0.9997	−0.0263
	4	0.0825	2.7812	0.9408	0.3388

. The parameters in the angular variable ζ distribution function and the corresponding R² coefficient are shown in

Table 3. Fitted parameters in angular variable ζ.

Distribution	m	${\mathit{w}}_{\mathit{m}}$	${\mathit{k}}_{\mathit{m}}$	$\sin \left({\mathit{\mu }}_{\mathit{m}}\right)$	$\cos \left({\mathit{\mu }}_{\mathit{m}}\right)$	R2
FVMF with M = 2	1	0.2544	0.9229	−0.0856	0.9963	0.9906
	2	0.7456	0.4296	0.9019	−0.4318

. From Figure 3a and Table 2, it can be seen that the FVMF with M = 4 was the best fitted model for the wind direction data. In addition, in order to estimate the parameters of the FVMF model, we redefined $\sin \left({\mu }_m\right)$ and $\cos \left({\mu }_m\right)$ instead of ${\mu }_m$. Here, the ${\mu }_m$ should be calculated by $\sin \left({\mu }_m\right)$ and $\cos \left({\mu }_m\right)$ to describe the mean wind direction. The FVMF with M = 4 reflected that the most dominant wind direction was 4.69 rad (268.6°). As displayed in Figure 3b and Table 3, the FVMF with M = 2 was suitable to describe the angular variable ζ distribution.

Figure 4 shows the joint probability distribution of actual wind direction and speed at a height of 70 m. Figure 5 shows the joint probability distribution at a height of 70 m, which was obtained with the angular–linear distribution model. According to the proposed joint probability distribution at a height 70 m, the wind speed distribution under 16 wind directions could be determined. The corresponding R² value can be seen in

Table 4. The R² value of the 70-m high wind speed distribution under 16 wind directions.

Wind Direction	N	NNE	NE	ENE	E	ESE	SE	SSE	S
R2	0.9831	0.9122	0.909	0.9341	0.9137	0.9354	0.9799	0.9997	0.9886
Wind Direction	SSW	SW	WSW	W	WNW	NW	NNW
R2	0.9853	0.9958	0.971	0.9848	0.997	0.985	0.9869

. As can be seen from Figure 4 and Figure 5 and Table 4, the proposed model had a great degree of fit to the sample data at a height 70 m. Moreover, in order to verify the validity of the selected model, the proposed model was used to construct the joint probability distribution at a height of 50 m. It should be borne in mind that the wind direction varied slightly from 70 m to 50 m.

Figure 6 shows the joint probability distribution at a height of 50 m. In the proposed joint distribution model, W–L distribution was selected to process the wind speed data, and the FVMF with M = 4 was the fitted model for the wind direction data. Moreover, the FVMF with M = 2 was chosen to describe the angular variable ζ distribution. The parameters in the joint probability distribution at a height of 50 m and the corresponding R² coefficient are shown in

Table 5. Fitted parameters in the joint probability distribution at a height of 50 m.

Distribution	m	${\mathit{w}}_{\mathit{m}}$	${\mathit{\alpha }}_{\mathit{m}}/{\mathit{\mu }}_{\mathit{m}}$	${\mathit{\beta }}_{\mathit{m}}\left(\mathit{m}/\mathit{s}\right) / {\mathit{\sigma }}_{\mathit{m}}$	${\mathit{k}}_{\mathit{m}}$	$\sin \left({\mathit{\mu }}_{\mathit{m}}\right)$	$\cos \left({\mathit{\mu }}_{\mathit{m}}\right)$	R2
W–L	1	0.6108	5.3568	1.6104	—	—	—	0.9933
	2	0.3892	1.8266	0.3853	—	—	—
FVMF with M = 4	1	0.1999	—	—	0.1999	5.2058	0.4319	0.9864
	2	0.1567	—	—	0.1567	8.3011	−0.1414
	3	0.5610	—	—	0.5610	1.5859	−0.9997
	4	0.0825	—	—	0.0825	2.7812	0.9408
FVMF with M = 2	1	0.4410	—	—	0.6048	−0.1081	0.9941	0.9903
	2	0.5590	—	—	0.6030	0.8675	−0.4974

. According to the proposed joint probability distribution at a height of 50 m, the wind speed distribution under 16 wind directions could be determined. The corresponding R² value can be seen in

Table 6. The R² value of the 50-m high wind speed distribution under 16 wind directions.

Wind Direction	N	NNE	NE	ENE	E	ESE	SE	SSE	S
R2	0.9868	0.9288	0.9033	0.9139	0.9337	0.9468	0.9729	0.9835	0.9835
Wind Direction	SSW	SW	WSW	W	WNW	NW	NNW
R2	0.9808	0.9955	0.9737	0.9857	0.9981	0.9882	0.9896

. It can be seen from the above analysis that the proposed model could well present wind characteristics at different heights.

The wind direction was recorded in 16 intercardinal directions. Each wind direction of the 16 intercardinal directions had a corresponding wind speed and a wind power density distribution. Due to the limitation of space, we have only given the results of two wind directions—wind direction 12 (E) and wind direction 13 (W)—in Figure 7. The figure was drawn by applying the joint probability density distribution at a height of 70 m. As can be seen from Figure 7, the most probable wind speeds for wind direction 12 and 13 were 6.5 m/s and 6.95 m/s, respectively; the wind speeds carrying maximum energy were 9.34 m/s and 9.63 m/s, respectively.

The FVMF with M = 4 was a good fitted model for the wind direction data and could be used to investigate the wind direction characteristics for four seasons. The four seasons were defined as follows: (a) spring: March, April, and May; (b) summer: June, July, and August; (c) Autumn: September, October, and November; (d) winter: December, January, and February. Figure 8 depicts the wind direction distributions over four seasons. It was drawn using the wind direction data of 70 m for three years. The parameters in the wind direction model over four seasons and the corresponding R² coefficient are shown in

Table 7. Fitted parameters in wind direction models over four seasons.

Season	Distribution	m	${\mathit{w}}_{\mathit{m}}$	${\mathit{k}}_{\mathit{m}}$	$\sin \left({\mathit{\mu }}_{\mathit{m}}\right)$	$\cos \left({\mathit{\mu }}_{\mathit{m}}\right)$	R2
Spring	FVMF with M = 4	1	0.1640	6.8674	0.3582	−0.9336	0.9871
		2	0.0309	7.4065	0.9851	0.1721
		3	0.7488	1.0602	−0.9580	0.2869
		4	0.0563	23.9653	−0.0576	0.9983
Summer	FVMF with M = 4	1	0.3752	1.0583	−0.9999	−0.0100	0.9873
		2	0.3126	3.4264	0.5993	−0.8005
		3	0.2492	4.6498	−0.1368	0.9906
		4	0.0630	6.4372	0.8626	0.5058
Autumn	FVMF with M = 4	1	0.4161	3.1634	−0.9855	−0.1696	0.9859
		2	0.0599	3.4360	0.9456	0.3254
		3	0.2950	4.5123	0.3778	−0.9259
		4	0.2290	6.5931	−0.2147	0.9767
Winter	FVMF with M = 4	1	0.0803	8.5141	0.4261	−0.9047	0.984
		2	0.0310	4.0085	0.9529	0.3031
		3	0.7966	1.5032	−0.9989	−0.0477
		4	0.0922	22.7153	−0.0708	0.9975

. As can be seen from Figure 8, the most dominant wind direction was 5 rad (286.6°) in the spring season, 6 rad (348.8°) in the summer season, 4.5 rad (260.2°) in the autumn season, and 4.7 rad (267.1°) in the winter season. There were more than two prevailing directions in the four seasons. It was concluded from Figure 8 that the wind direction had seasonal features.

7. Conclusions

This paper presents the analysis of wind speed and wind direction using joint probability distribution methods. The joint distribution was obtained with wind direction and speed distributions. By fitting the wind data at different heights, the proposed model could well present wind characteristics. Thus, the proposed model was used to investigate the wind direction and wind speed. The main findings of this paper are as follows:

(1) W–W, L–L, and W–L distribution with higher R² coefficient provide a good fitting for wind speed data at a height of 70 m. In this paper, the best fitting model for wind speed was W–L distribution, which could accurately represent wind speed characteristics.

(2) The FVMF with M = 4 was a good fitted model for the wind direction data at a height of 70 m and could accurately describe the wind direction distribution with multiple peaks. Moreover, the dominant wind direction was 4.69 rad (268.6°) in the three years.

(3) The joint probability distribution provided a good fit to the measured wind data at different heights, and the R² of the joint probability distribution was greater than 0.99, indicating that this model was able to explain more than 99% of wind data. With the joint probability distribution, the most probable wind speeds for wind direction 12 and 13 were 6.5 m/s and 6.95 m/s at a height of 70 m, respectively; the wind speed carrying maximum energy were 9.34 m/s and 9.63 m/s at a height of 70 m, respectively.

(4) The wind direction had seasonal features. The most dominant wind direction was 5 rad (286.6°) in the spring season, 6 rad (348.8°) in the summer season, 4.5 rad (260.2°) in the autumn season, and 4.7 rad (267.1°) in the winter season. Moreover, there were more than two prevailing directions in the four seasons.