# Raindrop Size Distribution Retrieval Using Joint Dual-Frequency and Dual-Polarization Microwave Links

Authors: Kun Song; Xichuan Liu; Taichang Gao; Binsheng He
Journal: Advances in Meteorology
 (2019)

---

## 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×103 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.


---

## Body

## 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 (ZHH) and differential reflectivity (ZDR) [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 ZHH, ZDR, and specific differential phase (KDP). However, KDP 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 ZHH, ZDR, and KDP, while the “constrained-gamma” method uses ZHH, ZDR, 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 KDP 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. Section2 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 as(1)γrain=4.343×∫0∞CextD,f,T⋅NDdD,where D (mm) is the raindrop equivalent volume diameter; CextD,f,T (m2) is the extinction cross section of the raindrop, which is related to the raindrop shape, frequency, and polarization of the microwave and temperature; and ND (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]:(2)ND=N0exp−ΛD,0<D<Dmax,where N0 (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 as(3)ND=N0Dμexp−ΛD,0<D<Dmax,where μ is a shape parameter. M-P distribution is a special case of gamma distribution when μ=0. 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 theD increases, which is shown in Table 1.Table 1
Raindrop shape vs. equivalent diameter.
D
<3 mm
3∼4 mm
4∼5.5 mm
>5.5 mm
Raindrop shape
Sphere
Flat ellipsoid
Ellipsoid of flat bottom
Ellipsoid of bottom with grooveThe axial ratio of the raindrop is defined to represent raindrop shape:(4)ratio=ba,where a (mm) is the semiminor axis and b (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(5)ba=1.03−0.062D.The raindrop extinction cross section (Cext) 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 Cext (m2) 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 Cext.According to the T-matrix method,Cext of the raindrop can be calculated by following equation [29]:(6)Cext,α=−1k12E0inc2Re∑n=1∞∑m=−nnamnpmn∗+bmnqmn∗,where k1 is the wave number in the embedding medium; E0inc is the electric field intensity of the incident wave; amn and bmn are the expansion coefficients of the incident field; pmn and qmn 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 Figure1. 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.Figure 1
Energy of each microwave link is attenuated by rainfall and other factors between the microwave transmitter and receiver, and the attenuation is recorded in the data acquisition system. Thus, rain-induced attenuation of three links can be obtained with the help of extraction algorithm of the rain-induced attenuation. h and v represent horizontal and vertical polarizations, respectively.Table 2
Frequencies and polarizations of the microwave links.
Link
Frequency (GHz)
Polarization
1
f
1
Horizontal polarization (h)
2
f
2
Horizontal polarization (h)
3
f
3
Vertical polarization (v)


#### 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):(7)γrain,1=4.343×∫0∞Cext,hD,f1,T⋅NDdD,γrain,2=4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,3=4.343×∫0∞Cext,vD,f1,T⋅NDdD,where γrain,i is the rain-induced attenuation of the ith link. Then, equation (7) can be modified to the following form:(8)γrain,1γrain,2=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,1γrain,3=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,vD,f1,T⋅NDdD.Substituting the equation of the gamma DSD, equation (8) can be rewritten as follows:(9)γrain,1γrain,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD,γrain,1γrain,3=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD.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

γ
rain
,
1, γrain,2, and γrain,3 can be obtained by microwave links mentioned with equation (7). Meanwhile, Cext 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(10)Fopt,1=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,2,Fopt,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,3.The vector form of equation (10) is Fopt=fx−κ where(11)fx=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD,x=μ,Λ,κ=γrain,1γrain,2γrain,1γrain,3.Therefore, the problem of solving the forward model is converted to be a minimization problem:(12)xopt=argminxFoptx.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 valuex=x0 and the error accuracy ε=10−5. Then, the objective function is Foptx0.Step 2.
SetμM=0.001 and k=0.Step 3.
CalculateFoptxk+1 and δxk=JkTJk+μMI−1JkTεk, where Jk is the value of the Jacobian matrix at point xk and I is the unit matrix.Step 4.
IfFoptxk+1>Foptxk, set xk+1=xk+δxk and μM=10μM and then proceed to Step 3. Otherwise, if δxk<ε, stop the iteration; if δxk≥ε, let xk+1=xk+δxk,μM=0.1μM and proceed to Step 3.Figure2 shows a flowchart about solving the numerical solutions of the DSD parameters by using the LM algorithm.Figure 2
Flowchart of inversion using the LM optimization algorithm.Λ
′ and μ′, the numerical solutions of two DSD parameters, can be obtained by using the method presented above. Then, N0′ can be calculated by reversing equation (1), with Λ′ and μ′ substituted into equation (1):(13)N0′=γrain,14.343×∫0∞Cext,hD,f1,T⋅Dμ′⋅exp−Λ′⋅DdD.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 Figure3.Figure 3
Flowchart of the retrieval using DFDP link.


## 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 as(1)γrain=4.343×∫0∞CextD,f,T⋅NDdD,where D (mm) is the raindrop equivalent volume diameter; CextD,f,T (m2) is the extinction cross section of the raindrop, which is related to the raindrop shape, frequency, and polarization of the microwave and temperature; and ND (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]:(2)ND=N0exp−ΛD,0<D<Dmax,where N0 (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 as(3)ND=N0Dμexp−ΛD,0<D<Dmax,where μ is a shape parameter. M-P distribution is a special case of gamma distribution when μ=0. 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 theD increases, which is shown in Table 1.Table 1
Raindrop shape vs. equivalent diameter.
D
<3 mm
3∼4 mm
4∼5.5 mm
>5.5 mm
Raindrop shape
Sphere
Flat ellipsoid
Ellipsoid of flat bottom
Ellipsoid of bottom with grooveThe axial ratio of the raindrop is defined to represent raindrop shape:(4)ratio=ba,where a (mm) is the semiminor axis and b (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(5)ba=1.03−0.062D.The raindrop extinction cross section (Cext) 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 Cext (m2) 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 Cext.According to the T-matrix method,Cext of the raindrop can be calculated by following equation [29]:(6)Cext,α=−1k12E0inc2Re∑n=1∞∑m=−nnamnpmn∗+bmnqmn∗,where k1 is the wave number in the embedding medium; E0inc is the electric field intensity of the incident wave; amn and bmn are the expansion coefficients of the incident field; pmn and qmn are the expansion coefficients of the scattering field; and α represents the polarization direction of the incident wave.


## 2.1.1. DSD

DSD function was commonly expressed to be an exponential distribution, named M-P distribution [1]:(2)ND=N0exp−ΛD,0<D<Dmax,where N0 (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 as(3)ND=N0Dμexp−ΛD,0<D<Dmax,where μ is a shape parameter. M-P distribution is a special case of gamma distribution when μ=0. 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 theD increases, which is shown in Table 1.Table 1
Raindrop shape vs. equivalent diameter.
D
<3 mm
3∼4 mm
4∼5.5 mm
>5.5 mm
Raindrop shape
Sphere
Flat ellipsoid
Ellipsoid of flat bottom
Ellipsoid of bottom with grooveThe axial ratio of the raindrop is defined to represent raindrop shape:(4)ratio=ba,where a (mm) is the semiminor axis and b (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(5)ba=1.03−0.062D.The raindrop extinction cross section (Cext) 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 Cext (m2) 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 Cext.According to the T-matrix method,Cext of the raindrop can be calculated by following equation [29]:(6)Cext,α=−1k12E0inc2Re∑n=1∞∑m=−nnamnpmn∗+bmnqmn∗,where k1 is the wave number in the embedding medium; E0inc is the electric field intensity of the incident wave; amn and bmn are the expansion coefficients of the incident field; pmn and qmn 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 Figure1. 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.Figure 1
Energy of each microwave link is attenuated by rainfall and other factors between the microwave transmitter and receiver, and the attenuation is recorded in the data acquisition system. Thus, rain-induced attenuation of three links can be obtained with the help of extraction algorithm of the rain-induced attenuation. h and v represent horizontal and vertical polarizations, respectively.Table 2
Frequencies and polarizations of the microwave links.
Link
Frequency (GHz)
Polarization
1
f
1
Horizontal polarization (h)
2
f
2
Horizontal polarization (h)
3
f
3
Vertical polarization (v)


### 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):(7)γrain,1=4.343×∫0∞Cext,hD,f1,T⋅NDdD,γrain,2=4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,3=4.343×∫0∞Cext,vD,f1,T⋅NDdD,where γrain,i is the rain-induced attenuation of the ith link. Then, equation (7) can be modified to the following form:(8)γrain,1γrain,2=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,1γrain,3=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,vD,f1,T⋅NDdD.Substituting the equation of the gamma DSD, equation (8) can be rewritten as follows:(9)γrain,1γrain,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD,γrain,1γrain,3=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD.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

γ
rain
,
1, γrain,2, and γrain,3 can be obtained by microwave links mentioned with equation (7). Meanwhile, Cext 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(10)Fopt,1=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,2,Fopt,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,3.The vector form of equation (10) is Fopt=fx−κ where(11)fx=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD,x=μ,Λ,κ=γrain,1γrain,2γrain,1γrain,3.Therefore, the problem of solving the forward model is converted to be a minimization problem:(12)xopt=argminxFoptx.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 valuex=x0 and the error accuracy ε=10−5. Then, the objective function is Foptx0.Step 2.
SetμM=0.001 and k=0.Step 3.
CalculateFoptxk+1 and δxk=JkTJk+μMI−1JkTεk, where Jk is the value of the Jacobian matrix at point xk and I is the unit matrix.Step 4.
IfFoptxk+1>Foptxk, set xk+1=xk+δxk and μM=10μM and then proceed to Step 3. Otherwise, if δxk<ε, stop the iteration; if δxk≥ε, let xk+1=xk+δxk,μM=0.1μM and proceed to Step 3.Figure2 shows a flowchart about solving the numerical solutions of the DSD parameters by using the LM algorithm.Figure 2
Flowchart of inversion using the LM optimization algorithm.Λ
′ and μ′, the numerical solutions of two DSD parameters, can be obtained by using the method presented above. Then, N0′ can be calculated by reversing equation (1), with Λ′ and μ′ substituted into equation (1):(13)N0′=γrain,14.343×∫0∞Cext,hD,f1,T⋅Dμ′⋅exp−Λ′⋅DdD.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 Figure3.Figure 3
Flowchart of the retrieval using DFDP link.


## 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):(7)γrain,1=4.343×∫0∞Cext,hD,f1,T⋅NDdD,γrain,2=4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,3=4.343×∫0∞Cext,vD,f1,T⋅NDdD,where γrain,i is the rain-induced attenuation of the ith link. Then, equation (7) can be modified to the following form:(8)γrain,1γrain,2=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,hD,f2,T⋅NDdD,γrain,1γrain,3=4.343×∫0∞Cext,hD,f1,T⋅NDdD4.343×∫0∞Cext,vD,f1,T⋅NDdD.Substituting the equation of the gamma DSD, equation (8) can be rewritten as follows:(9)γrain,1γrain,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD,γrain,1γrain,3=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD.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

γ
rain
,
1, γrain,2, and γrain,3 can be obtained by microwave links mentioned with equation (7). Meanwhile, Cext 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(10)Fopt,1=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,2,Fopt,2=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD−γrain,1γrain,3.The vector form of equation (10) is Fopt=fx−κ where(11)fx=∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f2,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,hD,f1,T⋅Dμ⋅exp−Λ⋅DdD∫0∞Cext,vD,f1,T⋅Dμ⋅exp−Λ⋅DdD,x=μ,Λ,κ=γrain,1γrain,2γrain,1γrain,3.Therefore, the problem of solving the forward model is converted to be a minimization problem:(12)xopt=argminxFoptx.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 valuex=x0 and the error accuracy ε=10−5. Then, the objective function is Foptx0.Step 2.
SetμM=0.001 and k=0.Step 3.
CalculateFoptxk+1 and δxk=JkTJk+μMI−1JkTεk, where Jk is the value of the Jacobian matrix at point xk and I is the unit matrix.Step 4.
IfFoptxk+1>Foptxk, set xk+1=xk+δxk and μM=10μM and then proceed to Step 3. Otherwise, if δxk<ε, stop the iteration; if δxk≥ε, let xk+1=xk+δxk,μM=0.1μM and proceed to Step 3.Figure2 shows a flowchart about solving the numerical solutions of the DSD parameters by using the LM algorithm.Figure 2
Flowchart of inversion using the LM optimization algorithm.Λ
′ and μ′, the numerical solutions of two DSD parameters, can be obtained by using the method presented above. Then, N0′ can be calculated by reversing equation (1), with Λ′ and μ′ substituted into equation (1):(13)N0′=γrain,14.343×∫0∞Cext,hD,f1,T⋅Dμ′⋅exp−Λ′⋅DdD.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 Figure3.Figure 3
Flowchart of the retrieval using DFDP link.


## 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 Table3. 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.Table 3
Characteristics of the links.
Link
f (GHz)
Polarization
Resolution (dB)
Length (km)
1
25
Horizontal (h)
0.01
1
2
38
Horizontal (h)
0.01
1
3
38
Vertical (v)
0.01
1Four 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:(14)MB=∑ildatar,i−datat,il,NMB=∑ildatar,i−datat,i∑ildatat,i,NME=∑ildatar,i−datat,i∑ildatat,i,RMSE=∑ildatar,i−datat,i2l,where datar is the retrieval value and datat is the true value.


### 3.1. Path-Averaged DSD Parameters

Figures4 and 5 show the retrieval DSD parameters vs. the true values of the first group. The retrieval values of Λ, μ, and N0 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.Figure 4
Retrieval results of the DSD parameters (first group); the blue line and the red line represent the true values and the retrieval values, respectively: (a) time sequence of Λ; (b) time sequence ofμ; (c) time sequence of N0.
(a)
(b)
(c)Figure 5
Retrieval results of the DSD parameters (first group); scatters are the true values vs. the retrieval values, and the solid lines represent the fitted regression lines from scatters.
(a)
(b)
(c)Table 4
Statistical parameters of the first group for the results.
MB
RMSE
NME (%)
NMB (%)
Λ (mm−1)
0.08
0.46
11
2
μ
0.19
1.12
28
6
N
0 (mm−1·m−3)
116.47
1.23 × 103
11
1The second group retrieval DSD parameters are shown in Figures6 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.Figure 6
Retrieval results of the DSD parameters (second group): the retrieval values with and without noise are distributed near the true values in all three DSD parameters.
(a)
(b)
(c)Figure 7
Retrieval results of the DSD parameters (second group): most of scatters are distributed in the line with the slope equal to one.
(a)
(b)
(c)Table 5
Statistical parameters of the second group for the results.
MB
RMSE
NME (%)
NMB (%)
Λ (mm−1)
−0.01
0.24
6
−1
μ
−0.01
0.62
14
−1
N
0 (mm−1·m−3)
−4.32
1.23 × 103
4
−1


### 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]:(15)R=π6∫0∞D3⋅ND⋅VDdD,where V(m/s) is the raindrop terminal velocity. It can be calculated by the following equation [34]:(16)V=0,D≤0.03mm,4.323D−0.03,0.03mm<D≤0.6mm,9.65−e−0.6D,D>0.6mm.The second method is based on the well-known power law relation betweenγrain and R:(17)γrain=k×Rα,where k and α are the function of frequency f (GHz) and the temperature. R can be estimated by the inverse transformation of equation (17) on the condition that the γrain is measured:(18)R=γraink1/α.At present, the specific attenuation model of rain-induced attenuation ratio proposed by International Telecommunication Union (ITU) is widely used. The values ofk and α (25 GHz-h, 38 GHz-h and 38 GHz-v) are given in Table 6 [35].Table 6
Values ofk and α (ITU model).
25 GHz-h
38 GHz-h
38 GHz-v
k
0.1571
0.4001
0.3844
α
0.9991
0.8816
0.8552Figures8 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.Figure 8
Retrieval results ofR (first group). R-True represents the true rainfall intensity. R-DFDP represents the rainfall intensity calculated by the first method based on the DSD while R-25 GHz-h, R-38 GHz-h, and R-38 GHz-v represent the rainfall intensity calculated by the second method of the ITU model in different frequencies and polarizations.Figure 9
Retrieval results ofR (second group).Table 7
Statistical parameters of the first group forR.
MB (mm/h)
RMSE (mm/h)
NME (%)
NMB (%)
R-DFDP
−0.01
0.14
1.0
−0.1
R-25 GHz-h
−1.15
1.29
8.4
−8.4
R-38 GHz-h
−0.52
1.13
6.5
−3.8
R-38 GHz-v
0.45
1.59
8.1
3.3Table 8
Statistical parameters of the second group forR.
MB (mm/h)
RMSE (mm/h)
NME (%)
NMB (%)
R-DFDP
0.01
0.27
0.4
0.1
R-25 GHz-h
−3.74
4.00
7.2
−7.2
R-38 GHz-h
4.80
8.49
10.3
9.3
R-38 GHz-v
11.60
16.01
22.7
22.5


## 3.1. Path-Averaged DSD Parameters

Figures4 and 5 show the retrieval DSD parameters vs. the true values of the first group. The retrieval values of Λ, μ, and N0 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.Figure 4
Retrieval results of the DSD parameters (first group); the blue line and the red line represent the true values and the retrieval values, respectively: (a) time sequence of Λ; (b) time sequence ofμ; (c) time sequence of N0.
(a)
(b)
(c)Figure 5
Retrieval results of the DSD parameters (first group); scatters are the true values vs. the retrieval values, and the solid lines represent the fitted regression lines from scatters.
(a)
(b)
(c)Table 4
Statistical parameters of the first group for the results.
MB
RMSE
NME (%)
NMB (%)
Λ (mm−1)
0.08
0.46
11
2
μ
0.19
1.12
28
6
N
0 (mm−1·m−3)
116.47
1.23 × 103
11
1The second group retrieval DSD parameters are shown in Figures6 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.Figure 6
Retrieval results of the DSD parameters (second group): the retrieval values with and without noise are distributed near the true values in all three DSD parameters.
(a)
(b)
(c)Figure 7
Retrieval results of the DSD parameters (second group): most of scatters are distributed in the line with the slope equal to one.
(a)
(b)
(c)Table 5
Statistical parameters of the second group for the results.
MB
RMSE
NME (%)
NMB (%)
Λ (mm−1)
−0.01
0.24
6
−1
μ
−0.01
0.62
14
−1
N
0 (mm−1·m−3)
−4.32
1.23 × 103
4
−1


## 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]:(15)R=π6∫0∞D3⋅ND⋅VDdD,where V(m/s) is the raindrop terminal velocity. It can be calculated by the following equation [34]:(16)V=0,D≤0.03mm,4.323D−0.03,0.03mm<D≤0.6mm,9.65−e−0.6D,D>0.6mm.The second method is based on the well-known power law relation betweenγrain and R:(17)γrain=k×Rα,where k and α are the function of frequency f (GHz) and the temperature. R can be estimated by the inverse transformation of equation (17) on the condition that the γrain is measured:(18)R=γraink1/α.At present, the specific attenuation model of rain-induced attenuation ratio proposed by International Telecommunication Union (ITU) is widely used. The values ofk and α (25 GHz-h, 38 GHz-h and 38 GHz-v) are given in Table 6 [35].Table 6
Values ofk and α (ITU model).
25 GHz-h
38 GHz-h
38 GHz-v
k
0.1571
0.4001
0.3844
α
0.9991
0.8816
0.8552Figures8 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.Figure 8
Retrieval results ofR (first group). R-True represents the true rainfall intensity. R-DFDP represents the rainfall intensity calculated by the first method based on the DSD while R-25 GHz-h, R-38 GHz-h, and R-38 GHz-v represent the rainfall intensity calculated by the second method of the ITU model in different frequencies and polarizations.Figure 9
Retrieval results ofR (second group).Table 7
Statistical parameters of the first group forR.
MB (mm/h)
RMSE (mm/h)
NME (%)
NMB (%)
R-DFDP
−0.01
0.14
1.0
−0.1
R-25 GHz-h
−1.15
1.29
8.4
−8.4
R-38 GHz-h
−0.52
1.13
6.5
−3.8
R-38 GHz-v
0.45
1.59
8.1
3.3Table 8
Statistical parameters of the second group forR.
MB (mm/h)
RMSE (mm/h)
NME (%)
NMB (%)
R-DFDP
0.01
0.27
0.4
0.1
R-25 GHz-h
−3.74
4.00
7.2
−7.2
R-38 GHz-h
4.80
8.49
10.3
9.3
R-38 GHz-v
11.60
16.01
22.7
22.5


## 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 retrieveN0, Λ, 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.



---

*Source: 7251870-2019-06-02.xml*