Development of ZJU High-Spectral-Resolution Lidar for Aerosol and Cloud: Extinction Retrieval

Abstract

The retrieval of the extinction coefficients of aerosols and clouds without assumptions is the most important advantage of the high-spectral-resolution lidar (HSRL). The standard method to retrieve the extinction coefficient from HSRL signals depends heavily on the signal-to-noise ratio (SNR). In this work, an iterative image reconstruction (IIR) method is proposed for the retrieval of the aerosol extinction coefficient based on HSRL data, this proposed method manages to minimize the difference between the reconstructed and raw signals based on reasonable estimates of the lidar ratio. To avoid the ill-posed solution, a regularization method is adopted to reconstruct the lidar signals in the IIR method. The results from Monte-Carlo (MC) simulations applying both standard and IIR methods are compared and these comparisons demonstrate that the extinction coefficient and the lidar ratio retrieved by the IIR method have smaller root mean square error (RMSE) and relative bias values than the standard method. A case study of measurements made by Zhejiang University (ZJU) HSRL is presented, and their results show that the IIR method not only obtains a finer structure of the aerosol layer under the condition of low SNR, but it is also able to retrieve more reasonable values of the lidar ratio.

1. Introduction

Accurate measurements of atmospheric aerosol properties are essential to the study of climate change, since aerosols affect the Earth’s climate system by scattering and absorbing solar radiation, as well as influencing the cloud properties [1]. The lidar technique is widely recognized as an important tool for remote sensing, which can provide the vertical distribution of aerosols and clouds [2]. However, the requirement for a priori assumptions of the lidar ratio in the retrieval of aerosol optical properties limits the accuracy of elastic backscatter lidars [3,4,5]. Raman Lidar permits the independent measurement of aerosol extinction coefficient and backscatter coefficient by detecting the elastic and Raman scattering signals, but it is mainly used during nighttime because of the large background radiation in daytime [6]. High-spectral-resolution lidar (HSRL) takes advantage of the different spectral broadening of Mie scattering and Cabannes-Brillouin scattering to discriminate aerosol returns from molecular returns with a spectral discriminator. Therefore, the aerosol backscatter coefficient and extinction coefficient can be retrieved directly in HSRL [6,7,8].

The intensity of lidar signals attenuate inversely as the square of the distance. The signal-to-noise (SNR) of the lidar signals is affected by the solar background light and the noise of the acquisition circuit [9,10]. The SNR decreases dramatically with the increase of the detection distance, and to retrieve the aerosol extinction coefficient accurately, the standard method in HSRL requires sufficient SNR [11]. The molecular channel in HSRL, which eliminates Mie scattering by the spectral discrimination filter, has a lower signal intensity. In addition, the standard method based on the numerical derivative makes the retrieval process extremely sensitive to the SNR, which greatly affects the accuracy of the aerosol extinction coefficient in HSRL [12].

In order to obtain aerosol optical properties accurately, various methods have been used to improve the SNR of the lidar signals, such as the application of wavelet transform (WT) filter, ensemble Kalman filter, empirical mode decomposition (EMD) approach, and adaptive filter [13,14,15,16]. These methods employ the low-pass filters at the lidar signals across the temporal and spatial axes, which is beneficial to the retrieval of optical properties. However, the low-pass filters improve the SNR of the lidar signal only to a certain extent. Since the retrieval of aerosol extinction coefficient is sensitive to noise, Pappalardo et al. [17] introduced different algorithms for aerosol extinction and backscatter retrieval in the European Aerosol Research Lidar Network (EARLINET). Esselborn et al. [18] used a first-order Savitzky-Golay (SG) filter to perform sliding window fitting on the optical thickness, which suppressed fluctuations in the aerosol extinction coefficient. However, these methods might introduce unnecessary biases when the observations of aerosols and clouds are not spatially uniform. Pornsawad et al. [12] performed an iterative regularization method to retrieve the aerosol extinction coefficient in Raman lidar, which produced a smoother extinction coefficient profile than the sliding linear fitting method under the condition of high SNR with a low temporal and spatial resolution. Marais et al. [19] proposed an approach to simultaneously denoise and invert aerosol backscatter and extinction coefficient in HSRL, which improved the retrieval under the condition of low SNR. Nevertheless, the proposed method is suitable for the photon-limited atmospheric lidar observations and the retrieval of aerosols is affected by the range of the clear-sky region.

