Remote sensing of raindrop size distribution using the coherent Doppler lidar

The coherent Doppler wind lidar (CDL) shows capability in precipitation detection. Retrieval of the raindrop size distribution (DSD) using CDL is still challenging work, as both accurate backscattering cross section at the working wavelength and reflectivity spectrum of raindrop are required. Firstly, the Mie theory and the vectorial complex ray model (VCRM) are applied to calculate backscattering cross section for small spheric raindrops and large oblate raindrops, respectively. Secondly, an iterative deconvolution method is proposed to separate the reflectivity spectrum of raindrop from the lidar power spectrum, which is a superposition of rain and aerosol components. An accompanying aerosol signal model considering the effect of temporal window, from the same height and time, is used to improve the accuracy and robustness of the iteration. In experiment, a co-located micro rain radar (MRR) is used for comparison. Good agreements are obtained despite tremendous differences in wavelength and scattering characteristics. As an example, at 600 m height, the R2 of linear fitting to the mean rain velocity and mean raindrop diameter between CDL and MRR are 0.96 and 0.93, respectively. © 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Raindrop size distribution is the fundamental microphysical property of precipitation. Its vertical profiles are helpful in revealing the formation and evolution processes of precipitation, such as collision, coalescence, melting, breakup, and evaporation [1][2][3][4][5]. Besides, in weather radar applications, the information on raindrop size distribution (DSD) can be used to calibrate rainfall rate (R) estimation which is derived from the radar reflectivity (Z) based on the equation commonly known as Z-R relation [6].
Ground-based disdrometers can provide only point measurements of DSD [7,8]. It has been suggested to use Doppler radar at vertical incidence to derive the DSD profiles via the relation between terminal fall velocity and drop size [9]. The K-band micro rain radar is then developed under the assumption of zero vertical wind [10]. Some researchers reported simultaneously measurement of echoes from Bragg scattering caused by inhomogeneities of the atmospheric refractive index, and from Rayleigh scattering by hydrometeors, using very high frequency (VHF) and ultra high frequency (UHF) radar wind profilers [11][12][13][14][15]. The information of ambient air motion makes it possible to remove the effect of vertical wind velocity and turbulence on the derived DSD profiles by a deconvolution process [16,17].
At the wavelength regime of lidar applications, echoes from aerosols (Mie scattering) and raindrops (geometrical scattering) can both be detected during precipitation events [18][19][20], which provide a potential to derive the DSD by lidars. Two-color lidar measurements have been used to estimate the DSD of drizzle [21][22][23]. A combination of Doppler lidar and cloud radar was also applied to estimate the DSD based on the different rain velocities measured by two instruments [24]. Aoki et al. introduced the method used in the radar wind profilers to the coherent Doppler lidar and derived the DSDs using an empirical backscattering efficiency function [25]. An accurate backscattering cross section is critical in retrieving DSDs. The main difference in precipitation detection between radar and lidar is their operation wavelength and corresponding scattering characteristics. For long-wave radars, Rayleigh scattering is applied, which implies a 6th power relation between the backscattering cross section and drop size. While in lidars, geometric optics approximation (external reflection) is often used, implying a dependency on the square of the raindrop size [24,26,27]. The application of theoretical backscattering cross section in DSD retrieval using coherent Doppler lidar and comparison of the results between lidar and radar have not been reported.
In this study, an all-fiber coherent Doppler lidar is applied for precipitation detection with vertical stare mode. The pulse width is set as 600 ns to achieve relatively narrow linewidth and thus accurate spectrum estimate. An aerosol signal model is built by convolution of a Gaussian-shaped turbulence spectrum and sinc-like window spectrum, improving the accuracy and robustness of the deconvolution process. As backscattering cross section is considered, it is calculated by using Mie theory [28] and vectorial complex ray model (VCRM) [29] for small spheric drops and large oblate drops, respectively. The DSD results are compared with co-located MRR data.

