Feasibility Study of Rain Rate Monitoring from Polarimetric GNSS Propagation Parameters

In this work, the feasibility of estimating rain rate based on polarimetric Global Navigation Satellite Systems (GNSS) signals is explored in theory. After analyzing the cause of polarimetric signals, three physical-mathematical relation models between co-polar phase shift (KHH, KVV), specific differential phase shift (KDP), and rain rate (R) are respectively investigated. These relation models are simulated based on four different empirical equations of nonspherical raindrops and simulated Gamma raindrop size distribution. They are also respectively analyzed based on realistic Gamma raindrop size distribution and maximum diameter of raindrops under three different rain types: stratiform rain, cumuliform rain, and mixed clouds rain. The sensitivity of phase shift with respect to some main influencing factors, such as shape of raindrops, frequency, as well as elevation angle, is also discussed, respectively. The numerical results in this study show that the results by scattering algorithms T-matrix are consistent with those from Rayleigh Scattering Approximation. It reveals that they all have the possibility to estimate rain rate using the KHH-R, KVV-R or KDP-R relation. It can also be found that the three models are all affected by shape of raindrops and frequency, while the elevation angle has no effect on KHH-R. Finally, higher frequency L1 or B1 and lower elevation angle are recommended and microscopic characteristics of raindrops, such as shape and size distribution, are deemed to be important and required for further consideration in future experiments. Since phase shift is not affected by attenuation and not biased by ground clutter cancellers, this method has considerable potential in precipitation monitoring, which provides new opportunities for atmospheric research.


Introduction
The Global Navigation Satellite Systems (GNSS), mainly including the United States' Global Position System (GPS), Russia's GLONASS, China's BeiDou, and European Union's GALILEO, have been rapidly developed in the past few decades [1].Along with its vigorous development of GNSS, the non-navigational applications of GNSS signals have attracted more and more attention, especially in the field of Atmosphere-Ocean and Space environmental remote sensing.Currently, some research has been put into operational systems, with some of these operational systems having been investigated by numerous experiments.
The zenith tropospheric delay (ZTD) along the GNSS signal path, which is one of the major error sources for GNSS positioning, can be used to detect precipitable water vapor (PWV) [2,3], and its three-dimensional distribution can be retrieved by ground-based GNSS station networks.The results of these research endeavors created a new science, GNSS meteorology.The GNSS Radio Occultation (RO) signals received from mountain-based or space-based (Low-Earth-Orbit (LEO)) receivers are used to retrieve atmospheric temperature profiles, atmospheric moisture profiles, and total electron content (TEC) [4][5][6].Using reflected signals of GNSS from different underlying surface, we can monitor the parameters of ocean and land surface such as sea surface height, ocean wind speed and direction, snow thickness, ice thickness, and soil moisture [7][8][9].
A concept of heavy rain detection based on polarimetric GNSS-RO signals was proposed by Cardellach et al. in 2010 [10], which provided a new application direction.They aimed to use the double-polarization GNSS receivers placed on Spanish Low Earth Orbiter (LEO) satellite PAZ to receive RO signals across a precipitation area.Then, they did some theoretical research and ground-based field campaign using the polarimetric parameter and polarimetric phase shift, to validate its feasibility [11,12].Yan et al. [13] had analyzed the sensitivity of cross-polarization discrimination (XPD) of GNSS signals with rain rate by simulation.Up to now, the method of rain rate detection using GNSS signals is still under research.
The goal of this study is to explore the possibility of other polarimetric GNSS propagation parameters for rain monitoring in theory.The polarimetric radio-propagation parameters also contain specific co-polar attenuation, cross-polar attenuation, co-polar phase shift, and specific differential phase shift.Since GNSS signals belong to the L-band [14], rain-induced attenuation is quite weak [15].Therefore, we mainly investigate the co-polar phase shift and specific differential phase shift of GNSS signals.The occurrence of phase shift is due to the nonspherical shapes of raindrops and the canting of raindrops [16].Phase shift has been widely utilized in polarimetric weather radar in order to estimate rain rate and analyze the characteristics of different hydrometeors [17,18].The phase shift in weather radar reflects the backward scattering features of raindrops, while in this study, that of GNSS signals comes from the results of forward scattering.
This paper is organized as follows.Section 2 details the relation models between three polarimetric propagation parameters and rain rate, namely K HH -R, K VV -R, and K DP -R, based on different shapes of raindrops and Gamma raindrop size distribution, respectively.Using simulated raindrop size distribution, the three relation models are simulated and the feasibility of estimating rain rate using these models is analyzed in Section 3; the effects of scattering algorithm, shape of raindrops, frequency, and elevation angle on the three models are also discussed systematically; furthermore, the three relation models are investigated under three different real rain conditions.Section 4 summarizes some meaningful conclusions.

