- Email: [email protected]
Enhancement of lidar backscatters signal-to-noise ratio using empirical mode decomposition method
Optics Communications 267 (2006) 137–144 www.elsevier.com/locate/optcom
Enhancement of lidar backscatters signal-to-noise ratio using empirical mode decomposition method Songhua Wu *, Zhishen Liu, Bingyi Liu Ocean Remote Sensing Laboratory of Ministry of Education of China, Ocean Remote Sensing Institute (ORSI), Ocean University of China, No. 5 YuShan Road, Qingdao 266003, China Received 21 September 2005; received in revised form 1 April 2006; accepted 30 May 2006
Abstract Lidar is being widely used to monitor meteorological parameters and atmospheric constituents. Applications include meteorology, environmental pollution, atmospheric dynamics and global climate change. Signal processing for lidar applications involve highly nonlinear models and consequently nonlinear filtering. In this paper, we applied a new method, empirical mode decomposition to the lidar signal processing. The denoising approach is done by removal of the proper intrinsic mode functions. The data from the simulation and measurements are analyzed to evaluate this method comparing with the traditional low-pass filter and the multi-pulse averaging. Results show that it is effective and superior to the band-pass filter and the averaging method. The denoising method also allows less averaging laser shots which is important for the real-time monitoring and for the low cost laser transmitter. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Lidar; Signal-to-noise ratio; Empirical mode decomposition
1. Introduction Lidar is an active remote sensing instrument that emits laser pulses towards the atmosphere or target and measures the backscattering. The signal-to-noise ratio (SNR) of the lidar backscattering often attenuates due to the noise and interferences, such as stochastic turbulences, background noise, dark current, electronics readout noise and atmospheric turbulence. In order to improve the precision of the measurement and subsequent analysis, lidar researchers usually use the multiple-pulse average (temporal average) or running-average (spatial average) approach to smooth the lidar signal [1–3]. The conventional average approach is a kind of lowpass filtering process based on the least square. System errors are assumed to have a Poisson distribution, and
*
Corresponding author. Tel.: +86 532820321750; 53282032909. E-mail address: [email protected] (S. Wu).
fax:
+86
0030-4018/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2006.05.069
average pffiffiffi over n pulses decrease the noise magnitude from 1= n. However, such a process has a large bandwidth and poor cut-off spectral property, which is inefficient to handle the non-stationary noise especially in the far distance backscattering. It is also difficult to distinctly define the statistical properties of the lidar signal without a priori knowledge. The Fourier smoothing technique is another widely used technology for reducing noise. A common assumption of this process is that the information of a signal can be separated from the noise by taking into account that the signal varies slowly in comparison to the noise. Since lidar signals represent spatially varying information, setting a particular cutoff frequency may result in signal distortion. Even wavelet analysis is essentially an adjustable window Fourier spectral analysis [4]. Because of the limited length of the basic wavelet function, it is difficult for wavelets to quantitatively define the energy–frequency–time distribution. It is also difficult for wavelets to describe local frequency properties because the wavelet transform is based on predetermined stepping process.
S. Wu et al. / Optics Communications 267 (2006) 137–144
Recently, a new nonlinear technique, empirical mode decomposition (EMD), has been developed by Huang et al. to adaptively represent non-stationary signals as sums of intrinsic mode functions (IMFs) [4,5]. IMFs are derived directly from the data and are not restricted by linearity or a priori conceptions. Instead the method allows the modes to be nonlinear while still requiring local orthogonality in a least square sense. In this work, the lidar data from the simulation and measurements are analyzed by the EMD method. The layout of this paper is as follows. We briefly introduce in Section 2 the lidar system which data of this paper come from, and illustrate in Section 3 the principle of EMD. In Section 4.1, the simulated lidar signal is used to estimate the denoising performance of the EMD-based method and the traditional low-pass filter. In Section 4.2, the practical applications to lidar measurements are evaluated by comparing the SNR and power density spectrum of the original and denoised lidar signal. We propose to use the power density spectrum to determine how many IMFs can be removed. The possibility of using less averaging time to obtain the comparative SNR is also discussed. 2. Lidar system The lidar signal for analysis is from our Cabanne-Mie Doppler wind lidar based on the iodine vapor filter as shown in Fig. 1 [6,7]. The laser transmitter is an injection-seeded lamp-pumped frequency-doubled Nd:YAG laser which can be tuned and locked to the iodine absorption spectrum [8]. The laser pulse energy is around 100 mJ with the repetition rate of 10 Hz at 532 nm. A scanner with two azimuth and elevation mirrors enables the hemispherical measurement. The backscatterings are collected by the
Schmidt–Cassegrain telescope and coupled into a fiber cable that delivers the signal to the frequency discriminator, iodine vapor filter with a temperature controller. The wind speed is retrieved from the ratio of backscattering signals from the frequency discrimination channel and the energy reference channel. The lidar signal from the energy reference channel with 100-shot averaged are plotted in Fig. 2. The Doppler shift can be retrieved from both the Rayleigh scattering and the Mie scattering. Because of this unique property, the lidar signal in this work covers the low altitude atmosphere that varies more rapidly and is more sensitive to the ground interference than the high altitude atmosphere. Therefore, signal processing requires rapid and efficient techniques. We use the lidar backscatterings from the energy reference channel to illustrate the denoising process.
5
1.2x10
Range-corrected Lidar Signal
138
5
1.0x10
4
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
Range (km)
Fig. 2. The range corrected lidar return signal.
Fig. 1. The diagram of Doppler lidar with iodine vapor filter.
12
S. Wu et al. / Optics Communications 267 (2006) 137–144
3. Empirical mode decomposition EMD method decomposes the time series data into a series of IMFs with a zero local mean. The local mean is extracted by calculating the mean of the envelope of the data. This mean is iteratively subtracted from the current data until the residual has a local mean. This residual is then the first IMF that contains the highest frequencies of the time series. This process is so-called ‘sifting’ process. The subsequent IMFs can be found by subtracting the first IMF from the original data and repeating the above ‘sifting’ process. In this way all the IMFs can be extracted and the last IMF usually is a monotonic trend. The detailed information of decomposition refers to Ref. [4]. We use a measured lidar signal (Fig. 2) to illustrate the typical IMFs of the lidar backscattering. There are 21 IMFs and a trend mode obtained by the EMD method. Comparing this with the traditional Fourier expansion, one can see the efficiency of the EMD: the expansion of a turbulence data set with only 22 terms. Because we concern the high-frequency and small-scale components in the denoising process, we only show the first 6 IMFs in Fig. 3. From the result, the data is separated into locally non-overlapping time scale components. As we know, the lidar return signal is a kind of typical non-stationary time series. In lidar inversion studies, the data is usually normalised to the transmitting range. The noise components are also amplified in this range calibration process so that the real signals at far distance may
2x10
4
2.0x10
IMF2
4
0.0
-2.0x10
4
2.0x10
4
IMF3
4
2.0x10
4
4
2.0x10
4
0.0
-2.0x10
4
2.0x10
4
IMF6
1
10
0.0
-2.0x10
IMF5
4.1. Simulation
0.0
-2.0x10
IMF4
4. Data analysis
We use the lidar backscatters simulation program to obtain the ‘‘true’’ signal and the noisy lidar signal. It’s necessary because the ‘‘ground truth’’ is difficult to find to evaluate the denoising performance in the actual measurements. In the simulation, the Rayleigh and a cloud backscatters are simulated as the ‘‘true’’ data without noise using the same system parameters as our Doppler lidar (Fig. 4). The Gaussian white noise is introduced to simulate the ‘‘measured’’ lidar signal as shown in Fig. 5.
4
0 -2x1 0
be submerged by the background noise. The time-varying (or distance-changing) characteristics can be clearly seen in each IMF of Fig. 3, which is difficult for the Fourierbased method. The necessary conditions to describe a nonlinear and non-stationary time series in terms of basis function are completeness, orthogonality, locality and adaptivity. Completeness is needed to guarantee the precision of the basis expansion. Orthogonality is necessary to limit the energy leakage. Locality is required due to the non-stationary nature of the time series considered. In a non-stationary time series, there is no global time scale. Adaptivity is also crucial for the nonlinear and non-stationary time series in order to adjust to the local changes in the data. Only by resolving those changes one can account for the underlying nonlinear and non-stationary physics at hand. The EMD was developed to address these conditions and has been found to be very useful. A practical problem of the spline fitting in the ‘sifting’ process should be handled with care. The cubic spline fitting can have swings near the ends of the data, which could propagate inward and corrupt the data during the iterative decomposition. Huang et al. have devised a patent numerical method to eliminate the end effects. To avoid the divergence of the spline fitting, we herein add zero extrema at the both ends of data.
Signal Intensity
IMF1
139
0.0
-2.0x10
0
10
-1
10
4
0
2
4
6
8
10
Range (km) Fig. 3. The first 6 IMFs of the lidar return signal.
12
-2
0
2
4
6
8
10
12
14
Range (km)
Fig. 4. The simulated ‘‘true’’ lidar signal without noise.
16
140
S. Wu et al. / Optics Communications 267 (2006) 137–144 1
Imf6 0
1
10
Signal Intensity
1
Imf5 0 0
10
1
Imf4 0 1
-1
10
Imf3 -2
0
2
4
6
8
10
12
14
0
16
Range (km)
1
Fig. 5. The simulated ‘‘measured’’ lidar signal with the Gaussian white noise.
Imf2 0 2
There are 13 IMFs extracted from the simulated lidar signal, and the first 6 IMFs are shown in Fig. 6. Similar to Fig. 3, IMFs are sorted from high to low frequency although they will rarely have a constant frequency. We firstly evaluate the energy distribution of all the IMFs by power spectral density (PSD) (as shown in Fig. 7). We assume that high-frequency IMFs contain only noise and turbulence. This is a conservative estimation because those IMFs may contain useful signal. Whereas
IMF1
2 1 0 -1 -2 0.5
IMF2 0.0 -0.5 0.5
IMF3
0.0 -0.5 0.5
IMF4 0.0 -0.5 0.5
IMF5 0.0 -0.5
IMF6
2 1 0 -1 -2 0
2
4
6
8
10
12
14
16
Range (km) Fig. 6. The IMFs of the simulated ‘‘measured’’ lidar signal.
Imf1 1 0 0.0
0.1
0.2
0.3
0.4
0.5
Frequency (Hz)
Fig. 7. The IMFs’ PSD of the simulated ‘‘measured’’ signal.
the energy of the high frequency IMFs are of little value in the total backscattering, it is still practicable to improve SNR by subtracting high-frequency modes from the data. In addition, high-order IMFs often contain small spatial (or time) scale fluctuations that are much less than those of the wind speed we concerned. Therefore, high-frequency modes can also be removed to get the proper spatial resolution. The power spectrum of IMFs clearly shows that this method serves as a series band-pass filters, which is consistent with the result by Flandrin that EMD is kind of filterbank [9]. Wu and Huang also deduced in 2004 that the energy-density function of IMFs from white noise is chisquared distributed [10]. We compare the EMD-based result with the result by a Butterworth low-pass filter in the equivalent low-frequency band. There are four IMFs (IMF1–IMF4) that removed form the data in the EMD denoising process. Therefore, the cut-off frequency (3 dB band) of the butterworth lowpass filter is 0.06 Hz that is equal to the centre frequency of IMF4. Figs. 8 and 9 are the results by the low-pass filter and the EMD-based method, respectively. The lidar return by the EMD-based method is smoother than that of the low-pass filtering. The difference between the noise level of two results is not sufficient to evaluate two methods. The local details of the results in Fig. 10 distinctly show that the result by the EMD-based method is more consistent with the true data in the region of the near distance (Fig. 10(a)) and the cloud area (Fig. 10(b)) where sudden-change structures dominate. In the case of the lowpass filter, the denoised data are distorted because the
S. Wu et al. / Optics Communications 267 (2006) 137–144
141
De4Imfs LowPassButt True
1
2x10 1
Signal Intensity
Signal Intensity
10
0
10
1
10
-1
10
-2
0
2
4
6
8
10
12
14
16
Range (km)
0.1
(a)
0.2
0.3
0.4
0.5
Range (km)
Fig. 8. The denoised result by the Butterworth low-pass filter. De4Imfs LowPassButt True 0
10
Signal Intensity
1
Signal Intensity
10
0
10
-1
10
4
(b) -2
0
2
4
6
8
10
12
14
5
6
Range (km)
16
Range (km)
Fig. 10. The local details of the denoised data: (a) the near distance; (b) the cloud area.
Fig. 9. The denoised result by the EMD-based method.
entire high-frequency components are removed without any discrimination. Fig. 11 shows the linear regression of the ‘‘measured’’ lidar data and the denoised data as the function of the ‘‘true’’ data. We can see improvements of the standard deviation (SD) and correlation coefficient (R). In the case of noisy ‘‘measured’’ data SD is 0.144 and R is 0.952. For the low-pass filter, SD is 0.067 and R is 0.988. The EMD-based method has the best linear fit that SD is 0.039 and R is 0.996. In Fourier smoothing filtering, a common assumption is that the information content of a signal can be separated from the noise if the signal varies slowly in comparison to the noise. Since lidar backscattering signals represent spatially varying frequency components, setting a particular cutoff frequency could not handle the non-stationary noise and the sudden-change structure. The low-pass filter usually sets a cutoff frequency to remove noise such that all the frequency components above the cutoff are zero (or decay to zero). The low-pass filter can smooth the noisy components but blurs the edges of the signal. These distortions may lead to the bigger error than the noise itself, espe-
cially for the system based on the differential algorithm, such as the direct-detection Doppler lidar and the differential absorption lidar. 4.2. Measurements In practical processing, how many IMFs to be removed are determined by the noise level and the range resolution of the lidar signal. The PSD is used to analyze the IMFs of the lidar signal and the denoised data. The PSD of the measured lidar data (Fig. 2) is shown in Fig. 12. The highest frequency of PSD, 0.2 Hz, corresponds to the spatial frequency of 5 m that is twice as the range resolution of 2.5 m of the lidar signal (Nyquist frequency). The centre of the 5th IMF’s PSD locates at the 0.02 Hz corresponding to 50 m spatial resolution that satisfies the resolution requirement of the line-of-sight velocity retrieval. The first 4 IMFs have higher spatial frequency than the 5th IMF (Fig. 12(b)). In other words, regardless of noise distribution there are five high frequency IMFs at most that can be subtracted in order to obtain the 50 m resolution.
142
S. Wu et al. / Optics Communications 267 (2006) 137–144 100
20
Lidar signal Linear fit
Imf6 10
"measured" lidar signal
0 0.00 20
10
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.10
0.15
0.20
Imf5 10 0 0.00 20 1
Imf4 10 0 0.00 20
0.1 0.1
1
(a)
10
100
Imf3 10 0
"true" lidar signal
0.00 20
Denoised signal by Low-pass filter
100
Imf2 10 Lidar signal Linear fit
0 0.00 20
Imf1 10
10
0 0.00
0.05
(a) 1
Frequency (Hz) 20
Imf6 10 0 0.1 0.1
1
(b)
10
0.00 20
100
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.00
0.01
0.02
0.03
0.04
0.05
0.00
0.01
0.02
0.03
0.04
0.05
Imf5 10
"true" lidar signal
0 0.00 20
Denoised signal by EMD-based method
100
Lidar signal Linear fit
0 0.00 20
10
Imf3 10 0 0.00 20
1
Imf2 10 0 20 0.1 0.1
(c)
Imf4 10
1
10
100
Imf1 10
"true" lidar signal
Fig. 11. Comparison between the ‘‘true’’ data and (a) the simulated ‘‘measured’’ data, (b) the denoised data by low-pass filter and (c) the denoised data by the EMD-based method.
The denoised signal is achieved by subtracting the first 5 IMFs from the original signal (Fig. 13). The denoised signal is more smooth and less of fluctuations than the original signal. The partial data of the original signal and the denoised data are replotted in Fig. 14 for easy viewing. The fluctuations with the large magnitude but the small-
0
(b)
Frequency (Hz)
Fig. 12. Power Density spectrum of IMFs: (a) the PSD of IMFs; (b) the PSD with zoomed frequency coordinates.
scale are sufficiently suppressed in far distance. At the same time, local structures of the lidar return are preserved. For example, strong backscatters at 8 km are not distorted after denoising.
Range-corrected Lidar Signal
S. Wu et al. / Optics Communications 267 (2006) 137–144 1.2x10
5
1.0x10
5
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
4
the noise from the dark current and the readout electronics. Nb and Ne are estimated by the square of standard deviation of the lidar return signal in far distance where the background noise and electronics noise dominate. The SNR of the original and denoised data are calculated and plotted in Fig. 6. As mentioned in Section 2, the wind speed is calculated from the ratio of backscattered returns. The measurement uncertainty of LOS wind speed, VLOS (m/s), is defined in the following equation: dvLOS ¼ 0
2
4
6
8
10
12
Range (km)
Fig. 13. The denoised result of the lidar signal.
4
Original Denoised
ð2Þ
SV is the measurement sensitivity expressed in (m/s)1. As shown in Fig. 15, the SNR of the lidar return is obviously improved by a factor of 2.3. Consequently, EMD denoising can improve at least twice the retrieve accuracy. We compared this method with the multi-pulse average. The 1000-shot averaged data (As shown in Fig. 16) successive to the 100-shot averaged lidar signal in Fig. 2 are compared. The PSD of 100-shot averaged, the 1000-shot
4
6x10
100
Denoised data Original data
4
5x10
10
4
4x10
3.0
3.2
(a)
3.4
3.6
3.8
4.0
Range (km)
1
Original Denoised
5
10
Range-corrected Lidar Signal
1 S V ðSNRÞ
SNR
Range-corrected Lidar Signal
7x10
143
0.1
0
2
4
6
8
10
12
Range (km)
4
10
Fig. 15. SNR of the original and denoised data.
3
10
5
1.2x10
1000 shots Averaged
2
10 10.0
10.2
(b)
10.4
10.6
10.8
11.0
Range (km)
Fig. 14. The partial original and denoised data. lidar return (a) from 3 km to 4 km, (b) from 10 km to 11 km.
The denoising result is evaluated by SNR defined as follows: Ns SNR ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ns þ Nb þ Ne
ð1Þ
where Ns is the laser backscattering signal, Nb is the noise introduced by photons from the sky background, and Ne is
Range-corrected Lidar Signal
5
1.0x10
4
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
Range (km)
Fig. 16. Lidar return signal with 1000-shot averaged.
12
S. Wu et al. / Optics Communications 267 (2006) 137–144
Power
144 10
7
10
6
10
5
10
4
10
3
10
2
10
1
10
0
10
-1
10
-2
10
-3
tem can be less power consuming, which are important for the real-time monitoring and the low-cost laser transmitter. EMD-based method is also an option for the lidar working on the fast scanning mode such as planar position indicator (PPI) of the Doppler lidar and the volume imaging lidar.
Original data 1000 shots averaged data Denoised
5. Conclusion
4
In conclusion, it is the first application of the Empirical Mode Decomposition to the analysis of the lidar data. EMD analysis is implemented to reduce the noise and keep the significance of the signal. Comparing with the traditional low-pass filter and multi-pulse averaging method, the EMD-based method shows the attractive characteristics on handling the time-varying and nonlinear components of the lidar signal. We also propose to use the PDS as the criterion for the denoising approach. Furthermore, EMD is possibly the only computational method of the real-time lidar signal processing. Thus, EMD denoising is also useful in the case of scanning lidars requiring fast measurement.
4
Acknowledgements
0.00
0.05
0.10
0.15
0.2 0
Frequency (Hz)
Fig. 17. The PSD of 100-shot averaged data, denoised data and 1000-shot averaged data.
5
1.2x10
1000 shots Averaged Denoised
Range-corrected Lidar Signal
5
1.0x10
8.0x10
6.0x10
This work is partly supported by the Natural and Science Foundation of China Nos. 40505003 and 40427001 under grant and the research foundation of Information Science and Engineering Department, OUC.
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
12
References
Range (km)
Fig. 18. 1000-Shot averaged lidar return signal and 100-shot averaged data denoised by EMD.
averaged and the denoised data are shown in Fig. 17. The power of the 1000-shot averaged data at the high frequency is smaller than that of the original data. The EMD denoising shows even better noise suppressing performance than that of the 1000-shot averaging. Comparing with the 1000-shot averaged data, the power of the denoised data at high frequency decreases by a factor of 10. The most attractive specific is the EMD-based approach only needs a small quantity of data for average when it achieves comparative performance and obtains the instantaneous atmospheric motion (Fig. 18). With this method, lidar researchers may reduce the measurement time and the sys-
[1] M.J.T. Milton, E.T. Woods, Appl. Opt. 26 (1987) 2598. [2] N. Menyuk, D.K. Killinger, C.R. Menyuk, Appl. Opt. 21 (1982) 3377. [3] J.M. Dias, E.S. Fonseca, D.P. Resendes, SPIE: Application of Lidar to Current Atmospheric 3757 (1999) 158. [4] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, Proc. Roy. Soc. Lond. A 454 (1998) 903. [5] Norden E. Huang, Man-Li C. Wu, Steven R. Long, Samuel S.P. Shen, Wendong Qu, Per Gloersen, Kuang L. Fan, Proc. Roy. Soc. Lond. A 459 (2003) 2317. [6] Liu Zhi-Shen, Chen Wei-Biao, Zhang Ting-Lu, John W. Hair, C.Y. She, Appl. Phys. B. 64 (1997) 561. [7] Zhi-Shen Liu, Dong Wu, Jin-Tao Liu, Kai-Lin Zhang, Wei-Biao Chen, Xiao-Quan Song, Johnathan W. Hair, Chiao-Yao She, Appl. Opt. 41 (33) (2002) 7079. [8] Song-Hua Wu, Da-Peng Sun, Zhi-Shen Liu, SPIE 4893 (2003) 348. [9] Patrick Flandrin, Gabriel Rilling, Paulo Goncalves, IEEE Signal Process. Lett. 11 (2004) 112. [10] Z. Wu, N.E. Huang, Proc. Roy. Soc. Lond. A 460 (2004) 1597.
Enhancement of lidar backscatters signal-to-noise ratio using empirical mode decomposition method Songhua Wu *, Zhishen Liu, Bingyi Liu Ocean Remote Sensing Laboratory of Ministry of Education of China, Ocean Remote Sensing Institute (ORSI), Ocean University of China, No. 5 YuShan Road, Qingdao 266003, China Received 21 September 2005; received in revised form 1 April 2006; accepted 30 May 2006
Abstract Lidar is being widely used to monitor meteorological parameters and atmospheric constituents. Applications include meteorology, environmental pollution, atmospheric dynamics and global climate change. Signal processing for lidar applications involve highly nonlinear models and consequently nonlinear filtering. In this paper, we applied a new method, empirical mode decomposition to the lidar signal processing. The denoising approach is done by removal of the proper intrinsic mode functions. The data from the simulation and measurements are analyzed to evaluate this method comparing with the traditional low-pass filter and the multi-pulse averaging. Results show that it is effective and superior to the band-pass filter and the averaging method. The denoising method also allows less averaging laser shots which is important for the real-time monitoring and for the low cost laser transmitter. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Lidar; Signal-to-noise ratio; Empirical mode decomposition
1. Introduction Lidar is an active remote sensing instrument that emits laser pulses towards the atmosphere or target and measures the backscattering. The signal-to-noise ratio (SNR) of the lidar backscattering often attenuates due to the noise and interferences, such as stochastic turbulences, background noise, dark current, electronics readout noise and atmospheric turbulence. In order to improve the precision of the measurement and subsequent analysis, lidar researchers usually use the multiple-pulse average (temporal average) or running-average (spatial average) approach to smooth the lidar signal [1–3]. The conventional average approach is a kind of lowpass filtering process based on the least square. System errors are assumed to have a Poisson distribution, and
*
Corresponding author. Tel.: +86 532820321750; 53282032909. E-mail address: [email protected] (S. Wu).
fax:
+86
0030-4018/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2006.05.069
average pffiffiffi over n pulses decrease the noise magnitude from 1= n. However, such a process has a large bandwidth and poor cut-off spectral property, which is inefficient to handle the non-stationary noise especially in the far distance backscattering. It is also difficult to distinctly define the statistical properties of the lidar signal without a priori knowledge. The Fourier smoothing technique is another widely used technology for reducing noise. A common assumption of this process is that the information of a signal can be separated from the noise by taking into account that the signal varies slowly in comparison to the noise. Since lidar signals represent spatially varying information, setting a particular cutoff frequency may result in signal distortion. Even wavelet analysis is essentially an adjustable window Fourier spectral analysis [4]. Because of the limited length of the basic wavelet function, it is difficult for wavelets to quantitatively define the energy–frequency–time distribution. It is also difficult for wavelets to describe local frequency properties because the wavelet transform is based on predetermined stepping process.
S. Wu et al. / Optics Communications 267 (2006) 137–144
Recently, a new nonlinear technique, empirical mode decomposition (EMD), has been developed by Huang et al. to adaptively represent non-stationary signals as sums of intrinsic mode functions (IMFs) [4,5]. IMFs are derived directly from the data and are not restricted by linearity or a priori conceptions. Instead the method allows the modes to be nonlinear while still requiring local orthogonality in a least square sense. In this work, the lidar data from the simulation and measurements are analyzed by the EMD method. The layout of this paper is as follows. We briefly introduce in Section 2 the lidar system which data of this paper come from, and illustrate in Section 3 the principle of EMD. In Section 4.1, the simulated lidar signal is used to estimate the denoising performance of the EMD-based method and the traditional low-pass filter. In Section 4.2, the practical applications to lidar measurements are evaluated by comparing the SNR and power density spectrum of the original and denoised lidar signal. We propose to use the power density spectrum to determine how many IMFs can be removed. The possibility of using less averaging time to obtain the comparative SNR is also discussed. 2. Lidar system The lidar signal for analysis is from our Cabanne-Mie Doppler wind lidar based on the iodine vapor filter as shown in Fig. 1 [6,7]. The laser transmitter is an injection-seeded lamp-pumped frequency-doubled Nd:YAG laser which can be tuned and locked to the iodine absorption spectrum [8]. The laser pulse energy is around 100 mJ with the repetition rate of 10 Hz at 532 nm. A scanner with two azimuth and elevation mirrors enables the hemispherical measurement. The backscatterings are collected by the
Schmidt–Cassegrain telescope and coupled into a fiber cable that delivers the signal to the frequency discriminator, iodine vapor filter with a temperature controller. The wind speed is retrieved from the ratio of backscattering signals from the frequency discrimination channel and the energy reference channel. The lidar signal from the energy reference channel with 100-shot averaged are plotted in Fig. 2. The Doppler shift can be retrieved from both the Rayleigh scattering and the Mie scattering. Because of this unique property, the lidar signal in this work covers the low altitude atmosphere that varies more rapidly and is more sensitive to the ground interference than the high altitude atmosphere. Therefore, signal processing requires rapid and efficient techniques. We use the lidar backscatterings from the energy reference channel to illustrate the denoising process.
5
1.2x10
Range-corrected Lidar Signal
138
5
1.0x10
4
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
Range (km)
Fig. 2. The range corrected lidar return signal.
Fig. 1. The diagram of Doppler lidar with iodine vapor filter.
12
S. Wu et al. / Optics Communications 267 (2006) 137–144
3. Empirical mode decomposition EMD method decomposes the time series data into a series of IMFs with a zero local mean. The local mean is extracted by calculating the mean of the envelope of the data. This mean is iteratively subtracted from the current data until the residual has a local mean. This residual is then the first IMF that contains the highest frequencies of the time series. This process is so-called ‘sifting’ process. The subsequent IMFs can be found by subtracting the first IMF from the original data and repeating the above ‘sifting’ process. In this way all the IMFs can be extracted and the last IMF usually is a monotonic trend. The detailed information of decomposition refers to Ref. [4]. We use a measured lidar signal (Fig. 2) to illustrate the typical IMFs of the lidar backscattering. There are 21 IMFs and a trend mode obtained by the EMD method. Comparing this with the traditional Fourier expansion, one can see the efficiency of the EMD: the expansion of a turbulence data set with only 22 terms. Because we concern the high-frequency and small-scale components in the denoising process, we only show the first 6 IMFs in Fig. 3. From the result, the data is separated into locally non-overlapping time scale components. As we know, the lidar return signal is a kind of typical non-stationary time series. In lidar inversion studies, the data is usually normalised to the transmitting range. The noise components are also amplified in this range calibration process so that the real signals at far distance may
2x10
4
2.0x10
IMF2
4
0.0
-2.0x10
4
2.0x10
4
IMF3
4
2.0x10
4
4
2.0x10
4
0.0
-2.0x10
4
2.0x10
4
IMF6
1
10
0.0
-2.0x10
IMF5
4.1. Simulation
0.0
-2.0x10
IMF4
4. Data analysis
We use the lidar backscatters simulation program to obtain the ‘‘true’’ signal and the noisy lidar signal. It’s necessary because the ‘‘ground truth’’ is difficult to find to evaluate the denoising performance in the actual measurements. In the simulation, the Rayleigh and a cloud backscatters are simulated as the ‘‘true’’ data without noise using the same system parameters as our Doppler lidar (Fig. 4). The Gaussian white noise is introduced to simulate the ‘‘measured’’ lidar signal as shown in Fig. 5.
4
0 -2x1 0
be submerged by the background noise. The time-varying (or distance-changing) characteristics can be clearly seen in each IMF of Fig. 3, which is difficult for the Fourierbased method. The necessary conditions to describe a nonlinear and non-stationary time series in terms of basis function are completeness, orthogonality, locality and adaptivity. Completeness is needed to guarantee the precision of the basis expansion. Orthogonality is necessary to limit the energy leakage. Locality is required due to the non-stationary nature of the time series considered. In a non-stationary time series, there is no global time scale. Adaptivity is also crucial for the nonlinear and non-stationary time series in order to adjust to the local changes in the data. Only by resolving those changes one can account for the underlying nonlinear and non-stationary physics at hand. The EMD was developed to address these conditions and has been found to be very useful. A practical problem of the spline fitting in the ‘sifting’ process should be handled with care. The cubic spline fitting can have swings near the ends of the data, which could propagate inward and corrupt the data during the iterative decomposition. Huang et al. have devised a patent numerical method to eliminate the end effects. To avoid the divergence of the spline fitting, we herein add zero extrema at the both ends of data.
Signal Intensity
IMF1
139
0.0
-2.0x10
0
10
-1
10
4
0
2
4
6
8
10
Range (km) Fig. 3. The first 6 IMFs of the lidar return signal.
12
-2
0
2
4
6
8
10
12
14
Range (km)
Fig. 4. The simulated ‘‘true’’ lidar signal without noise.
16
140
S. Wu et al. / Optics Communications 267 (2006) 137–144 1
Imf6 0
1
10
Signal Intensity
1
Imf5 0 0
10
1
Imf4 0 1
-1
10
Imf3 -2
0
2
4
6
8
10
12
14
0
16
Range (km)
1
Fig. 5. The simulated ‘‘measured’’ lidar signal with the Gaussian white noise.
Imf2 0 2
There are 13 IMFs extracted from the simulated lidar signal, and the first 6 IMFs are shown in Fig. 6. Similar to Fig. 3, IMFs are sorted from high to low frequency although they will rarely have a constant frequency. We firstly evaluate the energy distribution of all the IMFs by power spectral density (PSD) (as shown in Fig. 7). We assume that high-frequency IMFs contain only noise and turbulence. This is a conservative estimation because those IMFs may contain useful signal. Whereas
IMF1
2 1 0 -1 -2 0.5
IMF2 0.0 -0.5 0.5
IMF3
0.0 -0.5 0.5
IMF4 0.0 -0.5 0.5
IMF5 0.0 -0.5
IMF6
2 1 0 -1 -2 0
2
4
6
8
10
12
14
16
Range (km) Fig. 6. The IMFs of the simulated ‘‘measured’’ lidar signal.
Imf1 1 0 0.0
0.1
0.2
0.3
0.4
0.5
Frequency (Hz)
Fig. 7. The IMFs’ PSD of the simulated ‘‘measured’’ signal.
the energy of the high frequency IMFs are of little value in the total backscattering, it is still practicable to improve SNR by subtracting high-frequency modes from the data. In addition, high-order IMFs often contain small spatial (or time) scale fluctuations that are much less than those of the wind speed we concerned. Therefore, high-frequency modes can also be removed to get the proper spatial resolution. The power spectrum of IMFs clearly shows that this method serves as a series band-pass filters, which is consistent with the result by Flandrin that EMD is kind of filterbank [9]. Wu and Huang also deduced in 2004 that the energy-density function of IMFs from white noise is chisquared distributed [10]. We compare the EMD-based result with the result by a Butterworth low-pass filter in the equivalent low-frequency band. There are four IMFs (IMF1–IMF4) that removed form the data in the EMD denoising process. Therefore, the cut-off frequency (3 dB band) of the butterworth lowpass filter is 0.06 Hz that is equal to the centre frequency of IMF4. Figs. 8 and 9 are the results by the low-pass filter and the EMD-based method, respectively. The lidar return by the EMD-based method is smoother than that of the low-pass filtering. The difference between the noise level of two results is not sufficient to evaluate two methods. The local details of the results in Fig. 10 distinctly show that the result by the EMD-based method is more consistent with the true data in the region of the near distance (Fig. 10(a)) and the cloud area (Fig. 10(b)) where sudden-change structures dominate. In the case of the lowpass filter, the denoised data are distorted because the
S. Wu et al. / Optics Communications 267 (2006) 137–144
141
De4Imfs LowPassButt True
1
2x10 1
Signal Intensity
Signal Intensity
10
0
10
1
10
-1
10
-2
0
2
4
6
8
10
12
14
16
Range (km)
0.1
(a)
0.2
0.3
0.4
0.5
Range (km)
Fig. 8. The denoised result by the Butterworth low-pass filter. De4Imfs LowPassButt True 0
10
Signal Intensity
1
Signal Intensity
10
0
10
-1
10
4
(b) -2
0
2
4
6
8
10
12
14
5
6
Range (km)
16
Range (km)
Fig. 10. The local details of the denoised data: (a) the near distance; (b) the cloud area.
Fig. 9. The denoised result by the EMD-based method.
entire high-frequency components are removed without any discrimination. Fig. 11 shows the linear regression of the ‘‘measured’’ lidar data and the denoised data as the function of the ‘‘true’’ data. We can see improvements of the standard deviation (SD) and correlation coefficient (R). In the case of noisy ‘‘measured’’ data SD is 0.144 and R is 0.952. For the low-pass filter, SD is 0.067 and R is 0.988. The EMD-based method has the best linear fit that SD is 0.039 and R is 0.996. In Fourier smoothing filtering, a common assumption is that the information content of a signal can be separated from the noise if the signal varies slowly in comparison to the noise. Since lidar backscattering signals represent spatially varying frequency components, setting a particular cutoff frequency could not handle the non-stationary noise and the sudden-change structure. The low-pass filter usually sets a cutoff frequency to remove noise such that all the frequency components above the cutoff are zero (or decay to zero). The low-pass filter can smooth the noisy components but blurs the edges of the signal. These distortions may lead to the bigger error than the noise itself, espe-
cially for the system based on the differential algorithm, such as the direct-detection Doppler lidar and the differential absorption lidar. 4.2. Measurements In practical processing, how many IMFs to be removed are determined by the noise level and the range resolution of the lidar signal. The PSD is used to analyze the IMFs of the lidar signal and the denoised data. The PSD of the measured lidar data (Fig. 2) is shown in Fig. 12. The highest frequency of PSD, 0.2 Hz, corresponds to the spatial frequency of 5 m that is twice as the range resolution of 2.5 m of the lidar signal (Nyquist frequency). The centre of the 5th IMF’s PSD locates at the 0.02 Hz corresponding to 50 m spatial resolution that satisfies the resolution requirement of the line-of-sight velocity retrieval. The first 4 IMFs have higher spatial frequency than the 5th IMF (Fig. 12(b)). In other words, regardless of noise distribution there are five high frequency IMFs at most that can be subtracted in order to obtain the 50 m resolution.
142
S. Wu et al. / Optics Communications 267 (2006) 137–144 100
20
Lidar signal Linear fit
Imf6 10
"measured" lidar signal
0 0.00 20
10
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.05
0.10
0.15
0.20
0.10
0.15
0.20
Imf5 10 0 0.00 20 1
Imf4 10 0 0.00 20
0.1 0.1
1
(a)
10
100
Imf3 10 0
"true" lidar signal
0.00 20
Denoised signal by Low-pass filter
100
Imf2 10 Lidar signal Linear fit
0 0.00 20
Imf1 10
10
0 0.00
0.05
(a) 1
Frequency (Hz) 20
Imf6 10 0 0.1 0.1
1
(b)
10
0.00 20
100
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.01
0.02
0.03
0.04
0.05
0.00
0.01
0.02
0.03
0.04
0.05
0.00
0.01
0.02
0.03
0.04
0.05
Imf5 10
"true" lidar signal
0 0.00 20
Denoised signal by EMD-based method
100
Lidar signal Linear fit
0 0.00 20
10
Imf3 10 0 0.00 20
1
Imf2 10 0 20 0.1 0.1
(c)
Imf4 10
1
10
100
Imf1 10
"true" lidar signal
Fig. 11. Comparison between the ‘‘true’’ data and (a) the simulated ‘‘measured’’ data, (b) the denoised data by low-pass filter and (c) the denoised data by the EMD-based method.
The denoised signal is achieved by subtracting the first 5 IMFs from the original signal (Fig. 13). The denoised signal is more smooth and less of fluctuations than the original signal. The partial data of the original signal and the denoised data are replotted in Fig. 14 for easy viewing. The fluctuations with the large magnitude but the small-
0
(b)
Frequency (Hz)
Fig. 12. Power Density spectrum of IMFs: (a) the PSD of IMFs; (b) the PSD with zoomed frequency coordinates.
scale are sufficiently suppressed in far distance. At the same time, local structures of the lidar return are preserved. For example, strong backscatters at 8 km are not distorted after denoising.
Range-corrected Lidar Signal
S. Wu et al. / Optics Communications 267 (2006) 137–144 1.2x10
5
1.0x10
5
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
4
the noise from the dark current and the readout electronics. Nb and Ne are estimated by the square of standard deviation of the lidar return signal in far distance where the background noise and electronics noise dominate. The SNR of the original and denoised data are calculated and plotted in Fig. 6. As mentioned in Section 2, the wind speed is calculated from the ratio of backscattered returns. The measurement uncertainty of LOS wind speed, VLOS (m/s), is defined in the following equation: dvLOS ¼ 0
2
4
6
8
10
12
Range (km)
Fig. 13. The denoised result of the lidar signal.
4
Original Denoised
ð2Þ
SV is the measurement sensitivity expressed in (m/s)1. As shown in Fig. 15, the SNR of the lidar return is obviously improved by a factor of 2.3. Consequently, EMD denoising can improve at least twice the retrieve accuracy. We compared this method with the multi-pulse average. The 1000-shot averaged data (As shown in Fig. 16) successive to the 100-shot averaged lidar signal in Fig. 2 are compared. The PSD of 100-shot averaged, the 1000-shot
4
6x10
100
Denoised data Original data
4
5x10
10
4
4x10
3.0
3.2
(a)
3.4
3.6
3.8
4.0
Range (km)
1
Original Denoised
5
10
Range-corrected Lidar Signal
1 S V ðSNRÞ
SNR
Range-corrected Lidar Signal
7x10
143
0.1
0
2
4
6
8
10
12
Range (km)
4
10
Fig. 15. SNR of the original and denoised data.
3
10
5
1.2x10
1000 shots Averaged
2
10 10.0
10.2
(b)
10.4
10.6
10.8
11.0
Range (km)
Fig. 14. The partial original and denoised data. lidar return (a) from 3 km to 4 km, (b) from 10 km to 11 km.
The denoising result is evaluated by SNR defined as follows: Ns SNR ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ns þ Nb þ Ne
ð1Þ
where Ns is the laser backscattering signal, Nb is the noise introduced by photons from the sky background, and Ne is
Range-corrected Lidar Signal
5
1.0x10
4
8.0x10
4
6.0x10
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
Range (km)
Fig. 16. Lidar return signal with 1000-shot averaged.
12
S. Wu et al. / Optics Communications 267 (2006) 137–144
Power
144 10
7
10
6
10
5
10
4
10
3
10
2
10
1
10
0
10
-1
10
-2
10
-3
tem can be less power consuming, which are important for the real-time monitoring and the low-cost laser transmitter. EMD-based method is also an option for the lidar working on the fast scanning mode such as planar position indicator (PPI) of the Doppler lidar and the volume imaging lidar.
Original data 1000 shots averaged data Denoised
5. Conclusion
4
In conclusion, it is the first application of the Empirical Mode Decomposition to the analysis of the lidar data. EMD analysis is implemented to reduce the noise and keep the significance of the signal. Comparing with the traditional low-pass filter and multi-pulse averaging method, the EMD-based method shows the attractive characteristics on handling the time-varying and nonlinear components of the lidar signal. We also propose to use the PDS as the criterion for the denoising approach. Furthermore, EMD is possibly the only computational method of the real-time lidar signal processing. Thus, EMD denoising is also useful in the case of scanning lidars requiring fast measurement.
4
Acknowledgements
0.00
0.05
0.10
0.15
0.2 0
Frequency (Hz)
Fig. 17. The PSD of 100-shot averaged data, denoised data and 1000-shot averaged data.
5
1.2x10
1000 shots Averaged Denoised
Range-corrected Lidar Signal
5
1.0x10
8.0x10
6.0x10
This work is partly supported by the Natural and Science Foundation of China Nos. 40505003 and 40427001 under grant and the research foundation of Information Science and Engineering Department, OUC.
4
4.0x10
4
2.0x10
0.0 0
2
4
6
8
10
12
References
Range (km)
Fig. 18. 1000-Shot averaged lidar return signal and 100-shot averaged data denoised by EMD.
averaged and the denoised data are shown in Fig. 17. The power of the 1000-shot averaged data at the high frequency is smaller than that of the original data. The EMD denoising shows even better noise suppressing performance than that of the 1000-shot averaging. Comparing with the 1000-shot averaged data, the power of the denoised data at high frequency decreases by a factor of 10. The most attractive specific is the EMD-based approach only needs a small quantity of data for average when it achieves comparative performance and obtains the instantaneous atmospheric motion (Fig. 18). With this method, lidar researchers may reduce the measurement time and the sys-
[1] M.J.T. Milton, E.T. Woods, Appl. Opt. 26 (1987) 2598. [2] N. Menyuk, D.K. Killinger, C.R. Menyuk, Appl. Opt. 21 (1982) 3377. [3] J.M. Dias, E.S. Fonseca, D.P. Resendes, SPIE: Application of Lidar to Current Atmospheric 3757 (1999) 158. [4] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, Proc. Roy. Soc. Lond. A 454 (1998) 903. [5] Norden E. Huang, Man-Li C. Wu, Steven R. Long, Samuel S.P. Shen, Wendong Qu, Per Gloersen, Kuang L. Fan, Proc. Roy. Soc. Lond. A 459 (2003) 2317. [6] Liu Zhi-Shen, Chen Wei-Biao, Zhang Ting-Lu, John W. Hair, C.Y. She, Appl. Phys. B. 64 (1997) 561. [7] Zhi-Shen Liu, Dong Wu, Jin-Tao Liu, Kai-Lin Zhang, Wei-Biao Chen, Xiao-Quan Song, Johnathan W. Hair, Chiao-Yao She, Appl. Opt. 41 (33) (2002) 7079. [8] Song-Hua Wu, Da-Peng Sun, Zhi-Shen Liu, SPIE 4893 (2003) 348. [9] Patrick Flandrin, Gabriel Rilling, Paulo Goncalves, IEEE Signal Process. Lett. 11 (2004) 112. [10] Z. Wu, N.E. Huang, Proc. Roy. Soc. Lond. A 460 (2004) 1597.