Three-Dimensional Wind Measurement Based on Ultrasonic Sensor Array and Multiple Signal Classification

The wind power industry continues to experience rapid growth worldwide. However, the fluctuations in wind speed and direction complicate the wind turbine control process and hinder the integration of wind power into the electrical grid. To maximize wind utilization, we propose to precisely measure the wind in a three-dimensional (3D) space, thus facilitating the process of wind turbine control. Natural wind is regarded as a 3D vector, whose direction and magnitude correspond to the wind’s direction and speed. A semi-conical ultrasonic sensor array is proposed to simultaneously measure the wind speed and direction in a 3D space. As the ultrasonic signal transmitted between the sensors is influenced by the wind and environment noise, a Multiple Signal Classification algorithm is adopted to estimate the wind information from the received signal. The estimate’s accuracy is evaluated in terms of root mean square error and mean absolute error. The robustness of the proposed method is evaluated by the type A evaluation of standard uncertainty under a varying signal-to-noise ratio. Simulation results validate the accuracy and anti-noise performance of the proposed method, whose estimated wind speed and direction errors converge to zero when the SNR is over 15 dB.


Introduction
As wind energy is considered one of the most promising renewable energies, it is widely used for electric power generation [1,2]. However, despite the continuing increase in installed capacity of wind power worldwide, the wind power that is actually utilized is much less than would be expected, because of the unstable, intermittent, and highly volatile nature of wind. Furthermore, due to the characteristics of randomness and fluctuation, the electricity generated by the wind will be excessive or insufficient, making it difficult to integrate into the power grid [3,4]. On one hand, one possibility aiming at the redundant electricity is to simply abandon the fluctuating wind power [5,6]. However, that would come with not only a huge loss of electricity but also a mass waste of wind resources and wind power equipment. According to a report from the National Energy Administration of China [7], the national average rate of wind abandonment in 2018 is~7%, which is equivalent to 27.7 billion kWh, resulting in an enormous loss of about US $2.1 billion. On the other hand, a solution for overcoming the deficiency of wind power generation is to store the wind energy. Hybrid wind/compressed air energy storage (CAES) systems are used to transform the intermittent wind resources into a constantly available power supply in Germany, the USA, and even remote Arctic areas [8]. However, the underground geology may bring risks to practical applications of CAES [9]. In contrast, accurate wind measurement in a 3D space can provide data that can be used to control wind horizontal and vertical components, namely, the azimuth and pitch angles, respectively. We adopted the MUSIC algorithm for wind measurement due to two advantages: the first is that the MUSIC algorithm makes the measurement accuracy independent of the precision of the signal transmission time in the wind. In addition, the environment noise and the sampling noise are suppressed, which reduces the influence of noise on the accuracy of measurement results [35]. The second advantage is that the MUSIC algorithm can be combined with a variety of sensor arrays and applied to different scenarios. Except for the arc array mentioned by Li et al. [35], the MUSIC algorithm can also be used on a uniform circular array [38], which estimates the dual angle of azimuth and pitch angle of the incident signal. In addition, other array structures, such as the uniform linear (ULA) [39], L-shaped array [40], and conformal array [41], can be combined with the MUSIC algorithm for diverse scenarios and purposes. To measure these components simultaneously, especially adding the measurement of the pitch angle, we propose a novel semi-conical ultrasonic sensor array in 3D space, which is different from the planar array proposed by Li et al. [35]. The array consists of six ultrasonic sensors, where the transmitting sensor is on the peak of the cone, and the five receiving sensors are evenly positioned on the undersurface arc of the cone. According to the spatial relationship between the transmitter and the receiver, we simulate a 3D wind vector, and calculate the theoretic transmission time of ultrasonic signal. The theoretic transmission time is then used to determine the array manifold matrix, which is critical for the implementation of the MUSIC algorithm. After receiving the true ultrasonic signal, the MUSIC algorithm divides the covariance matrix of the received information into signal and noise subspaces. Under ideal conditions, the signal and noise subspaces are completely orthogonal, and can be used to calculate the spectral function. By searching for the peak value of the spectrum with the speed searching step of 0.1 m/s and the direction searching step of 1º, we can obtain the corresponding information of wind speed, pitch angle, and azimuth angle simultaneously. Simulation results indicate that the proposed method can measure the azimuth angle in (0 • , 360 • ), the pitch angle in (0 • , 90 • ). Furthermore, the proposed method demonstrates better performances than the state-of-the-art method in terms of convergence speed, estimation error, and variance.
The rest of the study is organized as follows. Section 2 describes the structure of the proposed sensor array and the wind measurement principle. Section 3 discusses the simulations and a comparison with the state-of-the-art method, which analyzes the accuracy and anti-noise performances of the methods. Section 4 includes our conclusions and suggestions for further research.