In this work, a novel method based on an iterative image reconstruction (IIR) technique is proposed to retrieve the aerosol extinction coefficient in HSRL with improving accuracy under low SNR conditions. For a photomultiplier tube (PMT) operating in the analog detection regime, a previous study had proved that the noise of the lidar signals is the compounded Poisson noise [20]. As the atmospheric conditions are relatively stable in a short period, lidar signals reflecting similar atmospheric conditions are related. Taking advantage of these properties of lidar signals, the IIR method reconstructs the “image” of the “feature” of the lidar signals on the basis of reasonable estimates of the lidar ratio by minimizing the difference between the reconstructed signals with the observations under an appropriate assumption of signal noises. The “feature” is identified with a scattering ratio profile greater than one standard deviation from the expected clear-sky [21]. Using an iterative regularization method to carry out the best optimization estimates of the lidar ratio, the IIR method transforms the ill-posed problem of the numerical derivative into a well-posed problem of the image reconstruction. The IIR method has the potential to retrieve aerosol optical properties reasonably and automatically as well as significantly improve the detection capability of the HSRL under the condition of low SNR.

The remainder of this paper is organized as follows. Section 2 provides an introduction of HSRL principles, the lidar noise distribution model and the methods of aerosol extinction coefficient retrieval, including the standard method and the IIR method. Section 3 shows the Monte-Carlo (MC) simulations for noise analysis. Section 4 presents a case study from the ambient measurements made by the Zhejiang University (ZJU) HSRL system to demonstrate the capabilities of the proposed approach. Section 5 discusses the results obtained. The paper ends with the conclusions in Section 6.

2. Methodology

2.1. The Basic HSRL Theory

The HSRL technique employs a spectral discrimination filter to discriminate aerosol returns from molecular returns and measures the aerosol extinction coefficient and backscatter coefficient independently. Mie scattering from aerosol particles has nearly the same spectral distribution as that of the incident laser pulse, while the molecular Cabannes-Brillouin scattering is broadened by a few GHz due to the effect of Doppler broadening [22,23]. The general HSRL system we considered in this study with only two channels for the sake of simplicity is shown in Figure 1. The laser beam is collimated and expanded by an expander and then passes through the atmosphere. The backscattered signal carrying atmospheric information is received by the telescope and divided into two channels by a beam splitter, which are referred to as the combined channel and molecular channel. With the combined channel, the backscattered lidar signals of both aerosols and molecules are measured. The molecular channel is equipped with a spectral discrimination filter to eliminate the aerosol returns and pass the wings of the molecular spectrum.

The models of the combined channel and molecular channel could be derived from the single-scatter lidar equation [18,24]. For ease of discussion in this study, all parameters will be expressed in a matrix form. Raw data consists of $N$ range bins (row axis) and $K$ profiles (column axis). The lidar signals can be presented as:

$$P_C=C_C\cdot ({\beta }_a+{\beta }_m)\cdot e^{-2\tau }+B_C$$

(1)

$$P_M=C_M\cdot (T_a\cdot {\beta }_a+T_m\cdot {\beta }_m)\cdot e^{-2\tau }+B_M$$

(2)

where the symbol $·$ indicates that the matrices are multiplied pointwise. $P\in {ℝ}^{N\times K}$ represents the received power; the subscripts C and M represent the combined channel and the molecular channel, respectively; $C\in {ℝ}^{N\times K}$ is the calibration matrix relative to the system constants, such as the optical efficiency, the sensitivity of detector, the overlap function, the detection range, etc.; $\beta \in {ℝ}^{N\times K}$ is the backscatter coefficient; the subscripts $a$ and $m$ represent aerosols and molecules, respectively; $T_m$ and $T_a$ are the transmittances of the spectral discrimination filter to molecular signal and aerosol signal, respectively; $B\in {ℝ}^{N\times K}$ represents the background noise; the $\tau (z)\in {ℝ}^{N\times K}$ describes the atmospheric optical thickness and can be expanded to:

$$\tau =Q[{\alpha }_a+{\alpha }_m]$$

(3)

where $Q\in {ℝ}^{N\times K}$ represents the integral matrix, which is a lower triangular matrix; $\alpha \in {ℝ}^{N\times K}$ is the extinction coefficient.

The only unknowns in this model are the aerosol backscatter coefficient ${\beta }_a$ and the aerosol extinction coefficient ${\alpha }_a$, and the rest are pre-generated calibration matrices. The ‘one-equation, two-unknowns’ problem could be solved by combining the combined channel and the molecular channel. The detailed discussion about how to retrieve the aerosol backscatter coefficient ${\beta }_a$ and the aerosol extinction coefficient ${\alpha }_a$ will be presented in Section 2.3 and Section 2.4

