A method based on iterative morphological filtering and multiple scattering for detecting layer boundaries and extinction coefficients with LIDAR

Layer boundaries detection with LIDAR is of great significance for the meteorological and environmental research. Apart from the background noise, multiple scattering can also seriously affect the detection results in LIDAR signal processing. To alleviate these issues, a novel approach was proposed based upon morphological filtering and multiple scattering correction with multiple iterations, which essentially acts as a weighted algorithm with multiple scattering factors in different filtering scales, and applies integral extinction coefficients as media to perform correction. Simulations on artificial signals and real LIDAR signals support this approach.


Introduction
Detecting and quantifying the spatial distributions and sizes of clouds and aerosol particles can be useful for the climate change investigation and aeronautical meteorological research [1][2][3][4]. Many equipments have been utilized to detect the optical properties of clouds and aerosol particles, such as nephelometer and sun-photometer, but LIDAR is a popular probe that can be used to scan the spacial physical properties in three dimensions, thereby reflecting clouds and aerosol particles distribution completely [5,6]. However, the original LIDAR signal and related retrieval parameters are affected by background noise and multiple scattering seriously. Thus, it is necessary to develop practical methods to mitigate these problems.
As a form of interference, Background noise exhibits non-linear and non-stationary. In the early research [7], windows and adaptive algorithm were proposed to denoise the LIDAR signal, which applies multiple point windows to obtain the layer base with a slope of the original signal and takes the complexities of the LIDAR signal into account. In addition, signal decomposition techniques were also introduced to eliminate signal noise, such as wavelet, empirical mode decomposition (EMD), and their variants [8][9][10][11] etc. Non-linear and non-stationary are considered, but low-pass filtering based on the allover LIDAR signal performs in these methods [8][9][10][11], thereby the boundary characteristics cannot be revealed completely. In 2010, a differential zero-crossing method [12] was demonstrated to process the LIDAR signal, which applies the signal intensity and cloud information to reduce the error. Layer boundary characteristics can be considered, but signal pretreatment must be carried out before layer detection; some useful information of layer boundaries may be missed accordingly. Recently, some segmentation methods [13][14][15] about simplifying the original signal were proposed, which conduct based on characteristic points and arrange the signal simplification and layer detection in one step. Useful signal messages can be retained and reveal the layer boundary characteristics, whereas determining simplification scale is a challenge in the segmentation process. Furthermore, in the case of complicated signal, filtering approach is more mature and comprehensive than most simplification approaches.
Multiple scattering is a kind of optical interference from the layers and photons interactions. Thus, it is generally recognized that the optical properties (optical depth, extinction coefficients, LIDAR ratio, etc.) can be used to estimate multiple scattering characteristics of the cloud. The experimental model on multiple scattering was proposed in 1976 [16]. Then, the importance of considering multiple scattering effect in layer detection was first emphasized by Platt et al. [17]. Moreover, multiple scattering factors were presented with simpler and more practical description of the multiple scattering degrees [17]. In 2013, an initial scheme for correcting the cloud height by multiple scattering factors was introduced [18]. The error is described and calculated by a geometry description of part extinction coefficients. However, it is a pity that the light fluctuation of extinction coefficients is ignored in the noncloud areas, and influences of multi-layer cirrus on multiple scattering factors are not considered.
Aiming at alleviating the noise and multiple scattering effects on layer detection, a novel scheme is introduced based on signal iterative morphological filtering and multiple scattering factors. Characteristic point detection and signal simplification perform in one step without choosing the simplifying window, and the iterative morphological filtering is employed to process the detected signal combining with the multiple scattering correction. Moreover, the overall extinction coefficients are applied to calculate the multiple scattering factors, which make the ultimate values more accurate. The layout of this paper is arranged as follows: In Sect. 2, the overall strategies of this method are shown in the flowchart. Detailed descriptions are demonstrated step-by-step. Consequently, the experiments with simulated and real signals are exhibited, respectively, in Sect. 3, and the performance evaluations of this approach are illustrated. Finally, the conclusions are shown in Sect. 4.