Instruments
The two remote sensing instruments: the CDL and the MRR are installed on the campus of the University of Science and Technology of China (31.83°N, 117.25°E). The key parameters are summarized in Table 1. Besides, a ground-based optical disdrometer (second-generation particle size and velocity, Parsivel-2 [30]) is deployed for comparative measurements. All instruments are collocated within a diameter of 15 m. The compact all-fiber CDL operates with an eye-safe wavelength of 1.5 µm. It has been applied for researches on the atmospheric boundary layer height [31], gravity waves [32], and turbulence [33]. During the precipitation experiment, the CDL is set to operate at vertical stare mode with a temporal resolution of 1 second. The detailed information about the system is presented in previous works [20,34,35].
The Micro Rain Radar is a vertically pointing frequency modulated continuous wave (FMCW) radar (MRR-2, METEK Germany). It operates with microwave at a frequency of 24.23 GHz with modulation of 0.5-15 MHz according to the height resolution. The wavelength is selected after trade-offs between the sensitivity to small raindrops and the rain attenuation effect [10]. In this study, the temporal and spatial resolutions of the MRR are set as 1 min and 100 m, respectively.

Principle
The directly detected information of rain by Doppler lidar or radar is the Doppler-shifted reflectivity spectrum, which is the superposition of the signal from individual raindrops at various falling velocities determined by the size. Considering the backscattering cross section as the weight, the reflectivity spectrum of rain in the velocity domain can be expressed as [12] where v and D are the Doppler velocity and the equivalent spherical diameter, respectively. N(D) is the raindrop size distribution to be solved in this work; S rain (v) is retrieved from the Doppler power spectrum measured by the CDL, and σ bk (D) is calculated by theoretical models, as introduced in the following sections. dD/dv is derived from the relationship between the terminal falling velocity and raindrop size. The relationship is given by an empirical equation [9] v where v and D are in units of m/s and mm, respectively. ρ and ρ 0 are the respective air densities at the observation height and sea level. The parameter c 0 in Eq. (1) is a calibration factor depending on the instrument parameters and path attenuation. It can be obtained by solving the lidar equation accurately. However, due to the strong attenuation of light by raindrops, this approach is difficult. In this work, we use the radar reflectivity Z outputted by the MRR as reference information to calibrate c 0 . Nevertheless, the shape of DSD and the derived parameters, such as the mean rain velocity and mean diameter of raindrops, are retrieved by CDL independently.

Backscattering cross section
The scattering properties of raindrops depend mainly on the wavelength-dependent refractive index and geometry. An accurate backscattering cross section is critical in retrieving DSDs as expressed in Eq. (1). The Mie theory [28] offers a rigorous solution to the scattering of electromagnetic waves from a homogenous sphere. But the sphere assumption is only valid for the raindrops of diameter less than 1.0 mm. Otherwise, the effect of non-spheric shape should be considered [36]. On the one hand, we use the MiePlot tool developed by Philip Laven [37] to calculate the backscattering cross section of spheric drops. Considering raindrops with large size parameters πD/λ, where λ is the working wavelength, MiePlot gives stable results. On the other hand, considering the large non-spheric raindrops, a ray-tracing model called VCRM developed by Kuan Fang Ren [29] is used. Comparing to numerical calculation methods, such as T-matrix and Finite-Difference Time-Domain (FDTD), VCRM is not limited to relatively small size parameters [38]. By introducing the property of the wavefront in the geometrical optics model, the VCRM can calculate precisely the interaction of light with an object of smooth surface. It has been verified elsewhere by comparison with other models and oblate droplet experiments [38,39].
The complex refractive index for liquid water is required in the calculation. At the wavelength of 1.5 µm, the value is 1.32 + 0.000135j where the imaginary part means absorption of the light [40]. Considering the shape model of non-spheric raindrops, we approximate large raindrops as oblate spheroids with an axis ratio represented by a 4th order polynomial [41,42] where a and b are the radii of the horizontal long axis, and the vertical short axis along the gravity direction, respectively. D is the equivalent spherical diameter in unit of centimeter. The axis ratio can reach 0.7 when the diameter D equals 5.0 mm. Then, the values of a and b can be calculated by using Eq. (3) and the volume equation Note that the backscattering cross section is calculated under the equilibrium condition where a balance between forces of surface tension and hydrostatic pressure is assumed. The complex oscillation phenomenon of raindrops and multiple scattering are not considered here [43,44].
As the backscattering efficiency (unitless) Q bk (D) is more often displayed, it is calculated and shown in Fig. 1. The relation between backscattering cross section σ bk (D) and backscattering efficiency Q bk (D) are expressed as The backscattering efficiency at an individual diameter will exhibit a highly variable behavior, due to the very large size parameters. While the actual scattering processes are associated with a large number of raindrops of nonuniform size. Therefore, an average operation is applied to the backscattering efficiency at each diameter by 100 drops with log-normal size distribution and a standard deviation of 1% of the central diameter, to remove the high-frequency oscillation [37,45]. The averaged results are shown in Fig. 1. The gray circles with the diameter range of 0.01-6 mm are calculated using Mie theory by considering the raindrops as water spheres. The backscattering efficiency of spheres decreases quickly from 0.1 mm due to the strong absorption of light within the drop, and reaches the geometric optical limit (P=0, external reflection) at the diameter of about 2 mm. For the actual oblate spheroids, the results are calculated using VCRM and shown as red squares from 1 to 6 mm. The backscattering efficiency of oblate spheroids increases with the diameter from about 2 mm. This is caused by the increased cross section area in the horizontal axis (perpendicular to the laser beam). We combine the data of Mie (0.1-1 mm) and the data of VCRM (1-6 mm), and fit the results by a smooth spline. The fitting line is used to retrieve the DSD in this study.