2.2. The Noise Model in HSRL

It is well known that the noise of single-shot lidar signals can be modeled by a Poisson probability mass function (PMF) [19,25]. Let $Y_i\in {ℝ}^{N\times K}$ be the Poisson noisy observations of lidar signals $P_i$. The Poisson PMF of $Y_i$ is defined as:

$$P(y_{n,k}|P_{i,n,k})={\displaystyle \prod _{n=1}^N{\displaystyle \prod _{k=1}^K\frac{{(P_{i,n,k})}^{y_{n,k}}}{y_{n,k}!}\exp (-P_{i,n,k})}}$$

(4)

where y_n,k is the possible observation of n-th row axis and k-th column.

The multiplication process of the PMT changes the primary distribution of incident photons and it can be modeled by a compounded Poisson PMF [20]. However, the compounded Poisson distribution lacks concise mathematical expression. To obtain sufficient SNR, the lidar signals generally would be processed by temporal averaging and the Gaussian PMF can be used as an approximation [20,25]. The Gaussian PMF of $Y_i$ is defined as:

$$P(y_{n,k}|P_{i,n,k})={\displaystyle \prod _{n=1}^N{\displaystyle \prod _{k=1}^K\frac{1}{\sqrt{2\pi }{\sigma }_{n,k}}\exp (-\frac{(y_{n,k}-P_{i,n,k}){)}^2}{2{\sigma }_{n,k}{}^2})}}$$

(5)

where ${\sigma }_{n,k}$ is the standard deviation of n-th row axis and k-th column of lidar signals.

Although the noise of lidar observations can be reduced by averaging, it only helps to a certain extent. Even in the case of tiny fluctuations of lidar signals, the aerosol extinction coefficient retrieved by the standard method would be severely affected [26]. In the data processing of noisy images, the image reconstruction technique is widely employed in many applications, such as medical imaging systems, confocal microscopy and interferometry [27,28]. It provides the possibility to retrieve optical properties, since the lidar observations could be regarded as a noisy two-dimensional image with temporal axis and spatial axis. Detailed instructions for the retrieval of optical properties based on the image reconstruction technique will be presented in Section 2.4.

2.3. Standard Method

Ideally, where there is no noise and the calibration matrices are fully known, the aerosol backscatter coefficient ${\beta }_a$ and the aerosol extinction coefficient ${\alpha }_a$ can be retrieved algebraically [29]. From the HSRL models Equations (1) and (2), the aerosol backscatter coefficient can be obtained from the observations $Y_C$ and $Y_M$ as:

$${\beta }_a=\frac{(T_m-K)\cdot {\beta }_m}{K-T_a}$$

(6)

$$K=\frac{(Y_M-B_M)}{(Y_C-B_C)}\cdot \frac{C_C}{C_M}$$

(7)

Substituting Equation (6) into Equation (2), we can obtain the optical thickness from the HSRL:

$$\tau =-\frac{1}{2}\cdot \ln (\frac{(K-T_a)\cdot (Y_M-B_M)}{(T_m-T_a){\beta }_{m}K})$$

(8)

Based on Equation (3), the aerosol extinction coefficient can be readily obtained using:

$${\alpha }_a=Q^T(\tau )-a_m$$

(9)

where the matrix $Q^T$ is the differential operator, which is the inverse matrix of the integral matrix $Q$. Then, the intensive parameter, lidar ratio $S_a$, can be simply derived by the ratio of extinction to backscatter:

$$S_a=\frac{{\alpha }_a}{{\beta }_a}$$

(10)

However, when the SNR of the lidar signals is insufficient, an averaging function is generally applied to the noisy observations $Y_i$ to reduce the noise variance. Furthermore, especially for the retrieval of the aerosol extinction coefficient, a denoising filter, such as the SG filter, should also be applied on the spatial axis of the optical thickness. This is equivalent to a linear regression on a specified number of data points within a moving window, and the slope of the fitted straight line in each point is taken as the value of the derivative [18]. In the standard method, as in this present study, the SG filter would be an optional denoising filter performed on the lidar signals and the optical thickness.

2.4. Iterative Image Reconstruction (IIR) Method