Shapes of Raindrops
Since it is the main cause of phase shift, considerations regarding the shapes of raindrops are essential for studying phase shift.Their final shapes are affected by hydrostatic, surface tension, and aerodynamic forces [19].Therefore, they are actually nonspherical.These shapes are often described by axis ratio r = a/b, where a and b represent the lengths of semiminor and semimajor axes, respectively.The relation between axis ratio and equivolumetric spherical diameter, also called shape-size relation, was researched by many experts.Pruppacher and Beard [20] developed an empirical linear equation (hereinafter named PB model) by tunnel measurements, which was used in early radar measurements.Based on different experimental data, some nonlinear equations, even fourth-order polynomials, were obtained by other researchers, such as Beard and Chuang (hereinafter named BC model) [21], Brandes et al. (hereinafter named Brandes model) [22], and Thurai et al.
(hereinafter named Thurai model) [23].The expressions of these shape models are listed in Table 1.
Here, D eq is the diameter of the equivolumetric spherical drop in mm.Among these shape models, since the BC model introduced aerodynamic pressure to the equilibrium condition, it is more precise and often used in many calculations or regarded as a comparative object [19,24,25].Therefore, it was selected for our following simulation research.

Raindrop Size Distribution
Phase shift of GNSS signals is also related to raindrop size distribution.There are many models to express raindrop size distribution, such as Laws-Parsons (LP) [26], Marshall-Palmer (MP), drizzle (JD), thunderstorm (JT) [27], and Gamma distribution [28].However, according to the analysis data, many naturally occurring raindrop size distributions can be represented by the three-parameter Gamma distribution [28] N where D is diameter of the equivolumetric spherical drop in millimeters; N(D)dD is the number of drops of size D to D + dD; N 0 is the intercept parameter in mm −1−µ •m −3 ; µ is the shape parameter in mm −1 ; and Λ is the slope parameter in mm −1 .The values of these parameters can be obtained from measured raindrop data or radar data.This Gamma distribution is widely accepted and used by radar meteorologists and other researchers to model natural raindrop size distributions [11,29,30].The rain rate R can be computed from [31] R mm h with the terminal fall velocity of the drop V ∞ (D) in meters per second.There are several empirical formulas to express V ∞ (D).A widely used one is given by Gunn and Kinzer [32] as

Calculation Models of Phase Shift Parameters
These polarimetric propagation parameters, co-polar phase shift (K HH , K VV ) and specific differential phase shift (K DP ), were considered.They are respectively given by [33] where λ is the wave length; Re indicates the real part; N(D) is the raindrop size distribution; and f HH and f VV are forward-scattering amplitudes at horizontal and vertical polarization state, respectively.Now, the scattering amplitudes for nonspherical drops can be computed from many algorithms [25].However, considering the frequencies of GNSS signals (L-band) and the shape-size features of raindrops, the Rayleigh Scattering Approximation and T-matrix are used in this simulation research [16,34,35].K HH , K VV , and K DP are all in deg•km −1 units.In Equations ( 4)-( 6), canting angles of raindrops are assumed to be 0 • .
Based on the above equations, the relation models between polarimetric propagation parameters K HH , K VV , K DP , and rain rate R, namely K HH -R, K VV -R, and K DP -R, can be respectively obtained.These relations will be simulated and discussed in the following section to study the feasibility of estimating rain rate.

