A Novel Lidar Signal-Denoising Algorithm Based on Sparrow Search Algorithm for Optimal Variational Modal Decomposition

Abstract

Atmospheric lidar is susceptible to the influence of light attenuation, sky background light, and detector dark currents during the detection process. This results in a large amount of noise in the lidar return signal. To reduce noise and extract a useful signal, a novel denoising method combined with variational modal decomposition (VMD), the sparrow search algorithm (SSA) and singular value decomposition (SVD) is proposed. The SSA is used to optimize the number of decomposition layers K and the quadratic penalty factor α values of the VMD algorithm. Some intrinsic mode function (IMF) components obtained from the VMD-SSA decomposition are grouped and reconstructed according to the interrelationship number selection criterion. Then, the reconstructed signal is further denoised by combining the strong noise-reduction ability of SVD to obtain a clean lidar return signal. To verify the effectiveness of the VMD-SSA-SVD method, the method is compared and analysed with wavelet packet decomposition, empirical modal decomposition (EMD), ensemble empirical modal decomposition (EEMD), and adaptive noise-complete ensemble empirical modal decomposition (CEEMD), and its noise-reduction effect is considerably improved over that of the other four methods. The method can eliminate the complex noise in the lidar return signal while retaining all the details of the signal. The signal is not distorted, the waveform is smoother, and far-field noise interference can be suppressed. The denoised signal is closer to the real signal with higher accuracy, which shows the feasibility and the practicality of the proposed method.

1. Introduction

Lidar is an important active measurement method because of its high temporal and spatial resolution. Its laser light source has several advantages, including good monochromaticity, strong coherence, and high collimation. Moreover, lidar can be used for fine detection of optical and microphysical properties of atmospheric aerosols, as well as temperature, relative humidity, and other parameters. Thus, lidar technology has undergone rapid development and has obtained outstanding results. Therefore, it has been widely used in atmospheric remote sensing and environmental monitoring. However, the return signal from lidar is easily affected by noise produced from the sky background and dark currents of the light detector. With the increase in detection range, the return signal strength becomes increasingly weaker and the far-field signal is easily drowned in the noise [1]. Therefore, it is important to achieve noise removal and extract the effective signal.

The level of signal processing is becoming increasingly mature, and there are more new methods for lidar return signal processing. In the traditional method of signal processing, the Fourier transform is used to extract the target signal by analysing the frequencies of the noise and the useful signal, but it cannot be used for local time-frequency analysis and is only applicable to linear and smooth signals. However, for lidar signals with nonlinear and nonstationary features, the Fourier transform does not achieve the expected results and easily causes distortion of the denoised signal [2].

As a new transform, the wavelet transform is built based on the traditional Fourier transform, which decomposes the signal into different frequency components, localizes the signal time-frequency, and extracts the signal by stretching the translation sequence for denoising. In 2009, Wang et al. applied the wavelet multiresolution thresholding method to the denoising analysis of lidar return signals [3]. In 2011, Mao et al. used the wavelet packet to finely extract the original lidar signal to make up for the shortcomings of the wavelet transform and effectively filter out the noise [4]. In 2016, Qin et al. used an adaptive method combining wavelet analysis and a neural network to improve the signal-to-noise ratio (SNR) of the lidar signal and obtained a better denoising effect [5].

In recent years, the empirical mode decomposition (EMD) method has been rapidly developed in the field of signal denoising, and scientists have performed much research. Huang et al. in 1998 adopted the Hilbert–Huang transform method combining EMD and Hilbert analysis to decompose the signal from small to large into a series of intrinsic mode functions (IMFs) as well as residuals to achieve noise reduction on seismic wave signals [6]. This method not only makes up for the shortcomings of wavelet transforms with good adaptivity but also allows the processing of nonstationary signals. In 2009, Zheng et al. applied the EMD method for filtering acquired lidar signals and produced important experimental results [7]. In 2009, Wu et al. proposed a method that can solve the problem of EMD constrained by both mode mixing and endpoint effect constraints to some extent, namely the ensemble empirical mode decomposition (EEMD) method [8]. Wang et al. applied this EEMD improved method to through-wall radar data denoising [9], and it can be clearly seen that the modal aliasing is improved and the denoising effect is better. In 2012, Li et al. used the EEMD method to process the lidar return signal [10], and the obtained signal had better smoothness and retained the useful signal at low frequencies. In 2020, Cheng et al. proposed a segmented singular value decomposition (SVD) denoising method based on EEMD and boosted wavelet transform [11], which is ideal for noise reduction in lidar return signals.

In 2011, on the basis of EEMD, Torresetal et al. proposed the complete ensemble empirical mode decomposition (CEEMD) method [12], which adds pairs of positive and negative opposite white noise signals to the original signal and then performs EMD decomposition for these two opposite signals. The method can eliminate the modal aliasing problem in EMD and solve the problem of residual auxiliary noise in EEMD reconstruction. In 2021, Ma et al. proposed CEEMD combined with improved wavelet thresholding for lidar signal denoising [13], and the experimental effect was remarkable with improved SNR. Although EMD methods have been continuously improved, they still have some defects, such as easy disturbance by noise in the obtainment of inaccurate IMF components, problems related to modal mixing and endpoint effects, the lack of a strict mathematical definition, and reliance on an empirical foundation.

