- Email: [email protected]
A novel noise reduction method for space-borne full waveforms based on empirical mode decomposition
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.com/locate/ijleo
Original research article
A novel noise reduction method for space-borne full waveforms based on empirical mode decomposition Zhijie Zhanga, Xiangfeng Liub,*, Rong Shub, Feng Xieb, Fengxiang Wangb, Zhihui Liub, Hanwei Zhanga, Zhenhua Wangc a
School of Surveying and Land Information Engineering, Henan Polytechnic University, Jiaozuo, China Shanghai Institute of Technical Physics of The Chinese Academy of Sciences, Shanghai, China c College of Information, Shanghai Ocean University, Shanghai, China b
A R T IC LE I N F O
ABS TRA CT
Keywords: Full waveform Noise reduction Empirical mode decomposition Laser altimetry GaoFen-7
In this project, a novel noise reduction method based on empirical mode decomposition (EMD), referred to as EMD-AdaptiveP, is proposed and investigated with the expectation of providing an effective enhancement for the full waveform of nonlinear and nonstationary signals obtained by space-borne laser altimetry and the forthcoming GaoFen-7 (GF-7) mapping satellite in China. This method combines EMD with a similar strategy of wavelet filtering to adapt to the characteristics of each waveform. The reconstruction of an effective waveform signal is implemented through reverse superimposition of its intrinsic mode functions (IMFs) and the residual, which is thresholded by the three-sigma principle of background noise of the waveform signal and an additional check mechanism in the lower or higher signal of noise levels, after the waveform is decomposed by EMD. This method is qualitatively and quantitatively examined with simulated waveform data, real ice, cloud and land elevation satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) waveforms and GF-7 proto-waveform data, in contrast with other IMF selection schemes, like EMD-soft, EMD-hard and EMD-Wavelet, or the traditional filtering methods, such as Gaussian, μ\λ and wavelet. The results show that the method of EMD-AdaptiveP (1) is sufficiently robust to denoise the waveform signal with different levels of signal-to-noise ratio, (2) has the adaptive ability to distinguish effective signal IMFs from noise IMFs in signal reconstruction, and (3) performs noise reduction more effectively than traditional methods.
1. Introduction Space-borne laser altimetry is an active remote sensing technique with high temporal‒spatial resolution and a high accuracy range, which can provide rapid acquisition of geospatial information in large regions and at high-precision benchmarks on a global scale, even in special areas (e.g., polar regions) [1]. Since the ice, cloud and land elevation satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) for Earth observation launched in 2003 [2], full waveform data have been widely used for the observation of polar ice sheet change; vertical distribution of vegetation, clouds and aerosols; estimation of biomass; sea surface height; and monitoring of climatic phenomena [1], water and land elevation [3–5]. GaoFen-7 (GF-7), the first high-resolution stereo mapping satellite in China, equipped with a double-line array camera and a laser altimeter, will be launched in 2019. Compared with discrete return systems, the full waveform, with the characteristics of waveform shapes, widths, intensities and skewness for the laser echoes,
⁎
Corresponding author. E-mail addresses: [email protected], [email protected] (X. Liu).
https://doi.org/10.1016/j.ijleo.2019.163581 Received 17 July 2019; Accepted 10 October 2019 0030-4026/ © 2019 Elsevier GmbH. All rights reserved.
Please cite this article as: Zhijie Zhang, et al., Optik - International Journal for Light and Electron Optics, https://doi.org/10.1016/j.ijleo.2019.163581
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
not only incorporates responses from ground level but also provides additional geometric and physical information about multiplescattering responses along the path [6]. However, full waveform data are always contaminated by various noises in the process of transmission and acquisition, such as ray radiation from the sun and sky; dark current noise; and thermal noise from the optical detection system [7]. Furthermore, the echo signal is very weak in its laser intensity with the increase of the detection range [8]. That would be the major impediment to full waveform decomposition and extraction of ground feature parameters (e.g., ranges, slope and roughness) [9,10]. The waveform of echo signals demonstrates both nonlinear and nonstationary characteristics, and it may be submerged in the noises. Therefore, noise reduction of the full waveform is an indispensable key step in improving the retrieval accuracy of the data and ensuring the reliability of elevation-related applications [11]. At present, fixed-width Gaussian filters [12], average filters [13], μ\λ filtering [14], Fourier low-pass filtering [15] and wavelet filtering [16] have been extensively used to denoise the full waveform. However, some drawbacks come with these methods; fixedwidth Gaussian filters, average filtering and μ\λ filtering face interference by improper filtering parameters, which lead to waveform distortion, such as peak amplitude shrinkage or pulse width increase. Fourier filtering only demonstrates the signals in a frequency domain or a time domain; the unsuitable cut-off frequency may cause signal distortion. Although wavelet filtering shows frequency and energy content in the time domain, it needs to choose a priori basis function with the original signal, based on experience, and a preset wavelet basis function cannot adapt to all kinds of waveforms. In short, the common problem with these methods is that they need to suitably select the filtering parameters. Inappropriate parameters will cause distortion of waveforms [7]. However, empirical mode decomposition (EMD), devised by Huang et al. [17], can filter and extract the trend term with nonlinear and nonstationary signals, which is a fully signal-dependent approach based on the local characteristic time scale that adaptively decomposes a given signal into several intrinsic mode functions (IMFs) and a residual. Researchers were already familiar with this practice from various forms of signal filtering, such as acoustic signals [18], phonocardiograms [19], pressure signals of loop airlift reactors [20], signals for fault diagnosis of rolling bearings [21] and signal filters attempted in the enhancement of waveform data [22]. The basic principle of these methods is to reconstruct the signal within previously filtered IMFs by direct removal of the first several IMFs with high frequency [23], combination with the traditional filters [24], thresholding for preset threshold functions [25] or fast Fourier transformation (FFT) [22] toward all IMFs to select the high frequency IMFs. The common problem in the above EMD-based filters is that there is no easy way to adaptively distinguish between the noise IMFs and effective signal IMFs and denoise the signal without losing useful information; this is still an open issue [26]. Therefore, to effectively reduce the various noises in full waveform and apply that reduction technique to the future GF-7 full waveform, a novel noise reduction method based on EMD, named as EMD-AdaptiveP, is proposed in this paper. It can adaptively distinguish effective signal IMFs from noise IMFs and select the threshold for the construction of a full waveform with different levels of signal-to-noise ratio (SNR). To test and verify the effectiveness of the proposed method in a full waveform signal, some experiments with simulation waveform, GLAS waveform and the GF-7 proto-waveform were carried out. Meanwhile, the EMD-AdaptiveP was compared to other different IMF selection schemes, such as EMD-soft, EMD-hard, EMD-Wavelet and so on. Furthermore, comparisons of this method with the traditional methods of Gaussian, μ\λ and wavelet were performed with real waveform data. The rest of this paper is organised as follows: In Section 2, we describe the noise reduction method based on EMD. In Section 3, this method is examined along with the simulation waveform and real waveform from laser altimetry, comparing it to other IMF selection and traditional noise reduction methods. Finally, in Section 4 the conclusions are drawn from the data gathered.
2. Proposed method For the full waveform, the noises are regarded as white Gaussian noise [10], including quantum noise, thermal noise, amplifier noise, dark current noise, background noise and so on to divide the IMFs of the original full waveforms into noise IMFs and effective signal IMFs and determine which IMFs contain most of the effective signal components, beginning from a coarse, and progressing to a fine scale. The reconstruction of the denoised signal is obtained by summing up the effective IMFs reversely until the denoised full waveform meets the preset three-sigma threshold. The reasonability and validity of this idea is discussed as regards the following aspects: the completeness of full waveforms after EMD decomposition, the estimation of involved noise parameters to set the threshold of stopping signal reconstruction and the conception of the overall procedure.
2.1. The completeness of full waveform EMD decomposition Compared to wavelet filtering, EMD is a fully data-driven method and does not necessitate any basis function as a prior, based on the sequential extraction of energy associated with various intrinsic time scales of the signal, starting from the finer temporal scales (high frequency modes) and moving toward coarser ones (low frequency modes). Each mode should possess an equal number of extrema and zero crossings. The detail of the EMD algorithm was described in references [27,28]. To validate the completeness of EMD in space-borne full waveform processing, a typical GLAS waveform was used. There are six IMFs decomposed by EMD, as shown in Fig. 1, and the first IMF has the highest frequency component, while the frequency decreases with the increment of the IMFs’ order. Fig. 2 shows the original full waveform and reconstruction error between the sum of all the IMFs and the residual, which displays the decomposition of EMD thoroughly. The reconstruction error demonstrates that the total sum of the IMFs and the residual matches well with the original signal and that the order of amplitude reached 10−17, which can be attributed to the rounding error of computers. The completeness lays out the foundation of the following IMFs’ selection analysis and waveform reconstruction. 2
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 1. IMFs of a typical GLAS full waveform decomposed by EMD.
Fig. 2. A typical GLAS full waveform and reconstruction error concerning all the IMFs and the residual.
2.2. Estimation of noise mean and standard deviation Determining how to design the threshold that would stop iterations of IMFs is key to properly filtering noise signals and preserving the characteristic waveform signal. Usually, the value of noise is unknown, while the background noise mixed in the returned waveform signal has a relatively stable standard deviation. The beginning and ending parts of a waveform signal are generally considered stable background noise [29]. In this paper, we adopt this method to estimate the mean value and the standard deviation of noise, calculated as Eqs. (1) and (2). N
μnoise =
σnoise =
∑i =1
Vi (1)
N
N
∑i =1
(Vi − μN )2 N−1
(2) 3
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 3. The workflow for the denoising strategy (checking constraints in the case of the lower and higher SNR levels shown in the red dotted frame).
Where Vi is the amplitude of the i th sample, N is the number of samples and μN is the mean value for the sampled noise. In practice, echo components for the unknown detecting range may appear at the beginning or ending of the waveform. These may lead to an overestimated noise level with a higher estimated mean and standard deviation. Therefore, we take the smaller value from the beginning noise ( μbeg ) and the ending noise ( μend ) and even take the corresponding standard deviation (σbeg , σend ) to be the noise threshold in the event that we miss the effective signal 2.3. Designing of the proposed filter for full waveform signal After a waveform is decomposed by EMD, the adaptive threshold for waveform reconstruction from the reverse superposition of its IMFs and the residual is based on the root mean squared error (RMSE) between the original waveform amplitude and the denoised waveform amplitude, which should be three times the standard deviation over the mean background noise. That means the differing amplitude between the pre- and post-denoising waveforms should fall within the requirements for categorisation as noise. Given that most of the significant structures of the signal are concentrated in the lower frequency modes (last IMFs) and decrease toward the high-frequency modes (first IMFs), the strategy behind this filtering method is to partially reconstruct the effective parts of the decomposed signal from the last several IMFs, up to the penultimate IMF, where the number P of the last IMF is adaptive to the noise levels of each waveform. Thus the method known as EMD-AdaptiveP differs from other methods by directly removing the first few high-frequency IMFs, whose collective denoising effects depend on the setting of their number, or by using soft and hard thresholds to compensate for the effective signals that may be removed by directly discarding the high-frequency IMFs. Fig. 3 shows the entire procedure of this strategy. In the experiments with simulated waveforms in the next section, when the noise is at medium level, the basic noise threshold can easily be met with better noise suppression effects. However, when the noise is 4
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
lower (e.g., 10db or 15db), results show that the reconstructed signal by summing the last one or two IMFs can also meet the threshold with a severe loss of effective signal parts. In this case, it is necessary to add a check mechanism for the waveforms with lower SNR levels. We propose that the difference between the maximum peak of the original waveform and the smoothed waveform should also be below the noise threshold. It would also need to be bigger than the noise mean plus 1.5 times the standard deviation to avoid an insufficient smoothing effect. For the high SNR level, all the IMFs are summed, and none of the IMFs that cannot denoise the full waveform are discarded. Therefore, the other constraint, stipulating that the extreme order (k = maxIMFs–1) of IMFs is the penultimate IMF, is applied. 3. Experimental results and analysis To test and verify the effectiveness of the proposed novel filter method, experiments using simulation waveforms, GLAS waveforms and GF-7 proto-waveforms are conducted in this section. Furthermore, the qualitative and quantitative comparisons and analysis among other EMD-based methods and traditional filtering methods are carried out. 3.1. Simulation waveform Referencing the GLAS samples of land type waveforms and the parameters of GF-7, the simulation waveforms were simulated and tested. For the transmitted pulse with the form of the Gaussian model, the mixed Gauss model was adopted when simulating the pure received waveform signal. In the experiment, random coarse numbers (1, 10), (6, 40) were set as the amplitude value and the peak width, respectively; a random integer (100, 500) was set as the peak location. The maximum number of peaks was three. Then the random (normally distributed) white Gaussian noise was superimposed on the pure signals. To simulate the effect of digitiser sampling, the return samples were rounded to integer amplitude values, with a length of 544 ns. 3.1.1. Full waveform reconstruction under different SNR levels To validate the effectiveness of this method, two groups of waveforms with fixed peak numbers under different SNR levels and varied peak numbers with a fixed SNR level denoting the complexity of the waveforms were tested. Fig. 4 and 5, respectively, show the recovery processes of the signals by the residuals and their IMFs. In Fig. 4(a)–(d), the simulated full waveforms have three peak numbers under different SNR levels of 10db, 20db, 30db and 40db. The results show that with the increase of SNR, the number of noise IMFs decreases. The discarded first few IMFs are 3 (Fig. 4(a)), 2 (Fig. 4(b)), 1 (Fig. 4(c)) and 1 (Fig. 4(d)), respectively. In Fig. 5(a)‒(d), the simulated full waveforms have varied peak numbers of 1, 2, 3 and 4 under the SNR level of 15db. The results show that with the increase of the peak numbers, the numbers of noise IMFs vary. The discarded first few IMFs are 3 (Fig. 5(a)), 2 (Fig. 5(b)), 2 (Fig. 5(c)) and 2 (Fig. 5(d)) respectively, which all maintain the peak shape well, even with two very close peaks (Fig. 5(b)). These demonstrate that the proposed EMD-AdaptiveP method can adaptively reconstruct the varying waveform signals with different SNR levels and peak numbers while preserving the detailed waveform shape effectively. 3.1.2. Evaluation and comparison of denoising performance To evaluate the performance of this method, the EMD-AdaptiveP was compared with five other IMF selection schemes, using a typical simulated full waveform by directly removing either the first or the first two high-frequency IMF methods (EMD-1IMF or EMD-2IMFs); to avoid the loss of the effective waveform signal, after the direct removal of the first two IMFs, the universal threshold (either EMD-soft or EMD-hard) was used to extract the effective parts, which were then added to the reconstructed signal; the EMDWavelet was then performing wavelet filtering of each IMF decomposed by EMD. In this experiment, depending on the characteristics of the waveform, DN4 was chosen as the wavelet basis, and the soft threshold method was used to reconstruct the waveform signal. The denoising effect was as shown in Fig. 6(a)‒(f). The results obviously illustrate that EMD-1IMF has an inadequate denoising effect. There are still lots of high-frequency noises in denoised waveforms. The denoising effects of the remaining four comparison methods are similar, except in the case of EMD-Wavelet (Fig. 6(e)), which caused the decline of the peak value. The proposed EMD-AdaptiveP smoothed the waveform nicely. To validate the effectiveness of this method and avoid the resulting randomisation for each individual waveform, 2000 waveforms with different SNR levels were examined in this experiment. The quantitative assessment indicator was the consistency ratio of waveform decomposition from the Levenberg-Marquardt method. The consistency ratio is defined by two constraints: 1) the fitting number of the Gaussian components is equal to that of the simulated ones; 2) the average error of peak position is less than 1 ns. The statistics of consistency ratio for the experiments were determined with the maximum, minimum, median, mean and standard deviation (Std.), as shown in Fig. 7 and Table 1. Table 1 shows that the EMD-AdaptiveP has the highest consistency ratio of 0.7045 with a Std. of 0.0715. From the perspective of mean value, the consistency ratio of EMD-AdaptiveP has improved by 33.21%, 5.44%, 6.27%, 5.18% and 4.56%, compared with the EMD-1IMF, EMD-2IMFs, EMD-soft, EMD-hard and EMD-Wavelet methods, respectively. Looking at the trend of the performances as shown in Fig. 7, the EMD-AdaptiveP has a much better effect than do the other methods, especially when the SNR levels are lower. The consistency shows a downward trend when the waveforms have a higher SNR (about 35db). The reason may be that there was insufficient denoising with the removal of one IMF or that removing multiple IMFs led to the loss of effective signals. However, the maximum mean of 0.6394 illustrates the strong robustness of this approach under varying SNR levels. 5
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 4. Reconstruction of simulated waveform with fixed peak numbers under different SNR levels of 10db, 20db, 30db and 40db, respectively, by reverse IMF superimposition (the blue solid line is the original noisy waveform; the red solid line is the sum of the residual and effective signal IMFs’ components).
3.2. ICESat/GLAS waveform This method was tested with real ICESat/GLAS waveforms and then compared with traditional filtering methods, such as Gaussian filtering, μ\λ filtering and wavelet filtering. In the experiment, some quantitative indices [29,30] such as mean squared error, (MSE), mean absolute error (MAE) and cross correlation (R2) between the original noisy waveform and denoised waveform were taken to assess the degree of deviation and similarity between the methods. The output SNR of the signal and the peak signal-to-noise ratio (PSNR) were also taken to assess whether more useful signals were contained in the denoised waveforms. These indices can be calculated with the following Eqs. (3)‒(7):
MSE =
1 N
N
MAE =
N
∑n=1 (Raws (n) − Smooths (n))2
∑n=1
(3)
|Raws (n) − Smooths (n)| (4)
N 6
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 5. Reconstruction of simulated full waveform, with varied peak numbers of 1, 2, 3 and 4, respectively under SNR, equals 15 after reverse IMF superimposition (the blue solid line is the original noisy waveform; the red solid line is the sum of the residual and effective signal IMFs’ components). N
⎛ ⎞ ∑n=1 |Raws (n)|2 SNR = 10*log ⎜ N ⎟ ⎜∑ |Raws (n) − Smooths (n)|2 ⎟ ⎝ n=1 ⎠ N
PSNR = 10*log
∑n=1 N
∑n=1
max 2
(Raws (n) − Smooths (n))2
(6) 2
N
R2 =
(5)
(Raws (n) − Raws (n))(Smooths − Smooths ) ⎤ ⎡∑ ⎣ n=1 ⎦ N
∑n=1
N
(Raws (n) − Raws (n))2∑
n=1
(Smooths − Smooths )2
(7)
Where Raws represents the original noisy waveform signal; Smooths represents the denoised waveform signal; and N is the length of the waveform sample 3.2.1. Comparison with traditional filtering methods For Gaussian filtering, the width 6 ns of the GLAS transmit pulse was used as the Gaussian kernel, and 13 bins were adopted as the filtering length; for the μ\λ filtering, the values of λ and μ were 0.6307 and -0.6372, respectively; for the wavelet filtering, the “Daubechies 4” basis with the level 4 composition and the soft threshold method was employed; for the EMD-AdaptiveP method, 50 samples of beginning and ending parts of waveforms were used to estimate the noise mean and standard deviation. In order to show 7
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 6. Denoised full waveform obtained from five IMF selection schemes: (a) EMD-1IMF, (b) EMD-2IMFs, (c) EMD-hard, (d) EMD-soft, (e) EMDWavelet, (f) EMD-AdaptiveP (the blue one is the noisy full waveform, the red one is the denoised waveform).
Fig. 7. Consistency ratio variation of simulated full waveform along with the increase of SNR levels. Table 1 The statistics of consistency ratio of different IMF selection schemes. Items
EMD-1IMF
EMD-2IMFs
EMD-soft
EMD-hard
EMD-Wavelet
EMD-AdaptiveP
Maximum Minimum Median Mean Std.
0.5145 0.3555 0.4975 0.4800 0.0568
0.6485 0.4320 0.6405 0.6064 0.0788
0.6490 0.4535 0.6235 0.6017 0.0671
0.6515 0.4290 0.6440 0.6079 0.0805
0.6750 0.4340 0.6380 0.6090 1.2202
0.7045 0.4985 0.6685 0.6394 0.0715
the denoising effect of these methods more clearly, only the parts of the effective waveform signal based on about 200 samples are displayed (the original waveform samples numbered 544) in Fig. 8(a)‒(d). The results show that each method has effectively denoised the waveform signal. However, there is an obvious peak decrease in the denoised signal conducted by Gaussian and μ\λ filtering, as marked with a blue circle in Fig. 8(a) and (b). In addition, there is a spike noise accompanying the 130 ns treated by wavelet filtering, as marked with a blue circle in Fig. 8(c). EMD-AdaptiveP, as shown in Fig. 8(d), not only smooths the waveform but also maintains the main characteristics of the waveform. According to the statistical results listed in Table 2, EMD-AdaptiveP has the highest output SNR (20.285325 dB) and PSNR (33.589603); the smallest MSE (0.000106) and MAE (0.007867); and the biggest R2 (0.99435), which demonstrates that the method can greatly enhance the output SNR and PSNR, outdoing traditional filters. 3.2.2. Comparison with EMD-based filtering methods In addition, there were some qualitative and quantitative analyses conducted with different EMD-based filtering methods for the same GLAS waveform. In this waveform, the first IMF was removed by EMD-AdaptiveP, which is the same as EMD-1IMF. Therefore, EMD-1IMF was omitted in the comparisons. Fig. 9(a)‒(d) illustrates the denoised waveforms conducted by EMD-2IMFs, EMD-hard, EMD-soft and EMD-AdaptiveP. The results clearly show that many useful signal components were discarded by EMD-2IMFs, as marked by the red rectangle boxes in Fig. 9(a). There were sudden changes or distortions due to EMD-hard or EMD-soft thresholding, 8
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 8. Denoised ICESat/GLAS waveform by different filtering methods (the red one is the original noisy waveform, the green one is the denoised waveform). Table 2 Quantitative evaluation of denoised effects by the traditional methods and by EMD-AdaptiveP. Items
MSE
MAE
SNR
PSNR
R2
Gaussian μ\λ Wavelet EMD-AdaptiveP
0.00033 0.00029 0.00078 0.00011
0.01158 0.01094 0.01739 0.00787
15.35926 15.93050 11.60523 20.28533
28.66354 29.23478 24.90950 33.58960
0.98458 0.98627 0.95724 0.99435
as marked by the red rectangle boxes in Fig. 9 (b) and (c). However, the EMD-AdaptiveP method performed more effectively than these methods, and the statistical results listed in Table 3 further demonstrate that EMD-AdaptiveP has the biggest output PSNR (33.58960) and the highest correlation (0.99435) with the original waveform, which gives it advantages over other EMD-based methods. 3.3. GF-7 proto-waveform Because the GF-7 satellite has not yet been launched, an echo waveform acquired by laser altimetry was developed in a simple surface test environment. Thus, this waveform is referred to as a proto-waveform, which may be different from a waveform in a real in-orbit environment. However, it does not affect the validation of the filtering method. As with the processing of ICESat/GLAS waveforms, when denoising the proto-waveform, different traditional filters and EMD-based denoising methods were used to test the effectiveness of the proposed EMD-AdaptiveP method. 3.3.1. Comparison with traditional filtering methods In the experiment, the GF-7 proto-waveform was obtained through the laser device for GF-7 laser altimetry. The sample length of the GF-7 proto-waveform is 400, and the intensity is the counts value as shown in Fig. 10. In the EMD-AdaptiveP method, 20 gates of the beginning and ending parts were taken to calculate noise mean and standard deviation. Fig. 10(a)‒(d) displays the denoised effect of these methods. Similarly, each method denoised the waveform effectively, but Gaussian and μ\λ filters led to the decrease of waveform amplitude (as marked by the blue circles), compared with the wavelet filter and EMD-AdaptiveP approach. The corresponding assessment indices were listed in Table 4. Results show that EMD-AdaptiveP has the highest SNR, with 32.06322 dB, and a 9
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 9. Denoised effects of ICESat/GLAS waveform with different EMD-based methods (the blue one is the original noisy full waveform, the orange one is the denoised waveform).
Table 3 Quantitative evaluation of denoised effects with different EMD-based filters. Items
MSE
MAE
SNR
PSNR
R2
EMD-2IMFs EMD-hard EMD-soft EMD-AdaptiveP
0.00175 0.00020 0.00044 0.00011
0.02575 0.01083 0.01584 0.00787
8.10677 17.62631 14.11089 20.28533
21.41105 30.93058 27.41516 33.58960
0.90367 0.98955 0.97635 0.99435
PSNR of 43.951, as well as the smallest MSE of 3.9193 and a MAE of 1.95796, and the biggest R2 of 0.99894. These findings demonstrate that this method offers more improvement over the others. 3.3.2. Comparison with EMD-based filtering methods For the same GF-7 proto-waveform, the EMD–based filtering methods were carried out. In the EMD-AdaptiveP method, the first IMF was removed. The denoised waveform was the same as with the EMD-1IMF method. The rest of the results are displayed in Fig. 11(a)‒(d). The data show that the EMD-2IMFs clearly broadened the original waveform and reduced the peak, as marked in Fig. 11(a). EMD-hard occurred when there was a sudden change (as marked by the red rectangle) at 200 ns and 250 ns, and EMD-soft deviated slightly from the original waveform, near 190 ns and 240 ns, as seen in Fig. 11(b) and (c). The EMD-AdaptiveP performed more effectively than the above methods. The quantitative results in Table 5 demonstrate that this method significantly enhanced the SNR of the original waveform. 4. Conclusion The full waveform of space-borne laser altimetry can provide additional geometric and physical information on the scattering substance along the path, which plays an important role in many fields such as atmosphere, ice, water, forestry, land elevation and so on. However, the waveform of the echo signal demonstrates that it is both nonlinear and nonstationary and may be submerged within various degrees of noises from the process of transmission and acquisition. Noise reduction is a key step in ensuring the retrieval accuracy of data. To effectively reduce the various noises in the full waveform, and to provide a reference for the processing of future GF-7 full 10
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 10. Comparison of denoised GF-7 proto-waveform by different filtering methods (the red one is the original noisy full waveform, the green one is the denoised waveform).
Table 4 Quantitative evaluation of denoising effects of traditional and EMD-AdaptiveP filters. Items
MSE
MAE
SNR
PSNR
R2
Gaussian μ\λ Wavelet EMD-AdaptiveP
8.06063 6.07467 4.47712 3.91930
2.21114 1.79445 1.95882 1.95796
28.93162 30.16009 31.48533 32.06322
40.81940 42.04786 43.37311 43.95100
0.99815 0.99884 0.99874 0.99889
waveform data, a novel noise reduction method named EMD-AdaptiveP was proposed in this paper, which has the adaptive ability to distinguish the effective signal IMFs from all the reconstructed noise IMFs, and to select the threshold for the construction of a full waveform signal, according to different noise grades. This method combines EMD with the strategy of wavelet-like signals to adapt to the characteristics of each waveform. The reconstruction of an effective waveform signal part is implemented through reverse superimposition of its IMFs and its residual, with a threshold of three-sigma principle of background noise for the signal and the additional check mechanism. Some experiments with the simulation waveform, the GLAS waveform, and the GF-7 proto-waveform were carried out to test and verify the effectiveness of this method. It was compared with other different IMF selection schemes, like EMD-soft, EMD-hard, EMD-Wavelet and so on and the traditional filtering methods of Gaussian, μ\λ and wavelet. The results show that the proposed EMD-AdaptiveP method (1) has the adaptive ability to denoise the waveform signal with different levels of SNR; (2) proves itself to be more suitable for distinguishing the effective signal IMFs from all the reconstructed noise IMFs; and (3) performs more effectively on noise reduction compared to the traditional methods of Gaussian, μ\λ and wavelet.
Acknowledgments This work was supported by the Project of Major Special Research Projects of the High Resolution Earth Observation System (Civil Part) (Project No. 06-Y20A17-9001-17/18), a Project funded by China Postdoctoral Science Foundation (Project No. 2018M632172) and a Capacity Building Project of Local Colleges and Universities of Shanghai Science and Technology Commission (Project No. 19050502100). The authors would like to thank the National Snow and Ice Data Center for providing ICESat data.
11
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 11. Denoised effects of GF-7 proto-waveform filtered by different EMD-based methods (the blue one is the original noisy waveform, the orange one is the denoised waveform). Table 5 Quantitative evaluation of denoising effects of different EMD-based filters. Items
MSE
MAE
SNR
PSNR
R2
EMD-2IMFs EMD-hard EMD-soft EMD-AdaptiveP
199.44777 4.59385 16.14926 3.91930
6.13336 2.02139 2.93238 1.95796
14.99702 31.37354 25.91379 32.06322
26.88480 43.26132 37.80157 43.95100
0.94428 0.99870 0.99577 0.99889
References [1] X. Wang, X. Cheng, P. Gong, H. Huang, Z. Li, X. Li, Earth science applications of ICESat/GLAS: a review, Int. J. Remote Sens. 32 (2011) 8837–8864. [2] S. Luthcke, C. Carabajal, D. Rowlands, D. Pavlis, Improvements in spaceborne laser altimeter data geolocation, Surv. Geophys. 22 (2001) 549–559. [3] H.J. Zwally, B. Schutz, W. Abdalati, J. Abshire, C. Bentley, A. Brenner, J. Bufton, J. Dezio, D. Hancock, D. Harding, ICESat’s laser measurements of polar ice, atmosphere, ocean, and land, J. Geodyn. 34 (2002) 405–445. [4] K.D. Fieber, I.J. Davenport, J.M. Ferryman, R.J. Gurney, J.P. Walker, J.M. Hacker, Analysis of full-waveform LiDAR data for classification of an orange orchard scene, ISPRS J. Photogramm. Remote. Sens. 82 (2013) 63–82. [5] H.D. Pritchard, S.R.M. Ligtenberg, H.A. Fricker, D.G. Vaughan, M.R.V.D. Broeke, L. Padman, Antarctic ice-sheet loss driven by basal melting of ice shelves, Nature 484 (2012) 502. [6] E. Khalefa, I.P. Smit, A. Nickless, S. Archibald, A. Comber, H. Balzter, Retrieval of savanna vegetation canopy height from ICESat-GLAS spaceborne LiDAR with terrain correction, IEEE Geosci. Remote. Sens. Lett. 10 (2013) 1439–1443. [7] Y. Cheng, J. Cao, Q. Hao, Y. Xiao, F. Zhang, W. Xia, K. Zhang, H. Yu, A novel de-noising method for improving the performance of full-waveform LiDAR using differential optical path, Remote Sens. 9 (2017) 1109. [8] F. Pirotti, Analysis of full-waveform LiDAR data for forestry applications: a review of investigations and methods, iForest-Biogeosci. For. 4 (2011) 100. [9] C.S. Gardner, Ranging performance of satellite laser altimeters, IEEE Trans. Geosci. Remote. Sens. 30 (1992) 1061–1072. [10] J.L. Bufton, Laser altimetry measurements from aircraft and spacecraft, Proceedings of the IEEE, (1989), pp. 463–477. [11] H. Duong, R. Lindenbergh, N. Pfeifer, G. Vosselman, ICESat full-waveform altimetry compared to airborne laser scanning altimetry over the Netherlands, IEEE Trans. Geosci. Remote. Sens. 47 (2009) 3365–3378. [12] A. Brenner, Derivation of Range and Range Distributions From Laser Pulse Waveform Analysis for Surface Elevations, Roughness, Slope, and Vegetation Heights. Algorithm Theoretical Basis Document V4. 1, http://www. csr. utexas. edu/glas/pdf/Atbd_20031224. pdf (2003). [13] H. Azami, K. Mohammadi, B. Bozorgtabar, An improved signal segmentation using moving average and Savitzky-Golay filter, J. Signal Inf. Process. 3 (2012) 39. [14] G. Taubin, A signal processing approach to fair surface design, Conference on Computer Graphics & Interactive Techniques (1995). [15] S.-L. Yu, J.-C. Gu, Removal of decaying DC in current and voltage signals using a modified Fourier filter algorithm, IEEE Trans. Power Deliv. 16 (2001) 372–379. [16] C.E. Parrish, Exploiting full-waveform lidar data and multiresolution wavelet analysis for vertical object detection and recognition, IEEE Int. Geosci. Remote
12
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Sens. Symposium (2008). [17] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis Proceedings of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences, (1998), pp. 903–995. [18] J.-T. Yu, S.-Y. Zhao, Q. Wang, De-nosing of acoustic emission signals based on empirical mode decomposition and wavelet transform, Harbin Gongye Daxue Xuebao(Journal of Harbin Institute of Technology) 43 (2011) 88–92. [19] H. Sun, W. Chen, J. Gong, An improved empirical mode decomposition-wavelet algorithm for phonocardiogram signal denoising and its application in the first and second heart sound extraction, 2013 6th International Conference on Biomedical Engineering and Informatics, IEEE (2013) 187–191. [20] L. Luo, Y. Yan, P. Xie, J. Sun, Y. Xu, J. Yuan, Hilbert–Huang transform, Hurst and chaotic analysis based flow regime identification methods for an airlift reactor, Chem. Eng. J. 181 (2012) 570–580. [21] Z. Peng, W.T. Peter, F. Chu, A comparison study of improved Hilbert–huang transform and wavelet transform: application to fault diagnosis for rolling bearing, Mech. Syst. Signal Process. 19 (2005) 974–988. [22] C. Dai, D. Liang, Y. Wang, S.U. Jinshan, The LiDAR echo denoising algorithm based on EMD and FFT, Chin. J. Electron Devices (2015). [23] D. Li, L. Xu, X. Li, L. Ma, Full-waveform LIDAR signal filtering based on empirical mode decomposition method, 2013 IEEE International Geoscience and Remote Sensing Symposium-IGARSS, IEEE (2013) 3399–3402. [24] A.-O. Boudraa, J.-C. Cexus, Denoising via empirical mode decomposition, Proc. IEEE ISCCSP 4 (2006) (2006). [25] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inf. Theory 41 (1995) 613–627. [26] M. Zhou, C. Li, L. Ma, H. Guan, Land cover classification from full-waveform Lidar data based on support vector machines, international archives of the photogrammetry, Remote Sens. Spatial Inf. Sci. 41 (2016). [27] Y. Lei, J. Lin, Z. He, M.J. Zuo, A review on empirical mode decomposition in fault diagnosis of rotating machinery, Mech. Syst. Signal Process. 35 (2013) 108–126. [28] Q. Xie, B. Xuan, S. Peng, J. Li, W. Xu, H. Han, Bandwidth empirical mode decomposition and its application, International Journal of Wavelets, Multiresolution Inf. Process. 6 (2008) 777–798. [29] X. Lai, M. Zheng, A denoising method for LiDAR full-waveform data, Math. Probl. Eng. 2015 (2015). [30] J. Chang, L. Zhu, H. Li, F. Xu, B. Liu, Z. Yang, Noise reduction in Lidar signal using correlation-based EMD combined with soft thresholding and roughness penalty, Opt. Commun. 407 (2018) 290–295.
13
Contents lists available at ScienceDirect
Optik journal homepage: www.elsevier.com/locate/ijleo
Original research article
A novel noise reduction method for space-borne full waveforms based on empirical mode decomposition Zhijie Zhanga, Xiangfeng Liub,*, Rong Shub, Feng Xieb, Fengxiang Wangb, Zhihui Liub, Hanwei Zhanga, Zhenhua Wangc a
School of Surveying and Land Information Engineering, Henan Polytechnic University, Jiaozuo, China Shanghai Institute of Technical Physics of The Chinese Academy of Sciences, Shanghai, China c College of Information, Shanghai Ocean University, Shanghai, China b
A R T IC LE I N F O
ABS TRA CT
Keywords: Full waveform Noise reduction Empirical mode decomposition Laser altimetry GaoFen-7
In this project, a novel noise reduction method based on empirical mode decomposition (EMD), referred to as EMD-AdaptiveP, is proposed and investigated with the expectation of providing an effective enhancement for the full waveform of nonlinear and nonstationary signals obtained by space-borne laser altimetry and the forthcoming GaoFen-7 (GF-7) mapping satellite in China. This method combines EMD with a similar strategy of wavelet filtering to adapt to the characteristics of each waveform. The reconstruction of an effective waveform signal is implemented through reverse superimposition of its intrinsic mode functions (IMFs) and the residual, which is thresholded by the three-sigma principle of background noise of the waveform signal and an additional check mechanism in the lower or higher signal of noise levels, after the waveform is decomposed by EMD. This method is qualitatively and quantitatively examined with simulated waveform data, real ice, cloud and land elevation satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) waveforms and GF-7 proto-waveform data, in contrast with other IMF selection schemes, like EMD-soft, EMD-hard and EMD-Wavelet, or the traditional filtering methods, such as Gaussian, μ\λ and wavelet. The results show that the method of EMD-AdaptiveP (1) is sufficiently robust to denoise the waveform signal with different levels of signal-to-noise ratio, (2) has the adaptive ability to distinguish effective signal IMFs from noise IMFs in signal reconstruction, and (3) performs noise reduction more effectively than traditional methods.
1. Introduction Space-borne laser altimetry is an active remote sensing technique with high temporal‒spatial resolution and a high accuracy range, which can provide rapid acquisition of geospatial information in large regions and at high-precision benchmarks on a global scale, even in special areas (e.g., polar regions) [1]. Since the ice, cloud and land elevation satellite (ICESat)/Geoscience Laser Altimeter System (GLAS) for Earth observation launched in 2003 [2], full waveform data have been widely used for the observation of polar ice sheet change; vertical distribution of vegetation, clouds and aerosols; estimation of biomass; sea surface height; and monitoring of climatic phenomena [1], water and land elevation [3–5]. GaoFen-7 (GF-7), the first high-resolution stereo mapping satellite in China, equipped with a double-line array camera and a laser altimeter, will be launched in 2019. Compared with discrete return systems, the full waveform, with the characteristics of waveform shapes, widths, intensities and skewness for the laser echoes,
⁎
Corresponding author. E-mail addresses: [email protected], [email protected] (X. Liu).
https://doi.org/10.1016/j.ijleo.2019.163581 Received 17 July 2019; Accepted 10 October 2019 0030-4026/ © 2019 Elsevier GmbH. All rights reserved.
Please cite this article as: Zhijie Zhang, et al., Optik - International Journal for Light and Electron Optics, https://doi.org/10.1016/j.ijleo.2019.163581
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
not only incorporates responses from ground level but also provides additional geometric and physical information about multiplescattering responses along the path [6]. However, full waveform data are always contaminated by various noises in the process of transmission and acquisition, such as ray radiation from the sun and sky; dark current noise; and thermal noise from the optical detection system [7]. Furthermore, the echo signal is very weak in its laser intensity with the increase of the detection range [8]. That would be the major impediment to full waveform decomposition and extraction of ground feature parameters (e.g., ranges, slope and roughness) [9,10]. The waveform of echo signals demonstrates both nonlinear and nonstationary characteristics, and it may be submerged in the noises. Therefore, noise reduction of the full waveform is an indispensable key step in improving the retrieval accuracy of the data and ensuring the reliability of elevation-related applications [11]. At present, fixed-width Gaussian filters [12], average filters [13], μ\λ filtering [14], Fourier low-pass filtering [15] and wavelet filtering [16] have been extensively used to denoise the full waveform. However, some drawbacks come with these methods; fixedwidth Gaussian filters, average filtering and μ\λ filtering face interference by improper filtering parameters, which lead to waveform distortion, such as peak amplitude shrinkage or pulse width increase. Fourier filtering only demonstrates the signals in a frequency domain or a time domain; the unsuitable cut-off frequency may cause signal distortion. Although wavelet filtering shows frequency and energy content in the time domain, it needs to choose a priori basis function with the original signal, based on experience, and a preset wavelet basis function cannot adapt to all kinds of waveforms. In short, the common problem with these methods is that they need to suitably select the filtering parameters. Inappropriate parameters will cause distortion of waveforms [7]. However, empirical mode decomposition (EMD), devised by Huang et al. [17], can filter and extract the trend term with nonlinear and nonstationary signals, which is a fully signal-dependent approach based on the local characteristic time scale that adaptively decomposes a given signal into several intrinsic mode functions (IMFs) and a residual. Researchers were already familiar with this practice from various forms of signal filtering, such as acoustic signals [18], phonocardiograms [19], pressure signals of loop airlift reactors [20], signals for fault diagnosis of rolling bearings [21] and signal filters attempted in the enhancement of waveform data [22]. The basic principle of these methods is to reconstruct the signal within previously filtered IMFs by direct removal of the first several IMFs with high frequency [23], combination with the traditional filters [24], thresholding for preset threshold functions [25] or fast Fourier transformation (FFT) [22] toward all IMFs to select the high frequency IMFs. The common problem in the above EMD-based filters is that there is no easy way to adaptively distinguish between the noise IMFs and effective signal IMFs and denoise the signal without losing useful information; this is still an open issue [26]. Therefore, to effectively reduce the various noises in full waveform and apply that reduction technique to the future GF-7 full waveform, a novel noise reduction method based on EMD, named as EMD-AdaptiveP, is proposed in this paper. It can adaptively distinguish effective signal IMFs from noise IMFs and select the threshold for the construction of a full waveform with different levels of signal-to-noise ratio (SNR). To test and verify the effectiveness of the proposed method in a full waveform signal, some experiments with simulation waveform, GLAS waveform and the GF-7 proto-waveform were carried out. Meanwhile, the EMD-AdaptiveP was compared to other different IMF selection schemes, such as EMD-soft, EMD-hard, EMD-Wavelet and so on. Furthermore, comparisons of this method with the traditional methods of Gaussian, μ\λ and wavelet were performed with real waveform data. The rest of this paper is organised as follows: In Section 2, we describe the noise reduction method based on EMD. In Section 3, this method is examined along with the simulation waveform and real waveform from laser altimetry, comparing it to other IMF selection and traditional noise reduction methods. Finally, in Section 4 the conclusions are drawn from the data gathered.
2. Proposed method For the full waveform, the noises are regarded as white Gaussian noise [10], including quantum noise, thermal noise, amplifier noise, dark current noise, background noise and so on to divide the IMFs of the original full waveforms into noise IMFs and effective signal IMFs and determine which IMFs contain most of the effective signal components, beginning from a coarse, and progressing to a fine scale. The reconstruction of the denoised signal is obtained by summing up the effective IMFs reversely until the denoised full waveform meets the preset three-sigma threshold. The reasonability and validity of this idea is discussed as regards the following aspects: the completeness of full waveforms after EMD decomposition, the estimation of involved noise parameters to set the threshold of stopping signal reconstruction and the conception of the overall procedure.
2.1. The completeness of full waveform EMD decomposition Compared to wavelet filtering, EMD is a fully data-driven method and does not necessitate any basis function as a prior, based on the sequential extraction of energy associated with various intrinsic time scales of the signal, starting from the finer temporal scales (high frequency modes) and moving toward coarser ones (low frequency modes). Each mode should possess an equal number of extrema and zero crossings. The detail of the EMD algorithm was described in references [27,28]. To validate the completeness of EMD in space-borne full waveform processing, a typical GLAS waveform was used. There are six IMFs decomposed by EMD, as shown in Fig. 1, and the first IMF has the highest frequency component, while the frequency decreases with the increment of the IMFs’ order. Fig. 2 shows the original full waveform and reconstruction error between the sum of all the IMFs and the residual, which displays the decomposition of EMD thoroughly. The reconstruction error demonstrates that the total sum of the IMFs and the residual matches well with the original signal and that the order of amplitude reached 10−17, which can be attributed to the rounding error of computers. The completeness lays out the foundation of the following IMFs’ selection analysis and waveform reconstruction. 2
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 1. IMFs of a typical GLAS full waveform decomposed by EMD.
Fig. 2. A typical GLAS full waveform and reconstruction error concerning all the IMFs and the residual.
2.2. Estimation of noise mean and standard deviation Determining how to design the threshold that would stop iterations of IMFs is key to properly filtering noise signals and preserving the characteristic waveform signal. Usually, the value of noise is unknown, while the background noise mixed in the returned waveform signal has a relatively stable standard deviation. The beginning and ending parts of a waveform signal are generally considered stable background noise [29]. In this paper, we adopt this method to estimate the mean value and the standard deviation of noise, calculated as Eqs. (1) and (2). N
μnoise =
σnoise =
∑i =1
Vi (1)
N
N
∑i =1
(Vi − μN )2 N−1
(2) 3
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 3. The workflow for the denoising strategy (checking constraints in the case of the lower and higher SNR levels shown in the red dotted frame).
Where Vi is the amplitude of the i th sample, N is the number of samples and μN is the mean value for the sampled noise. In practice, echo components for the unknown detecting range may appear at the beginning or ending of the waveform. These may lead to an overestimated noise level with a higher estimated mean and standard deviation. Therefore, we take the smaller value from the beginning noise ( μbeg ) and the ending noise ( μend ) and even take the corresponding standard deviation (σbeg , σend ) to be the noise threshold in the event that we miss the effective signal 2.3. Designing of the proposed filter for full waveform signal After a waveform is decomposed by EMD, the adaptive threshold for waveform reconstruction from the reverse superposition of its IMFs and the residual is based on the root mean squared error (RMSE) between the original waveform amplitude and the denoised waveform amplitude, which should be three times the standard deviation over the mean background noise. That means the differing amplitude between the pre- and post-denoising waveforms should fall within the requirements for categorisation as noise. Given that most of the significant structures of the signal are concentrated in the lower frequency modes (last IMFs) and decrease toward the high-frequency modes (first IMFs), the strategy behind this filtering method is to partially reconstruct the effective parts of the decomposed signal from the last several IMFs, up to the penultimate IMF, where the number P of the last IMF is adaptive to the noise levels of each waveform. Thus the method known as EMD-AdaptiveP differs from other methods by directly removing the first few high-frequency IMFs, whose collective denoising effects depend on the setting of their number, or by using soft and hard thresholds to compensate for the effective signals that may be removed by directly discarding the high-frequency IMFs. Fig. 3 shows the entire procedure of this strategy. In the experiments with simulated waveforms in the next section, when the noise is at medium level, the basic noise threshold can easily be met with better noise suppression effects. However, when the noise is 4
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
lower (e.g., 10db or 15db), results show that the reconstructed signal by summing the last one or two IMFs can also meet the threshold with a severe loss of effective signal parts. In this case, it is necessary to add a check mechanism for the waveforms with lower SNR levels. We propose that the difference between the maximum peak of the original waveform and the smoothed waveform should also be below the noise threshold. It would also need to be bigger than the noise mean plus 1.5 times the standard deviation to avoid an insufficient smoothing effect. For the high SNR level, all the IMFs are summed, and none of the IMFs that cannot denoise the full waveform are discarded. Therefore, the other constraint, stipulating that the extreme order (k = maxIMFs–1) of IMFs is the penultimate IMF, is applied. 3. Experimental results and analysis To test and verify the effectiveness of the proposed novel filter method, experiments using simulation waveforms, GLAS waveforms and GF-7 proto-waveforms are conducted in this section. Furthermore, the qualitative and quantitative comparisons and analysis among other EMD-based methods and traditional filtering methods are carried out. 3.1. Simulation waveform Referencing the GLAS samples of land type waveforms and the parameters of GF-7, the simulation waveforms were simulated and tested. For the transmitted pulse with the form of the Gaussian model, the mixed Gauss model was adopted when simulating the pure received waveform signal. In the experiment, random coarse numbers (1, 10), (6, 40) were set as the amplitude value and the peak width, respectively; a random integer (100, 500) was set as the peak location. The maximum number of peaks was three. Then the random (normally distributed) white Gaussian noise was superimposed on the pure signals. To simulate the effect of digitiser sampling, the return samples were rounded to integer amplitude values, with a length of 544 ns. 3.1.1. Full waveform reconstruction under different SNR levels To validate the effectiveness of this method, two groups of waveforms with fixed peak numbers under different SNR levels and varied peak numbers with a fixed SNR level denoting the complexity of the waveforms were tested. Fig. 4 and 5, respectively, show the recovery processes of the signals by the residuals and their IMFs. In Fig. 4(a)–(d), the simulated full waveforms have three peak numbers under different SNR levels of 10db, 20db, 30db and 40db. The results show that with the increase of SNR, the number of noise IMFs decreases. The discarded first few IMFs are 3 (Fig. 4(a)), 2 (Fig. 4(b)), 1 (Fig. 4(c)) and 1 (Fig. 4(d)), respectively. In Fig. 5(a)‒(d), the simulated full waveforms have varied peak numbers of 1, 2, 3 and 4 under the SNR level of 15db. The results show that with the increase of the peak numbers, the numbers of noise IMFs vary. The discarded first few IMFs are 3 (Fig. 5(a)), 2 (Fig. 5(b)), 2 (Fig. 5(c)) and 2 (Fig. 5(d)) respectively, which all maintain the peak shape well, even with two very close peaks (Fig. 5(b)). These demonstrate that the proposed EMD-AdaptiveP method can adaptively reconstruct the varying waveform signals with different SNR levels and peak numbers while preserving the detailed waveform shape effectively. 3.1.2. Evaluation and comparison of denoising performance To evaluate the performance of this method, the EMD-AdaptiveP was compared with five other IMF selection schemes, using a typical simulated full waveform by directly removing either the first or the first two high-frequency IMF methods (EMD-1IMF or EMD-2IMFs); to avoid the loss of the effective waveform signal, after the direct removal of the first two IMFs, the universal threshold (either EMD-soft or EMD-hard) was used to extract the effective parts, which were then added to the reconstructed signal; the EMDWavelet was then performing wavelet filtering of each IMF decomposed by EMD. In this experiment, depending on the characteristics of the waveform, DN4 was chosen as the wavelet basis, and the soft threshold method was used to reconstruct the waveform signal. The denoising effect was as shown in Fig. 6(a)‒(f). The results obviously illustrate that EMD-1IMF has an inadequate denoising effect. There are still lots of high-frequency noises in denoised waveforms. The denoising effects of the remaining four comparison methods are similar, except in the case of EMD-Wavelet (Fig. 6(e)), which caused the decline of the peak value. The proposed EMD-AdaptiveP smoothed the waveform nicely. To validate the effectiveness of this method and avoid the resulting randomisation for each individual waveform, 2000 waveforms with different SNR levels were examined in this experiment. The quantitative assessment indicator was the consistency ratio of waveform decomposition from the Levenberg-Marquardt method. The consistency ratio is defined by two constraints: 1) the fitting number of the Gaussian components is equal to that of the simulated ones; 2) the average error of peak position is less than 1 ns. The statistics of consistency ratio for the experiments were determined with the maximum, minimum, median, mean and standard deviation (Std.), as shown in Fig. 7 and Table 1. Table 1 shows that the EMD-AdaptiveP has the highest consistency ratio of 0.7045 with a Std. of 0.0715. From the perspective of mean value, the consistency ratio of EMD-AdaptiveP has improved by 33.21%, 5.44%, 6.27%, 5.18% and 4.56%, compared with the EMD-1IMF, EMD-2IMFs, EMD-soft, EMD-hard and EMD-Wavelet methods, respectively. Looking at the trend of the performances as shown in Fig. 7, the EMD-AdaptiveP has a much better effect than do the other methods, especially when the SNR levels are lower. The consistency shows a downward trend when the waveforms have a higher SNR (about 35db). The reason may be that there was insufficient denoising with the removal of one IMF or that removing multiple IMFs led to the loss of effective signals. However, the maximum mean of 0.6394 illustrates the strong robustness of this approach under varying SNR levels. 5
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 4. Reconstruction of simulated waveform with fixed peak numbers under different SNR levels of 10db, 20db, 30db and 40db, respectively, by reverse IMF superimposition (the blue solid line is the original noisy waveform; the red solid line is the sum of the residual and effective signal IMFs’ components).
3.2. ICESat/GLAS waveform This method was tested with real ICESat/GLAS waveforms and then compared with traditional filtering methods, such as Gaussian filtering, μ\λ filtering and wavelet filtering. In the experiment, some quantitative indices [29,30] such as mean squared error, (MSE), mean absolute error (MAE) and cross correlation (R2) between the original noisy waveform and denoised waveform were taken to assess the degree of deviation and similarity between the methods. The output SNR of the signal and the peak signal-to-noise ratio (PSNR) were also taken to assess whether more useful signals were contained in the denoised waveforms. These indices can be calculated with the following Eqs. (3)‒(7):
MSE =
1 N
N
MAE =
N
∑n=1 (Raws (n) − Smooths (n))2
∑n=1
(3)
|Raws (n) − Smooths (n)| (4)
N 6
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 5. Reconstruction of simulated full waveform, with varied peak numbers of 1, 2, 3 and 4, respectively under SNR, equals 15 after reverse IMF superimposition (the blue solid line is the original noisy waveform; the red solid line is the sum of the residual and effective signal IMFs’ components). N
⎛ ⎞ ∑n=1 |Raws (n)|2 SNR = 10*log ⎜ N ⎟ ⎜∑ |Raws (n) − Smooths (n)|2 ⎟ ⎝ n=1 ⎠ N
PSNR = 10*log
∑n=1 N
∑n=1
max 2
(Raws (n) − Smooths (n))2
(6) 2
N
R2 =
(5)
(Raws (n) − Raws (n))(Smooths − Smooths ) ⎤ ⎡∑ ⎣ n=1 ⎦ N
∑n=1
N
(Raws (n) − Raws (n))2∑
n=1
(Smooths − Smooths )2
(7)
Where Raws represents the original noisy waveform signal; Smooths represents the denoised waveform signal; and N is the length of the waveform sample 3.2.1. Comparison with traditional filtering methods For Gaussian filtering, the width 6 ns of the GLAS transmit pulse was used as the Gaussian kernel, and 13 bins were adopted as the filtering length; for the μ\λ filtering, the values of λ and μ were 0.6307 and -0.6372, respectively; for the wavelet filtering, the “Daubechies 4” basis with the level 4 composition and the soft threshold method was employed; for the EMD-AdaptiveP method, 50 samples of beginning and ending parts of waveforms were used to estimate the noise mean and standard deviation. In order to show 7
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 6. Denoised full waveform obtained from five IMF selection schemes: (a) EMD-1IMF, (b) EMD-2IMFs, (c) EMD-hard, (d) EMD-soft, (e) EMDWavelet, (f) EMD-AdaptiveP (the blue one is the noisy full waveform, the red one is the denoised waveform).
Fig. 7. Consistency ratio variation of simulated full waveform along with the increase of SNR levels. Table 1 The statistics of consistency ratio of different IMF selection schemes. Items
EMD-1IMF
EMD-2IMFs
EMD-soft
EMD-hard
EMD-Wavelet
EMD-AdaptiveP
Maximum Minimum Median Mean Std.
0.5145 0.3555 0.4975 0.4800 0.0568
0.6485 0.4320 0.6405 0.6064 0.0788
0.6490 0.4535 0.6235 0.6017 0.0671
0.6515 0.4290 0.6440 0.6079 0.0805
0.6750 0.4340 0.6380 0.6090 1.2202
0.7045 0.4985 0.6685 0.6394 0.0715
the denoising effect of these methods more clearly, only the parts of the effective waveform signal based on about 200 samples are displayed (the original waveform samples numbered 544) in Fig. 8(a)‒(d). The results show that each method has effectively denoised the waveform signal. However, there is an obvious peak decrease in the denoised signal conducted by Gaussian and μ\λ filtering, as marked with a blue circle in Fig. 8(a) and (b). In addition, there is a spike noise accompanying the 130 ns treated by wavelet filtering, as marked with a blue circle in Fig. 8(c). EMD-AdaptiveP, as shown in Fig. 8(d), not only smooths the waveform but also maintains the main characteristics of the waveform. According to the statistical results listed in Table 2, EMD-AdaptiveP has the highest output SNR (20.285325 dB) and PSNR (33.589603); the smallest MSE (0.000106) and MAE (0.007867); and the biggest R2 (0.99435), which demonstrates that the method can greatly enhance the output SNR and PSNR, outdoing traditional filters. 3.2.2. Comparison with EMD-based filtering methods In addition, there were some qualitative and quantitative analyses conducted with different EMD-based filtering methods for the same GLAS waveform. In this waveform, the first IMF was removed by EMD-AdaptiveP, which is the same as EMD-1IMF. Therefore, EMD-1IMF was omitted in the comparisons. Fig. 9(a)‒(d) illustrates the denoised waveforms conducted by EMD-2IMFs, EMD-hard, EMD-soft and EMD-AdaptiveP. The results clearly show that many useful signal components were discarded by EMD-2IMFs, as marked by the red rectangle boxes in Fig. 9(a). There were sudden changes or distortions due to EMD-hard or EMD-soft thresholding, 8
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 8. Denoised ICESat/GLAS waveform by different filtering methods (the red one is the original noisy waveform, the green one is the denoised waveform). Table 2 Quantitative evaluation of denoised effects by the traditional methods and by EMD-AdaptiveP. Items
MSE
MAE
SNR
PSNR
R2
Gaussian μ\λ Wavelet EMD-AdaptiveP
0.00033 0.00029 0.00078 0.00011
0.01158 0.01094 0.01739 0.00787
15.35926 15.93050 11.60523 20.28533
28.66354 29.23478 24.90950 33.58960
0.98458 0.98627 0.95724 0.99435
as marked by the red rectangle boxes in Fig. 9 (b) and (c). However, the EMD-AdaptiveP method performed more effectively than these methods, and the statistical results listed in Table 3 further demonstrate that EMD-AdaptiveP has the biggest output PSNR (33.58960) and the highest correlation (0.99435) with the original waveform, which gives it advantages over other EMD-based methods. 3.3. GF-7 proto-waveform Because the GF-7 satellite has not yet been launched, an echo waveform acquired by laser altimetry was developed in a simple surface test environment. Thus, this waveform is referred to as a proto-waveform, which may be different from a waveform in a real in-orbit environment. However, it does not affect the validation of the filtering method. As with the processing of ICESat/GLAS waveforms, when denoising the proto-waveform, different traditional filters and EMD-based denoising methods were used to test the effectiveness of the proposed EMD-AdaptiveP method. 3.3.1. Comparison with traditional filtering methods In the experiment, the GF-7 proto-waveform was obtained through the laser device for GF-7 laser altimetry. The sample length of the GF-7 proto-waveform is 400, and the intensity is the counts value as shown in Fig. 10. In the EMD-AdaptiveP method, 20 gates of the beginning and ending parts were taken to calculate noise mean and standard deviation. Fig. 10(a)‒(d) displays the denoised effect of these methods. Similarly, each method denoised the waveform effectively, but Gaussian and μ\λ filters led to the decrease of waveform amplitude (as marked by the blue circles), compared with the wavelet filter and EMD-AdaptiveP approach. The corresponding assessment indices were listed in Table 4. Results show that EMD-AdaptiveP has the highest SNR, with 32.06322 dB, and a 9
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 9. Denoised effects of ICESat/GLAS waveform with different EMD-based methods (the blue one is the original noisy full waveform, the orange one is the denoised waveform).
Table 3 Quantitative evaluation of denoised effects with different EMD-based filters. Items
MSE
MAE
SNR
PSNR
R2
EMD-2IMFs EMD-hard EMD-soft EMD-AdaptiveP
0.00175 0.00020 0.00044 0.00011
0.02575 0.01083 0.01584 0.00787
8.10677 17.62631 14.11089 20.28533
21.41105 30.93058 27.41516 33.58960
0.90367 0.98955 0.97635 0.99435
PSNR of 43.951, as well as the smallest MSE of 3.9193 and a MAE of 1.95796, and the biggest R2 of 0.99894. These findings demonstrate that this method offers more improvement over the others. 3.3.2. Comparison with EMD-based filtering methods For the same GF-7 proto-waveform, the EMD–based filtering methods were carried out. In the EMD-AdaptiveP method, the first IMF was removed. The denoised waveform was the same as with the EMD-1IMF method. The rest of the results are displayed in Fig. 11(a)‒(d). The data show that the EMD-2IMFs clearly broadened the original waveform and reduced the peak, as marked in Fig. 11(a). EMD-hard occurred when there was a sudden change (as marked by the red rectangle) at 200 ns and 250 ns, and EMD-soft deviated slightly from the original waveform, near 190 ns and 240 ns, as seen in Fig. 11(b) and (c). The EMD-AdaptiveP performed more effectively than the above methods. The quantitative results in Table 5 demonstrate that this method significantly enhanced the SNR of the original waveform. 4. Conclusion The full waveform of space-borne laser altimetry can provide additional geometric and physical information on the scattering substance along the path, which plays an important role in many fields such as atmosphere, ice, water, forestry, land elevation and so on. However, the waveform of the echo signal demonstrates that it is both nonlinear and nonstationary and may be submerged within various degrees of noises from the process of transmission and acquisition. Noise reduction is a key step in ensuring the retrieval accuracy of data. To effectively reduce the various noises in the full waveform, and to provide a reference for the processing of future GF-7 full 10
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 10. Comparison of denoised GF-7 proto-waveform by different filtering methods (the red one is the original noisy full waveform, the green one is the denoised waveform).
Table 4 Quantitative evaluation of denoising effects of traditional and EMD-AdaptiveP filters. Items
MSE
MAE
SNR
PSNR
R2
Gaussian μ\λ Wavelet EMD-AdaptiveP
8.06063 6.07467 4.47712 3.91930
2.21114 1.79445 1.95882 1.95796
28.93162 30.16009 31.48533 32.06322
40.81940 42.04786 43.37311 43.95100
0.99815 0.99884 0.99874 0.99889
waveform data, a novel noise reduction method named EMD-AdaptiveP was proposed in this paper, which has the adaptive ability to distinguish the effective signal IMFs from all the reconstructed noise IMFs, and to select the threshold for the construction of a full waveform signal, according to different noise grades. This method combines EMD with the strategy of wavelet-like signals to adapt to the characteristics of each waveform. The reconstruction of an effective waveform signal part is implemented through reverse superimposition of its IMFs and its residual, with a threshold of three-sigma principle of background noise for the signal and the additional check mechanism. Some experiments with the simulation waveform, the GLAS waveform, and the GF-7 proto-waveform were carried out to test and verify the effectiveness of this method. It was compared with other different IMF selection schemes, like EMD-soft, EMD-hard, EMD-Wavelet and so on and the traditional filtering methods of Gaussian, μ\λ and wavelet. The results show that the proposed EMD-AdaptiveP method (1) has the adaptive ability to denoise the waveform signal with different levels of SNR; (2) proves itself to be more suitable for distinguishing the effective signal IMFs from all the reconstructed noise IMFs; and (3) performs more effectively on noise reduction compared to the traditional methods of Gaussian, μ\λ and wavelet.
Acknowledgments This work was supported by the Project of Major Special Research Projects of the High Resolution Earth Observation System (Civil Part) (Project No. 06-Y20A17-9001-17/18), a Project funded by China Postdoctoral Science Foundation (Project No. 2018M632172) and a Capacity Building Project of Local Colleges and Universities of Shanghai Science and Technology Commission (Project No. 19050502100). The authors would like to thank the National Snow and Ice Data Center for providing ICESat data.
11
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Fig. 11. Denoised effects of GF-7 proto-waveform filtered by different EMD-based methods (the blue one is the original noisy waveform, the orange one is the denoised waveform). Table 5 Quantitative evaluation of denoising effects of different EMD-based filters. Items
MSE
MAE
SNR
PSNR
R2
EMD-2IMFs EMD-hard EMD-soft EMD-AdaptiveP
199.44777 4.59385 16.14926 3.91930
6.13336 2.02139 2.93238 1.95796
14.99702 31.37354 25.91379 32.06322
26.88480 43.26132 37.80157 43.95100
0.94428 0.99870 0.99577 0.99889
References [1] X. Wang, X. Cheng, P. Gong, H. Huang, Z. Li, X. Li, Earth science applications of ICESat/GLAS: a review, Int. J. Remote Sens. 32 (2011) 8837–8864. [2] S. Luthcke, C. Carabajal, D. Rowlands, D. Pavlis, Improvements in spaceborne laser altimeter data geolocation, Surv. Geophys. 22 (2001) 549–559. [3] H.J. Zwally, B. Schutz, W. Abdalati, J. Abshire, C. Bentley, A. Brenner, J. Bufton, J. Dezio, D. Hancock, D. Harding, ICESat’s laser measurements of polar ice, atmosphere, ocean, and land, J. Geodyn. 34 (2002) 405–445. [4] K.D. Fieber, I.J. Davenport, J.M. Ferryman, R.J. Gurney, J.P. Walker, J.M. Hacker, Analysis of full-waveform LiDAR data for classification of an orange orchard scene, ISPRS J. Photogramm. Remote. Sens. 82 (2013) 63–82. [5] H.D. Pritchard, S.R.M. Ligtenberg, H.A. Fricker, D.G. Vaughan, M.R.V.D. Broeke, L. Padman, Antarctic ice-sheet loss driven by basal melting of ice shelves, Nature 484 (2012) 502. [6] E. Khalefa, I.P. Smit, A. Nickless, S. Archibald, A. Comber, H. Balzter, Retrieval of savanna vegetation canopy height from ICESat-GLAS spaceborne LiDAR with terrain correction, IEEE Geosci. Remote. Sens. Lett. 10 (2013) 1439–1443. [7] Y. Cheng, J. Cao, Q. Hao, Y. Xiao, F. Zhang, W. Xia, K. Zhang, H. Yu, A novel de-noising method for improving the performance of full-waveform LiDAR using differential optical path, Remote Sens. 9 (2017) 1109. [8] F. Pirotti, Analysis of full-waveform LiDAR data for forestry applications: a review of investigations and methods, iForest-Biogeosci. For. 4 (2011) 100. [9] C.S. Gardner, Ranging performance of satellite laser altimeters, IEEE Trans. Geosci. Remote. Sens. 30 (1992) 1061–1072. [10] J.L. Bufton, Laser altimetry measurements from aircraft and spacecraft, Proceedings of the IEEE, (1989), pp. 463–477. [11] H. Duong, R. Lindenbergh, N. Pfeifer, G. Vosselman, ICESat full-waveform altimetry compared to airborne laser scanning altimetry over the Netherlands, IEEE Trans. Geosci. Remote. Sens. 47 (2009) 3365–3378. [12] A. Brenner, Derivation of Range and Range Distributions From Laser Pulse Waveform Analysis for Surface Elevations, Roughness, Slope, and Vegetation Heights. Algorithm Theoretical Basis Document V4. 1, http://www. csr. utexas. edu/glas/pdf/Atbd_20031224. pdf (2003). [13] H. Azami, K. Mohammadi, B. Bozorgtabar, An improved signal segmentation using moving average and Savitzky-Golay filter, J. Signal Inf. Process. 3 (2012) 39. [14] G. Taubin, A signal processing approach to fair surface design, Conference on Computer Graphics & Interactive Techniques (1995). [15] S.-L. Yu, J.-C. Gu, Removal of decaying DC in current and voltage signals using a modified Fourier filter algorithm, IEEE Trans. Power Deliv. 16 (2001) 372–379. [16] C.E. Parrish, Exploiting full-waveform lidar data and multiresolution wavelet analysis for vertical object detection and recognition, IEEE Int. Geosci. Remote
12
Optik - International Journal for Light and Electron Optics xxx (xxxx) xxxx
Z. Zhang, et al.
Sens. Symposium (2008). [17] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis Proceedings of the Royal Society of London, Series A: Mathematical, Physical and Engineering Sciences, (1998), pp. 903–995. [18] J.-T. Yu, S.-Y. Zhao, Q. Wang, De-nosing of acoustic emission signals based on empirical mode decomposition and wavelet transform, Harbin Gongye Daxue Xuebao(Journal of Harbin Institute of Technology) 43 (2011) 88–92. [19] H. Sun, W. Chen, J. Gong, An improved empirical mode decomposition-wavelet algorithm for phonocardiogram signal denoising and its application in the first and second heart sound extraction, 2013 6th International Conference on Biomedical Engineering and Informatics, IEEE (2013) 187–191. [20] L. Luo, Y. Yan, P. Xie, J. Sun, Y. Xu, J. Yuan, Hilbert–Huang transform, Hurst and chaotic analysis based flow regime identification methods for an airlift reactor, Chem. Eng. J. 181 (2012) 570–580. [21] Z. Peng, W.T. Peter, F. Chu, A comparison study of improved Hilbert–huang transform and wavelet transform: application to fault diagnosis for rolling bearing, Mech. Syst. Signal Process. 19 (2005) 974–988. [22] C. Dai, D. Liang, Y. Wang, S.U. Jinshan, The LiDAR echo denoising algorithm based on EMD and FFT, Chin. J. Electron Devices (2015). [23] D. Li, L. Xu, X. Li, L. Ma, Full-waveform LIDAR signal filtering based on empirical mode decomposition method, 2013 IEEE International Geoscience and Remote Sensing Symposium-IGARSS, IEEE (2013) 3399–3402. [24] A.-O. Boudraa, J.-C. Cexus, Denoising via empirical mode decomposition, Proc. IEEE ISCCSP 4 (2006) (2006). [25] D.L. Donoho, De-noising by soft-thresholding, IEEE Trans. Inf. Theory 41 (1995) 613–627. [26] M. Zhou, C. Li, L. Ma, H. Guan, Land cover classification from full-waveform Lidar data based on support vector machines, international archives of the photogrammetry, Remote Sens. Spatial Inf. Sci. 41 (2016). [27] Y. Lei, J. Lin, Z. He, M.J. Zuo, A review on empirical mode decomposition in fault diagnosis of rotating machinery, Mech. Syst. Signal Process. 35 (2013) 108–126. [28] Q. Xie, B. Xuan, S. Peng, J. Li, W. Xu, H. Han, Bandwidth empirical mode decomposition and its application, International Journal of Wavelets, Multiresolution Inf. Process. 6 (2008) 777–798. [29] X. Lai, M. Zheng, A denoising method for LiDAR full-waveform data, Math. Probl. Eng. 2015 (2015). [30] J. Chang, L. Zhu, H. Li, F. Xu, B. Liu, Z. Yang, Noise reduction in Lidar signal using correlation-based EMD combined with soft thresholding and roughness penalty, Opt. Commun. 407 (2018) 290–295.
13