The Ultrasonic Sensor Array
Wind measurements can be based on several types of ultrasonic sensor array structures, of which the arc ultrasonic sensor array and the ULA structure are the most prominent [35,39]. When combined with the MUSIC algorithm, the arc array structure takes precedence over the ULA structure in terms of estimate variance [39]. However, these experiments were performed in the 2D space. To measure the natural wind in the 3D space, we propose a novel semi-conical structure sensor array based on the arc array structure. The sensor array consists of six ultrasonic sensors, numbered 0-5, wherein Sensor 0 is used for transmitting ultrasonic signals, and Sensors 1-5 are receiving sensors. Figure 1 depicts the semi-conical structure in the 3D space. The transmitting ultrasonic sensor (Sensor 0) is located on the vertex of the cone, with height H. The five receiving ultrasonic sensors (Sensors 1-5) are evenly positioned on the undersurface arc of the cone. To clearly display the positional relationship of the five receiving sensors on the horizontal plane, we plot the undersurface of the cone in Figure 2. As illustrated in Figures 1 and 2, the five receiving sensors are evenly located on the semi arc with the center O and the radius R, where the angle between every two adjacent receiving sensors is α. The distance D between the transmitting sensor and each receiving sensor can be deduced by D =

Premise Assumptions
The transmission of signal through wireless channel is complicated, which makes it difficult to completely describe the physical environment and establish a rigorous mathematical model. As mentioned in the literature [35], the arrival signals of the wavefront in the array system can be considered a spherical wave. Therefore, the propagation time of the signal to each sensor is not only related to the direction, but also to the distance of the signal travels. In this case, the wind speed and wind direction information can be estimated with the MUSIC algorithm on the grounds that the output data matrix of the sensor array is only related to the wind speed and wind direction.
Although the wind speed and direction can be regarded as continuous signals, we discretize them by the sampling frequency of 100 MHz, resulting in the short sampling period of 10 ns.

Premise Assumptions
The transmission of signal through wireless channel is complicated, which makes it difficult to completely describe the physical environment and establish a rigorous mathematical model. As mentioned in the literature [35], the arrival signals of the wavefront in the array system can be considered a spherical wave. Therefore, the propagation time of the signal to each sensor is not only related to the direction, but also to the distance of the signal travels. In this case, the wind speed and wind direction information can be estimated with the MUSIC algorithm on the grounds that the output data matrix of the sensor array is only related to the wind speed and wind direction.
Although the wind speed and direction can be regarded as continuous signals, we discretize them by the sampling frequency of 100 MHz, resulting in the short sampling period of 10 ns.

Premise Assumptions
The transmission of signal through wireless channel is complicated, which makes it difficult to completely describe the physical environment and establish a rigorous mathematical model. As mentioned in the literature [35], the arrival signals of the wavefront in the array system can be considered a spherical wave. Therefore, the propagation time of the signal to each sensor is not only related to the direction, but also to the distance of the signal travels. In this case, the wind speed and wind direction information can be estimated with the MUSIC algorithm on the grounds that the output data matrix of the sensor array is only related to the wind speed and wind direction. Although the wind speed and direction can be regarded as continuous signals, we discretize them by the sampling frequency of 100 MHz, resulting in the short sampling period of 10 ns. According to the validity of Taylor's hypothesis of frozen turbulence [42], we assume that the amplitude of the wind speed and direction are approximately constant throughout a sampling period. The propagation speed of the signal on each path is only affected by the component of the wind vector parallel to the propagation path as in transit-time Ultrasonic Flowmeter [43,44]. The key to adopting the MUSIC algorithm is to build the array manifold matrix, which implies the information of wind speed and direction. Ultimately, the measured wind information can be regarded as the average ones over the sampling period.