In 2014, Dragomiretskiy et al. proposed variational modal decomposition (VMD), which has obvious advantages for processing nonlinear and nonstationary signals [14]. In principle, the EMD method is based on time domain recursive decomposition, whereas the VMD method is based on frequency and completely nonrecursive decomposition, so the VMD method can better solve the problems of the EMD method. The VMD method uses an iterative approach to search for the optimal solution of the variational model to determine the bandwidth and centre frequency of the IMF components. The VMD method has good robustness and is capable of adaptively dissecting the signal frequency domain to effectively separate IMFs, and it has been applied in some fields, such as dam signal monitoring and part fault diagnosis [15]. In 2018, Xu et al. proposed an efficient denoising method combining the whale optimization algorithm (WOA) and VMD and successfully applied it to the noise reduction of lidar return signals, producing a substantially better noise-reduction effect than wavelet analysis and other methods [16]. In 2019, Luo et al. proposed a joint noise reduction for ball mill vibration signals by using VMD combined with SVD; this reduction can remove the noise in the mill cylinder signal and improve the reliability of the subsequent analysis and processing of the signal [17].

In this paper, by using the sparrow search algorithm (SSA), a novel VMD-SSA-SVD denoising algorithm is proposed in which the SSA is used to optimize the VMD parameters, and then the SVD algorithm is employed to denoise the initially decomposed signal. To verify the feasibility and practicality, the VMD-SSA-SVD denoising algorithm is applied to the Mie scattering lidar signal and compared with other methods. The experimental results show that it can eliminate complex noise in the lidar return signal while preserving the complete details of the signal.

2. Methodology

2.1. Lidar Equation

By using the interaction between the laser light and atmospheric particles or molecules to produce scattering and attenuation phenomena, the lidar return equation can be determined from the energy of the backscattering signal at the detection distance,

$$P\left(r\right)=P_{0}C{r}^{-2}\beta \left(r\right)\times e^{\left. [-2{{\displaystyle \int }}_0^r\alpha \left(r^{\prime }\right)d{r}^{\prime }\right]},$$

(1)

where $P\left(r\right)$ is the lidar return signal power detected at distance r, $P_0$ is the laser transmit power, and C is the system calibration constant, including the optical loss of the transmitting and receiving systems, the effective receiving area of the receiving system, and the rest of the system constants. $\beta \left(r\right)={\beta }_{\mathrm{M}}\left(r\right)$ + ${\beta }_{\mathrm{A}}\left(r\right)$ represents the backscattering coefficient, which is composed of the atmospheric molecule backscattering coefficient ${\beta }_{\mathrm{M}}\left(r\right)$ and aerosol backscattering coefficient ${\beta }_{\mathrm{A}}\left(r\right)$. The extinction coefficient is denoted $\mathsf{\alpha }\left(r\right)={\alpha }_{\mathrm{M}}\left(r\right)$ + ${\alpha }_{\mathrm{A}}\left(r\right)$ and is composed of the atmospheric molecule extinction coefficient of ${\alpha }_{\mathrm{M}}\left(r\right)$ and aerosol extinction coefficient ${\alpha }_{\mathrm{A}}\left(r\right)$.

$\beta \left(r\right)$ and $\alpha \left(r\right)$ are the two unknown quantities, and usually, two unknown quantities in an equation cannot be solved algebraically. Therefore, common algorithms, such as the Collis slope method [18], Klett’s method [19] and Fernald’s method [20], are used. It is necessary to invert the aerosol extinction coefficient by using the lidar ratio, namely assuming the ratio of the aerosol extinction coefficient to the backscattering coefficient.

2.2. VMD Algorithm

Compared with the EMD-like method, VMD abandons the complex processing of circular screening and uses the advantages of its own Wiener filtering characteristics, such as high decomposition accuracy, fast computation, strong theoretical support, and high noise robustness, to achieve adaptive ground signal decomposition. In addition, it can overcome the defects of EMD-type methods such as the modal aliasing phenomenon and endpoint effect, has strong resistance to modal aliasing, and can obtain an obvious filtering effect.

The VMD algorithm determines the centre frequency and bandwidth of each intrinsic mode component by iteratively searching for the optimal solution of the variational model. By iteratively searching, it can determine the centre frequency and bandwidth of each of the intrinsic mode components and effectively dissect the signal from low to high frequencies. Essentially, it is the process of solving the variational problem [21].

With regard to the construction of the variational problem, to find k modal functions $u_k\left(t\right)$, the sum of the estimated bandwidths of all their modes is minimized, and the constraint is that the input signal f is equal to the sum of the individual modes. The finite bandwidth of each mode is evaluated by the following steps.

The analytic signal of each mode ${\mu }_k$ is calculated by the Hilbert transform to obtain the one-sided spectrum of each mode function

$$\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right),$$

(2)

where $\delta \left(t\right)$ is denoted as the shock function, t is the time, and $\left\{u_k\left(t\right)\right\}=\left\{u_1\left(t\right),u_2\left(t\right),\dots u_k\left(t\right)\right\}$ denotes the k modal components obtained by VMD.

The centre frequency is modulated according to exponential mixing to convert the spectrum of each mode to the corresponding fundamental band

$$\left. [\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right)\right]e^{-j{\omega }_{k}t},$$

(3)

where $\left\{{\omega }_k\right\}=\left\{{\omega }_1,{\omega }_2,\dots {\omega }_k\right\}$ represents the centre frequency of each ${\mu }_k\left(t\right)$.