Retrieval of reflective spectrum
The shape of the power spectrum under precipitation conditions is the sum of the backscattering spectra from the aerosols and raindrops. However, the observed rain spectrum is broader than the true reflectivity spectrum, due to the effect of turbulence, wind shear, and window effects. Therefore, the observed spectrum under precipitation conditions can be expressed as [46,47] where v is the Doppler velocity (downward +); S air (v) is the turbulent spectrum from aerosols; S rain (v) is the reflectivity spectrum of rain given by Eq. (1); W(v) is the spectrum of the preset temporal window, which can be calculated theoretically. The asterisk represents the convolution operation. The turbulent spectrum of aerosol can be approximated as Gaussian distribution [12], given by where I air , v 0 , σ air are the intensity, mean vertical wind velocity, and spectrum width, respectively. Thus, Eq. (6) can be rewritten as where S ′ air (v) = S air (v) * W(v). In some previous works, the window effect is neglected and S ′ air (v) is treated as S air (v), using a Gaussian shape approximation [25,48]. This approximation introduces obvious error, as a short-time window function leads to side-lobes in spectral domain that can't be ignored. In other works, the window spectrum is firstly removed by a direct deconvolution before retrieving the raindrop scattering spectrum. Two times of deconvolution are needed [46]. Here, we built an accompanying aerosol signal model S ′ air (v), which is the convolution of the Gaussian-shaped turbulent spectrum and the window spectrum, to improve the accuracy and robustness when retrieving the rain reflectivity spectrum S rain (v).
Two main methods have been reported to derive the DSD from Doppler spectra: the parametric fitting method and the iterative deconvolution method [17,25]. In the parametric fitting method, the DSD is assumed to have a functional form, such as exponential distribution or Gamma distribution. The parameters are determined by nonlinear least square fitting. This method is sensitive to the initial values. In the iterative deconvolution method, no particular DSD shape is assumed. Here, we use the iterative method to retrieve the DSD. The following steps are performed.
1) Determine the aerosol signal S ′ air by the least square fitting method, as shown by the orange dash line in Fig. 2. 2) Calculate the initial estimate of the reflectivity spectrum (the star dash line in Fig. 2) by S rain,i=0 = S − S ′ air . Limit the available value of S rain,i=0 to the velocity range of 0∼10 m/s, to avoid non-realistic raindrops. What's more, the values that lower than the noise level are set as zero.
3) Calculate the spectrum S i (the green line in Fig. 2) using Eq. (8). Calculate the new estimates of the reflectivity spectrum S rain,i+1 = S rain,i × (S/S i ). A 10-points Gaussian moving average is pre-performed to S/S i , to reduce the amplified noise during the iterations.
4) The iteration continues from step 3 until the change of the reflectivity spectrum is negligible or the preset maximum iteration number is reached. The maximum iteration number is proportional to the data number used in the moving average process in Step 3. Too many iterations will amplify the noise. Aoki et al. set the number of iterations to 5. In our calculation, 6 iterations give better results. Figure 2 gives two examples of iterative deconvolution. A Rectangle window is used in the data processing preset in the CDL. The sidelobes of the aerosol peak extend to the velocity range Fig. 2. Two examples of iterative deconvolution. (a) The aerosol peak is higher than the rain peak and the vertical wind velocity is near zero, (b) the rain peak is higher than the aerosol peak, and an upward wind of −1.17 m/s exists. of the rain signal. By using an aerosol signal model in which the window effect is combined, a more accurate reflectivity spectrum is obtained. Although a longer pulse width will cause a low inherent spatial resolution, the resulting narrower spectrum width is predictable to improve the inversion accuracy. Figure 3 shows the measurement results of the CDL and MRR on June 13, 2020. The rain rate and visibility on the ground level are recorded by the Parsivel-2 and a co-located visibility sensor (Vaisala PWD50), respectively. Except for several minutes short pause, the precipitation continues from 00:00 to 18:00 with accumulated rainfall of 64 mm. The wideband carrier-tonoise ratio (CNR) is defined as the ratio of total signal power to noise power over the entire spectral bandwidth [35]. The CNR with values larger than 0 dB as shown in Fig. 3(a) mainly corresponds to clouds. The CDL is set to clean the telescope lens automatically every 5 minutes, to mitigate the effect of wavefront distortion on the receiving efficiency. Therefore, the CNR shows periodic stripes during precipitation. The spectrum width and skewness have been applied to precipitation identification in previous works [20,33]. The spectrum width will be broadened during precipitation conditions due to additional signal peaks from raindrops. Similarly, the skewness will depart from zero if the rain signal exists. Compared to the spectrum width, the skewness is more sensitive in identifying the precipitation but also requires relatively higher CNR. As shown in Fig. 3(c), the positive values of skewness mean that the aerosol peak is larger than the rain peak, and vice versa. Figure 3(d) shows the attenuated radar reflectivity recorded by the MRR, where high values corresponding to high rain rates.

