Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Coherent Doppler wind lidar signal denoising adopting variational mode decomposition based on honey badger algorithm

Open Access Open Access

Abstract

Coherent Doppler wind lidar (CDWL) is used to measure wind velocity distribution by using laser pulses. However, the echo signal is easily affected by atmospheric turbulence, which could decrease the effective detection range of CDWL. In this paper, a variation modal decomposition based on honey badger algorithm (VMD-HBA) is proposed and demonstrated. Compared with conventional VMD-based methods, the proposed method utilizes a newly developed HBA to obtain the optimal VMD parameters by iterating the spectrum fitness function. In addition, the Correlation Euclidean distance is applied to identify the relevant mode and used to reconstruct the signal. The simulation results show that the denoising performance of VMD-HBA is superior to other available denoising methods. Experimentally, this combined method was successfully realized to process the actual lidar echo signal. Under harsh detection conditions, the effective detection range of the homemade CDWL system is extended from 13.41 km to 20.61 km.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Coherent Doppler wind lidar (CDWL) is a lidar which uses Doppler effect to obtain target motion information. With the strong clear sky detection ability, it has been used in meteorological research [13], aviation safety [4], and other fields. Nevertheless, the intensity of lidar echo signal decays rapidly with the extension of detection range according to the lidar equation [5]. In addition, there are background light interference, random noise, and other factors in the detection process, especially in harsh detection conditions, which would drown the weak echo signal in noise [6,7]. Therefore, it is important to expand the detection range of CDWL by seeking a denoising algorithm to extract the backscattered signal from noise effectively.

Noise in the lidar echo signal is a typical non-stationary signal which has the characteristics of continuous amplitude and random phase in the time domain. With the continuous development of the signal processing field, such as wavelet transform (WT), empirical mode decomposition (EMD), and other denoising methods have emerged for the processing of such signals. R. Hussein et al. used a WT method based on histogram threshold function to denoise the partial discharged signal. Compared with ordinary wavelet transform methods, it has unique advantages in signal-to-noise ratio (SNR) and noise elimination. [8]. However, WT cannot adaptively find the optimal combination for different problems. In contrast, EMD technology makes up for the defects of WT [9]. Based on the traditional EMD, Su et al. proposed a new frequency conversion resolution decomposition method (VFEMD), which extended the original EMD to the high-resolution scale and overcame the limitation that the resolution only depends on the length of the decomposed signal [10]. Nonetheless, EMD and its variants also have some defects, such as mode aliasing, over decomposition, and end effects.

Variational modal decomposition (VMD), proposed by Dragomiretskiy and Zosso in 2014, is a new time-frequency analysis method [11]. It can decompose multi-component signals into multiple single-component AM and FM signals at one time, therefore can deal with nonlinear and non-stationary signals effectively. Using its unique Wiener filtering characteristics, VMD can validly make up for the problems of mode aliasing and the boundary effect of EMD. It shows the advantages of better complex data decomposition accuracy and better anti-noise interference. With its excellent characteristics, VMD is already used in fault diagnosis, biomedicine, signal denoising, and other fields [1215]. However, two crucial parameters need to be manually selected: the decomposition mode number K and the quadratic penalty $\alpha$, as well as their selection could seriously affect the final decomposition effect. Moreover, this parameter combination can only get a set of relatively optimal results through many artificial attempts, thus this process will inevitably consume a lot of time and computer memory. This problem is not solved until some novel methods were proposed in recent years. In order to determine related parameters of VMD adaptively, many researchers combined VMD with swarm intelligence optimization algorithms, such as particle swarm optimization algorithm (PSO), whale optimization algorithm (WOA), grasshopper optimization algorithm (GOA), etc. [1618]. In this way, the optimal solution of parameter combination can be found adaptively, meanwhile a better denoising effect can be obtained.

In this paper, we proposed an algorithm of VMD combined with HBA and used to expand the detection range of CDWL. Firstly, the proposed method utilizes HBA to optimize the spectrum fitness function of signal in each range gate. Then the signal is denoised by VMD with the optimal parameters $[K,\;\alpha ]$. Finally, the signal is reconstructed by using the relevant mode, which is selected by calculating the Correlation Euclidean distance between each mode and the input signal. The effectiveness of our method has been demonstrated through multiple verifications by simulation and experiment.

2. Principle of proposed VMD-HBA denoising method

2.1. Brief description of VMD

Based on classical Wiener filtering, the VMD algorithm decomposes an input real signal into a series of band-limited intrinsic mode functions (BLIMFs), and each mode has a center frequency [11]. In order to evaluate the bandwidth of each BLIMF to minimize it, and the sum of all decomposed modes needs to be equal to the original signal, the following constrained variational problem is constructed.