Gaussian smoothing of the demodulated signal enables us to estimate the bandwidth of each modal component, that is, to calculate the square of the gradient diophantine. From this, a constrained variational model can be constructed and expressed as

$$\begin{gathered} \underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{min}\left\{{\displaystyle \sum }_{k=1}^K{\partial }_t\parallel \left. [\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\parallel }_2^2\right\} \\ s.t.\hspace{1em}\hspace{1em}{{\displaystyle \sum }}_{k=1}^K{u}_k=f. \end{gathered}$$

(4)

With regard to the solution of the variational fractional problem, as a sufficiently large positive number, the quadratic penalty factor can guarantee the reconstruction accuracy of the input signal, and the Lagrangian penalty operator can make the reconstructed accuracy more precise. Therefore, by introducing the quadratic penalty factor $\alpha$ and the Lagrangian penalty operator $\lambda$, the constrained variational problem of Equation (3) can be transformed into an unconstrained variational problem. The combination of the two can be written as the generalized Lagrangian expression as

$$\begin{gathered} L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha {{\displaystyle \sum }}_k\parallel {\partial }_t\left. [\delta \left(t\right)+\frac{j}{\pi t}*u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\parallel }_2^2+\parallel f\left(t\right)-{{\displaystyle \sum }}_i{u}_i\left(t\right){\parallel }_2^2 \\ +\lambda \left(t\right),f\left(t\right)-{{\displaystyle \sum }}_k{u}_k\left(t\right). \end{gathered}$$

(5)

To obtain the “saddle point” of the generalized Lagrangian expression L, iterative updating of ${\mu }_k, {\omega }_k$,$\lambda$, i.e., alternating direction multiplication, is performed to obtain the optimal solution of the constrained variational model under the condition that it is within the required accuracy. To update ${\mu }_k^{n+1}$, the following equations can describe the corresponding minimum problem:

$$u_k^{n+1}=\underset{u_k\in x}{argmin}\left\{\alpha \parallel {\partial }_t\left. [\delta \left(t\right)+\frac{j}{\pi t}*u_k\left(t\right)\right]e^{-j{\omega }_{k}t}{\parallel }_2^2+\parallel f\left(t\right)-{{\displaystyle \sum }}_i{u}_i\left(t\right)+\frac{\lambda \left(t\right)}{2}{\parallel }_2^2\right\}.$$

(6)

Based on the Parseval/Plancherel Fourier isometric transform, which converts Equation (5) to the frequency domain, the above minimum problem can be solved by

$${\stackrel{\wedge }{u}}_k^{n+1}\left(\omega \right)=\frac{f\left(\omega \right)-{{\displaystyle \sum }}_{i\ne k}{\stackrel{\wedge }{u}}_i\left(\omega \right)+\frac{\stackrel{\wedge }{\lambda \left(\omega \right)}}{2}}{1+2\alpha {\left(\omega -{\omega }_k\right)}^2},$$

(7)

where ${\stackrel{\wedge }{u}}_k^{n+1}$ represents the current signal Wiener filtering and $\left\{u_k\left(t\right)\right\}$ is the real part obtained after the Fourier transform of $\left\{{\stackrel{\wedge }{u}}_k\left(\omega \right)\right\}$.

By following the same steps, the centre frequency ${\omega }_k$ is updated as above,

$${\omega }_k^{n+1}=\frac{{{\displaystyle \int }}_0^{\infty }\omega {\left|\stackrel{\wedge }{u_k\left(w\right)}\right|}^{2}d\omega }{{{\displaystyle \int }}_0^{\infty }{\left|\stackrel{\wedge }{u_k\left(w\right)}\right|}^{2}d\omega }.$$

(8)

To reconstruct the signal accurately and efficiently, the Lagrangian multipliers are updated as

$${\stackrel{\wedge }{\lambda }}^{n+1}\left(\omega \right)={\stackrel{\wedge }{\lambda }}^n\left(\omega \right)+\tau \left(\stackrel{\wedge }{f}\left(\omega \right)-{{\displaystyle \sum }}_k{\stackrel{\wedge }{u}}_k^{n+1}\left(\omega \right)\right) .$$

(9)

The whole iterative process is modified by using Equations (7) and (8) to determine whether the constraint is satisfied. If it is satisfied, the loop is stopped to obtain a finite number of IMF components. If it is not satisfied, the loop is returned to Equation (6) and is repeated until the convergence stopping condition is satisfied. The stopping condition is

$${{\displaystyle \sum }}_k\parallel {\stackrel{\wedge }{u}}_k^{n+1}-{\stackrel{\wedge }{u}}_k^n{\parallel }_2^2/\parallel {\stackrel{\wedge }{u}}_k^n\parallel <\epsilon ,$$

(10)

where $\epsilon$ is the convergence determination accuracy,$\epsilon$> 0.

Figure 1 shows a flow chart of the specific VMD.

From the process of variational model solution, it is clear that the way to effectively set the parameters of VMD in the process of decomposing lidar signal, retain all the information of the original sequence as much as possible, and avoid modal mixing is crucial to the effect of decomposition and the accuracy of the final model prediction results. The modal number $K$, the quadratic penalty factor $\alpha$, the iteration step $\tau$, and the convergence tolerance $\epsilon$ are all parameters that need to be set in advance for the VMD decomposition, and these parameters must be varied adaptively for different signals to achieve the best decomposition results. In fact, it is found that $\tau$ and $\epsilon$ have small impact on the decomposition results, and the default values are usually adopted directly in the standard VMD; whereas $K$ and $\alpha$ have the greatest impact on the decomposition results.