It should be noticed that the aerosol backscatter coefficient could be used to constrain the retrieval of the aerosol extinction coefficient, since there is a linear relation (i.e., lidar ratio) between the aerosol backscatter coefficient and extinction coefficient. The range of the lidar ratio is about 0-100 typically. Once the correct estimate of the lidar ratio is obtained, the aerosol extinction coefficient could be algebraically computed by the product of the aerosol backscatter coefficient and the lidar ratio. In this work, the IIR method is proposed to reconstruct the “image” of the “feature” of the lidar signals on the basis of a reasonable estimates of the lidar ratio under appropriate assumption of signals noise. In this study, the Gaussian PMF is utilized when the observations are from the analog mode lidar measurements. The characteristics of the proposed method can be summarized into two points: (1) considering the relation of the lidar signals and utilizing the feature detection method to obtain the “feature” of the “image”, (2) transforming the ill-posed problem of the numerical derivative into a well-posed one of the estimates of the lidar ratio. This “feature” is identified with the Feature Detection Algorithm from a previous work to represent the data associated with aerosols and clouds [21]. The region where the scattering ratio profile is greater than one standard deviation from the expected clear-sky is determined as the “feature.” The reconstruction of the “image” of the “feature” must fit with the observations in HSRL to ensure the accuracy of the optical property retrieval. This process would be achieved by using an iterative optimization method on the basis of an adaption of the sparse Poisson image reconstruction algorithm (SPIRAL) [27,30].

The retrieval of the backscatter coefficient in the IIR method is the same as the standard method. To estimate the lidar ratio, the molecular channel in the lidar signal from Equation (2) can be expressed as a function of the lidar ratio:

$$g(S_a)=C_G\cdot \exp (-2Q[{\beta }_a\cdot S_a+{\beta }_m\cdot S_m])+B_M$$

(11)

$$C_G=C_M(T_a\cdot {\beta }_a+T_m\cdot {\beta }_m)$$

(12)

For the present study the molecular channel noise is considered as a Gaussian noise for general cases [20]. It would also be feasible to adopt the Poisson noise in a similar way, and detailed information can be seen in [19]. The image reconstruction technique is used to estimate the lidar ratio by minimizing the negative log-likelihood of Equation (5). The negative log of Equation (5) can be expressed as:

$$l(S_a,Y_M)=\{\begin{array}{l}{1}_N^T[\frac{{(Y_M-g(S_a))}^2}{2{\sigma }_{}{}^2}-\ln \frac{1}{\sqrt{2\pi }{\sigma }_{}}]1_K^{},feature\\ 0\begin{array}{cccc}\begin{array}{cccc}\begin{array}{cc}\begin{array}{ccc}& & \end{array}& \end{array}& & & \end{array}& & & ,\end{array}clear-sky\end{array}$$

(13)

where the $\sigma \in {ℝ}^{N\times K}$ is the matrix of the standard deviation and the elements of vectors 1_N and 1_K are ones. Throughout this paper, Equation (13) is referred to as the loss function. A penalty function is added to the loss function to regularize the estimates of the lidar ratio $S_a$. Since the lidar ratio of an aerosol layer is relatively stable, it is reasonable to consider that the total variation (TV) seminorm is an adequate penalty function to constrain the smoothness of the lidar ratio $S_a$. The specific TV seminorm used in this paper is:

$$pen(S_a)={\Vert S_a\Vert }_{TV}={\displaystyle \sum _{\begin{array}{l}n=1,\\ k=1\end{array}}^{\begin{array}{l}N-1,\\ K\end{array}}{F}_{n,k}\left|S_a{}_{,n,k}-S_a{}_{,n+1,k}\right|}+{\displaystyle \sum _{\begin{array}{l}n=1,\\ k=1\end{array}}^{\begin{array}{l}N,\\ K-1\end{array}}{F}_{n,k}\left|S_a{}_{,n,k}-S_a{}_{,n,k+1}\right|}$$

(14)

where $F\in {ℝ}^{N\times K}$ is 1 for “feature” and 0 for clear-sky; the subscripts are indices to the rows and columns of $S_a$ and $F$. Then, the estimates of the lidar ratio is turned into a regularization problem:

$$\stackrel{-}{S_a}=\underset{S_a\in U}{\arg \min }\{l(S_a,Y_M)+\lambda \cdot {\Vert S_a\Vert }_{TV}\}$$

(15)

where $\lambda$ is referred to as the regularization parameter; $U$ is the constrain set of the lidar ratio $S_a$.

This form of regularization problem could be solved using SPIRAL with the constraint of TV seminorm [27]. This iterative optimization algorithm is well established and numerically stable. The principle of the SPIRAL will not be repeated in this work. Instead, an adaption of the SPIRAL is used to implement the data process, where the gradient matrix of the loss function is expressed as:

$$\begin{array}{l}\nabla l(S_a,Y_M)=1_N^T[2(Y_M-\mathrm{g}(S_a))\frac{\partial g(S_a)}{\partial S_a}]1_K^{}\\ =2\frac{(Y_M-\mathrm{g}(S_a\left)\right){\beta }_a}{{\sigma }^2}\cdot Q^T[g(S_a)]\end{array}$$

(16)

The IIR method still posed some problems: the initial value of the lidar ratio $S_{initial}$ when the first iteration starts running and the determination of the regularization parameter $\lambda$. For the first iteration, the average lidar ratio of the aerosol layer obtained from the standard method is employed as an initial guess. Regularization parameter $\lambda$ can be found by manually tuning parameters and adjusting while it is more practical to choose it automatically. Various methods have been presented to choose the regularization parameter $\lambda$ [12,27,30,31]. In this work, the regularization factor $\lambda$ is selected using cross-validation [30]. The observations $Y_i$ is separated into two groups $Y^{training}$ and $Y^{validation}$ via hold-out, where $Y_i=Y^{training}+Y^{validation}$. The regularization factor is selected from a set of non-negative numbers, for example $\lambda$ = {10⁻², 10^−1.8, 10^−1.6, …, 10} [31]. For every $\lambda$, the IIR method is used with training set $Y^{training}$ and then the negative log-likelihood of the validation set $Y^{validation}$ is calculated. The regularization factor $\lambda$ that minimizes the negative log-likelihood of the $Y^{validation}$ is chosen. Figure 2 shows the flow diagram for the retrieval of the optical properties in the IIR method, where $W_j$ is the retrieval window of profile $j$ and wins is the half width of the window. The different region of clear-sky area we utilized may change the results of the “feature” and the retrieval based on the feature detection is easy to reproduce once the “feature” is determined. Therefore, combining the feature detection algorithm with the SPIRAL of Gaussian model, the IIR method has the potential to obtain the optical properties of the “feature” reasonably and automatically. The lidar ratio of each profile could be retrieved through the corresponding retrieval window. Hereafter, the final result is obtained by averaging the lidar ratio in every $W_j$. More details about SPIRAL can be seen in [32].

3. Error Analysis Based on Monte-Carlo Simulations

3.1. The Specifications of MC Simulations

In order to verify the feasibility of the IIR method, MC simulations are employed, which allow us to compare the results obtained by different methods with the true values input into the simulations. For this section, the synthetic dataset was created on the basis of real observations of cirrus clouds between 8.5–10 km with 109 profiles and used to estimate the performance of different methods. The measured backscatter coefficient of the cirrus clouds is employed and the extinction coefficient of the cirrus clouds is created with the simulated lidar ratio. The lidar ratio of the cirrus clouds is modeled in Gaussian with a value of 26.28 ± 3.63 sr, while it is set to $8\pi /3$ sr for the clear-sky region [33]. The modeled inputs are presented in Figure 3.

The lidar signals are generated from forward MC simulations based on the HSRL model mentioned in Section 2. Then, a Gaussian distributed random noise with standard deviation equal to the square root of the signals is incorporated into the lidar signals. For the simulation, the scene represents a real dataset with low SNR where the resolution is 7.5 m by 10 s to compare the performance of two methods. The system specifications employed by MC simulations are listed in

Table 1. System specifications employed by Monte-Carlo (MC) simulations.

Specifications	Value
Laser wavelength (nm)	532
Laser energy (mJ)	200
Diameter of primary telescope (mm)	280
Interference filter bandpass (nm)	0.3
Range resolution (m)	7.5
Temporal resolution (s)	10
Molecular transmittance of spectral discrimination filter $T_m$	0.19
Aerosol transmittance of spectral discrimination filter $T_a$ (e-12)	2.52
Type of noise	Gaussian

, the $T_m$ is a constant for simplicity, and the $T_a$ is obtained by the measured length and temperature of the iodine molecular absorption filter.

3.2. Comparison Between the Standard Method and IIR Method

In this subsection, the performance of the standard and the IIR method is compared against each other. The root mean square error (RMSE) is used to indicate the performance of different methods, along with the relative bias. The RMSE can be defined as:

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^K{(M_{i,j}^{true}-M_{i,j}^{retrieval})}^2}}}{N\times K}}$$

(17)