Experiment and analysis
The DSDs obtained from CDL with 1-s resolution are then averaged both in the time and space domain to match the data resolution of MRR. To validate the performance, the DSDs are compared with that of MRR as shown in Fig. 4. The datasets are separated into three regimes of rain rate, namely, light rain (<1 mm/h), moderate rain (1-10 mm/h), and heavy rain (>10 mm/h), using the ground-level rain rate recorded by Parsivle-2. The lowest range gate of MRR is excluded from the comparison due to the near-field effect, as advised by the manufacture. Thus, the results at three heights from 200 m are plotted. In light rain, the difference of DSD between CDL and MRR becomes large from the diameter of 2 mm, as shown in Fig. 4(a1) and Fig. 4(b1). This may be caused by the low number concentration of large drops at light rain, and less sensing volume of CDL compared to the MRR. In moderate rain, the average DSD of CDL shows a great agreement with that of MRR at the diameter range of 1-5 mm. In heavy rain, the difference is relatively large due to small samples. At each condition, the number concentration of small drops with diameter less than 1 mm shows significant differences. The MRR gives the largest concentration at small drops, while the Parsivel-2 shows the lowest concentration. It has been reported that the MRR may overestimate the number of small drops due to spectrum aliasing and vertical wind [49,50]. The Parsivel-2 has also been pointed out to underestimate the number of small drops [30]. The CDL gives more reasonable values.
Actually, the diameter-dependent backscattering cross section function causes large errors at small drops. As analyzed in [10], for a certain instrument noise level s, the error of retrieved DSD can be expressed as s/σ bk (D). For the Rayleigh approximation in radar, σ bk ∝ D 6 holds. The error increases rapidly with a speed of inverse of 6th order of D as D approaches zero. While in the lidar, the error follows s/D 2 relation for the geometric approximation. Therefore, the error of small drops retrieved by CDL is much smaller than that obtained from MRR relative to large drops.
For further comparisons and statistics between the CDL and MRR, the reflectivity-weighted rain velocityv in the sense of light scattering, and the mass-weighted mean diameter D m are calculatedv The reflectivity spectra in the sense of light scattering are firstly calculated and compared between the CDL and MRR, as shown in Fig. 5. Although the influence of the overestimated concentration of small raindrops is insignificant on the integrated precipitation parameters such as rain rate (RR), liquid water content (LWC), and radar reflectivity (Z), it may greatly affect the calculated reflectivity spectrum in the sense of light scattering. As shown in Fig. 5(b) and Fig. 5(c), the reflectivity spectra of MRR show very high values at small velocity range due to the overestimation of small drop number. Therefore, a lower limit of the velocity (or diameter) is needed in the integration of Eq. (9) and Eq. (10) to get reasonable values. The position of the local minimum is used to determine the lower limit v min , below which the data of MRR are regarded as unreliable and not included in the comparison. Besides, in MRR retrieval, zero vertical wind is assumed due to its inability in wind detection. Any non-zero value will introduce a bias in the retrieved parameters. We use the wind velocity measured by the CDL to correct the MRR product. However, improper vertical wind correction may further deteriorate the quality of estimated parameters [10], due to the different sensing volumes between the CDL and MRR. Therefore, if the wind correction gives worse results, then the original results are used.
The scatterplot of mean rain velocityv and mean diameter D m at three heights between the CDL and MRR are shown in Fig. 6. The corresponding rain rates are displayed by the color of the scatters. The linear fitting results present good agreements at each height. The coefficient of determination (R 2 ) reaches at least 0.93. The slightly better consistency at higher height may The values of mean rain velocity are also given. v min is the lower limit in calculating the mean rain velocity. be due to the decrease of turbulence intensity with height increasing. The values of mean rain velocity and diameter show positive relation to the rain rate. The drizzle below 0.1 mm/h is dominated by small raindrops below 1 mm. While the positive relation becomes indistinct when the rain rate exceeds 1 mm/h, which means that it's the drop concentration but not the size that determines the rain rate. Although there exist large differences in the wavelength and scattering characteristics between the CDL and MRR, the derived DSD parameters show good agreement, except for the small drops. According to the error analysis, the CDL gives more reasonable values. However, there are still some limitations in using CDL to detect rainfall and other integrated parameters independently. To obtain the absolute value of DSD, the lidar equation should be solved exactly, which is difficult due to the following reasons. The first one is the strong attenuation of the transmitted laser by raindrops on the path, especially at high rain rates. The second one is the sharp reduction of coupling efficiency due to wavefront distortion caused by accumulated water on the telescope lens, which can be seen clearly from Fig. 3(a). Besides, the detection range of CDL is limited by thick clouds near the ground, as shown in Fig. 3(a) from 3:00 to about 4:30. It's hard for the laser to penetrate the thick cloud. When there exists a cloud signal, the rain signal could be submerged totally, because the backscattering efficiency of small particles with a diameter of 10-100 µm is much larger than that of the raindrops, as shown in Fig. 1. Consequently, the received signal is dominated by cloud particles.

Conclusion
In this paper, the DSD was retrieved based on the power spectrum of a 1.5 µm CDL using the theoretical backscattering cross section and an iterative deconvolution method. The accuracy and robustness were improved after considering the window effect of the accompanying aerosol signal model. The DSD retrieved from CDL was compared with that from a co-located MRR. The results showed good agreement. The MRR was found to overestimate the concentration of small drops, while the CDL gives more reasonable values, which is demonstrated using a ground-level Parsivel-2.
The radar reflectivity from MRR was used as reference information to determine the absolute value of DSD for CDL. The shape of the DSD and the derived parameters, such as the mean rain velocity and mean diameter of raindrops, were retrieved by CDL independently. For future work, we will try to solve the lidar equation during precipitation conditions. We also plan to investigate the evolution of vertical profiles of DSD and the influence of ambient wind using CDL and other instruments.

Disclosures. The authors declare no conflicts of interest.
Data Availability. Data underlying the results presented in this paper can be obtained from the authors upon reasonable request.