According to the analysis, it can be concluded that if $K$ is too small, under-decomposition will occur (i.e., the signal cannot be completely decomposed and some important information in the original signal will be lost) and if $K$ is too large, over-decomposition will occur (i.e., the decomposition result will have false components or the signal of the same frequency will be decomposed into different components). If $\alpha$ is too small, the bandwidth of the modal components will be too large, and some components will be included in other modal components. If $\alpha$ is too large, the bandwidth of the modal component will be too small, and important information in the original signal will be lost.

Therefore, when VMD decomposition of lidar signals is performed, an appropriate method is needed to select a suitable combination of [$K$,$\alpha$], which can suppress the phenomenon of modal mixing in VMD decomposition under noisy environment to a certain extent, and speed up the decomposition speed and improve the accuracy.

2.3. Sparrow Search Algorithm (SSA)

At present, there is no fixed method for VMD parameter selection; if the parameters are selected empirically, the results may be greatly biased. The SSA is a swarm intelligence optimization algorithm proposed to simulate the foraging and antipredator behaviour of sparrows and is based on their predator-avoiding behaviour during foraging [22], which has some advantages, such as the advantages of not easily falling into local optimal solutions, strong search traversal, fast convergence, short optimization search time, and more robustness.

In this paper, in order to achieve adaptive selection of the parameters of the VMD method, SSA is used to search for the optimal number of decomposition layers $K$ and penalty factor $\alpha$ to make up for the deficiencies of VMD to a certain extent. Considering that the SSA algorithm needs to determine an adaptation function when searching for parameters, the local minimal first value is used as the adaptation value in the optimization search process. The minimization of the local minimum first value is the final goal of the search. The accuracy of the SSA optimized VMD for lidar signal decomposition should be ensured; otherwise, it will indirectly affect the subsequent selection of the sensitive IMF components containing the most abundant information of the original lidar signal, and then affect the retrieval accuracy of the extinction coefficient profile, thus reducing the denoising effect of the lidar signal.

The whole sparrow population is divided into two fixed proportion subpopulations of discoverers and joiners according to the different division of labor, by which the scouts act as a certain proportion of individuals and send warning messages when there is danger. They update their positions, respectively, according to their respective rules.

Figure 2 shows the flow chart of the VMD-SSA algorithm.

2.4. SVD Algorithm

If the conventional SVD calculation method in matrix theory is used to decompose larger matrices, the computational workload for the experimenter is increased, and computational errors occur in the calculation process. After a long theoretical development period, a QR decomposition method in SVD was proposed by the mathematician Golub. The SVD principle is applicable to the decomposition of matrices. For a one-dimensional signal sequence, it is necessary to transform the one-dimensional signal sequence into a two-dimensional matrix. The two-dimensional matrix is decomposed into vector spaces by the SVD method, and then the noise-free components are selected for reconstruction.

The principle of SVD noise reduction is as follows: first, the original noisy signal is constructed into a Hankel matrix H of order p × q; then the SVD is used to decompose the matrix H into singular values to obtain the singular value matrix of the signal ${\sigma }_{i }$. Finally, the singular values corresponding to the noise are set to zero, and then the signal is reconstructed by using the SVD inverse operation to obtain the noise-reduced signal.

The real original signal x(t) is generally composed of two parts: the real signal s(t) and the noisy signal n(t). The energy of the real signal s(t) and the noisy signal n(t) are not correlated with each other. The energy of the real signal is relatively concentrated and the energy of the noisy signal is more dispersed, so the SVD can be regarded as dividing the Hankel matrix H into two uncorrelated spaces, i.e., the larger singular value corresponds to the reconstruction of the real signal, and the smaller singular value corresponds to the reconstruction of the noisy signal.

The SVD decomposition of the Hankel matrix H yields

$$H=U\sum V^T=U_s{{\displaystyle \sum }}_s{V}_s^T=U_n{{\displaystyle \sum }}_n{V}_n^T={{\displaystyle \sum }}_{i=1}^l{\sigma }_i{u}_i{v}_i^T ,$$

(11)

where ${\sigma }_i$ is the singular value corresponding to the demarcation between the real signal and the noisy signal. In the SVD noise-reduction process, the most critical step is to select the component matrices that match the characteristics of the original signal for reconstruction. Selecting fewer matrix components in the reconstruction process leads to serious distortion of the reconstructed signal, whereas selecting more matrix components results in too much noise in the reconstructed signal. In this paper, we use the singular value difference spectrum method to select the appropriate matrix components to reconstruct the noise-reduction signal.

2.5. VMD-SSA-SVD Noise-Reduction Method