where $M\in \{{\alpha }_a,S_a\}\in {ℝ}_{}^{N\times K}$; the superscript “true” and “retrieval” represent the true value and the retrieval value, respectively.

Hereafter, we only focus on the extinction coefficient retrieval of aerosols and clouds (i.e., features) identified by a feature detection method since the accurate measurements of the clear-sky region is beyond the capabilities of the HSRL. The SG filter is used as the low-pass filter in the standard method for both temporal and spatial axes. The SG filters of both temporal axis and spatial axis are first-order polynomials, whose window sizes are nine profiles and nine altitude bins, respectively. The moving SG window for optical thickness is 71 altitude bins. These SG filter parameters are selected based on the minimization of the extinction coefficient RMSE. The window of the IIR method is nine profiles as a balance of efficiency and accuracy.

Figure 4 shows the results obtained by the standard method (left panel) and the IIR method (right panel). The true values are shown in Figure 3. Figure 4a,c and e are the results obtained by the standard method, which represent the extinction coefficient, the absolute errors of the extinction coefficient of the true values and the retrieval values, and the lidar ratio, respectively. Figure 4b,d and f are obtained by the IIR method. The dashed lines in f are selected for further analysis, whose profile numbers are 19, 48, 78, and 99, respectively. The extinction coefficient obtained by the standard method is smoother than that obtained by the IIR method. However, the absolute errors of the extinction coefficient from the IIR method (Figure 4d) are smaller than those obtained by the standard method (Figure 4c). Comparing to the input values in Figure 3b, it could be noticed that the structure of the extinction coefficient from the IIR method is closer to the input values than that from the standard method. Furthermore, the comparisons of the true values in Figure 3c and the results obtained by two methods indicate that the IIR method is able to achieve better results of the lidar ratio by the image reconstruction technique. The standard method would require more accumulated profiles and altitude bins (i.e., a higher SNR but a lower spatial and temporal resolution) to reduce the variance of the noise to a sufficient level until the slope of the optical thickness are non-negative (see Equation (9)).

The reason why the IIR method does not suffer from the same problem (negative slope of the optical depth) as the standard method is because it does not try to retrieve the aerosol extinction coefficient algebraically. The IIR method estimates the lidar ratio with the loss function Equation (13) in conjunction with the HSRL model to “fit” with the noisy observations. The lidar ratio is always limited to the bounds of 0–100, as it should be physically. Figure 5a,c show four profiles of the extinction coefficient and the lidar ratio obtained by different methods, respectively. The locations of four profiles are shown in Figure 4d. The lines in black are the true values, while the blue and red lines represent the results of the “feature” obtained by the standard and the IIR method, respectively. What we concerned most is the region containing a “feature,” i.e., dense aerosols and clouds. Therefore, a feature detection method is utilized to exclude the clear-sky region [21]. It is obvious that the structures of the extinction coefficient and the lidar ratio obtained by the IIR method are closer to the true values than those obtained by the standard method. The relative biases of the retrieval of the extinction coefficient and the lidar ratio are presented in Figure 5b,d, respectively. It could be found that the extinction coefficient and the lidar ratio obtained by the IIR method have smaller biases than those from the standard method. The relative bias of the extinction obtained by the IIR method is mostly within 20%, while that obtained by the standard method is much bigger. From the

Table 2. Error analysis of optical properties. RMSE: root mean square error.

Extinction	RMSE	Relative bias
Standard method	0.015	0.285
IIR method	0.009	0.170
IIR method (Poisson)	0.010	0.176
Lidar Ratio	RMSE	Relative bias
Standard method	2.451	0.226
IIR method	0.870	0.085
IIR method (Poisson)	0.950	0.105

, the mean relative bias of the extinction coefficient obtained by the IIR method is 17.0%, while it is up to 28.5% with the standard method. The maximum relative bias of the lidar ratio obtained by the standard method is up to 100%. It indicates that the standard method is not suitable for the further aerosol analysis under the condition of low SNR, which is sensitive to the extinction coefficient errors, such as the feature classification, the retrieval of microphysical properties, etc. [34]. The relative biases of profile 48 (the second profile) from the IIR method at the start of “feature” may be caused by the system constant and the initial value of the lidar ratio. The mean relative bias of the lidar ratio obtained by the IIR method is 8.5%, while it is up to 22.6% with the standard method. The green lines in Figure 5 represent the results obtained by the IIR method based on the Poisson model. The comparisons between the red lines and the green lines indicate that the results of the Gaussian model are similar to those of the Poisson model and the relative bias of the Gaussian model is smaller than that of the Poisson model.

