Determining Ultrasound Arrival Time by HHT and Kurtosis in Wind Speed Measurement

Abstract

The determination of ultrasonic echo signal onset time is the core of performing the time difference method to calculate wind speed. However, in practical cases, background noise makes precise determination extremely difficult. This paper carries out research on the accurate determination of onset time, exploring the advantages of an improved method based on the combination of Hilbert-Huang Transform (HHT) and high-order statistics (kurtosis). Performing Hilbert-Huang Transform to the received wave is aimed at determining a rough arrival time, around which a fixed size of data is extracted as initial sample to avoid a false pick. Then the fourth-order kurtosis of a smaller sample, extracted successively by a moving window from the initial sample, is calculated. The minimum point corresponds to the initial onset time. This approach was tested on a real ultrasonic echo signal dataset, acquired in a wind tunnel with an ultrasonic anemometer. The proposed method showed satisfying results in both ideal cases and low signal-to-noise ratio (SNR) environment, compared with traditional onset time determination approaches, including Akaike Information Criterion (AIC-picker), Short-term Average over Long-term Average (STA/LTA), and Teager-Kaiser energy operator (TKEO). The experimental results acquired by the HHT-kurtosis method demonstrated that the proposed method possesses a high accuracy.

1. Introduction

During the process of wind tunnel experiment, wind speed is one of the most critical parameters that needs to be accurately measured [1,2]. Precise determination of wind speed can be of great significance in many fields relating to weather forecasting, power generating, farming, engineering, and manufacturing [3,4]. For instance, in the field of wind power utilization, wind speed determination is vital for selecting sites and wind machines [5]. The traditional method generally adopts mechanical anemometers, thermal anemometers, or pitot tube anemometers. However, the anemometers mentioned above have some limitation. Mechanical anemometers can be easily damaged, especially for the rotating component inside. Besides, this component needs an initial speed to rotate due to its inertia, hence it is almost impossible to measure low-speed wind. Thermal anemometers have a limited measurement range, and its performance is heavily dependent on environmental factors, including temperature and noise. Due to the low accuracy of the pressure sensor, it is exorbitantly consuming to design a pitot tube anemometer to precisely detect the pressure difference of low-speed wind. In contrast, ultrasonic anemometer is a reliable method with a wide measurement range and no start-up wind speed limitation while being applied in the field of wind speed determination.

Recent mainstream approaches utilizing ultrasonic anemometers can be categorized into vortex-shedding method [6] and time-of-flight (TOF) method [7]. Measuring the vortex frequency, vortex-shedding method is based on the proportional relationship between ultrasound and wind speed. However, this approach is difficult to conduct and temperature influence must be considered. The idea of the TOF method is to implement a sensor array and measure a set of ultrasound traveling times between transmitters and receivers [8], and therefore both wind speed and wind direction can be determined. The mainstream application of TOF estimation is time difference method [9,10], whose theory is to calculate the travel time difference between tailwind ultrasound and headwind one, following the Equation (1),

$$Vcos\theta =\frac{l}{2}(\frac{1}{t_1}-\frac{1}{t_2}),$$

(1)

where $l$ is the distance between a pair of emitter and receiver, and $V$ is the wind speed, having an angle of $\theta$ by the east-west direction. As shown in Figure 1,$t_1$, and $t_2$ are the travel time of ultrasound in tailwinds and headwinds respectively.

As shown in Equation (1), in order to acquire a trustworthy wind speed, picking a precise arrival time is of vital importance. Previously, the automatic onset time determination was carried out by a picking algorithm, which generally falls into three categories. The simplest technique is to use an amplitude threshold picker, which applies an empirically determined value as a threshold to filter the received wave [11]. However, the determination of this value can be puzzling and this method is not applicable in a high-noise environment. An improved method in this category is to apply a short-term-average over long-term-average (STA/LTA) filter to avoid false pick [12], but window length determination remains difficult. The second type of automatic onset time picker is to apply a running extraction window along the time-axis and successively calculate the characteristic function value corresponding to each picked sample. As assumed, at the time of signal arrival, the value of this characteristic function changes significantly. In this category, prevailing methods include Akaike information criterion (AIC) picker [13,14], Teager-Kaiser energy operator (TKEO) [15], and kurtosis-based picker [16]. The final type is to calculate the matching rate between the received wave and a reference waveform [17,18], based on the assumption of the existence of a reference waveform similar to the received wave. However, the calculated onset time corresponding to the maximum matching rate can be easily influenced by the shape of echo signal in the time domain.