We propose a noise-reduction method based on VMD-SSA-SVD and apply it to the processing of atmospheric lidar signals disturbed by noise. First, the SSA algorithm outperforms the existing algorithms in terms of search accuracy, convergence speed, stability and avoidance of local optima. The algorithm has good robustness, thus avoiding the intervention of human subjective factors. The SSA algorithm is applied to the parameter optimization of VMD, the optimal combination of influencing parameters is automatically filtered out [$K$,$\alpha$], and its decomposition yields several IMF components. Then, the IMF components and the original signal are subjected to the mutual relationship number calculation, the sensitive IMF components are screened according to the mutual relationship number selection criterion, and the sensitive IMF components are retained and reconstructed. Meanwhile the noise-containing components are eliminated. Because the sensitive IMF components containing the signal that was retained above are still affected by small noise, this paper uses SVD with good noise-cancellation ability to perform secondary denoising of the reconstructed signal, after its parameter optimization, to maximize the retention of useful information, eliminate useless information, improve the SNR and the accuracy of the noise-reduction algorithm, and perform successive segmentation SVD processing on each sensitive IMF component (i.e., convert the sensitive IMF components that are converted into Hankel matrices). Then successive segmented SVD is performed in turn. When the successive segmented SVD algorithm processes the lidar signal similar to the one containing DC components, the algorithm reconstructs the signal by selecting the second maximum mutation point in the singular value difference spectrum, and the reconstructed signal is the noise-reduction signal. In addition, four algorithms, namely wavelet packet, EMD, EEMD, and CEEMDAN, are used in the system to verify the effectiveness of this algorithm by comparing the noise-reduction effects of the five algorithms. The specific steps are as follows.

(1)

The signal is acquired and the optimal combination for penalty parameters α and modal numbers $K$ is found by using the VMD-SSA algorithm.

(2)

The signal undergoes VMD according to the parameters of (1).

(3)

The interrelationship number of each IMF component and the original signal is calculated, the interrelationship number threshold is determined, and the component greater than or equal to the threshold is selected as the sensitive modal component.

(4

The Hankel matrix of the sensitive modal component signal is obtained, SVD is performed, signal segmentation SVD noise-reduction processing is used, the effective order of the SVD is determined, and noise-reduction processing on the sensitive modal component through SVD is performed.

(5)

The modal components of the SVD are reconstructed to obtain the final noise-reduction signal.

Figure 3 shows a block diagram of the algorithm used in this paper.

2.5.1. Optimizing VMD Based on the SSA

In VMD, the decomposition mode number $K$ and the quadratic penalty factor $\alpha$ are important determinants of the decomposition results and have a great impact on the accuracy of the decomposition. Therefore, in practice, due to the complexity and variability of the signal, the artificial manual optimization of $K$ and $\alpha$ is unreliable and causes large deviation problems. Moreover, the selection of the two parameters is random and irregular. If the decomposition modal number $K$ is too large, it will be over decomposed; if the value of $K$ is too small, the decomposition will be incomplete and information will be lost. The secondary penalty factor $\alpha$ affects the bandwidth of the modal components, and the signal extraction results are also related to the bandwidth of different scales. Therefore, applying the sparrow search algorithm (SSA) to the optimization of the VMD parameters is considered to find the optimal combination of the two parameters.

The proposed SSA-based optimization VMD algorithm is detailed as follows.

(1)

In the case of the original signal input, f, the range of values of VMD parameters is determined, and the model parameters of SSA are set and initialized, which include the maximum number of iterations and the initial population number. The fitness function is the local minimal envelope entropy. In this paper, $K$ takes the integer in the interval [2, 15], $\alpha$ takes the integer in the interval [1000, 10,000], and the maximum number of iterations is set to 15. The number of populations is set to 30.

(2)

VMD of the signal and calculation of the initial fitness of each sparrow are next. Based on the optimization principle of the sparrow algorithm, the minimum fitness value is selected by updating the positions of the discoverer, the follower, and the position of the sparrow that is aware of the danger while continuously comparing the fitness values of each position.

(3)

Iterative loops are performed until the global optimal parameters and the best fitness values are determined or the maximum number of iterations is reached. The values are saved.

(4)

The previously obtained optimal parameters [$K$,$\alpha$] are used to perform the VMD of the signal.

(5)

The interrelationship values of each modal component with the original signal are calculated.

(6)

The effective components for noise reduction are selected and reconstructed.

2.5.2. Selection of Effective Components Based on the Number of Interrelationships

To extract the signal feature information more accurately, the decomposed effective IMF components are selected for reconstruction to achieve the elimination of noise in the original signal. The interrelation number method is used to select the effective IMF, which measures the amount of correlation information between two components and calculates the degree of interdependence that can be quantitatively expressed. This will be more accurate compared with the traditional correlation coefficient selection method. To achieve the selection of effective IMF components, we find the number of interrelationships between each IMF component and the original signal after decomposition, calculate the threshold value by the formula, and compare each correlation coefficient with its threshold value to select the component greater than the threshold value as the effective component. The number of interrelationships and the threshold value are calculated as

$${\rho }_{xy}=\frac{{{\displaystyle \sum }}_{n=1}^N\left. [\left(x_n-\overline{x}\right)\left(y_n-\overline{y}\right)\right]}{{\left. [{{\displaystyle \sum }}_{n=1}^N{\left(x_n-\overline{x}\right)}^2{{\displaystyle \sum }}_{n=1}^N{\left(y_n-\overline{y}\right)}^2\right]}^{\frac{1}{2}}}$$

(12)

$$\mu =\frac{\max \left({\rho }_{xy}\right)}{10\times \max \left({\rho }_{xy}\right)-3} ,$$

(13)

where $\max \left({\rho }_{xy}\right)$ is the maximum correlation coefficient, $x_{n }$ is the IMF component, $y_{n }$ is the original signal, and $\overline{x}$ and $\overline{y }$ are the mean values of the IMF component and the original signal, respectively.

2.5.3. SVD Segmentation Noise-Reduction Algorithm

To ensure the accuracy of the signal, the received signal is a long sequence signal. SVD increases the processing time when processing signals with large sample lengths. The nonstationary and nonlinear characteristics cause the denoised signal to lose detailed information. The short-time Fourier transform divides the nonlinear and nonstationary signal into several sequences of a certain length, and the more segments there are, the more it can be approximated as a smooth signal. Drawing on the short-time Fourier transform method, the nonstationary lidar signal is segmented into a stationary signal by adding a rectangular window of appropriate width. Therefore, the SVD algorithm is used to analyse the signal of each segment separately. Then the appropriate signal is selected to reconstruct each segment and connect each segment in turn to recover the noise-reduced signal. Figure 4 shows the flow chart of signal-segmentation SVD noise-reduction processing.

3. Results and Discussion

3.1. Evaluation Index of the Noise-Reduction Effect

To evaluate the noise-reduction effect and verify the effectiveness of noise reduction, two reference indicators are used: signal-to-noise ratio (SNR) and mean square error (MSE) [23]. SNR refers to the power ratio between the signal and noise, reflecting the quality of the signal. The higher the value is, the better the denoising effect. MSE refers to the variance between the original signal and the denoised signal, reflecting the similarity between the two signals before and after denoising, and the smaller the value is, the higher the similarity. The SNR and MSE are defined by the following equations:

$$SNR=10\mathit{log}\left(\frac{{{\displaystyle \sum }}_{i=1}^N{x}_i{}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-\widetilde{x_i}\right)}^2}\right)$$