Numerical Simulations Based on Simulated Raindrop Size Distribution
In this simulation, two kinds of scattering algorithms, Rayleigh Scattering Approximation and T-matrix, are selected to compare with each other.In order to use these algorithms, many factors, such as frequency and angle of incident electromagnetic waves, environmental temperature, refractive index of water, shapes of raindrops, as well as raindrop size distribution, are necessary.The L1 (1.57542 GHz) frequency of GPS is used in this paper to serve as an example.For the coordinate system in these algorithms, incident angle can be approximately equal to the elevation angle and is assumed to be 0 • here [25,34].The environmental temperature of 20 • C is used.The refractive index of water is calculated from Ray [36].The BC model is chosen to express the shapes of raindrops.
Since raindrop size distribution varies largely in different geographical locations and precipitation types, we consider the ranges of three parameters (N 0 , µ, Λ) of the Gamma size distribution to simulate natural raindrops size distributions.In considering previous research [19,29,30,37,38], the parameters of Gamma size distribution vary across a wide range and can be defined by 300 ≤ N 0 ≤ 300, 000 with the additional constraint that R ≤ 100 mm•h −1 .Figure 1 illustrates the simulated results of K HH -R, K VV -R, and K DP -R relations at the L1 (1575.42MHz) frequency of GPS signals based on simulated raindrop size distribution using Rayleigh Scattering Approximation (blue asterisk) and T-matrix (blue circles), respectively.It can be easily seen that the three polarimetric propagation parameters generally show upward trends with the increasing rain rate, and even their trends are similar.In these conditions, the value of K HH can be as large as 7 deg•km −1 , that of K VV can be more than 6.5 deg•km −1 , and the maximum value of K DP is less than 1 deg•km −1 .These results are also consistent with previous research by Ajaji [39].In contrast, the result of the T-matrix exhibited a very consistent correlation with that of the Rayleigh Scattering Approximation.This can be observed from the correlation coefficients for K HH , K VV , and K DP using the two kinds of scattering algorithms, which are all 1.0; and their root mean square errors are 0.0076, 0.0066, and 0.0011, respectively, all of which are rather small.Here, the number of statistical samples is 139,569.Therefore, the effect of scattering algorithms is quite weak, as such the T-matrix is used hereafter.
In Figure 1, the solid red line represents mean values of polarimetric propagation parameters for every constant rain rate.It is noted that for different scattering algorithms, the features of mean values of the three parameters are approximately the same, respectively.The parameters of the Gamma distribution used to derive the red line will be selected in our following simulation.
Figure 1 indicates that the variation of raindrop size distribution can lead to differences in the phase shift, so it is important to get the information of local raindrop size distribution in advance.Certainly, after getting the parameters of Gamma distribution, these polarimetric propagation parameters have the potential to estimate rain rate.

Sensitivity Analysis of Phase Shift with Respect to Shape of Raindrops
In order to analyze the effect of the shape of raindrops, four different shape models, PB model, BC model, Brandes model, and Thurai model, whose expressions were listed in Table 1, are used.Other parameters in this simulation have been depicted in Section 3.1.Figure 2 presents the polarimetric propagation parameters against rain rates for various shapes of raindrops at the L1

Sensitivity Analysis of Phase Shift with Respect to Shape of Raindrops
In order to analyze the effect of the shape of raindrops, four different shape models, PB model, BC model, Brandes model, and Thurai model, whose expressions were listed in Table 1, are used.
Other parameters in this simulation have been depicted in Section 3.1.Figure 2 presents the polarimetric propagation parameters against rain rates for various shapes of raindrops at the L1 (1575.42MHz) frequency of GPS signals.It can be observed that for K HH (Figure 2a) and K VV (Figure 2b) differences in the various shapes of raindrops between the two propagation parameters are not obvious, especially when rain rate is less than 50 mm•h −1 .However, when the rain rate becomes higher, K HH is largest in the PB model and smallest at Brandes model, while for K VV it shows the opposite character.Thus, the K DP is largest in the PB model and smallest in the Brandes model.This may be due to the fact that the range of Brandes model's equivolumetric spherical diameter only contains smaller raindrops (shown in Table 1), which induce a less specific differential phase shift.
From Figure 2, the effect of the shape of raindrops is weak when using K HH or K VV to estimate rain rate, while it should be considered when using K DP .2a) and KVV (Figure 2b) differences in the various shapes of raindrops between the two propagation parameters are not obvious, especially when rain rate is less than 50 mm•h −1 .However, when the rain rate becomes higher, KHH is largest in the PB model and smallest at Brandes model, while for KVV it shows the opposite character.Thus, the KDP is largest in the PB model and smallest in the Brandes model.This may be due to the fact that the range of Brandes model's equivolumetric spherical diameter only contains smaller raindrops (shown in Table 1), which induce a less specific differential phase shift.From Figure 2, the effect of the shape of raindrops is weak when using KHH or KVV to estimate rain rate, while it should be considered when using KDP.