$$\mathop {\min }\limits_{_{\{{{u_k}} \},\{{{\omega_k}} \}}} \left\{ {\sum\limits_{k = 1}^K {||{{\partial_t}[{({\delta (t) + {j / {\pi t}}} )\ast {u_k}(t)} ]{e^{ - j{\omega_k}t}}} ||_2^2} } \right\}$$
$$\sum\limits_{k = 1}^K {{u_k}} = f$$
where $\{{{u_k}} \}= \{{{u_1},{u_2},\;\ldots \;,{u_K}} \}$ is the kth mode, $\{{{\omega_k}} \}= \{{{\omega_1},{\omega_2},\;\ldots \;,{\omega_K}} \}$ is the center frequency of the each mode, K is the decomposition mode number, * is convolution, $({\delta (t) + {j / {\pi t}}} )\ast {u_k}(t)$ represents the analytical signal obtained by Hilbert transform, $||{{\partial_t}[{({\delta (t) + {j / {\pi t}}} )\ast {u_k}(t)} ]{e^{ - j{\omega_k}t}}} ||_2^2$ represents the bandwidth of the kth mode, and Eq. (2) represents the constraints of Eq. (1).

To solve the constrained variational problem, the quadratic penalty and Lagrange multiplier are introduced. The quadratic penalty can ensure the reconstruction accuracy of signal, and the Lagrange multiplier can increase the strictness of constraints. In this way, the above constrained variational problem can be transformed into an unconstrained variational problem. The augmented Lagrange expression is as follows:

$$\begin{array}{l} L({\{{{u_k}} \},\{{{\omega_k}} \},\lambda } )= \alpha \sum\limits_{k = 1}^K {\left\|{{\partial_t}[{({\delta (t) + {j / {\pi t}}} )\ast {u_k}(t)} ]{e^{ - j{\omega_k}t}}} \right\|_2^2} \\ + \left\|{f(t) - \sum\limits_{k = 1}^K {{u_k}(t)} } \right\|_2^2 + \left\langle {\lambda (t),f(t) - \sum\limits_{k = 1}^K {{u_k}(t)} } \right\rangle \end{array}$$
where $\alpha$ is the quadratic penalty, and $\lambda$ is the Lagrange multiplier.

The alternating direction multiplier method (ADMM) [19] is used to continuously iterate $u_k^{n + 1}, \omega _k^{n + 1}, \lambda _k^{n + 1}$ to find the saddle point of the augmented Lagrange expression as an optimal solution of the above constrained variational problem. Applying the principle of Parseval/ Plancherel Fourier isometric transform in the ADMM optimization concept for VMD, the minimum solution problem is transformed to the spectral domain, and the following results can be obtained:

$$\hat{u}_k^{n + 1}(\omega ) = \frac{{\hat{f}(\omega ) - \sum\limits_{i \ne k} {{{\hat{u}}_i}} (\omega ) + \frac{{\hat{\lambda }(\omega )}}{2}}}{{1 + 2\alpha {{({\omega - {\omega_k}} )}^2}}}$$
$$\omega _k^{n + 1} = \frac{{\int_0^{\textrm{ + }\infty } \omega {{|{{{\hat{u}}_k}(\omega )} |}^2}d\omega }}{{\int_0^{\textrm{ + }\infty } {{{|{{{\hat{u}}_k}(\omega )} |}^2}} d\omega }}$$
$${\hat{\lambda }^{n + 1}}(\omega ) = {\hat{\lambda }^n}(\omega ) + \tau \left( {\hat{f}(\omega ) - \sum\limits_{k = 1}^K {\hat{u}_k^{n + 1}(\omega )} } \right)$$

The algorithm iterates continuously according to Eq. (4)–(6) until the convergence condition or the maximum number of iterations are met, and jumps out of the iterative process. The convergence conditions shall meet

$$\sum\limits_{k = 1}^K {\left\|{\hat{u}_k^{n + 1} - \hat{u}_k^n} \right\|_2^2} /\left\|{\hat{u}_k^n} \right\|_2^2 < \varepsilon$$
where ${\hat{u}_k}$ is the spectrum of each mode, ${\omega _k}$ is the center frequency of each mode, $\hat{f}(\omega )$ is the spectrum of the signal, $\tau$ is the noise tolerance, $\varepsilon$ is the convergence criterion tolerance, generally ${10^{ - 6}}$, and n is the number of iterations.

The resulting $\hat{u}_k^{n + 1}$ separates the signal from the noise. According to the decomposition characteristics of VMD, noise generally exists in higher-order BLIMFs. Therefore, it is necessary to select the relevant modes for signal reconstruction by appropriate method.

2.2. Principle of HBA combined with VMD