Principles and methods
The back-scatter power P(r) received by the LIDAR detector can be presented as follows: where C is LIDAR constant, z is the range, and z b represents the layer base. b a z ð Þ and b m z ð Þ indicate the aerosol and molecule back-scattering coefficients, respectively. r a z ð Þ and r m z ð Þ are the aerosol and molecule extinction coefficients, respectively. e z ð Þ indicates the noise. The multiple scattering factor g z ð Þ is a significant parameter representing the degree of cloud multiple scattering. Below the cloud, the aerosol and molecular transmittance are ignored (their transmittances are set to 1). According to characteristics of the atmosphere layer as well as Eq (1), a novel layer detection method is proposed, which contains four key steps: preliminary layer detection, iterative morphological filtering, multiple scattering factors calculation, and layer boundaries correction. The flowchart of this scheme is demonstrated in Fig. 1

Preliminary layer detection
Preliminary layer detection is a process of boundaries determining and signal simplification; the purpose of this algorithm is to search the possible boundaries of layers, and to find a similar curve with fewer points that can describe the original signal characteristics. Unlike general methods applying low-pass filter and zero-crossing, this technique only focuses on the feature-point detection in principle. It can exclude the non-feature parts of the signal and retain the possible layer boundary points that would be eliminated by low-pass filter. Detail procedures are demonstrated as the following: Starting with a LIDAR signal that consists of n points, maxima, and minima of the signal are identified, termed p 1 ; p 2 ; . . .; p n 1 (maxima) and v 1 ; v 2 ; . . .; v n 2 (minima), respectively (Fig. 2). LIDAR signal is generally non-linear and non-stationary; thereby extremes (maxima and minima) at local waveform can carry the noise as well as the information of layers. Through investigating the relation of these extremes, the characteristics of different layers and noise are obtained. Then, the possible boundaries of layers can be discriminated by an appropriate threshold.
The range-corrected signal Pz 2 is applied instead of original signal P in layer detection. The relation of maxima and minima, hereafter referred to as the peak-base ratio (PBR), is considered to be: where P p i z 2 p i and P v j z 2 v i represent the range-corrected signals of the maxima p i and minima v i , respectively (p i and v i are the ith maxima and minima). The PBR may be an approximation, because the so-called peak or base can be effected by the noise, but it can be applied to investigate any possible layers of the LIDAR signal. In the study of PBR, Morille [11] assumed that above 7.5 km only, cloud layers contribute to the distribution of PBR; whereas under 7.5 km, both cloud and aerosol layers are components of the distribution of PBR. According to the probability density function (PDF) of PBR, Morille set a threshold value T to be 4 to separate cloud layers from aerosols for the particle layers over the full the LIDAR range, which can be translated to: where P b z 2 b and P p z 2 p are the range-corrected signals of the base and peak of a layer, respectively. In this study, 4 is selected as the threshold T based upon the previous research. Therefore, according to Eqs. (2) and (3), a novel function of layer detection is determined as: where z b and z t are the base and top height of the layer, respectively.

Iterative morphological filtering
After determining the possible boundaries of layers, morphological filtering factors are applied to measure the morphological gaps of two peak (or valley) in non-layer areas, referred to as D 1 ; D 2 ; . . .; D n ( Fig. 3a). Then, three important morphological filtering parameters can be expressed as: where D min and D max are the minimum and maximum morphological filtering scales, respectively, which determine the scale range of morphological filter and the interval D 0 that can be employed to implement iterative calculation. Then, the original LIDAR signal filtering performs with the spherical factor of the scale D min , and the smoothed signal P 0 is acquired. Similarly, a series of morphological filters can be used to process the original LIDAR signal with spherical morphological factors of the scales D min þ D 0 , D min þ 2D 0 ,…, D max ; denoised signals can be termed as P 1 , P 2 ,…, P nÀ1 . Through rolling the morphological factors along the waveform, the signal can be smoothed in different scales, which is illustrated in Fig. 3b.