Signal Model
In general, when M narrowband signals are incident on an array containing N array elements, the received signal can be expressed by the following complex envelope form, where u i (t) represents the amplitude of the signal. ϕ(t) denotes the phase of signal. ω 0 is the angular frequency of signal and ω 0 = 2π f , where f is the frequency of signal. According to the narrowband assumption, we can obtain: Therefore, it can be obtained that the expression of all signals received by the k th array element in the array is as follows; is the noise of k th array element at time t, τ ki denotes the transmission time when the i th signal reaches the k th array element, and g ki is the gain of i th signal on the k th array element. Based on the assumption that all array elements are isotropic and the coupling between them can be negligible, the g ki can be averaged to 1. Combined with Equation (2), the signal received by k th array element can be expressed as follows: Because there are five receiving sensors and only one signal resource in this paper, that is, N = 5, M = 1, we use S(t) = s(t) to denote the original signal transmitted by the Sensor 0. The received signal x k (t) at the k th sensor can be expressed as follows, where k = 1, 2, · · · 5, and τ k denotes the transmitting time. n k (t) denotes the noise on the k th receiving sensor, which is assumed to be an independently and identically distributed (i.i.d.) AWGN.
To speed up signal processing, we vectorize the received signal as follows.
Equation (4) can be further rewritten in the matrix form: where N(t) is the noise matrix and A is the array manifold matrix described by Equation (6): It is obvious that the key to obtaining array manifold matrix A is to obtain the transmission time of the ultrasonic signal to each receiving sensor. The transmission time can be calculated by analyzing the structure of the sensor array and the influence of natural wind on the signal propagation. As a result, the wind measurement problem is transformed into an identification problem of wind speed V, azimuth angle θ, and pitch angle β from the array manifold matrix A. The influence of the natural wind on the ultrasonic signal is thus reduced to the signal transmitting time τ in matrix A. In the case of no wind, τ should be identical at all five receiving sensors, which is where c = 340 m/s is the propagation speed of ultrasonic wave. When the transmission of the ultrasonic signal is influenced by the wind, τ i would be different for each receiving sensor, resulting in a more complex calculation. To simplify the calculation, we decompose the natural wind speed V into a vertical wind speed V t and a horizontal wind speed V h by vector decomposition, accordingly Similarly, the original propagation speed c of ultrasonic wave is decomposed into c h and c t , which are as follows, However, under the influence of wind, the actual speed of the ultrasonic signal is changed. We use V i to denote the ultrasonic signal speed received at the i th sensor. The value of V i is calculated based on the principle of vector decomposition and the positional relationship of the sensors in the 3D space. Accordingly, the horizontal components of V i are as follows, The vertical components V ti are identical for all the receiving sensors, which can be calculated by Therefore, we can infer V i by synthesizing V hi and V ti as follows, Accordingly, the transmission time τ i of the signal at the i th receiving sensor can be represented as follows, Sensors 2020, 20, 523 7 of 16 By combining Equations (8)-(13), we can obtain the relationship between the wind and the transmission time τ of the ultrasonic signal. However, it is not easy to retrieve the wind information V, θ, and β directly from the transmission time. Therefore, the MUSIC algorithm is used to identify the wind speed V, azimuth angle θ, and pitch angle β, as implied in the matrix.

