#### Abstract

Estimation of raindrop size distribution (DSD) is essential in many meteorological and hydrologic fields. This paper proposes a method for retrieving path-averaged DSD parameters using joint dual-frequency and dual-polarization microwave links of the telecommunication system. Detailed analyses of the rain-induced attenuation calculation are performed based on the T-matrix method. A forward model is established for describing the relation between the DSD and the rain-induced attenuation. Then, the method is proposed to retrieve propagation path DSD parameters based on Levenberg–Marquardt optimization algorithm. The numerical simulation for path-averaged DSD retrieval shows that the RMSEs of three gamma DSD parameters are 0.34 mm^{−1}, 0.81, and 3.21×10^{3} m^{−3}·mm^{−1}, respectively, in rainfall intensity above 30 mm/h. Meanwhile, the method can retrieve the rainfall intensity without the influence of variational DSD. Theoretical analyses and numerical simulations confirm that the method for retrieving path-averaged DSD parameters is promising. The method can complement existing DSD monitoring systems such as the disdrometer and provide high-resolution rainfall measurements with widely distributed microwave links without additional cost.

#### 1. Introduction

The raindrop size distribution (DSD) is the distribution of the number of raindrops according to their diameter (*D*). It was first introduced from the research weather radar [1], and further research has been carried on since then [2–4]. A capability to retrieve DSD parameters was shown to be crucial for studying the precipitation process, the weather radar measurements, the numerical weather prediction, and the weather modification [5–7].

Seliga and Bringi retrieved the concentration parameter and the drop median volume diameter of a two-parameter exponential DSD by using measurements of radar reflectivity (*Z*_{HH}) and differential reflectivity (*Z*_{DR}) [8, 9]. However, Ulbrich suggested that DSDs are more typically represented by gamma distribution, a three-parameter exponential DSD, measured by radar for short time periods commensurate [10]. Zhang et al. proposed a method [11] for retrieving three parameters of gamma DSD by using polarimetric radar. The results of the method show improvement over the existing techniques. Meanwhile, Zhang et al. addressed the measurement errors on the relation between the shape and slope parameters and reduced the bias and standard error [12]. Gorgucci et al. [13] propose a method for retrieving the DSD parameters from *Z*_{HH}, *Z*_{DR}, and specific differential phase (*K*_{DP}). However, *K*_{DP} is calculated from differential propagation phase with some noise. The method can produce large uncertainty, which causes some retrieval errors in space and time. Brandes et al. made a comparison of two polarimetric radar DSD retrieval algorithms [14]: the “beta” (*β*) method and the “constrained-gamma” method. The *β* method retrieves the DSD from *Z*_{HH}, *Z*_{DR}, and *K*_{DP}, while the “constrained-gamma” method uses *Z*_{HH}, *Z*_{DR}, and a constrained relation between the shape and slope parameters. The research shows that the “constrained-gamma” method provides reasonable retrieval results and the *β* method is sensitive to *K*_{DP} estimation errors.

Microwave links of the telecommunication system, designed for message passing between the signal towers, are known to be attenuated by atmospheric conditions in general. The specific rain-induced attenuation of a microwave link in the communication system at the frequencies of tens of GHz is attenuated seriously by the rain. It has been shown to be a practical method for estimating path-averaged rainfall intensity by using the measurement of the microwave link rain-induced attenuation [15–17]. This method can complement existing meteorological monitoring systems such as rain gauges network and weather radars, providing high-resolution rainfall measurements with widely distributed microwave links. Besides, this method can be used anytime with minimum supervision and without additional cost.

Rincon and Lang [18] and Berne and Schleiss [19] retrieved three parameters of DSD using the ratio of attenuations of dedicated two microwave links. In order to retrieve the parameters, a linear relationship between the scale parameter and the shape parameter is established based on DSD measurement data. However, the fitted relationship is not suitable for elsewhere. To solve it, a method is presented for retrieving DSD parameters using joint dual-frequency and dual-polarization (DFDP) microwave links without establishing the fitted relationship. The paper proposes a tentative idea for retrieving DSD in theory, and the method would be put into application with 5G communication in future.