(14)

$$MSE=\frac{1}{N}{{\displaystyle \sum }}_{i=1}^N{\left(x_i-\widetilde{x_i}\right)}^2,$$

(15)

where $x_i$ is the original signal, ${\tilde{x}}_i$ is the denoised signal, and N is the signal length.

3.2. Simulation Signal Verification

To verify the effectiveness of the VMD-SSA-SVD denoising method, four other methods, including wavelet thresholding, EMD, EEMD, and CEEMD, were used for comparison. Some simulation experiments were conducted by using a typical simulation signal similar to the lidar signal, i.e., the Bump test signal proposed by Donoho and Johnstone as the verification signal [24]. In a realistic environment, noise is often a complex phenomenon from many different sources. The lidar return signal is mainly susceptible to interference from thermal noise, scattered particle noise, background radiation noise, and generation of composite noise. It is assumed that the real noise is regarded as the sum of many random variables with different probability distributions, and each random variable is independent and obeys a Gaussian distribution according to the central limit theorem, which is mainly normalized to Gaussian white noise. The purpose of adding noise with a specified SNR to the simulation experiment is to quantify and evaluate the effect of the algorithm.

Figure 5 and Figure 6 show the simulated noise-reduction effects of the noisy Bump and Block signals with five different methods, respectively. In this paper, Gaussian white noise with SNR = 10 is added to both the Block test signal and the Bump test signal, and it can be seen that these five noise-reduction methods have a certain suppression effect on the noise. The wavelet packet noise-reduction method uses the db1 wavelet basis function for the three-layer decomposition. The EEMD method processes the noisy Bump signal, adding white noise with a standard deviation of 0.2 and averaging 100 times. The SSA searches for the optimal VMD parameter combination of the Bump signal. The convergence curve of parameter optimization shows that the optimal combination of parameters for $K$ is 3, and is 1613 for $\alpha$, and the number of SVD segments is 80. The five methods all improve the SNR of the signal. The effect of EMD and CEEMD processing is not satisfactory compared with the unprocessed signal, which exhibits a serious distortion phenomenon. The signal denoised by EEMD has many “burrs” and almost does not filter out most of the noise, so the processed SNR is not high. The VMD-SSA-SVD method has the best noise-reduction effect, which retains the signal details and almost restores the original Bump signal. The SNR after noise reduction is higher than that of the other methods, and the MSE is smaller than that of the other methods. It should be pointed out that because the Block test signal and Bump test signal have no unit, the amplitude of Figure 6 has no unit. Moreover, the point is the length of the selected data, because the test signal is periodic in nature, the data length can be determined according to signal integrity.

Figure 7 and Figure 8 show the denoising performance of the above methods for the noisy Bump and Block signals when the SNR varies from −5 dB to 10 dB. As seen, the VMD-SSA-SVD method achieves the expected best results for different input SNRs. Even in the worst case, the SNR of the Bump signal processed by VMD-SSA-SVD is still guaranteed to be as high as 7.93 dB when the input SNR is −5 dB, and the SNR of the Block signal is as high as 7.85 dB. In addition, we also compared the MSEs obtained by each method for different input SNRs. The results show that the VMD-SSA-SVD method can achieve the minimum error when reconstructing Bump and Block signals.

To demonstrate the applicability of the proposed VMD-SSA-SVD method and evaluate the effectiveness of the five noise-reduction methods above,

Table 1. The SNR and MSE after noise reduction by five methods when superimposed with noise with different SNRs.

Signal	SNR	Wavelet Threshold		EMD		EEMD		CEEMD		VMD-SSA-SVD
		SNR	MSE	SNR	MSE	SNR	MSE	SNR	MSE	SNR	MSE