Wind Measurement Based on MUSIC Algorithm
The process of the MUSIC algorithm is to first calculate the covariance matrix R X of the received signal X(t); next, perform eigenvalue decomposition on the obtained covariance matrix R X , resulting in a signal subspace U S and noise subspace U N ; and finally, based on the orthogonality of the array manifold matrix A and the noise subspace U N , the corresponding V, θ, and β can be obtained by performing spectral peak searching.
Specifically, the covariance matrix R X of the signal matrix X(t) can be obtained by Equation (14): where H denotes the conjugate transpose. R S and R N are the correlation matrices of signal and noise, respectively. σ 2 is the variance of noise, and I is the unit matrix.
As R X is symmetric, we can conduct eigenvalue decomposition on it, resulting in Therefore, σ 2 is the smallest eigenvalue of R X , which corresponds to N − M eigenvectors. These N − M eigenvectors comprise the noise subspace U N . The other M eigenvalues are relevant to the signals, whose eigenvectors constitute U S , the signal subspace. R X is a Hermitian matrix, thus the signal subspace U S and noise subspace U N are orthogonal, namely, U H S U N = 0. We can further infer that However, referring to Equation (14), we can also get R X U N = AR S A H U N + σ 2 IU N . Therefore, AR S A H U N = 0. We left multiply the In practical applications, R X cannot be directly obtained, while only the sampled covarianceR X can be calculated byR where L is the number of samplings. Theoretically, when L → ∞ ,R X and R X are consistent. In our experiment, we set L = 1000 to provide the best possible approximation of R X . Similarly, we can also perform eigenvalue decomposition onR X : where UŜ indicates the signal subspace consisting of the signal eigenvectors and UN represents the noise subspace. However, due to the deviation betweenR X and R X , A H UN 0. To retrieve the wind information V, θ, and β from the matrix A, the MUSIC algorithm constructs a spectral estimation formula: The maximum value of P music , namely the spectrum peak, corresponds to the smallest value of A H UN. According to Equations (4), (6), and (18), we calculate the theoretical maximum value of P music by searching through the measuring range of V, θ, and β. The corresponding values of V, θ, and β that generate the spectral peak are the measured results of the wind speed and direction information. Therefore, the wind speed V, azimuth angle θ, and pitch angle β of the natural wind in the 3D space are simultaneously obtained.

Simulations
Simulations are conducted using the ultrasonic sensor array shown in Figure 1, where the radius R and the height H are both set to 1 m, and the angle α between every two adjacent receiving sensors is 30 • . The propagation speed c and frequency f of the transmitted ultrasonic signal are 340 m/s and 40 kHz, respectively. The snapshot number L is 1000. With respect to the spectral search, we set the azimuth angle searching range to (0, 360 • ), the pitch angle range to (0, 90 • ), and the searching step length to 1 • . The searching range of the wind speed is (0, 60) m/s with the searching step length of 0.1 m/s for the MUSIC algorithm.
Based on the above experimental parameters, we performed the simulations. The noise on the receiving sensor was simulated by an AWGN with zero-mean, whose variance is the average power of noise. Figure 3 depicts the spectral searching result of SNR at 20 dB. As illustrated in Figure 3, the darkest red color is concentrated at the position corresponding to the peak value of the spectral function. In this case, the wind speed V = 35.7 m/s, wind azimuth angle θ = 146 • , and wind pitch angle β = 82 • . Figure 3b-d shows the slice results in the direction of wind speed, wind azimuth angle, and wind pitch angle, which demonstrate the estimation results of these components. These results are consistent with those in Figure 3a. of noise. Figure 3 depicts the spectral searching result of SNR at 20 dB. As illustrated in Figure 3, the darkest red color is concentrated at the position corresponding to the peak value of the spectral function. In this case, the wind speed = 35.7 m/s, wind azimuth angle = 146°, and wind pitch angle = 82°. Figure 3b-d shows the slice results in the direction of wind speed, wind azimuth angle, and wind pitch angle, which demonstrate the estimation results of these components. These results are consistent with those in Figure 3a.  To evaluate the stability of the proposed method, we performed 100 independent simulations with SNR ranging from 0 to 40 dB. The measurement deviation was quantitatively analyzed using the type A evaluation of standard uncertainty in Figure 4. As can be seen in Figure 4, the measurement uncertainty decreases with higher SNR. Specifically, the measurement uncertainty of wind speed V is 0.69 in the case of SNR at 0 dB, which sharply decreases to 0.09 at SNR of 5 dB, then drops to zero when the SNR is higher than or equal to 15 dB. The measurement uncertainty of pitch angle β shows a similar pattern, which begins at 2.35 at SNR = 0 dB, then falls to 0.42 at SNR = 5 dB and reaches zero at SNR ≥ 15 dB. However, the trend of the measurement uncertainty of wind azimuth angle θ is slightly different from those of the wind speed and pitch angle. The biggest uncertainty value of it is 2.41 when SNR = 0 dB, which then steeply falls to and remains at zero, since SNR is over 5 dB.
drops to zero when the SNR is higher than or equal to 15 dB. The measurement uncertainty of pitch angle shows a similar pattern, which begins at 2.35 at SNR = 0 dB, then falls to 0.42 at SNR = 5 dB and reaches zero at SNR ≥ 15 dB. However, the trend of the measurement uncertainty of wind azimuth angle is slightly different from those of the wind speed and pitch angle. The biggest uncertainty value of it is 2.41 when SNR = 0 dB, which then steeply falls to and remains at zero, since SNR is over 5 dB.  The estimation errors are measured by Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) in Figures 5 and 6, respectively. Both of these figures illustrate quite consistent and similar patterns. The average errors of the wind speed and pitch angle converge to zero when SNR ≥ 15 dB, and these trends are consistent with those of the measurement uncertainty indicated by Figure 4a,b. Similarly, the RMSE and MAE of wind azimuth angle achieved zero, as SNR at 5 dB, whose variation under different SNRs is in line with what is shown in Figure 4c. Therefore, we The estimation errors are measured by Root Mean Square Error (RMSE) and Mean Absolute Error (MAE) in Figures 5 and 6, respectively. Both of these figures illustrate quite consistent and similar patterns. The average errors of the wind speed V and pitch angle β converge to zero when SNR ≥ 15 dB, and these trends are consistent with those of the measurement uncertainty indicated by Figure 4a,b. Similarly, the RMSE and MAE of wind azimuth angle θ achieved zero, as SNR at 5 dB, whose variation under different SNRs is in line with what is shown in Figure 4c. Therefore, we can conclude that the proposed method is more reliable and more accurate with higher SNR, which is greater than 15 dB, to be precise. can conclude that the proposed method is more reliable and more accurate with higher SNR, which is greater than 15 dB, to be precise.  can conclude that the proposed method is more reliable and more accurate with higher SNR, which is greater than 15 dB, to be precise.