Sensitivity Analysis of Phase Shift with Respect to Frequency
Frequency not only has an effect on the refractive index of water, but the free-space wavenumber in calculating forward-scattering amplitudes, finally it also influences the phase shift.To have a view of the effect of different frequencies, GPS L1/L2 (1575.42MHz/1227.60MHz) and Beidou B1/B2 (1561.098MHz/1207.14MHz) are utilized here.Polarimetric propagation parameters against rain rates for frequencies of GPS and BeiDou signals are shown in Figure 3.It can be seen that polarimetric propagation parameters KHH, KVV, and KDP are all increasing with the growing frequency.Obviously, the results of L1 frequency are quite close to those of B1 frequency.This is no surprise because their frequencies have little difference.Also, it is similar for the results of L2 and those of B2.

Sensitivity Analysis of Phase Shift with Respect to Frequency
Frequency not only has an effect on the refractive index of water, but the free-space wavenumber in calculating forward-scattering amplitudes, finally it also influences the phase shift.To have a view of the effect of different frequencies, GPS L1/L2 (1575.42MHz/1227.60MHz) and Beidou B1/B2 (1561.098MHz/1207.14MHz) are utilized here.Polarimetric propagation parameters against rain rates for frequencies of GPS and BeiDou signals are shown in Figure 3.It can be seen that polarimetric propagation parameters K HH , K VV , and K DP are all increasing with the growing frequency.Obviously, the results of L1 frequency are quite close to those of B1 frequency.This is no surprise because their frequencies have little difference.Also, it is similar for the results of L2 and those of B2.
Figure 3 reveals that the effect of frequency is obvious, and in order to get larger K HH , K VV , and K DP when estimating rain rate, higher frequencies L1 and B1 are suggested.
Atmosphere 2016, 7, 159 7 of 12 Figure 3 reveals that the effect of frequency is obvious, and in order to get larger KHH, KVV, and KDP when estimating rain rate, higher frequencies L1 and B1 are suggested.

Sensitivity Analysis of Phase Shift with Respect to Elevation Angle
In the T-matrix scattering algorithm, elevation angle is an important input factor.Thus, it is necessary to take the elevation angle of GNSS signals into consideration.The elevation angle is set to 0°, 10°, 20°, 30°, and 40°.Other input factors have been mentioned in Section 3.1.The variations of polarimetric propagation parameters against rain rates for various elevation angles at the L1 (1575.42MHz) frequency of GPS signals have been calculated and presented in Figure 4.The elevation angle has no effect on KHH, which can be seen in Figure 4a.However, as can be observed from Figure 4b,c, its effect can be noted.KVV can be seen to have increased with the increase of elevation angle while KDP has shown an opposite trend.The maximum difference for KDP using 0° and 40° can reach up to 0.2324 deg•km −1 , which is not negligible, and therefore cannot be ignored.
Figure 4 tells us that the elevation angle should be known when using KVV and KDP to get rain rate, while it is not necessary for KHH.

Sensitivity Analysis of Phase Shift with Respect to Elevation Angle
In the T-matrix scattering algorithm, elevation angle is an important input factor.Thus, it is necessary to take the elevation angle of GNSS signals into consideration.The elevation angle is set to 0 • , 10 • , 20 • , 30 • , and 40 • .Other input factors have been mentioned in Section 3.1.The variations of polarimetric propagation parameters against rain rates for various elevation angles at the L1 (1575.42MHz) frequency of GPS signals have been calculated and presented in Figure 4.The elevation angle has no effect on K HH , which can be seen in Figure 4a.However, as can be observed from Figure 4b,c, its effect can be noted.K VV can be seen to have increased with the increase of elevation angle while K DP has shown an opposite trend.The maximum difference for K DP using 0 • and 40 • can reach up to 0.2324 deg•km −1 , which is not negligible, and therefore cannot be ignored.
Figure 4 tells us that the elevation angle should be known when using K VV and K DP to get rain rate, while it is not necessary for K HH .