Block	−5 dB	−4.9435	10.36	−5.711	12.363	−4.9089	10.278	−5.7292	12.415	7.8549	0.5439
	0 dB	−0.1372	3.4257	−1.4703	4.6565	0.0014	3.3181	−1.5215	4.7116	12.261	0.1979
	5 dB	4.9576	1.0599	1.5734	2.3104	5.2345	0.9944	2.2105	1.9952	16.071	0.082
	10 dB	18.691	0.0448	3.3386	1.5387	10.922	0.2684	9.9553	0.3353	20.682	0.0283
Bump	−5 dB	−4.9435	27.471	−5.843	33.793	−4.8836	27.094	−5.808	33.521	7.9391	1.4145
	0 dB	−0.1372	9.0832	−0.1268	9.0616	−0.016	8.8336	−0.3599	9.5612	11.47	0.627
	5 dB	4.9576	2.8104	3.7765	3.6887	4.9525	2.8137	4.2177	3.3324	15.133	0.269
	10 dB	11.98	0.5579	9.3555	1.0209	9.8018	0.9211	7.2264	1.668	17.292	0.164

shows the SNRs and MSEs of the Block signal and the Bump signal after noise reduction by the five methods under different input SNRs with superimposed Gaussian white noise of −5 dB, 0 dB, 5 dB, and 10 dB. The corresponding best SNR and MSE are shown in bold. The table shows that the VMD-SSA-SVD method outperforms the other four methods for the “Bumps” and “Blocks” signals in the listed input SNR range.

3.3. Comparison of Atmospheric Lidar-Measured Signal-Denoising Effect

The actual signals used in this paper were all detected by a small Mie scattering lidar developed at North Minzu University (106°06′ E, 38°29′ N).

Table 2. System parameters of the Mie scattering lidar.

Parameters	Parameter Index
Laser	Nd:YAG laser
Wavelength	1064 nm, 532 nm, 355 nm
Single pulse energy (1064 mm/532 mm/355 mm)	350 mJ/170 mJ/80 mJ
Impulse frequency	1–10 Hz
Pulse width (1064 nm)	≤10 ns
Beam diameter (1064 nm)	~9 mm
Divergence angle (1064 nm)	≤0.5 mrad
Telescope	Schmidt–Cassegrain type
Diameter	300 mm
Field of view angle	0.4 mrad
Horizontal angle	0–360°
Pitching angle	0–90°
Detection range	100–13,000 m (Night), 100–10,000 m (Daytime)

shows the system parameters of the Mie scattering lidar. There are 10,000 data points for each acquired raw signal, and the first 1000 data points are background noise that can be used to calculate the average noise of the data acquisition system and are subtracted during data preprocessing. Therefore, the starting point of the calculation was set to 1000.

Figure 9 shows the aerosol extinction coefficient profile and the corresponding denoising effect at a wavelength of 532 nm inverted by Klett’s method. The original signal Figure 9a contains considerable noise, and the useful signal is disturbed by the noise after 3 km. As the signal becomes very weak after 5 km, the noise interference also increases, and the signal is completely submerged in noise after 5 km. The extraction of the aerosol extinction coefficient profiles is improved by all five noise-reduction methods. Figure 9b shows the extinction coefficient profile denoised by the EMD method. The noise-reduction effect is not obvious, and the distortion is very serious and there is even a large abrupt change below 2 km. Figure 9c shows the extinction coefficient profile denoised by the EEMD method. The useful signal is still drowned in the noise, the noise-reduction effect is not obvious, and there is still serious noise interference after 4 km. Figure 9d shows the result using the CEEMD method. Although the curve is exceptionally smooth, more signal components are lost after 8 km, and the signal details are lost. Figure 9e shows the denoised profile using the wavelet packet method. Although the signal after 4 km is clearer than the original signal, the signal waveform is more obvious, the noise-reduction effect is somewhat improved, and most of the noise is filtered out; the overall profile is not smooth and noise is present. Figure 9f shows the result processed by the VMD-SSA-SVD method. Using the mutual relationship number as the IMF component screening condition and adding the SVD noise reduction, the noise can be removed while retaining more details of the extinction coefficient profile; the curve is smoother, the profile will not be distorted, and the noise interference of the far-field profile can be suppressed. In fact, for the far-field signal, the VMD-SSA-SVD method has a better noise-reduction effect and signal detail retention than the signals processed by the other four noise-reduction methods. Therefore, the proposed VMD-SSA-SVD noise-reduction method has higher accuracy than the above four noise-reduction methods. Moreover, not only is the noise suppression of the far-field signal very obvious but the denoised signal is closer to the real signal.

To more intuitively evaluate the effect of these five methods on the processing of this group of signals, the corresponding SNRs and MSEs were calculated and are listed in

Table 3. The SNR and MSE of signals denoised by the five methods are shown in Figure 9.

	Wavelet Threshold	EMD	EEMD	CEEMD	VMD-SSA-SVD
SNR	16.3059	8.4232	11.9542	16.2208	21.1368
MSE	5.299 × 10−7	1.412 × 10−6	3.139 × 10−7	2.345 × 10−7	7.56 × 10−8

. It can be seen that the proposed method has a higher SNR and a lower MSE than the first four and has the best noise-reduction effect.