Comparison with State-Of-The-Art Method
The measurements of natural wind based on the MUSIC algorithm are widely performed in the 2D space. One of the representative methods for wind speed and direction measurement in the 2D space (WSDM2D) was proposed in [35], based on an arc ultrasonic sensor array. To compare the proposed method with WSDM2D, we set the pitch angle to 90° to simulate the wind on the horizontal plane. The wind estimation problem was thus relegated to the 2D space in this case. We randomly simulated 100 groups of wind speed and direction to perform our estimates, where the former was uniformly distributed within 0-60 m/s, and the latter was uniformly distributed within 0 ∼ 360°. Figure 7 illustrates this distribution of the wind speed and direction data. Monte Carlo simulations of 100 runs were carried out to compare the two methods with SNR ranging from 0 dB to 40 dB. Similar to the results in Section 3.1, the errors in the wind estimates of both the proposed method and WSDM2D converged to zero when SNR ≥ 15 dB. To perform a comprehensive comparison of the accuracy and anti-noise performance of the proposed method with WSDM2D, we calculated RMSE and MAE of speed and direction at SNR varying from 0 to 20 dB, the results of which are presented in the Table 1. Three findings can be deduced from Table 1: 1. The trends of RMSEs and MAEs are consistent, which all decrease with higher SNR.

Comparison with State-Of-The-Art Method
The measurements of natural wind based on the MUSIC algorithm are widely performed in the 2D space. One of the representative methods for wind speed and direction measurement in the 2D space (WSDM2D) was proposed in [35], based on an arc ultrasonic sensor array. To compare the proposed method with WSDM2D, we set the pitch angle β to 90 • to simulate the wind on the horizontal plane. The wind estimation problem was thus relegated to the 2D space in this case. We randomly simulated 100 groups of wind speed and direction to perform our estimates, where the former was uniformly distributed within 0-60 m/s, and the latter was uniformly distributed within 0 ∼ 360 • . Figure 7 illustrates this distribution of the wind speed and direction data.