Simulation Calculation of Phase Shift under Real Rainfall
In Figure 1 and the above results, it can be seen that the effect of parameters on raindrop size distribution is obvious.In real rainfall, these parameters vary with position and time [28][29][30][37][38][39][40], which is an important factor in the model.In future experiments, we can use local statistical raindrop size distribution data or an empirical model.However, this may cause errors for rain rate estimation.Therefore, if possible, a disdrometer is needed in order to obtain detailed information regarding raindrop size distribution during a rain event.Thus, the errors can be reduced.Here, we take the results of raindrop size distribution obtained from a second-generation one-dimensional (1D) laser-optical disdrometer manufactured by OTT-Germany, which are used to calculate phase shift.The values of three parameters ( ) of the Gamma size distribution and the maximum diameters of raindrops during three different rain types, stratiform rain, cumuliform rain, and mixed clouds rain, can be seen in [41].The results of phase shift under these three rain types are shown in Figure 5.

Simulation Calculation of Phase Shift under Real Rainfall
In Figure 1 and the above results, it can be seen that the effect of parameters on raindrop size distribution is obvious.In real rainfall, these parameters vary with position and time [28][29][30][37][38][39][40], which is an important factor in the model.In future experiments, we can use local statistical raindrop size distribution data or an empirical model.However, this may cause errors for rain rate estimation.Therefore, if possible, a disdrometer is needed in order to obtain detailed information regarding raindrop size distribution during a rain event.Thus, the errors can be reduced.Here, we take the results of raindrop size distribution obtained from a second-generation one-dimensional (1D) laser-optical disdrometer manufactured by OTT-Germany, which are used to calculate phase shift.The values of three parameters (N 0 , µ, Λ) of the Gamma size distribution and the maximum diameters of raindrops during three different rain types, stratiform rain, cumuliform rain, and mixed clouds rain, can be seen in [41].The results of phase shift under these three rain types are shown in Figure 5.In Figure 5, since the parameters of raindrop size distribution are variable during the rain process, the phase shift presents different values.The rain rates computed from most of the samples are less than 10 mm•h −1 , and the corresponding phase shift is smaller.While the phase shift becomes larger at high rain rates.Also, it can be easily seen that no matter the type of the (real) rain event, the polarimetric propagation parameters KHH, KVV, and KDP all have a notable positive correlation with rain rate.Compared with Figure 1, the values of KHH, KVV, and KDP are all included in the ranges of the results in Figure 1, respectively.Obviously, their values are not the same under different rain types.Especially for KVV and KDP, it seems that KVV is the smallest and KDP is the largest under cumuliform rain.This is caused by the fact that more large-diameter raindrops occur in cumuliform rain and the maximum diameter of raindrops is 7.5-9.5 mm, while in stratiform rain and mixed clouds rain they are 3.25-4.25mm and 4.25-4.75mm, respectively [41].
The results calculated from these cases of real rain conditions demonstrate that the polarimetric propagation parameters KHH, KVV, and KDP are still sensitive to rain rate and have the possibility to estimate rain rate using realistic raindrop size distribution.In order to get more accurate phase shift processing, the detailed information of raindrop size distribution is required in experiments.