Honey badger algorithm is a new meta-heuristic algorithm for solving optimization problems first proposed by Fatma A. Hashim et al. in 2021 [20]. It features a simple structure, fewer adjustable parameters, faster convergence speed, stronger global search capability, and a strong ability to jump out of local optimization. The process of using HBA to find the optimal parameter combination $[K,\;\alpha ]$ of VMD is as follows:

  • (i) Initialization phase

    Initialize the population size and individual position of honey badger according to the following formula:

    $${x_i} = {l_i} + {r_1} \times ({{u_i} - {l_i}} )$$
    where ${r_1}$ is a random number within $({0,1} )$, ${x_i}$ is the ith individual position of N candidate individuals, ${u_i}$ and ${l_i}$ are respectively upper and lower bounds of the search domain. After combining with VMD, ${x_i}$ is a two-dimensional matrix, and the two dimensions correspond to K and $\alpha$ of VMD, respectively.

    To optimize the $[K,\;\alpha ]$ of VMD through HBA, it is necessary to set up an objective function. The principle of the function will be described in Section 4.1 in conjunction with the principle of CDWL.

    The intensity is related to concentration strength of the target and the distance between it and ith honey badger. If the intensity is high, the movement speed is fast, and vice versa.

    $${I_i} = {r_2} \times \frac{S}{{4\pi d_i^2}},\;S = {({{x_i} - {x_{i + 1}}} )^2},\;{d_i} = {x_{target}} - {x_i}$$
    where ${I_i}$ represents smell intensity of target, S represents the concentrated strength, ${d_i}$ denotes the distance between the target and the ith honey badger individual, ${x_{target}}$ is the target optimal position so far. In other words, ${x_{target}}$ is a set of K and $\alpha$ combinations of VMD that can minimize the objective function.

    The density factor $\delta$ controls the randomness of time variation to ensure a smooth transition between the digging phases and the honey phase. It decreases as the number of iterations increases, and the update process is as follows:

    $$\delta = C \times \exp \left( {\frac{{ - n}}{{{n_{\max }}}}} \right)$$
    where ${n_{\max }}$ is the maximum number of iterations, C is a constant $\mathrm{\ \ge }1$ (default = 2).

  • (ii) Digging phase

    In the digging phase, a honey badger performs a route similar to Cardioid shape, which can be simulated as follows:

    $${x_{new}} = {x_{target}} \times ({1 + F \times \gamma \times {I_i}} )+ F \times {r_3} \times \delta \times {d_i} \times |{\cos ({2\pi {r_4}} )\times [{1 - \cos ({2\pi {r_5}} )} ]} |$$
    where ${x_{new}}$ refer to the updated location of honey badger, $\gamma \ge 1$ (default = 6) is the ability of honey badger to obtain food, ${r_3}$, ${r_4}$ and ${r_5}$ are three different random numbers between 0 and 1, F is the flag to change the search direction, which is used to jump out of the local optimum. It can be expressed as:
    $$F = \left\{ {\begin{array}{cc} 1&{\textrm{ }{r_6} \le 0.5}\\ { - 1}&{\textrm{ }{r_6}\mathrm{\ > 0}\textrm{.5}\textrm{ }} \end{array}} \right.\;\;\;\;{r_6}\;\textrm{is}\;\textrm{a}\;\textrm{random}\;\textrm{number}\;\textrm{between}\;0\;\textrm{and}\;1$$

  • (iii) Honey phase

    The situation that a honey badger follows the honey guide bird to reach the beehive can be simulated as follows:

    $${x_{new}} = {x_{target}} + F \times {r_7} \times \delta \times {d_i}\;\;\;\;\;\;{r_7}\;\textrm{is}\;\textrm{a}\;\textrm{random}\;\textrm{number}\;\textrm{between}\;0\;\textrm{and}\;1$$

    Iterate the above process repeatedly until $n = {n_{\max }}$ is satisfied. In the iterative process, the optimal target is updated each time according to the fitness function, and the boundary judgment is performed. The final optimal combination is the VMD optimal $[K,\;\alpha ]$ found adaptively through the HBA.

2.3. Identification of the relevant mode with Correlation Euclidean distance

The Correlation Euclidean distance $E({f,{u_i}} )$ considers both correlation coefficient $R({f,{u_i}} )$ and Euclidean distance ${D_o}({f,{u_i}} )$. After testing different weights, they can best reflect the similarity between the two when their weights are 0.5 respectively. Some transformations were performed to make the two trends consistent with the order of magnitude as:

$$\begin{array}{l} E({f,{u_i}} )= 0.5 \cdot R({f,{u_i}} )+ 0.5 \cdot ({10/{D_o}({f,{u_i}} )} )\\ \;\;\;\;\;\;\;\;\;\;\;\; = \textrm{0}\textrm{.5} \cdot \frac{{\sum\limits_{n = 1}^N {({f(n )- \bar{f}} )} ({{u_i}(n )- {{\overline u }_i}} )}}{{\sqrt {\sum\limits_{k = 1}^N {{{({f(n )- \bar{f}} )}^2}} \sum\limits_{k = 1}^N {{{({{u_i}(n )- {{\overline u }_i}} )}^2}} } }} + 5/\sqrt {\sum\limits_{n = 1}^N {{{({f(n )- {u_i}(n )} )}^2}} } \end{array}$$
where f is the original time domain signal, $\bar{f}$ is the mean value of signal, N is the number of signal points, ${u_i}$ is the ith decomposed mode and ${\bar{u}_i}$ is the mean value of decomposed modes.

If the maximum value of $E({f,{u_i}} )$ appears at $BLIM{F_n}$, the reconstructed signal $\tilde{f}(t)$ can be expressed as:

$$\tilde{f}(t) = BLIM{F_n}$$

2.4. Procedure of proposed method

The flowchart of the VMD-HBA algorithm is shown in Fig. 1. The detailed steps of this method can be described as:

  • 1: Signal slice: the echo signal $f(t )$ is divided into different range gate signals ${f_k}(t )$ according to the range resolution, and $k = 1,2,\;\ldots \;,{N_{gate}}$, ${N_{gate}}$ is the number of range gates;
  • 2: Algorithm initialization: The range of the VMD parameters to be optimized is set. Here, the range of the mode number K is taken as $[2,\;15]$, the quadratic penalty factor $\alpha$ is assigned in the range $[1000,\;10000]$. Besides, the HBA model is initialized, including the maximum iteration number and the number of search agents. Here, the maximum number of iterations ${t_{\max }}$ is 30, and the initial population number N is 50;
  • 3: Core steps of the algorithm: the input signal ${f_k}(t )$ is decomposed by using VMD with different combinations of $[K,\;\alpha ]$, and the spectrum fitness function for each decomposition is calculated. The minimum value of the function and the corresponding $[K,\;\alpha ]$ are saved. Then apply optimal $[K,\;\alpha ]$ and Eq. (11) or Eq. (13) to update other combinations. Iterates continuously until the exit condition is met, and save the final optimal $[K,\;\alpha ]$;
  • 4: Signal decomposition: signal ${f_k}(t )$ decomposed by VMD with the optimal $[K,\;\alpha ]$;
  • 5: Signal reconstruction: the Correlation Euclidean distance between each decomposed mode and the original signal is calculated, meanwhile, the relevant mode is identified for signal reconstruction;
  • 6: Signal splice: the denoised signals in all range gates are spliced to obtain the complete denoised signal $F(t )$.

 figure: Fig. 1.

Fig. 1. Flowchart of the proposed method (VMD-HBA)

Download Full Size | PDF

3. Results of lidar signal denoising simulation

3.1. Evaluation index

In order to compare the performance of the proposed method with other denoising methods, the output spectrum detectability $(S{D_{out}})$ is adopted and defined as follows. In addition, the definition of parameters in the formula are consistent with those in the spectrum fitness function in Section 4.1.

$$\left\{ \begin{array}{l} S{D_{out}} = 10\log 10\left( {\frac{{\max ({Pf} )}}{{\sqrt {{{{{\sum\limits_{k \ne (m,n)} {({Pf(k) - mean} )} }^2}} / {({{N_{Pf}} - {N_{peak}}} )}}} }}} \right)({\textrm{dB}} )\\ mean = \frac{1}{{{N_{Pf}} - {N_{peak}}}}\sum\limits_{k \ne (m,n)} {Pf(k)} \end{array} \right.$$
$S{D_{out}}$ is utilized to evaluate the spectrum denoising performance of various algorithms. It is similar to the definition of output peak signal-to-noise ratio $(PSN{R_{out}})$ of the time-domain signal [21]. The variation trend of the $S{D_{out}}$ with the input SNR is shown in Fig. 2.

 figure: Fig. 2.

Fig. 2. Variation of the output SD with the input SNR

Download Full Size | PDF

3.2. Simulation under the constant wind speed in a single range gate

The covariance matrix based lidar signal simulation model is applied to simulate an ideal lidar echo signal [22,23], which is a signal in a single range gate. The synthetic signal consists of an ideal echo signal supplemented with 3 dB white Gaussian noise. Besides, the constant wind speed is set to 40 m/s.

In order to prove the superiority of the proposed method, we compared its performance with different denoising methods, namely, NeighCoeff-db6 WT (WT-db6), EMD combined with detrended fluctuation analysis (EMD-DFA), VMD method based on the NeighCoeff-db6 WT (WT-VMD), and VMD method based on the decomposition level of EMD (EMD-VMD).

By processing the noisy synthetic signal, the simulation results are compared in Fig. 3. Among these methods, VMD-HBA shows the best denoising performance and the highest $S{D_{out}}$. In addition, the spectrum signal peak is narrow and high, which is more convenient to identify in the process of wind speed inversion. Therefore, it can be concluded that our method has advantages in processing the lidar echo signal of a single range gate.

 figure: Fig. 3.

Fig. 3. Comparison of signal spectrum after denoising by different methods

Download Full Size | PDF

The denoising performance by the abovementioned methods with the input SNRs varying from – 9 dB to 6 dB is shown in Fig. 4. The VMD-HBA method obtained the best results, as expected, under different input SNRs. Even in the case of input SNR of -9 dB, VMD-HBA still guaranteed a high $S{D_{out}}$ up to 23.42 dB. The superiority and reliability of the proposed method were further proved by these simulations. But in a real detection situation, the echo signal is composed of multiple range gates. In order to further validate our method, the following is necessary to simulate the detection process of multiple range gates.

 figure: Fig. 4.

Fig. 4. Denoising performance with different input SNRs.

Download Full Size | PDF

3.3. Simulation under a sinusoidal wind speed in multiple range gates

In this section, the covariance matrix based lidar signal simulation model in different range gates is combined into a long-distance detection model to simulate lidar detection in real situations. Firstly, a total of 23 range gates are adopted. The input SNR decreases from 7 dB to – 15 dB, and each range gate decreases by 1 dB in turn. Secondly, one cycle sinusoidal wind speed model is selected as the standard wind speed of 23 range gates. It is used to verify the wind speed of the signal after denoising by VMD-HBA. Finally, 2048 points are utilized for Fourier transform to minimize the wind speed inversion error caused by spectrum resolution. Additionally, the number of sampling points in each gate is 512, corresponding to 153.6 m range resolution.

The variation of wind speed of the original signal, noisy signal, and denoised signal with detection range is displayed in Fig. 5. As shown, the wind speed of the noisy signal begins to deviate at the 14th gate. This is because the noise drowns the signal peak, resulting in errors in inversing the wind speed by using the peak recognition algorithm. The wind speed denoised by VMD-HBA shows superior denoising performance at gates 14-22, although the noise cannot be removed successfully at the 23th gate with a very low SNR $(\textrm{SNR} ={-} 15\;\textrm{dB})$.

 figure: Fig. 5.

Fig. 5. (a) Sinusoidal wind speed. (b) Wind speed of the noise signal. (c) Wind speed of the signal after denoising by VMD-HBA.

Download Full Size | PDF

The standard deviation of wind speed is defined here to measure the accuracy of wind speed after denoising, as shown in the following equation:

$${V_{sd}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{({V_{di}} - {V_{oi}})}^2}} }}{n}}$$
where ${V_d}$ is the denoised wind speed, ${V_o}$ is the standard wind speed, and n is the number of range gates corresponding to the effective detection range.