Multiple scattering factors calculation
The values of multiple scattering factors depend on several factors, particularly the cloud particle phase function, the cloud optical depth, and the cloud range etc [19]. Multiple scattering factors are derived through the following steps: (1) segment the simplified signal and establish the nonlinear equations to compute the LIDAR ratio. (2) Compute the multiple scattering factors with the calculated LIDAR ratio.
In segmentation, three feature points: the layer base z b , the layer top z t , and a point in a low-altitude region (the altitude lower than z b ) z c are employed to divide the simplified signal. z b and z t have been obtained in the preliminary layer detection; z c is the z corresponding to the minimum value of the function P z ð Þz 2 b m z ð Þ . The extinction coefficients at z b , z t , z c , and the LIDAR ratio for cloud are named as x 1 , x 2 , x 3 , and s a , respectively; then, non-linear equations are established as follows: where S z ð Þ is range-corrected signal computed as S z ð Þ ¼ P z ð Þz 2 and s m ¼ 8p=3 is the LIDAR ratio for molecular and assumed according to the Rayleigh scattering theory. By solving the equations, the LIDAR ratio can be obtained.
In Eq (9), the first and fifth equations determine the extinction coefficients of z b according to ref [20]. The second and third equations are constructed from Fernald formulation [21], with z ¼ z c and z ¼ z b . Chen's method [19] forms the fourth equation, which uses range-corrected signals to describe the laser transmission as well as the relation between the optical depth and extinction coefficients.
Opt Rev (2016) 23:646-656 649 After deriving the LIDAR ratio, the extinction coefficients can be expressed as [20]: where s 0 a is the LIDAR ratio outside the layers that is assumed to be 50. Inside the layer, the LIDAR ratio s a is computed using Eq (9). This technique proposes an assumed boundary z c to determine extinction coefficients r and segments the atmosphere according to layer detection.
The multiple scattering factor g is proposed to represent the degree of atmosphere multiple scattering, which varies with the thickness and density of the layer. Outside the layer, g is assumed to be 1. Inside, the multiple scattering factor g can be computed as [20 where m expresses the proportion of the single scattering and multiple scattering of light transmission in the atmosphere, which can reflect the statistics microscopic state of laser transmission. The Monte Carlo (MC) approach is used to compute g according to the ''Appendix''. Table 1 presents the parameters of simulations using the MC method according to the LIDAR model. Through applying the above method to P 0 , P 1 ,…, P nÀ1 , a series of multiple scattering factors and extinction coefficients can be acquired, termed as g 0 , g 1 ,…, g nÀ1 and r 0 , r 1 ,…, r nÀ1 , respectively.

The layer boundaries detection
After obtaining a series of multiple scattering factors and extinction coefficients via allover detection arrange, the corrected extinction coefficients can be calculated, according to: r 0 i z ð Þ is the ith revised extinction coefficients, by which a series of revised returned LIDAR signals can be calculated based on: where C is LIDAR constant, and b 0 i z ð Þ represents the ith revised aerosol back-scattering coefficients.
(2)-(5) to P 0 0 , P 0 1 ,…, P 0 nÀ1 , a series of revised layer boundaries can be obtained, named as z 1 b , z 2 b ,…, z n b and z 2 t , z 3 t ,…, z n t , respectively, which can outline the effect of noise with different degree of useful signal retention. Applying these corrected layer boundaries, the top and base values via the morphological scales can be acquired. Then, the final corrected top and base of layer are derived by the weighted least square method, which can be expressed as: Searching the extremes of Eqs. (14) and (15), the ultimate values of layer boundaries can be expressed as: where the processed layer boundaries are weighted according to the morphological filtering scales. This technique can make the corrected top and base accurate in the case of signal-noise balance.
3 Results and discussions

Experiment with simulated signals
An extensive experiment is performed to verify the proposed method using simulated signals. Three original artificial signals are produced based upon an approximate standard atmosphere model and different layers (one, two and three layers), which are simulated by linear function and quadratic function, respectively (Fig. 4). The layer boundaries of ''one layer'' are 4 and 6 km, those of ''two layers'' as well as ''three layers'' are 4, 4.8, 4.8, and 6, and 4, 4.5, 4.5, 5.2, 5.2, and 6, respectively, with a range resolution of 0.01 km. All of the layer intervals are set to be zero, which is beneficial for testing the denoising effects of different approaches on layer intervals. After simulating three original signals, a series of noises are mixed with them for creating the tested signals with different SNRs of 30DB, 40DB, and 50DB. Thus, nine artificial signals of different layer numbers and SNRs are formed. Figure 5 shows the results of one layer signal mixed with noises. It can be seen that in SNR = 30, the layers boundaries are very blur and many tiny peaks disturb the layer detection in the non-layer areas. ''One layer'' detection results using two methods are demonstrated in Table 2, where it can be observed that the results of our method are closer to 4 and 6 km than Xiong's. Especially, when SNR = 40DB, the error of our method is smaller 0.1 km than the one of Xiong's. It illustrates that the proposed method can improve the detection effect significantly in the case of fewer layers.
Because, layer detection is effected by complicated multiple layers and noise at the same time, the calculation errors in Tables 3 and 4 are larger than the ones in Table 2. It becomes difficult to detect the intersecting layer boundaries, such as the first-layer base and the secondlayer top. Using Xiong's approach intersecting layer boundaries are obtained with different values and have some larger deviations, while employing the proposed method the ones can be stabilized at single values regardless of the SNR degree. For example, in Table 3, when SNR = 30, the first-layer top 4.9 is larger than the secondlayer base 4.87 in the right column, which is unrealistic, but in the left column, the first-layer top 4.79 nearly equals the second-layer base 4.80. Moreover, the proposed approach can also effectively restrain the background noise for multiple layers. In Table 4, when SNR = 30, at 4.5 km, it can be found that the results of our method are very close to 4.5 km, the deviations are only 0.02 and 0.02 km, while the ones of Xiong's approach have large errors of 0.11 and 0.08 km.     Table 5, the proposed method is used to process the real signals. The obtained layer base heights are shown in Table 6, where the standard values detected by LUMSS are employed to verify two methods using LIDAR. One can observe that the detection heights using the proposed approach is closer to the standard value, compared with Xiong's method, and the mean of deviations with our method is less 45.6 m than the one with Xiong's method (the uncertainty is 10 m). Through the above four steps, the characteristics of LIDAR signals can be selected to determine the layer boundaries. Filter, iteration, and global correction precisely retain the features of trend and multiple scattering, and thus lead the LIDAR results to approach the outcomes of LUMSS.
The lowest blue line, the middle green line, the middle red line, the highest light, and the blue line represent the 1st, 2nd, 3rd, 4th, and 5th detected boundaries, respectively.
An actual case can be introduced to verify our method on multiple layers. LUMSS is also used to test the performance. With the LIDAR, two layers can be exhibited clearly with a range resolution of 0.01 km in Fig. 6. After the signal correction and iteration, final-detected boundaries can be obtained as 7.72, 8.70, 8.76, and 9.72 km. In Fig. 7, detection results via the number of iterations are showed, from which one can observe that the noise vibration acting on the first boundary (approach to 7.5 km) is weakened to a value that cannot affect the layer detection through 53 iterative filtering, and other boundaries also have the similar performances. It practically illustrates that iterative morphological filtering has a good robustness making detection results stabilized at a value with the iterative number increasing. Furthermore, at the boundaries of approximate 8.7 and 9.7 km, the fluctuations of detection results are weak, which displays good performances of the intervals detection of layers.
Detection results of real signal have been exhibited in Table 6 using the two methods. It clearly illustrates that the new approach outperforms Xiong's method in the gaps detection of layers (8.75 and 8.77 are much closer to 8.800 Furthermore, using Xiong's method, the detected layer tops can be higher than the ones using the proposed method, and closer to the standard values. It shows that the iterative morphological filtering restrains the noise better in the plat area of signals, and our method can determine the layer top at a lower altitude in the case of low-SNR.

Summary
A novel segmentation method based upon multiple scattering factors and iterative morphological filtering has been proposed to detect layer boundaries. This algorithm consists of four keys steps, which are the preliminary layer detection, iterative morphological filtering, multiple scattering factors calculation, and layer boundaries correction. Through the above four procedures, the waveform characteristics of the raw signal can be demonstrated in different scales, which is beneficial for combining with the multiple scattering factors to compute and revise the layer boundaries. Furthermore, as important parameters calculated, the integral extinction coefficients are employed to correct the preliminary detection results with the MC method. It is proved to be effective according to the simulations and experiments. The second-layer top 9.72 9.88 9.755