Conclusions
This paper has researched the theoretical feasibility of estimating rain rate based on the polarimetric propagation parameters of GNSS signals.Three theoretical models were set up that between co-polar phase shift KHH, KVV, specific differential phase shift KDP, and rain rate R, respectively.The models, KHH-R, KVV-R, and KDP-R, were simulated based upon assuming specific ranges of three parameters of Gamma raindrop size distribution and empirical nonlinear shape-size In Figure 5, since the parameters of raindrop size distribution are variable during the rain process, the phase shift presents different values.The rain rates computed from most of the samples are less than 10 mm•h −1 , and the corresponding phase shift is smaller.While the phase shift becomes larger at high rain rates.Also, it can be easily seen that no matter the type of the (real) rain event, the polarimetric propagation parameters K HH , K VV , and K DP all have a notable positive correlation with rain rate.Compared with Figure 1, the values of K HH , K VV , and K DP are all included in the ranges of the results in Figure 1, respectively.Obviously, their values are not the same under different rain types.Especially for K VV and K DP , it seems that K VV is the smallest and K DP is the largest under cumuliform rain.This is caused by the fact that more large-diameter raindrops occur in cumuliform rain and the maximum diameter of raindrops is 7.5-9.5 mm, while in stratiform rain and mixed clouds rain they are 3.25-4.25mm and 4.25-4.75mm, respectively [41].
The results calculated from these cases of real rain conditions demonstrate that the polarimetric propagation parameters K HH , K VV , and K DP are still sensitive to rain rate and have the possibility to estimate rain rate using realistic raindrop size distribution.In order to get more accurate phase shift processing, the detailed information of raindrop size distribution is required in experiments.

Conclusions
This paper has researched the theoretical feasibility of estimating rain rate based on the polarimetric propagation parameters of GNSS signals.Three theoretical models were set up that between co-polar phase shift K HH , K VV , specific differential phase shift K DP , and rain rate R, respectively.The models, K HH -R, K VV -R, and K DP -R, were simulated based upon assuming specific ranges of three parameters of Gamma raindrop size distribution and empirical nonlinear shape-size relation of raindrops.The results indicate that they all have the possibility to monitor rain rate.We also compared the results by T-matrix with those by Rayleigh Scattering Approximation, finding that they are very close to each other.
Additionally, some influencing factors, such as the shape of raindrops, frequency, and elevation angle have been considered and their sensitivity has been studied.It was found that the shape of raindrops and frequency of GNSS signals both have effect on the K HH -R, K VV -R, and K DP -R relationships, while the elevation angle only affects the K VV -R and K DP -R relationships.In order to get a larger phase shift, the use of a lower elevation angle and higher frequency L1 or B1 are suggested for future experiments.Moreover, the phase shift has been discussed based on realistic Gamma distribution under three different rain types, and the results also show the feasibility of this method.To reduce errors, the microscopic characteristics of rain, such as the shape of raindrops and raindrop size distribution should be obtained during a precipitation event.
It is quite new to use phase shift to estimate rain rate, which extends the application direction of GNSS signals and provides a new method for the area, real time, and passive measurement of rain rate.Certainly, there are some issues that need to be studied further.However, there is no need to worry about the detectability of phase shift because there are multi-frequency receivers with GPS/Beidou compatibility, which can measure the carrier phase of GPS and Beidou satellites precisely to millimeters [42][43][44].However, since an antenna is required to receive horizontal and vertical components of GNSS signals respectively, a dual-polarized antenna should be designed.Then, field campaigns are necessary to examine its validity.

1 )Figure 1 .
Figure 1.Results of polarimetric propagation parameters vs. rain rate at the L1 (1575.42MHz) frequency of GPS signals based on simulated raindrop size distribution.(a) K HH -R using Rayleigh Scattering Approximation; (b) K HH -R using T-matrix; (c) K VV -R using Rayleigh Scattering Approximation; (d) K VV -R using T-matrix; (e) K DP -R using Rayleigh Scattering Approximation; (f) K DP -R using T-matrix.

Figure 2 .
Figure 2. Polarimetric propagation parameters vs. rain rate for various shapes of raindrops at the L1 (1575.42MHz) frequency of GPS signals.(a) K HH -R; (b) K VV -R; (c) K DP -R.

Figure 3 .
Figure 3. Polarimetric propagation parameters vs. rain rate for frequencies of GPS and BeiDou signals.(a) K HH -R; (b) K VV -R; (c) K DP -R.

igure 4 .
Polarimetric propagation parameters vs. rain rate for various elevation angles at the L1 (1575.42MHz) frequency of GPS signals.(a) K HH -R; (b) K VV -R; (c) K DP -R.

Table 1 .
Expressions for different raindrop shape-size relations.
1575.42 MHz) frequency of GPS signals.It can be observed that for KHH (Figure (