In order to get a time-frequency spectrum, many signal processing methods have been advanced [19], including Wavelet Transform (WT) [20], Short Time Fourier Transform (STFT) [21], and Hilbert-Huang Transform (HHT) [22]. HHT is a superior algorithm among them, because this transformation does not follow the uncertainty principal, which enables acquiring a high resolution in both time domain and frequency domain. In this paper, we propose an improved method based on the combination of HHT and kurtosis characteristic function. First, Empirical Mode Decomposition (EMD) was applied to the received wave to get a series of intrinsic mode functions (IMFs). Next, Hilbert Transform (HT) was utilized on each of the IMFs to get an instantaneous frequency spectrum. A rough onset time was then determined based on the acquired frequency spectrum. Finally, a moving extraction window was used to calculate a series of fourth-order kurtosis value, and the minimum point corresponded to the calculated onset time. To verify the accuracy and stability of this method, it was applied to a dataset measured by a three-dimensional ultrasonic anemometer. The experimental results showed satisfying advantages compared to other three picking algorithms, including AIC, TKEO, and STA/LTA.

2. Materials and Methods

In this section, the theory of the proposed method is introduced, which is based on the combination of HHT and kurtosis characteristic function. Then, we briefly explain how the other three comparative methods are carried out, including AIC, TKEO, and STA/LTA.

2.1. Overview of HHT-Kurtosis Method

As suggested in the method name, HHT-kurtosis, the proposed method includes two steps, Hilbert-Huang Transform and kurtosis characteristic function calculation. As shown in Figure 2, the first step is to perform Hilbert-Huang Transform on the wave to determine a rough arrival time, around which a fixed size of data is extracted as the initial sample. Then in the second step, the fourth-order kurtosis of a smaller sample extracted successively by a moving window from the initial sample is calculated. The minimum point corresponds to the onset of the ultrasonic echo signal.

2.1.1. Hilbert-Huang Transform

Hilbert Transform, though effective in preserving time localities, requires input to be linear and stable [23]. Huang et al. proposed an improved method, Hilbert-Huang Transform, combining Empirical Mode Decomposition (EMD) and Hilbert Transform. EMD preprocessing makes input signals to be linear and stable, and therefore meets the input requirement of Hilbert Transform [24]:

Hilbert-Huang Transform can be divided into three steps, as is shown in Figure 3.

Step 1: perform EMD to the original ultrasonic echo signal, and get a series of IMFs and a residual.

Step 2: perform Hilbert Transform (HT) to each of IMFs to get the transient frequency spectrum.

Step 3: determine a rough arrival time by setting a threshold to the spectrum.

2.1.2. Kurtosis Characteristic Function Calculation

In order to pick precise onset time of ultrasonic further, higher order statistics method, which is an effective approach for processing non-Gaussian signals, is utilized to deal with this task. In this paper, kurtosis that is fourth-order cumulants was chosen as characteristic function because it is insensitive to the Gaussian noise and can achieve a higher accuracy with a smaller order and less calculation time [25,26].

Kurtosis, defined as the standardized fourth moment about the mean, is a statistical value describing how a given set of data is distributed [27]. For a dataset of N samples, represented as $x=\left\{x_1,x_2,\cdots ,x_N\right\}$, with a mean of $\overline{x}$, the discrete form of kurtosis can be calculated as shown in Equation (2),

$$K=\frac{E\left[{\left(X-\mu \right)}^4\right]}{E{\left[{\left(X-\mu \right)}^2\right]}^2}=\frac{\frac{1}{N}{{\displaystyle \sum }}_{i=1}^{N+1}{\left(x_i-\overline{x}\right)}^4}{{\left[\frac{1}{N}{{\displaystyle \sum }}_{i=1}^{N+1}{\left(x_i-\overline{x}\right)}^2\right]}^2}.$$