Results from MC simulations of the IIR method are compared with the standard method. From the error analysis in Table 2, the RMSE of the lidar ratio and the extinction coefficient obtained by the standard method are 2.451 and 0.015, respectively. The RMSE of the lidar ratio and the extinction coefficient obtained by the IIR method are 0.870 and 0.009, respectively. Both are much smaller than those from the standard method. Therefore, it is obvious that the estimates of the extinction coefficient and the lidar ratio based on the standard method are inaccurate relative to the IIR method. The more accurate optical properties from the IIR method are more suitable for the further aerosol analysis. The error analysis from the IIR method based on the Poisson model is also provided in Table 2. The results from two noisy models are similar and the errors of the IIR method are slightly smaller than those from the IIR method based on the Poisson model.

Figure 6a shows the comparisons of the lidar ratio of the “feature” obtained by the standard method, the IIR method, and the IIR method without feature detection, respectively. The profile number is 78. It is obvious that the results from the IIR method are better and closer to the true value. Figure 6b illustrates the results obtained by the IIR method without feature detection have larger relative bias compared with the IIR method. As the results of the IIR method are determined by all pixels of the “image,” the results of the “feature” are severely affected by the clear-sky region when using a Gaussian model (seen in Equation (13)). The results indicate that the weighted TV constraint from the feature detection method can improve the results and the feature detection method is utilized to eliminate the influence of the clear-sky area in the results of the “feature.”

4. Results

A case study of ambient measurements is presented in this section. This case study will also demonstrate that the optical properties retrieved by the IIR method are more stable compared to the standard method. In this case study, observations of the ZJU HSRL at 532 nm were used [21]. We used daytime scenes with a high vertical resolution of 7.5 m and a temporal resolution of 120 s to compare the IIR method and the standard method, since the noisy observations from daytime reflect the performance of different methods under the condition of lower SNR. Figure 7a–c represent the combined channel signals, the molecular channel signals, and the backscatter coefficient measured from ZJU HSRL (30°15′ N, 120°7′ E) between 08:25–15:05 UTC+8, respectively. The detection range is above 1 km to skip the incomplete overlap range between the divergence angle of the laser beam and the field-of-view (FOV) of the receiver, to avoid the possible errors of the correction of incomplete overlap function. The retrieval of the extinction coefficient is extremely sensitive to the correction of the overlap function in HSRL, as the optical thickness is directly affected by the overlap function (seen in Section 2.3, the calibration matrix $C$ includes overlap function) [35].

Figure 8a,b show the retrieval of the aerosol extinction coefficient and the lidar ratio based on the standard method using an SG filter. The sizes on the temporal axis and spatial axis are five profiles and nine altitude bins, respectively, while the size of the moving window for the optical thickness is 31. Figure 8c,d correspond to the aerosol extinction coefficient and the lidar ratio on the basis of the IIR method. The window for retrieval is selected as nine profiles. From Figure 8, it is clear that the IIR method is able to retrieve the smoother optical properties than the standard method. Furthermore, the optical properties of the residual aerosol with low aerosol loading around 14:00 UTC+8 could not be retrieved by the standard method due to the low SNR. However, reasonable results for the aerosol extinction coefficient and the lidar ratio can still be obtained by the IIR method. Furthermore, the extinction coefficient obtained by the IIR method has a similar structure with the observations in Figure 7 comparing that from the standard method during 09:00–14:00 UTC+8.

The reconstructed lidar signal and raw signal presented in Figure 9a demonstrate the effectiveness of the algorithm further. The IIR method reconstructs the “image” of lidar signals on the basis of reasonable estimates of the lidar ratio by minimizing the difference between the reconstructed signals with the observations under an appropriate assumption of Gaussian noise model. With the help of the pre-estimated backscatter coefficient, the IIR method obtains a better result. Figure 9b shows two single profiles of the lidar ratio retrieved by the standard method and the IIR method at 08:40 UTC+8, respectively. The average values of the lidar ratio obtained by the standard method and IIR method are 38.60 ± 12.18 sr and 44.90 ± 6.80 sr, respectively. Nevertheless, the lidar ratio obtained by the standard method is drastically changeable via the altitude, while the lidar ratio obtained by the IIR method is much more stable.

5. Discussion