According to Eq. (17), the standard deviation of denoised wind speed is 0.38 m/s, and the detection range is increased from 2.15 km to 3.38 km, with an increase of 57%.

The spectrum of noisy signal and denoised signal in two cases are compared in Fig. 6. At the 4th gate with higher SNR $(\textrm{SNR} = 4\;\textrm{dB})$, the denoised peak is smooth and sharp as well as the surrounding noise suppression effect is obvious. At the 22th gate with extremely low SNR $(\textrm{SNR} ={-} 14\;\textrm{dB})$, the correct signal peak has been submerged by noise. As expected, the signal spectrum denoised by VMD-HBA recovers the signal peak successfully. Although partial noise on both sides of the peak cannot be removed, it does not affect the peak identification in wind speed inversion. Therefore, our method can expand the detection range of lidar effectively. But the noise in the real situation is complex, so the following is necessary for the processing of lidar echo signal in the actual environment.

 figure: Fig. 6.

Fig. 6. Spectrum comparison of noisy signal and denoised signal: (a) Noisy signal in the 4th range gate. (b) Denoised signal in the 4th range gate. (c) Noisy signal in the 22th range gate. (d) Denoised signal in the 22th range gate.

Download Full Size | PDF

4. Results of real lidar signal denoising experiment

4.1. CDWL system and the corresponding spectrum fitness function