To further verify the denoising effect of the VMD-SSA-SVD algorithm on the complex return signal, we used the return signal containing one cloud layer. Figure 10 shows the profiles of the extinction coefficient inverted from the return signal containing one layer of obvious clouds. There is an obvious cloud layer between 3 km and 4 km, in which the distribution of the extinction coefficient is smooth. At approximately 5 km, the signal intensity and SNR decrease due to the attenuation of the laser beam energy, the trend of the extinction coefficient becomes increasingly blurred, and the useful signal is almost completely submerged in the noise. Figure 10b shows the extinction coefficient profile denoised by the EMD method, and the signal is not only severely distorted but also has abnormal mutation at 3 km. Figure 10c shows the extinction coefficient profile denoised by the EEMD method. The noise is hardly filtered out, especially at the peak of the extinction coefficient of the cloud layer, and the value of the extinction coefficient has a large distortion compared with the unprocessed value. Figure 10d shows the extinction coefficient profile obtained by the CEEMD method. Although the curve is very smooth, the signal is lost after 10 km, the distortion is very serious, and even if there are some abrupt changes, the noise elimination effect is far inferior. Figure 10e shows the extinction coefficient profile denoised by the wavelet packet method. Its noise-elimination effect and signal fidelity are relatively better than those above. Figure 10f shows the extinction coefficient profile processed by the VMD-SSA-SVD method. It is clear that it retains the details of the intact signal, can effectively identify the changing characteristics of the extinction coefficient, and does not lose the spatial resolution and signal fidelity to obtain the noise-elimination effect.

Table 4. The SNR and MSE of signals denoised by the five methods are shown in Figure 10.

	Wavelet Threshold	EMD	EEMD	CEEMD	VMD-SSA-SVD
SNR	17.3105	7.8042	9.9699	14.7930	22.3783
MSE	2.699 × 10−7	5.302 × 10−6	3.239 × 10−6	3.749 × 10−6	2.16 × 10−7

lists the SNRs and MSEs of the denoised profiles in Figure 10. The VMD-SSA-SVD method has a higher SNR and lower MSE than the other four methods, so its noise-reduction effect is the best among the five methods.

Figure 11 shows the extinction coefficient profile inverted from the return signal containing clouds; compared to Figure 10, this set of signals contains two distinct cloud layers. Figure 11a shows the original extinction coefficient profile containing a large amount of noise. There is a clear cloud layer between 3 km and 4 km with an extinction coefficient value of 0.049 ${\mathrm{km}}^{-1}$ and another clear cloud layer between 4 km and 6 km with an extinction coefficient of 0.038 ${\mathrm{km}}^{-1}$. The useful signal at approximately 6 km is almost completely submerged in noise. Figure 11b shows the extinction coefficient profile denoised by the EMD method. The curve is severely distorted and only shows the cloud layer between 3 km and 4 km, the signal of the cloud layer between 4 km and 6 km is lost, there are abnormal abrupt changes, especially at the peak of the extinction coefficient of the cloud layer, and the value of the extinction coefficient has a large change compared with the original value. Figure 11c shows the extinction coefficient profile processed by the EEMD method. The noise is almost not filtered out. Still, at the peak of the extinction coefficient of the two cloud layers, the extinction coefficient value produces some distortion. Figure 11d shows the extinction coefficient profile denoised by the CEEMD method. Although the curve is very smooth, it does not retain good signal details. There is also distortion at the peak of the extinction coefficients of the two cloud layers. Figure 11e shows the extinction coefficient profile obtained by the wavelet packet method, which has a relatively better noise-reduction effect and signal fidelity. Figure 11f shows the extinction coefficient profile processed by the VMD-SSA-SVD method, which not only retains the details of the signal completely but also does not produce distortion at the peak of the extinction coefficient. This indicates that this method has the best noise-elimination effect and higher accuracy.

Table 5. The SNR and MSE of signals denoised by the five methods are shown in Figure 11.

	Wavelet Threshold	EMD	EEMD	CEEMD	VMD-SSA-SVD
SNR	20.5950	8.6486	11.0205	17.5588	23.6226
MSE	5.285 × 10−6	1.161 × 10−5	3.904 × 10−5	4.235 × 10−6	7.23 × 10−7

lists the SNRs and MSEs of Figure 11. It is clear that the VMD-SSA-SVD method outperforms the other four in terms of both SNRs and MSEs.

4. Conclusions

From the above simulated and measured data, the denoising effect of EMD, EEMD, and CEEMD is the most unsatisfactory, and the processed signal has obvious distortion and large, abrupt changes, losing the original components of the signal. The wavelet packet method is better than the above three methods in terms of the denoising effect and signal fidelity, but the signal details are not fully restored. Moreover, there are subtle distortions. In this paper, a lidar signal-denoising method based on VMD-SSA-SVD is proposed. This method addresses the problem of optimal selection of VMD layers and penalty parameters based on the sparrow search algorithm and further denoises the initially decomposed signal by the successive segmentation SVD algorithm. In this paper, wavelet packet, EMD, EEMD, CEEMD and other common denoising methods are considered for comparison. In addition, the advantages and disadvantages of each noise-reduction method are analysed, and the comparative analysis shows that the proposed VMD-SSA-SVD processed signal has the highest SNR and the smallest MSE. This method can effectively improve the SNR of the lidar return signal while preserving the complete characteristics of the signal. Additionally, the denoising effect and signal fidelity of this method are better than those of the four other methods, which demonstrates its effectiveness, superior performance, and strong adaptability in lidar return signal denoising.