The complex condition of the lidar ratio varying with altitude is considered in MC simulations to evaluate the effectiveness of the proposed method. The results from MC simulations indicate that the structure of the extinction coefficient from the IIR method is closer to the input values than the extinction coefficient from the standard method. Furthermore, smaller values of the RMSE and bias of optical properties are obtained by the IIR method. The results demonstrate that the influence of noise in the estimates of optical properties is suppressed by the IIR method, and the results are more accurate than those from the standard method. Additionally, the error in the initial value of the lidar ratio may cause errors at the start region but the results converge and close to the real values rapidly as the retrieval range increases. The comparison between the IIR method with two models illustrate that the Gaussian model performs well on the extinction retrieval when the observations are from the analog mode lidar measurements. The feature detection method is utilized to eliminate the influence of the clear-sky area in the results of the “feature” and promote the retrieval results.

The results from the measurements of the ZJU HSRL, applying both standard and IIR methods, are compared. The IIR method reconstructs the “image” of lidar signals on the basis of reasonable estimates of the lidar ratio by minimizing the difference between the reconstructed signals with the observations, and the comparisons indicate that the reconstructed signal agrees with the raw signal. The average values of the lidar ratio obtained by the standard method and the IIR method are consistent and reasonable [36,37]. The lidar ratio from the standard method oscillates significantly with the height, and it is difficult to obtain reasonable values for the lidar ratio in the region of low aerosol loading with low SNR. However, the results from the IIR method demonstrate that the IIR method not only obtains a finer structure of the aerosol layer under the condition of low SNR, but is also able to retrieve more reasonable values of the lidar ratio, as the aerosol below the planetary boundary layer is relatively stable and generally well-mixed [38,39].

Compared to the standard method, the method proposed in this study performs well for the retrieval of the optical properties under low SNR conditions and obtains a finer structure of the aerosol layer. However, the limitation of the method is inevitable. First of all, the proposed method should be applied in more complicated scenes, such as when the boundary layer is not well-mixed, as well as the comparison of the results between the lidar and the in situ observations to verify the effectiveness of the method in future studies. Secondly, the method requires a prior initial value of the lidar ratio and the system constant, which plays important roles in the retrieval progress. The average lidar ratio of the aerosol layer obtained from the standard method might be an option for an initial guess, but a better method should be developed to decrease the error resulted from the initial values in future work. In addition, the study of which model is the best candidate to infer the lidar ratio could be conducted in the future, and the IIR method can be applied to improve observations with low SNR from airborne or satellite lidar.

6. Conclusions

An iterative image reconstruction method is proposed for the retrieval of aerosol extinction coefficients based on HSRL data. The IIR method reconstructs the “image” of “feature” of lidar signals on the basis of reasonable estimates of the lidar ratio by minimizing the difference between the reconstructed signals with the observations under an appropriate assumption of signal noises. Compared with the previous studies, the relation of the lidar signals is considered and the “feature” of the “image” is utilized to retrieve the lidar ratio and extinction coefficient in the IIR method. As the numerical derivative makes the retrieval process extremely sensitive to the SNR, the proposed method transforms the ill-posed problem of the extinction coefficient retrieval into a well-posed one by the iterative regularization method, which obtains more stable results. The feature detection method is utilized to eliminate the influence of the clear-sky area to the results of the “feature.” Combining the feature detection algorithm with the SPIRAL of Gaussian model, the IIR method has the potential to retrieve aerosol optical properties reasonably and automatically.

The proposed method was compared with the standard method by both synthetic data and real experimental data. The comparison based on the MC simulations between the standard method and the IIR method demonstrates that our proposed method yields the aerosol extinction coefficient and the lidar ratio with smaller RMSE and relative bias. The comparison between the IIR method and the IIR method based on the Poisson model illustrates that the Gaussian model performs well on the extinction retrieval when the observations are from the analog mode lidar measurements. The RMSE of the lidar ratio and aerosol extinction coefficient obtained by the IIR method are 0.870 and 0.009, respectively. The mean relative bias of the lidar ratio obtained by the IIR method is 8.5%, while it is up to 22.6% for the standard method. Finally, a case study from the ambient measurements by the ZJU HSRL system is presented. The results show that the IIR method not only obtain a finer structure of the aerosol layer under the condition of low SNR, but also is able to retrieve more reasonable values of the lidar ratio.

We will focus on the application of the method in more complicated scenes and in airborne or satellite lidar in our future work. The improvements of the choice of the initial value of lidar ratio and system constant should be studied in the future studies. Furthermore, the determination the best candidate of the model to retrieve the optical properties should be studied in future work.