Comparison with State-Of-The-Art Method
The measurements of natural wind based on the MUSIC algorithm are widely performed in the 2D space. One of the representative methods for wind speed and direction measurement in the 2D space (WSDM2D) was proposed in [35], based on an arc ultrasonic sensor array. To compare the proposed method with WSDM2D, we set the pitch angle to 90° to simulate the wind on the horizontal plane. The wind estimation problem was thus relegated to the 2D space in this case. We randomly simulated 100 groups of wind speed and direction to perform our estimates, where the former was uniformly distributed within 0-60 m/s, and the latter was uniformly distributed within 0 ∼ 360°. Figure 7 illustrates this distribution of the wind speed and direction data. Monte Carlo simulations of 100 runs were carried out to compare the two methods with SNR ranging from 0 dB to 40 dB. Similar to the results in Section 3.1, the errors in the wind estimates of both the proposed method and WSDM2D converged to zero when SNR ≥ 15 dB. To perform a comprehensive comparison of the accuracy and anti-noise performance of the proposed method with WSDM2D, we calculated RMSE and MAE of speed and direction at SNR varying from 0 to 20 dB, the results of which are presented in the Table 1. Three findings can be deduced from Table 1: 1. The trends of RMSEs and MAEs are consistent, which all decrease with higher SNR. Monte Carlo simulations of 100 runs were carried out to compare the two methods with SNR ranging from 0 dB to 40 dB. Similar to the results in Section 3.1, the errors in the wind estimates of both the proposed method and WSDM2D converged to zero when SNR ≥ 15 dB. To perform a comprehensive comparison of the accuracy and anti-noise performance of the proposed method with WSDM2D, we calculated RMSE and MAE of speed and direction at SNR varying from 0 to 20 dB, the results of which are presented in the Table 1. Three findings can be deduced from Table 1: 1.
The trends of RMSEs and MAEs are consistent, which all decrease with higher SNR.

2.
Accordingly, the biggest RMSE and MAE of wind speed and direction occur at SNR of 0 dB. The biggest wind speed RMSE of the proposed method is slightly bigger than that of WSDM2D, while the biggest direction RMSE is smaller. 3.
The speed and direction RMSEs of proposed method tend to zeros as SNR of 5 dB, outperforming the WSDM2D method, which converges to zero since SNR = 15 dB. Wind speed and direction MAEs of two methods have similar patterns with those of RMSEs.     In conclusion, the proposed method is more accurate and has better anti-noise performance than WSDM2D in wind speed and direction measurements, even in the 2D case.

Conclusions
We propose to measure the natural wind in the 3D space using an ultrasonic sensor array in the context of fluctuations in wind speed and production. More precisely, by building a novel semi-conical ultrasonic sensor array and using the MUSIC algorithm, we measure the wind speed , azimuth angle , and pitch angle simultaneously and accurately. The measurement deviation is quantified by the type A evaluation of standard uncertainty, and the estimation errors are evaluated In conclusion, the proposed method is more accurate and has better anti-noise performance than WSDM2D in wind speed and direction measurements, even in the 2D case.

Conclusions
We propose to measure the natural wind in the 3D space using an ultrasonic sensor array in the context of fluctuations in wind speed and production. More precisely, by building a novel semi-conical ultrasonic sensor array and using the MUSIC algorithm, we measure the wind speed V, azimuth angle θ, and pitch angle β simultaneously and accurately. The measurement deviation is quantified by the type A evaluation of standard uncertainty, and the estimation errors are evaluated using RMSE and MAE, respectively. The measurement errors of wind direction are within 1º, and those of wind speed within 0.1 m/s when the SNR is greater than 10 dB. Simulations show that the proposed method has better performance and higher accuracy when the SNR is greater than 15 dB, where errors are nearly zero. The measurement uncertainty and the estimation error illustrate similar patterns under the influence of varying SNRs, which further proves the consistency of the proposed method. Furthermore, we compared the accuracy of the proposed method with the state-of-the-art method in the 2D space. Monte Carlo simulations show that our method has higher accuracy and stronger noise immunity. Note that the structure of the sensor array is not limited to the one proposed in this study, and the accuracy of the measurement has room for improvement. To further verify the validity of the theory, we will perform practical experiments in the next step.