This paper is organized in four sections. Section 2 introduces the rain-induced attenuation calculation based on nonspherical raindrops and presents the DFDP method for retrieving path-averaged DSD parameters by using the DFDP microwave links. Section 3 proposes the simulation experiments for testing the DFDP method. Section 4 summarizes the material presented in the paper and discusses the results.

#### 2. Materials and Methods

##### 2.1. Calculation of the Rain-Induced Attenuation Based on Nonspherical Raindrops

The rain-induced attenuation is calculated based on a theoretical relation between the electromagnetic radiation and the rainfall. Microwave link, which propagates through the rainy area between the signal transmitter and receiver, is attenuated by absorption and scattering of raindrops. The rain-induced attenuation in theory [20, 21] is defined aswhere (mm) is the raindrop equivalent volume diameter; (m^{2}) is the extinction cross section of the raindrop, which is related to the raindrop shape, frequency, and polarization of the microwave and temperature; and (mm^{−1}·m^{−3}) is the DSD distribution function.

###### 2.1.1. DSD

DSD function was commonly expressed to be an exponential distribution, named M-P distribution [1]:where (mm^{−1}·m^{−3}) is a number concentration parameter and (mm^{−1}) is a slope parameter. However, some researchers indicate that natural rain DSD contains fewer of both large and small raindrops [10, 22, 23]. Thus, Ulbrich suggested gamma distribution for representing M-P distribution aswhere is a shape parameter. M-P distribution is a special case of gamma distribution when . Gamma DSD is capable of describing a boarder variation in raindrops. Meanwhile, Chinese researchers [24] point that gamma DSD can also typically represent Chinese regional DSD. Therefore, the DFDP method is designed for retrieving gamma DSD.

###### 2.1.2. Raindrop Shape Models and the Extinction Cross Section

Previous studies have indicated that raindrop shape would change as the *D* increases, which is shown in Table 1.

The axial ratio of the raindrop is defined to represent raindrop shape:where (mm) is the semiminor axis and (mm) is the semimajor axis of the raindrop. The shape is spherical when ratio = 1. However, raindrop shape is not always spherical with *D* increasing according to Table 1. Therefore, an accurate raindrop shape model is necessary. The Pruppacher–Beard (PB) model [25], Pruppacher–Pitter (PP) model [26], and Goddard model [27] are widely adopted for describing raindrop shape. In this paper, the PB model is adopted to calculate the rain-induced attenuation for the reason that three models have similar extinction efficiency factors from 15 GHz to 30 GHz (0 mm∼7 mm). The equation of PB model is written by

The raindrop extinction cross section (*C*_{ext}) should be calculated firstly in order to obtain the rain-induced attenuation based on equation (1). The T-matrix method [28], also called expanding boundary conditions method (EBCM), is selected to calculate the *C*_{ext} (m^{2}) in the study because it is used to obtain an analytical solution of Maxwell equation for particles with arbitrary shapes, which is suitable for calculation of the raindrop *C*_{ext}.

According to the T-matrix method, *C*_{ext} of the raindrop can be calculated by following equation [29]:where is the wave number in the embedding medium; is the electric field intensity of the incident wave; and are the expansion coefficients of the incident field; and are the expansion coefficients of the scattering field; and represents the polarization direction of the incident wave.

##### 2.2. DFDP Method for Retrieving Path-Averaged DSD Parameters

To estimate the parameters of DSD, the DFDP microwave measurement system is established. The system is shown in Figure 1. Three mutually independent microwave links propagate in the same path. The frequency and the polarization of each link are shown in Table 2. Note that three links with different frequencies or polarizations, respectively, can also work in the proposed method. DFDP links are here to illustrate the method for retrieving the DSD.

###### 2.2.1. Forward Model

The rain-induced attenuation of each link is different because of different frequencies and polarizations. Therefore, the following equations can be obtained according to equation (1):where is the rain-induced attenuation of the *i*^{th} link. Then, equation (7) can be modified to the following form:

Substituting the equation of the gamma DSD, equation (8) can be rewritten as follows:

Thus, equation (9) is used as the forward model to describe the relation between rain-induced attenuation of three links and DSD of the propagation path.

###### 2.2.2. Retrieval Method for Path-Averaged DSD Parameters

, , and can be obtained by microwave links mentioned with equation (7). Meanwhile, *C*_{ext} of the raindrop can be calculated by equation (6). Thus, the forward model can be seen as a nonlinear equation with two unknowns of path-averaged gamma DSD.

In order to solve the equation, firstly set the objective function to be

The vector form of equation (10) is where

Therefore, the problem of solving the forward model is converted to be a minimization problem:

To solve equation (12), Levenberg–Marquardt (LM) optimization algorithm [30] is adopted in the study. The LM algorithm is based on the Newton method, the Gauss–Newton method, and the steepest descent method. The detailed steps of solving equation (12) are as follows:

*Step 1. *Set the initial value and the error accuracy . Then, the objective function is .

*Step 2. *Set and .

*Step 3. *Calculate and , where is the value of the Jacobian matrix at point and is the unit matrix.

*Step 4. *If , set and and then proceed to Step 3. Otherwise, if , stop the iteration; if , let and proceed to Step 3.

Figure 2 shows a flowchart about solving the numerical solutions of the DSD parameters by using the LM algorithm.

and , the numerical solutions of two DSD parameters, can be obtained by using the method presented above. Then, can be calculated by reversing equation (1), with and substituted into equation (1):

Thus, all three parameters of the path-averaged gamma DSD are obtained. The flowchart of the DFDP method for retrieving the path-averaged DSD parameters is shown in Figure 3.

#### 3. Results and Discussion

Three links propagating the same path are rare in reality. However, many dual-frequency or dual-polarization links that the two links propagate the same path are spread in communication. Moreover, we can select the nearest link as the third link to combine the three links, approximatively, for retrieving the DSD. In this way, we designed the simulation experiments, an ideal scenario of three links in the same path, to verify the method. Meanwhile, the simulation results can convince the communication operator to provide the data of the links. The link characteristics, similar to the real links, are given in Table 3. The environmental temperature is 25°C, and the atmospheric pressure is 1013 hPa. DSD parameters, measured by using OTT disdrometer in June and July 2013, in Nanjing, China, are adopted to be true values for testing the retrieval [31]. The DSDs measured by using the OTT disdrometer are divided into two groups according to the rainfall intensity: the rainfall intensity of the first group is below 30 mm/h, and second one is above 30 mm/h.

Four statistical indicators are used to evaluate the DFDP method [32]: the mean bias (MB), the normalized mean bias (NMB), the normalized mean error (NME), and the root mean square error (RMSE). NMB and NME reflect the relative bias between the retrieval values and the true values. The retrieval results are considered to be good if the NMB and NME values are both less than 50%. MB and RMSE reflect the bias between the retrieval and the true value. The indicators are calculated as follows:where is the retrieval value and is the true value.

##### 3.1. Path-Averaged DSD Parameters

Figures 4 and 5 show the retrieval DSD parameters vs*.* the true values of the first group. The retrieval values of Λ, *μ*, and *N*_{0} are close to the true values from Figures 4(a)–4(c), which indicates good correlation with the true values. Figure 5 shows the scatter plot between the true and the retrieval values, and the solid line is the fitted regression line resulting from the scatters. According to the figure, the scatters are distributed nearby the straight line with the slope nearly equal to one, which indicates a satisfactory retrieval result in theory. The results indicate that the DFDP method can retrieve the path-averaged DSD parameters. Table 4 shows a quantitative analysis of the first group retrieval results. The NMEs of three parameters are 11%, 28%, and 11% and the NMBs are 2%, 6%, and 1% as can be seen from the table. The NMBs and NMEs of them are less than 30%, which indicates the satisfied retrieval accuracy. Meanwhile, the values of MBs and RMSEs also indicate that the method have a good retrieval accuracy in the rainfall intensity below 30 mm/h.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