(2)

At the time of wave arrival, the kurtosis will increase significantly from a relatively low value [28]. Thus, the minimum kurtosis corresponds to the calculated onset time.

The steps of calculating fourth-order kurtosis characteristic function are shown in Figure 4a:

Step 1: Extract a dataset of n values, represented as $x=\left\{x_1,x_2,\cdots ,x_n\right\}$, around the rough arrival time determined by Hilbert-Huang Transform.

Step 2: Slide a moving extraction window of N samples along x-axis, and record its first-order kurtosis values as $K_1\left(i\right), i\in \left[1,n-N+1\right]$.

Step 3: Apply a succession of transformation to the first-order kurtosis characteristic function as shown in steps (4)–(6).

Step 4: Perform the first transformation and calculate the second-order kurtosis characteristic function as shown in Equation (3),

$$\{\begin{array}{c}{F}_2\left(1\right)=F_1\left(1\right)\\ F_2\left(i+1\right)=F_2\left(i\right)+\delta \left(i\right)\times \left[F_1\left(i+1\right)-F_1\left(i\right)\right]\end{array},$$

(3)

where

$$\delta \left(i\right)=\{\begin{array}{c}1, F_1\left(\mathrm{i}+1\right)-F_1\left(\mathrm{i}\right)\ge 0\\ 0, F_1\left(\mathrm{i}+1\right)-F_1\left(\mathrm{i}\right)<0\end{array},$$

(4)

Step 5: Perform the second transformation and calculate the third-order kurtosis characteristic function as shown in Equation (5),

$$F_3\left(i\right)=F_2\left(i\right)-[\frac{F_2\left(N\right)-F_2\left(1\right)}{N}\times \left(i-1\right)+F_2\left(1\right)]$$

(5)

Step 6: Perform the final transformation and calculate the fourth-order kurtosis characteristic function as shown in Equation (6),

$$F_4\left(i\right)=F_3\left(i\right)-M_i$$

(6)

where $M_i$ is the local maximum of $F_3$.

As shown in Figure 4b, four curves record $F_1, F_2, F_3, \mathrm{and} F_4$ of the extracted initial sample, and the greatest minimum of $F_4$ corresponds to the arrival time of the wave.

2.2. Introduction of Other Comparative Methods

The AIC method is applied by calculating the local minimum of AIC function,

$$\mathrm{AIC}\left(\mathrm{k}\right)=klg\left(Var\left(x\left[1,k\right]\right)\right)+\left(l-k-1\right)lg\left(Var\left(x\left[k+1,l\right]\right)\right),$$

(7)

In this equation, the value k that can make AIC(k) reach the maximum is the calculated onset time [29].

The STA/LTA method is based on an empirically determined threshold. At onset point, the change of STA (average power in a short moving window) is much greater than LTA (average power in a long moving window), so the value n that makes $\frac{\mathrm{STA}}{\mathrm{LTA}}\ge \mathrm{threshold}$ is the calculated onset time [30].

To carry out TKEO method we first need to perform EMD and then calculate the maximum value of TKEO operator,

$$\psi \left[y\left(n\right)\right]=y{\left[n-1\right]}^2-y\left[n\right]\times y\left[n-2\right],$$

(8)

In this equation, the value n that can make $\psi \left[y\left(n\right)\right]$ maximum is the calculated onset time [31,32].

3. Experiments and Results

In this section, how actual experiments are carried out in a wind tunnel is explained. Then we perform four methods mentioned in the last section to the received ultrasonic signal and determine the onset time. Results show the great performance of the proposed method.

3.1. Experimental Device

To verify the performance of the proposed method further, actual experiments are carried out in a wind tunnel. All experiments are based on an experimental ultrasonic anemometer consisting of three pairs of ultrasonic sensors, which excites ultrasound and receives it from the corresponding sensor, whose structure is shown in Figure 5a. All ultrasonic sensors were produced by Airmar with a center frequency of 200 kHz. For gaining sufficient ultrasonic signals, experiments were performed under diverse wind velocities in the Low-Speed Wind Tunnel at the China Aerodynamics Research and Development Center, as is shown in Figure 5b. Consequently, a total of 300 sets of ultrasonic signals were sampled at a sampling rate of 1 MHz.

3.2. Accuracy Test under Different Signal-to-Noise Ratios

To test the reliability of the proposed the method, we conducted a series of experiments under different signal-to-noise ratios (SNRs). These experiments utilized received waves obtained by fixing the distance of one pair of transducers, whose standard arrival time was measured in the noiseless environment. As mentioned above, there are 300 sets of signals acquired from the proposed device. And these signals are added White Gaussian Noise (WGN) with different SNRs for testing the noise robustness of the proposed method. In order to show the satisfactory performances of our method, STA/LTA algorithm, AIC algorithm, and TKEO algorithm were performed under the same conditions to compare with the proposed method, respectively. The results are illustrated in Figure 6, Figure 7 and Figure 8, in which the green dots represent the distributions between measured arrival time by different algorithms and corresponding actual arrival time, and the orange lines termed trend line are a set of dots where the measured value equals the actual value.

As is shown in Figure 6, Figure 7 and Figure 8, under different SNRs, the measured results generated from the proposed method are close to the trend line, which denotes the deviations between the actual arrival time and the measured arrival time were quite small. The average deviation is within $5 \mathsf{\mu }\mathrm{s}$. The experimental results demonstrated the noise robustness of the proposed method. In contrast, the results from STA/LTA algorithm deviated from the trend line overall, and this method could generate erroneous measurement results possibly; AIC algorithm has a tendency to engender false pick with a fixed deviation; and the measured results from TKEO algorithm deviated from the trend line out of the range of $\pm 10 \mathsf{\mu }\mathrm{s}$, which is much worse than ours.

Table 1. Mean Square Error and R-square determined by the four methods.

	MSE(s2)			R-Square
	15 dB	10 dB	5 dB	15 dB	10 dB	5 dB
STA/LTA	1.42 × 10−9	1.99 × 10−9	2.53 × 10−9	0.7191	0.4680	0.2808
AIC	2.78 × 10−9	3.00 × 10−9	3.08 × 10−9	−1.0122	−1.0825	−0.7674
TKEO	2.36 × 10−10	2.45 × 10−10	2.85 × 10−10	0.9372	0.9279	0.9180
Our method	4.30 × 10−11	4.99 × 10−11	7.42 × 10−11	0.9871	0.9849	0.9751

indicates the Mean Square Error (MSE) and R-square determined by the four methods, which also demonstrates the great performance of the proposed method in determining arrival time of ultrasonic signals.

4. Discussion

As shown in Figure 6c, Figure 7c, and Figure 8c, the measured results obtained by AIC approach directly deviated from the trend line with a fixed distance obviously. After checking erroneous results, we find that there are two minima in the AIC function. The sample arrival point was considered to divide the signal into two different stationary processes. Similarly, the ending point of ultrasonic echo signal played the same role in dividing the signal into two different stationary processes. Therefore, the AIC function was likely to exist two minima corresponding to the onset point and end point of ultrasonic signals, respectively. If the minimum result was determined by the ending point, the method was more likely to generate false pick with fixed deviation.

Therefore, in the proposed method, we utilized HHT to determine an approximate range for kurtosis algorithm, which can avoid the fixed deviation, as is shown above. The experimental results demonstrate the great noise robustness and performances in arrival time determination of ultrasonic signals via the proposed method.

5. Conclusions

The purpose of this study is to develop an effective method for detecting the arrival time of ultrasonic signal, which is of vital importance during the wind speed determination. The main improvement regarding this study relates to the combination of Hilbert-Huang Transform and kurtosis characteristic function calculation. HHT was applied to the received wave to determine a rough arrival time. Then, the fourth-order kurtosis was calculated, and its minimum point corresponded to the precise arrival time of the ultrasonic echo signal.

The results obtained through this hybrid method were satisfying. To verify the advantages of our proposed method, a comparison among the results acquired by using the proposed method, AIC, STA/LTA, and TKEO was carried out. It turns out that, results acquired through the proposed method indicate the high determination accuracy of this method, under both the ideal case and real cases.