In order to illustrate the effectivity of VMD-HBA, the real lidar echo signals are processed. The experimental data were collected by the CDWL system developed by our group. The schematic diagram of lidar system is shown in Fig. 7, which is mainly composed of seed laser, slave laser, optical transceiver system, and data processing system. The CDWL delivers a 1645 nm single frequency laser with a pulse energy of 8mJ and a pulse repetition frequency of 300 Hz. The accumulated number is 150 and the beam elevation angle is $5\mathrm{^\circ }$.

 figure: Fig. 7.

Fig. 7. The schematic diagram of the CDWL system. DAB, digital acquisition board; ADC, analog to digital converter; AOM, acousto-optic modulator; PBS, polarizing beam splitter; QWP, quarter waveplate.

Download Full Size | PDF

The lidar echo signal acquisition experiment was carried out in Beijing $(116.31\mathrm{^\circ }E,\;36.65\mathrm{^\circ }N)$ on November 15, 2021. An acquisition board (Agilent U5303A) with a sampling frequency of 500 MHz was used to collect time-domain data. The spectrums of each pulse are accumulated and averaged, and then inversed back to the time domain and compressed into one signal. The power spectrum of the time-domain signal is obtained by Fourier transform, and the peak recognition algorithm is used to calculate the wind speed.

CDWL utilizes the Doppler shift of the echo signal to measure wind speed. In the range gate with a far detection range, the spectral noise will submerge the main peak of the signal, as shown in Fig. 8(a). This will undoubtedly lead to wrong wind speed inversion results.

 figure: Fig. 8.

Fig. 8. Explanation of the spectrum fitness function: (a) Signal spectrum in range gate with very low SNR; (b) Identification method of signal peak.

Download Full Size | PDF

In order to expand the detection range of CDWL, the spectrum of signal in each range gate should be optimized. Obviously, this can be solved during the iterative optimization of $[K,\;\alpha ]$ of VMD by HBA. Moreover, the objective function of the iteration in HBA should be set as a function related to the echo signal spectrum within each range gate. The corresponding objective function, namely, spectrum fitness function is proposed as follows:

$$\left\{ \begin{array}{l} fitness = 10\log 10\left( {\frac{{\max (Pf)}}{{\frac{1}{{{N_{Pf}} - {N_{peak}}}}\sum\limits_{k \ne (m,n)} {Pf(k)} }}} \right)({\textrm{dB}} )\\ \textrm{ }peak = Pf(m,n) \end{array} \right.$$
where $Pf$ is the spectrum of the echo signal in a range gate, ${N_{Pf}}$ is the number of signal spectrum points, ${N_{peak}}$ is the number of spectrum points corresponding to the signal peak, and $peak$ is the signal peak.

The identification rule of signal peak is: first find the corresponding point of the maximum value of signal spectrum, and continuously find the point with monotonic decreasing amplitude on the left and right until the requirements are not met. It is marked as m on the left and n on the right, as shown in Fig. 8 (b).

This spectrum fitness function represents the quality of the signal spectrum in each range gate. The larger the function value, the better the spectrum, and vice versa. It is used as the objective function of HBA iterative optimization.

4.2. Analysis of denoising results

As we know, the detection range of the CDWL system is greatly affected by atmospheric turbulence, which introduces lots of noises in bad weather. In good weather conditions, our CDWL system measured the wind speed with a detection range of 25 km [1]. While with strong atmospheric turbulence, the detection range is reduced. Now we show the detection range of our system can be expanded by using the algorithm of VMD combined with HBA.

To analyze the performance of the proposed method, the wind detection range with different accumulated numbers (denoted as N) and range resolution (denoted as $\mathrm{\Delta }r$) are compared for both noisy signal and denoised signal. Hypothetically, the wind speed measured under the initial parameters of the system $(N = 150,\;\mathrm{\Delta }r = 90\;\textrm{m})$ is used as the reference wind speed.

Under strong turbulence conditions, the wind detection range of our system before and after denoise are displayed in Fig. 9. It indicates that the effective wind detection range reduces as the accumulated number decreases, while the effective detection range extends evidently after the signal is denoised by the VMD-HBA. It is worth noting that the huge fluctuation of wind speed is empirically considered as the boundary of the wind measurement range. For the noisy signal, the effective detection range with $N = 150,\;60\;\textrm{and}\;30$ were about 13.41 km, 12.78 km, and 11.52 km, respectively, while for the denoised signal, it extends to about 20.61 km, 16.47 km, and 15.30 km.

 figure: Fig. 9.

Fig. 9. Wind detection range of noisy signal and denoised signal under different accumulated numbers: (a) $N = 150$, (b) $N = 60$, (c) $N = 30$.

Download Full Size | PDF

In addition to the influence of accumulated number, range resolution also has a great influence on the detection range of CDWL. The detection range with $N = 30\;\textrm{and}\;\mathrm{\Delta }r = 60\;\textrm{m}$ of noisy signal and denoised signal is displayed in Fig. 10. Under this condition, the effective detection range of the system is only 8.82 km. Nevertheless, the detection range after VMD-HBA denoising can still reach 14.40 km, even exceeding the detection range of the system under the initial parameters.

 figure: Fig. 10.

Fig. 10. Wind detection range of noisy signal and denoised signal under the condition of $N = 30,\;\mathrm{\Delta }r = 60\;\textrm{m}$

Download Full Size | PDF

For the long-range CDWL, it is important to measure the wind velocity with a long detection range and a high accuracy. Therefore, denoising results of three different accumulated numbers at 90 m range resolution are compared with the reference wind speed. It can be seen from Fig. 11 that within the detection range of the reference wind speed, the wind speed difference after denoising increases as the accumulated number decreases. More importantly, the wind speed standard deviation with $N = 30$ is only 0.48 m/s. The standard deviations of wind speed under different conditions calculated by Eq. (17) are summarized in Table 1.

 figure: Fig. 11.

Fig. 11. Difference between the wind speed of denoised signal and reference wind speed at different distances under three accumulated numbers: (a) $N = 150$, (b) $N = 60$, (c) $N = 30$.

Download Full Size | PDF

Tables Icon

Table 1. The wind speed standard deviation of noisy signals and denoised signals

The time resolution of CDWL is determined by its accumulated number. For a CDWL system with a 300 Hz repetition frequency, the time resolution is 0.5 s when 150 pulses are accumulated. According to this paper, for the $N = 30,\;\mathrm{\Delta }r = 60\;\textrm{m}$ condition, the effective detection range after denoising can still reach 14.40 km. Remarkably, it is even higher than 13.41 km under the initial conditions of the system. This proves that our method can improve the time resolution by 5 times, the range resolution by 1.5 times, and obtain a further detection range on the premise of ensuring the accuracy of wind speed.

5. Conclusion

In this paper, an algorithm of VMD combined with HBA is proposed. Firstly, HBA is used to optimize the spectrum fitness function of a single range gate signal, so as to search the optimal parameters $[K,\alpha ]$ of VMD. Secondly, the signal is denoised by VMD with the optimal parameters. Finally, the Correlation Euclidean distance is used to select the relevant mode for signal reconstruction. Wt-db6, EMD-DFA, WT-VMD, EMD-VMD, and VMD-HBA are used for denoising comparison on the same simulation model. It is verified that the signal denoised by VMD-HBA has the highest $S{D_{out}}$ under different input SNRs among these methods. In addition, the feasibility of using VMD-HBA to expand the detection range of lidar is tested under a sinusoidal wind speed simulation detection model. Experimentally, our method is applied to a homemade CDWL system to compare the effective detection range of noisy signal and denoised signal under different conditions. After denoising by VMD-HBA, the detection range with $N = 150\;\textrm{and}\;\mathrm{\Delta }r = 90\;\textrm{m}$ is extended from 13.41 km to 20.61 km. For some harsh conditions, the result means that the detection range increase can reach 53.7%. Moreover, the detection range with $N = 30\;\textrm{and}\;\mathrm{\Delta }r = 90\;\textrm{m}$ is extended from 11.52 km to 15.30 km, and the standard deviation of wind speed is only 0.48 m/s. Notably, for the $N = 30,\;\mathrm{\Delta }r = 60\;\textrm{m}$ condition, the detection range after denoising can still reach 14.40 km. Consequently, our denoising method can extract the signal from the noises, which increases the detection range. On the premise of ensuring the wind velocity measurement accuracy, the time resolution and the range resolution are improved. For 15 km detection, it only needs an accumulated number of 30, also the range resolution can reach 60 m. This means a 5 times improvement in time resolution and a 1.5 times improvement in range resolution.

Funding

National Natural Science Foundation of China (61627821).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. K. Wang, C. Gao, Z. Lin, Q. Wang, and C. Chen, “1645 nm Coherent Doppler wind lidar with a single-frequency Er:YAG laser,” Opt. Express 28(10), 14694–14704 (2020). [CrossRef]  

2. J. Yuan, K. Wu, T. Wei, L. Wang, Z. Shu, Y. Yang, and H. Xia, “Cloud Seeding Evidenced by Coherent Doppler Wind Lidar,” Remote Sens. 13(19), 3815 (2021). [CrossRef]  

3. J. Yuan, H. Xia, T. Wei, L. Wang, and Y. Wu, “Identifying cloud, precipitation, windshear, and turbulence by deep analysis of the power spectrum of coherent Doppler wind lidar,” Opt. Express 28(25), 37406–37418 (2020). [CrossRef]  

4. S. M. Hannon, K. S. Barr, D. K. Jacob, M. W. Phillips, U. N. Singh, and K. Mizutani, “Application of pulsed Doppler lidar in the airport terminal area,” Proc. SPIE 5653, 186–197 (2005). [CrossRef]  

5. C. F. Abari, X. Chu, R. M. Hardesty, and J. Mann, “A reconfigurable all-fiber polarization-diversity coherent Doppler lidar: principles and numerical simulations,” Appl. Opt. 54(30), 8999–9009 (2015). [CrossRef]  

6. J. Su, Y. Wang, and D. Liang, “Long range detection of line-array multi-pulsed coding lidar by combining the Accumulation coherence and Subpixel-energy detection method,” Opt. Express 23(12), 15174–15185 (2015). [CrossRef]  

7. H. T. Fang and D. S. Huang, “Noise reduction in lidar signal based on discrete wavelet transform,” Opt. Commun. 233(1-3), 67–76 (2004). [CrossRef]  

8. R. Hussein, K. B. Shaban, and A. H. El-Hag, “Wavelet Transform With Histogram-Based Threshold Estimation for Online Partial Discharge Signal Denoising,” IEEE Trans Instrum Meas 64(12), 3601–3614 (2015). [CrossRef]  

9. N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proc. Math. Phys. Eng. Sci. 454(1971), 903–995 (1998). [CrossRef]  

10. J. Su, Y. Wang, X. Yang, and X. Wang, “Enhancement of Weak Lidar Signal Based on Variable Frequency Resolution EMD,” IEEE Photonics Technol. Lett. 28(24), 2882–2885 (2016). [CrossRef]  

11. K. Dragomiretskiy and D. Zosso, “Variational mode decomposition,” IEEE Trans. Signal Process. 62(3), 531–544 (2014). [CrossRef]  

12. B. Hwa, W. B. Fei, and Z. Lu, “Application of variational mode decomposition optimized with improved whale optimization algorithm in bearing failure diagnosis,” ALEX ENG J 60(5), 4689–4699 (2021). [CrossRef]  

13. A. Mert, “ECG feature extraction based on the bandwidth properties of variational mode decomposition,” Physiol. Meas. 37(4), 530–543 (2016). [CrossRef]  

14. J. Lu, J. Yue, L. Zhu, D. Wang, and G. Li, “An improved variational mode decomposition method based on the optimization of salp swarm algorithm used for denoising of natural gas pipeline leakage signal,” Measurement 185, 110107 (2021). [CrossRef]  

15. J. Quan and L. Shang, “An Ensemble Model of Wind Speed Forecasting Based on Variational Mode Decomposition and Bare-Bones Fireworks Algorithm,” Math. Probl. Eng. 2021, 1–16 (2021). [CrossRef]  

16. J. Lu, X. Qu, D. Wang, J. Yue, and G. Li, “Signal filtering method of variational mode decomposition and Euclidean distance based on optimizing parameters of classification particle swarm optimization algorithm,” Trans. Inst. Meas. Control. 43(9), 2018–2029 (2021). [CrossRef]  

17. H. Li, J. Chang, X. Fan, Z. Liu, and B. Liu, “Efficient Lidar Signal Denoising Algorithm Using Variational Mode Decomposition Combined with a Whale Optimization Algorithm,” Remote Sens. 11(2), 126 (2019). [CrossRef]  

18. T. Hua, K. Dai, X. Zhang, Z. Yao, H. Wang, K. Xie, T. Feng, and H. Zhang, “Optimal VMD-based signal denoising for laser radar via Hausdorff distance and wavelet transform,” IEEE Access 7, 167997–168010 (2019). [CrossRef]  

19. E. Ghadimi, A. Teixeira, I. Shames, and M. Johansson, “Optimal parameter selection for the alternating direction method of multipliers (ADMM): quadratic problems,” IEEE Trans. Automat. Contr. 60(3), 644–658 (2015). [CrossRef]  

20. F. A. Hashim, E. H. Houssein, K. Hussain, M. S. Mabrouk, and W. Al-Atabany, “Honey Badger Algorithm: New metaheuristic algorithm for solving optimization problems,” Math Comput Simul 192, 84–110 (2022). [CrossRef]  

21. B. Qi, G. Yang, D. Guo, and C. Wang, “EMD and VMD-GWO parallel optimization algorithm to overcome Lidar ranging limitations,” Opt. Express 29(2), 2855–2873 (2021). [CrossRef]  

22. R. Frehlich and M. Yadlowsky, “Performance of mean-frequency estimators for Doppler radar and lidar,” J Atmos Ocean Technol 11(5), 1217–1230 (1994). [CrossRef]  

23. R. Frehlich, “Cramer-Rao bound for Gaussian random processes and applications to radar processing of atmospheric signals,” IEEE Trans Geosci Remote Sens 31(6), 1123–1131 (1993). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1.
Fig. 1. Flowchart of the proposed method (VMD-HBA)
Fig. 2.
Fig. 2. Variation of the output SD with the input SNR
Fig. 3.
Fig. 3. Comparison of signal spectrum after denoising by different methods
Fig. 4.
Fig. 4. Denoising performance with different input SNRs.
Fig. 5.
Fig. 5. (a) Sinusoidal wind speed. (b) Wind speed of the noise signal. (c) Wind speed of the signal after denoising by VMD-HBA.
Fig. 6.
Fig. 6. Spectrum comparison of noisy signal and denoised signal: (a) Noisy signal in the 4th range gate. (b) Denoised signal in the 4th range gate. (c) Noisy signal in the 22th range gate. (d) Denoised signal in the 22th range gate.
Fig. 7.
Fig. 7. The schematic diagram of the CDWL system. DAB, digital acquisition board; ADC, analog to digital converter; AOM, acousto-optic modulator; PBS, polarizing beam splitter; QWP, quarter waveplate.
Fig. 8.
Fig. 8. Explanation of the spectrum fitness function: (a) Signal spectrum in range gate with very low SNR; (b) Identification method of signal peak.
Fig. 9.
Fig. 9. Wind detection range of noisy signal and denoised signal under different accumulated numbers: (a) $N = 150$, (b) $N = 60$, (c) $N = 30$.
Fig. 10.
Fig. 10. Wind detection range of noisy signal and denoised signal under the condition of $N = 30,\;\mathrm{\Delta }r = 60\;\textrm{m}$
Fig. 11.
Fig. 11. Difference between the wind speed of denoised signal and reference wind speed at different distances under three accumulated numbers: (a) $N = 150$, (b) $N = 60$, (c) $N = 30$.

Tables (1)

Tables Icon

Table 1. The wind speed standard deviation of noisy signals and denoised signals

Equations (18)

Equations on this page are rendered with MathJax. Learn more.

min { u k } , { ω k } { k = 1 K | | t [ ( δ ( t ) + j / π t ) u k ( t ) ] e j ω k t | | 2 2 }
k = 1 K u k = f
L ( { u k } , { ω k } , λ ) = α k = 1 K t [ ( δ ( t ) + j / π t ) u k ( t ) ] e j ω k t 2 2 + f ( t ) k = 1 K u k ( t ) 2 2 + λ ( t ) , f ( t ) k = 1 K u k ( t )
u ^ k n + 1 ( ω ) = f ^ ( ω ) i k u ^ i ( ω ) + λ ^ ( ω ) 2 1 + 2 α ( ω ω k ) 2
ω k n + 1 = 0  +  ω | u ^ k ( ω ) | 2 d ω 0  +  | u ^ k ( ω ) | 2 d ω
λ ^ n + 1 ( ω ) = λ ^ n ( ω ) + τ ( f ^ ( ω ) k = 1 K u ^ k n + 1 ( ω ) )
k = 1 K u ^ k n + 1 u ^ k n 2 2 / u ^ k n 2 2 < ε
x i = l i + r 1 × ( u i l i )
I i = r 2 × S 4 π d i 2 , S = ( x i x i + 1 ) 2 , d i = x t a r g e t x i
δ = C × exp ( n n max )
x n e w = x t a r g e t × ( 1 + F × γ × I i ) + F × r 3 × δ × d i × | cos ( 2 π r 4 ) × [ 1 cos ( 2 π r 5 ) ] |
F = { 1   r 6 0.5 1   r 6   > 0 .5   r 6 is a random number between 0 and 1
x n e w = x t a r g e t + F × r 7 × δ × d i r 7 is a random number between 0 and 1
E ( f , u i ) = 0.5 R ( f , u i ) + 0.5 ( 10 / D o ( f , u i ) ) = 0 .5 n = 1 N ( f ( n ) f ¯ ) ( u i ( n ) u ¯ i ) k = 1 N ( f ( n ) f ¯ ) 2 k = 1 N ( u i ( n ) u ¯ i ) 2 + 5 / n = 1 N ( f ( n ) u i ( n ) ) 2
f ~ ( t ) = B L I M F n
{ S D o u t = 10 log 10 ( max ( P f ) k ( m , n ) ( P f ( k ) m e a n ) 2 / ( N P f N p e a k ) ) ( dB ) m e a n = 1 N P f N p e a k k ( m , n ) P f ( k )
V s d = i = 1 n ( V d i V o i ) 2 n
{ f i t n e s s = 10 log 10 ( max ( P f ) 1 N P f N p e a k k ( m , n ) P f ( k ) ) ( dB )   p e a k = P f ( m , n )
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.