The second group retrieval DSD parameters are shown in Figures 6 and 7. Meanwhile, Table 5 shows a quantitative analysis of the second group. The results are similar to the first group as can be seen from the figure and the table. However, the better result of the second group, compared with the first group, may indicate that the method has a higher retrieval accuracy in heavy rainfall. Maybe, the DFDP method is more stable in heavy rain for the reason that the value of the rain-induced attenuation is larger.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

##### 3.2. Rainfall Intensity

The path-averaged rainfall intensity (*R*, mm/h) can be calculated by two methods based on DFDP links. The first calculation is based on the retrieval results of DSD [33]:where (m/s) is the raindrop terminal velocity. It can be calculated by the following equation [34]:

The second method is based on the well-known power law relation between *γ*_{rain} and *R*:where and are the function of frequency *f* (GHz) and the temperature. can be estimated by the inverse transformation of equation (17) on the condition that the *γ*_{rain} is measured:

At present, the specific attenuation model of rain-induced attenuation ratio proposed by International Telecommunication Union (ITU) is widely used. The values of *k* and *α* (25 GHz-h, 38 GHz-h and 38 GHz-v) are given in Table 6 [35].

Figures 8 and 9 show the retrieval results of *R* about two groups by two methods, compared with the true values. The retrieval *R* is close to the true value in two groups from Figures 8 and 9. The results of both groups by two methods have good correlation with the true rainfall intensity. However, *R*-DFDP, based on the retrieval of the DSD, is closer to *R*-*True* than the retrieval *R* by the second method (*R*-25 GHz-h,*R*-38 GHz-h and *R*-38 GHz-v), which may indicate that variational DSD influences the accuracy of rainfall intensity estimation and the first method is more accurate than the second one, consequently. Moreover, as can be seen from the two figures, *R*-38 GHz-h and *R*-38 GHz-v are higher than *R*-*True* while *R*-25 GHz-h is less than *R*-*True* sometimes, especially when rainfall intensity is greater than 15 mm/h (more generally in the second group), meaning that the *R*-38 GHz-h and *R*-38 GHz-v are influenced badly by DSD. Tables 7 and 8 show the quantitative analysis of the two groups. The performance of the first method is slightly better than that of the second method because all the evaluation parameters of the first one are better than the second in two groups. Further, as can be seen from Table 8, *R*-38 GHz-v deviate from *R*-*True* seriously, indicating that the rain-induced attenuation of 38 GHz-v is influenced by variational DSD, badly. The results suggest that the first method has a better estimation accuracy of the rainfall intensity.

#### 4. Conclusions

In this paper, the DFDP method is proposed for retrieving path-averaged DSD parameters by using microwave DFDP links of the telecommunication system. The method can complement existing DSD monitoring systems such as the weather radars and provide high-resolution rainfall measurements with widely distributed microwave links without additional cost.

Detailed analyses of the rain-induced attenuation calculation based on PB raindrop model are performed. A forward model is established for describing the relation between the DSD and the rain-induced attenuation of the DFDP microwave links. Then, the model is converted into a minimization problem, and it is solved to get numerical solutions by using LM optimization algorithm. Thus, path-averaged DSD parameters can be obtained. Meanwhile, the rainfall intensity can be calculated based on the retrieval DSD parameters.

Theoretical analyses and numerical simulations confirm that the DFDP method for retrieving path-averaged DSD parameters is promising. The results of the simulation experiment, designed for verifying the DFDP method, indicate that the method can retrieve *N*_{0}, Λ, and *μ*, though there are some outliers. Meanwhile, the method is more suitable for the heavy rain according to the retrieval results of the two groups with different rainfall intensity. It is because that the rain-induced attenuation in heavy rain is bigger. In addition, rainfall intensity can be calculated with the acceptable accuracy based on the retrieval DSD parameters according to Figures 8 and 9 and Tables 7 and 8. Moreover, the method has less influences of variational DSD on retrieving the rainfall intensity than the method by using the power law relation.

The DFDP method is a theoretical method for the DSD retrieval by using the DFDP links of the communication system. More experiments need to be conducted with the telecommunication company collaboration in order to put the method into application.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This research was funded by the National Natural Science Foundation of China (41505135 and 41475020) and the Jiangsu Natural Science Foundation, China (BK20150708).