Retrieval of Raindrop Size Distribution Using Dual-Polarized Microwave Signals from LEO Satellites: A Feasibility Study through Simulations

In this paper, a novel approach for raindrop size distribution retrieval using dual-polarized microwave signals from low Earth orbit satellites is proposed. The feasibility of this approach is studied through modelling and simulating the retrieval system which includes multiple ground receivers equipped with signal-to-noise ratio estimators and a low Earth orbit satellite communicating with the receivers using both vertically and horizontally polarized signals. Our analysis suggests that the dual-polarized links offer the opportunity to estimate two independent raindrop size distribution parameters. To achieve that, the vertical and horizontal polarization attenuations need to be measured at low elevation angles where the difference between them is more distinct. Two synthetic rain fields are generated to test the performance of the retrieval. Simulation results suggest that the specific attenuations for both link types can be retrieved through a least-squares algorithm. They also confirm that the specific attenuation ratio of vertically to horizontally polarized signals can be used to retrieve the slope and intercept parameters of raindrop size distribution.


Introduction
Using the attenuation of electromagnetic signals to measure rainfall has received much attention recently. This approach is strongly supported by the theoretical basis for calculating attenuation due to the water content along a signal path, which has been established and validated by numerous propagation tests [1,2]. In 1991, it was suggested that custom-designed microwave links could be used for the reconstruction of rainfall fields [3]. The usage of the backhaul links of commercial wireless communication networks (CWCNs) for rainfall monitoring was first proposed by Messer et al. in 2006 [4] and then by Leijnse et al. in 2007 [5]. Since then, this technique has been improved rapidly through extensive studies (e.g., [6][7][8][9]). The Earth-space links of commercial geostationary Earth orbit (GEO) satellites have also been studied to serve the same purpose. It has been demonstrated that the GEO system is capable of path-integrated rain rate estimation [10][11][12], and that the estimation has a high accuracy for heavy rain events [13].
The recent development in large constellations of low Earth orbit (LEO) satellites for global broadband services hints at a new opportunity for rainfall measurements. Several companies are pursuing large constellations of LEO satellites to provide global broadband access, e.g., Starlink's system with more than 10,000 satellites and Telesat's system with 300 satellites [14]. Once the satellite broadband service is available to the public, it can be foreseen that millions of user terminals (ground receivers) will be deployed across the globe. A dynamic network of satellite-to-ground links will be formed on an unprecedented scale. It will become a very large sensor network that could provide off-the-shelf signal attenuation data for the observations of the atmosphere. Meanwhile, because of the increasing demand for multimedia applications and thus the requirement of high data rates, satellites broadband services need to use Ku (12~18 GHz), Ka (26.5~40 GHz), and/or Q/V (40~50 GHz) frequency bands [15]. Signals at these frequencies are severely attenuated by gas, cloud and rain in the atmosphere [16]. They are especially sensitive to liquid water, which indicates that measuring precipitation through them can be a valid approach. The opportunistic use of satellite communication links can provide global real-time observations at a very low cost.
One of the key advantages of using a LEO satellite system for rainfall estimation is that it can offer three-dimensional (3-D) rain field tomography. First proposed by Huang et al. in 2016 [17], this approach exploits the similarity in the revolutionary motion between the LEO satellites and the signal source in computed tomography (CT) [18] and performs tomographic reconstruction of rain attenuation fields. Various investigations have been carried out recently to further examine this approach. For example, it has been suggested that 3-D attenuation fields can be retrieved by using the estimated signal-tonoise ratio (SNR) at the ground receivers and an iterative method to estimate the sky noise [19]. Furthermore, some signal processing strategies have been explored to optimize the performance of the rain field retrieval [20][21][22].
The attenuation of microwave links also provides the opportunity for raindrop size distribution (DSD) measurements. Knowledge of the DSD is fundamental to quantitative precipitation estimation and realistic DSD representation is very important to atmospheric research. Studies have suggested that the DSD is best modeled by a gamma distribution [23], which is in the following form: in which D is the equivolumetric raindrop diameter in mm, N 0 is the intercept parameter in mm −1 ·m −3 , µ is a dimensionless shape parameter and Λ is the slope parameter in mm −1 . One straightforward approach to estimate the three parameters is through three different attenuation measurements. For instance, previous studies have shown that attenuation measurements of multi-frequency links can be used to retrieve the parameters [24,25]. Taking advantage of the oblateness of large raindrops, the approach of using attenuation measured by multi-polarization links has also been investigated [26,27]. Furthermore, the multi-polarization and multi-frequency approaches can be combined together to improve the performance of the DSD retrieval [28]. On the other hand, it has been suggested that not all three parameters are always needed for the representation of DSD. Some studies use fixed µ models (e.g., [25,29]), while others identify a µ − Λ relationship (e.g., [30,31]) so that the number of independent parameters is also reduced to two. Both approaches have their merits and other factors such as rain types and locations determine which approach performs better.
To date, existing studies use ground links at fixed locations to retrieve the pathaveraged DSD parameters. In this paper, we will explore the possibility of exploiting the tomographic capability of moving signal links offered by LEO satellites. Realistic system models are proposed for the satellite communication system in which both vertically and horizontally polarized links are in operation. Synthetic DSD fields are adopted to suit the purpose of the study, and we use the Pruppacher-Beard raindrop shape model [32] and the T-matrix method [33] to calculate the attenuation caused by electromagnetic scattering. Our calculation shows that the elevation angle of the link does not affect the specific attenuation very much for horizontally polarized signals. While for vertical polarization, the increase in specific attenuation for the same rain field could be more than 10% when the elevation angle changes from 40 to 90 degrees. The two attenuation fields retrieved by the two differently polarized links are examined for the purpose of DSD estimation. It is found that a linear relationship can be identified between the slope parameter and the retrieved specific attenuation ratio of vertically to horizontally polarized signals. Using this relationship, we examine a fixed µ model and investigate the system's ability to retrieve the slope (Λ) and intercept (N 0 ) parameters. Simulation results show that for areas of high Sensors 2021, 21, 6389 3 of 14 specific attenuations, the retrieval is more accurate, thereby both the retrieved slope and intercept parameters have a close agreement with their true values.
The main contributions of this study are summarized as follows.
• A theoretical model is built for the novel approach for raindrop size distribution retrieval using dual-polarized microwave signals from LEO satellites; • The feasibility of the approach is investigated through simulations of synthetic rain events and realistic satellite communication systems; • It is confirmed that the specific attenuation ratio of vertically to horizontally polarized signals can be used to retrieve the slope and intercept parameters of DSD.
The rest of the paper is organized as follows. Section 2 describes the satellite signal model, and explains how the SNR at the receivers is related to the path-integrated attenuation. It also establishes the DSD parametric model and examines the theoretical basis of DSD retrieval. In Section 3, simulation results of the synthetic DSD fields, the retrieved attenuation fields and DSD parameters are presented and discussed. Section 4 concludes the paper.

Satellite Signal Model
We use the same satellite signal model presented in previous studies on rain field [19] and cloud field retrieval [34]. The model consists of a LEO satellite passing directly over the area of interest and multiple ground receivers equipped with electronically steerable antennas. It should be noted that the receivers designed for satellite broadband will have a minimum elevation angle requirement in order to satisfy the minimum signal quality requirement for communication. When the minimum elevation angle is met, the satelliteto-ground link becomes operational and the SNRs at the ground receivers are estimated for the purpose of measuring the path-integrated attenuation. For one overpass, suppose that k denotes the elevation angle in degrees, and ρ(k) is the SNR at the receivers in the form of [19,20] ρ where A I (k) and A F (k) are the path-integrated rain attenuation and the free space path loss, in decibels, respectively. Parameter C(k) is an unknown baseline value for a receiver, which is mainly determined by the power of the transmitted signal, the antenna gain and the noise temperature of the receiving system. Noise figure F n (k) is a function of the sky temperature and the noise temperature of the receiving system. We assume that the noise temperature of the receiving system is known so F n (k) is a function of the unknown sky noise. A detailed description of C(k), A F (k) and F n (k) can be found in [19] (pp. 5437-5438). It is assumed that each receiver is operating on a dual-polarization mode, i.e., both vertical and horizontal polarizations are used at the same time. Hence, the system offers two estimated SNRs, namelyρ V (k) andρ H (k) at every receiver, for measuring path-integrated attenuations.

Specific Attenuation and DSD
The relationship between the specific attenuation (a in dB/km) and the DSD of (1) is given by [35] where Q(D) is the extinction cross section of the raindrop (in m 2 ), which is determined by the raindrop shape, temperature, and frequency, direction and polarization of the microwave signal. Because large raindrops are in an oblate shape, their extinction cross sections will be different for vertically and horizontally polarized microwave signals with the same wavelength. Figure 1 shows a sketch of the raindrop and signals with different polarizations, where Q H (D) and Q V (D) are the extinction cross sections for horizontally and vertically polarized links, respectively. As a result, the specific attenuations in the rain field, namely a V and a H , will also be different.
where is the semininor axis and is the semimajor axis of the raindrop. This model yields a very good approximation of raindrop shapes according to their diameters up to 7 mm [36], which enables us to utilize the T-matrix method [33,37] for calculating their extinction cross sections. The T-matrix method proposed by Waterman [38] is a powerful tool for accurately computing electromagnetic scattering by nonspherical particles. It is based on directly solving Maxwell's equations and yet its efficiency enables very fast computation. The fact that raindrops are approximately spheroids further simplifies the Tmatrix calculation.

Retrieval Method
The retrieval method aims to estimate the DSD parameters by using the attenuation measurements of both vertically and horizontally polarized signal links. For the LEO satellite communication systems, the elevation angle is constantly changing during an overpass. Assuming the raindrops have zero canting angles, we examine the extinction cross section for various drop sizes for 30 GHz signals at different elevation angles. Table 1 lists the values of ( ) when the elevation angle is at 40, 55, 70 and 90 degrees for raindrop sizes from 0.5 to 6.0 mm in diameter. Table 2 lists the ( ) to ( ) ratio to show the difference between the two. As the elevation angle goes higher and as the drop size becomes smaller, the gap between ( ) and ( ) diminishes. This indicates that for The exact shape of a raindrop is usually complicated as it is determined by the surface tension of the water and the air flow around it as it falls. However, there are simplified models available to suit the requirement of calculating Q(D). In this paper, we adopt the widely used Pruppacher-Beard model [32] to characterize the raindrop shape. The axial ratio of the raindrop is given by where x is the semininor axis and y is the semimajor axis of the raindrop. This model yields a very good approximation of raindrop shapes according to their diameters up to 7 mm [36], which enables us to utilize the T-matrix method [33,37] for calculating their extinction cross sections. The T-matrix method proposed by Waterman [38] is a powerful tool for accurately computing electromagnetic scattering by nonspherical particles. It is based on directly solving Maxwell's equations and yet its efficiency enables very fast computation. The fact that raindrops are approximately spheroids further simplifies the T-matrix calculation.

Retrieval Method
The retrieval method aims to estimate the DSD parameters by using the attenuation measurements of both vertically and horizontally polarized signal links. For the LEO satellite communication systems, the elevation angle is constantly changing during an overpass. Assuming the raindrops have zero canting angles, we examine the extinction cross section for various drop sizes for 30 GHz signals at different elevation angles. Table 1 lists the values of Q H (D) when the elevation angle is at 40, 55, 70 and 90 degrees for raindrop sizes from 0.5 to 6.0 mm in diameter. Table 2 lists the Q V (D) to Q H (D) ratio to show the difference between the two. As the elevation angle goes higher and as the drop size becomes smaller, the gap between Q V (D) and Q H (D) diminishes. This indicates that for better measurement of the DSD, we need to examine the horizontal and vertical polarization attenuations at lower elevation angles where their difference is more distinct.
With SNRs at multiple ground receivers estimated for both horizontally and vertically polarized links elevated at 40 • (ρ H (40) andρ V (40)), the specific attenuations a H (40) and a V (40) for a targeted area can be retrieved through a tomographic algorithm. This retrieval process involves dividing the area of interest into voxels, solving the specific attenuations of each voxel using algorithms such as least-squares [39] and iteratively updating the noise figures of the sky noises (see [19] (p. 5439) for details). With the two specific attenuations ready, from Equations (1) and (3), we obtain in which parameter N 0 is eliminated. By assuming that the shape parameter µ is a fixed known value, Equation (5) provides a one-to-one mapping between a V /a H and Λ. When Λ is successfully estimated, N 0 can then be solved using Equations (1) and (3).

Simulations and Results
Our proposed retrieval method was examined through a series of simulations. We use synthetic DSD fields combined with a realistic LEO satellite communication system to simulate the estimated SNRs. The retrieval results of the specific attenuations and the DSD parameters are compared with their true values.

Synthetic DSD Field
A vertical DSD field is generated synthetically for our simulations. The field is considered to remain static during the overpass time of the LEO satellite (approximately 289 s in our simulations). It is 6 km in length and 4 km in height, and divided into 12 × 8 = 96 voxels, with each voxels measuring 500 m by 500 m. The DSD profile for each voxel is then generated as follows. Firstly, the specific attenuation for each voxel for horizontally polarized 30 GHz signals with 40 • elevation angle (a H (40)) is generated according to a localized frontal rainfall event [19] simulated by the Weather Research and Forecasting (WRF) model. This event is a subset of a larger convection-permitting simulation covering the southern part of the Great Barrier Reef and adjacent coastline. As shown in Figure 2a, the entire rain field is contained within the 96 voxels. The specific attenuation varies from voxel to voxel. For instance, voxel A has the highest specific attenuation (a H (40) = 3.23 dB/km) and some voxels on the edge of the field have zero dB/km specific attenuations. Secondly, the shape parameter µ for all voxels is assumed to be 1, which is in line with the WDM6 Scheme of WRF. Thirdly, the field is divided into three layers (middle layer 2 km in height, top and bottom layers 1 km in height) and the slope parameters (Λ 1 , Λ 2 and Λ 3 ) are assigned for each layer (see Figure 2a). This is based on the assumption that for a particular rain event, the slope parameter is more spatially homogeneous than the intercept parameter [40]. According to the probability density function of Λ in recent studies [41], Λ 1 and Λ 3 are set to be 2.9 mm −1 and 4.1 mm −1 for the bottom and top layer, respectively. Initially assigned to 3.5 mm −1 , Λ 2 for the middle layer is able to take different values in the following simulations so that we can test the performance of the DSD retrieval. Finally, the intercept parameter N 0 for each voxel is calculated from a H (40) and Λ using Equations (1) and (3). The canting angles of raindrops follow a Gaussian distribution with a mean of 0 • and a standard deviation of 2 • [42,43]. In the calculation a MATLAB implementation of the T-Matrix method [44] is applied and the integral is computed numerically from D min = 0.01 mm to D max = 7.0 mm in 0.01 mm increments. Simulation results suggest that N 0 for voxels with nonzero specific attenuations is between 6.75 × 10 2 and 5.46 × 10 4 mm −1 ·m −3 when Λ 2 is 3.5 mm −1 . With N 0 and Λ ready, the specific attenuation for any elevation angle and for either polarizations can then be calculated. Figure 2b shows the true field of a V (40) for comparison, in which the specific attenuation for voxel A is 2.99 dB/km, 7.4% less than a H (40).
the southern part of the Great Barrier Reef and adjacent coastline. As shown in Figure 2a, the entire rain field is contained within the 96 voxels. The specific attenuation varies from voxel to voxel. For instance, voxel A has the highest specific attenuation ( (40) = 3.23 dB/km) and some voxels on the edge of the field have zero dB/km specific attenuations. Secondly, the shape parameter for all voxels is assumed to be 1, which is in line with the WDM6 Scheme of WRF. Thirdly, the field is divided into three layers (middle layer 2 km in height, top and bottom layers 1 km in height) and the slope parameters ( , and ) are assigned for each layer (see Figure 2a). This is based on the assumption that for a particular rain event, the slope parameter is more spatially homogeneous than the intercept parameter [40]. According to the probability density function of in recent studies [41], and are set to be 2.9 mm −1 and 4.1 mm −1 for the bottom and top layer, respectively. Initially assigned to 3.5 mm −1 , for the middle layer is able to take different values in the following simulations so that we can test the performance of the DSD retrieval. Finally, the intercept parameter for each voxel is calculated from (40) and using Equations (1) and (3). The canting angles of raindrops follow a Gaussian distribution with a mean of 0° and a standard deviation of 2° [42,43]. In the calculation a MATLAB implementation of the T-Matrix method [44] is applied and the integral is computed numerically from Dmin = 0.01 mm to Dmax = 7.0 mm in 0.01 mm increments. Simulation results suggest that for voxels with nonzero specific attenuations is between 6.75 × 10 and 5.46 × 10 mm −1 ⋅m −3 when is 3.5 mm −1 . With and ready, the specific attenuation for any elevation angle and for either polarizations can then be calculated. Figure 2b shows the true field of (40) for comparison, in which the specific attenuation for voxel A is 2.99 dB/km, 7.4% less than (40).

Retrieval of Specific Attenuation
In the simulations, we use a LEO satellite with a circular Keplerian motion trajectory of an orbital height of 1000 km and an inclination angle of 96°. The trajectory is on the same vertical plane as the rain field. There are 14 receivers evenly placed at the ground level, from −3.25 km to 3.25 km, with any two adjacent receivers being 0.5 km apart. The SNRs for the horizontally and vertically polarized links are generated using exactly the same parametric configurations as in previous studies [19,20]. It is worth noting that the estimation error of the SNRs is taken into account by using the Cramer-Rao lower bound (CRLB) [45]. This introduces an error of approximately 0.01 dB in the SNR measurements ( ) and ( ). We assume that there are 50 SNR measurements from 50 different elevation angles available, evenly distributed in the overpass of the satellite.

Retrieval of Specific Attenuation
In the simulations, we use a LEO satellite with a circular Keplerian motion trajectory of an orbital height of 1000 km and an inclination angle of 96 • . The trajectory is on the same vertical plane as the rain field. There are 14 receivers evenly placed at the ground level, from −3.25 km to 3.25 km, with any two adjacent receivers being 0.5 km apart. The SNRs for the horizontally and vertically polarized links are generated using exactly the same parametric configurations as in previous studies [19,20]. It is worth noting that the estimation error of the SNRs is taken into account by using the Cramer-Rao lower bound (CRLB) [45]. This introduces an error of approximately 0.01 dB in the SNR measurementsρ V (k) and ρ H (k). We assume that there are 50 SNR measurements from 50 different elevation angles available, evenly distributed in the overpass of the satellite.
We use a differential approach [20] and an iterative process [19] to estimate the specific attenuations of all voxels as well as the noise figure F n . Ideally, we want to retrieveâ H (40) andâ V (40) only fromρ H (40) andρ V (40). However, results suggest that links at 40 • only are not sufficient for the least-squares algorithm to solve the unknowns. On the other hand, for horizontally polarized links, the change in attenuation caused by the change in the elevation angle is relatively small within the targeted range of Λ. Table 3 lists the a H (k) to a V (40) ratio for different elevation angles and different values of Λ calculated by the T-Matrix method. It shows that the change in specific attenuation as the elevation angle increases never exceeds approximately 1%. In other words, considering a H being unchanged for all elevation angles will not introduce a significant error in the retrieval algorithm. Consequently, we just utilize the measurements ofρ H from all elevation angles to estimateâ H , which will be regarded asâ H (40). Same as previous studies [19], we use an iterative process (summarized in Algorithm 1) to update the noise figure for a more accurate estimation ofâ H (40). Figure 3a shows the estimatedâ H (40) for the entire field in one simulation. It can be seen that the retrieved field shows very close agreement with the true field in Figure 2a. In 100 independent simulations, we recordâ H (40) for the 10 voxels in the 4th row from the bottom (marked in Figure 3a) with nonzero true specific attenuations. As shown in Figure 3b, the true specific attenuations are marked by the red crosses and the average retrieved specific attenuations over 100 simulations are marked by the red circles. The average relative error, which is calculated by averaging over 100 simulations the ratio of the absolute error to the true specific attenuation, is shown by the blue triangles. It can be concluded that retrieval of specific attenuations is relatively more accurate for voxels with high specific attenuations. More specifically, the relative error for the middle four voxels is below 2%. Algorithm 1. The iterative process used to update the noise figure.
Input:ρ H (k), k = 1, 2, . . . , M; Output: α H for all voxels initialize F n = 0; while not the last iteration do Use least-squares to solve α H , Calculate A I (k) using α H , Update F n using α H and A I . For vertically polarized links, the change in a V with the elevation angle will be too large to ignore. Table 4 shows the value of a V (k)/a V (40) for different values of Λ. It can be seen that a V (54.34) already exceeds a V (40) by over 4% when Λ is 2.8 mm −1 . As a result, onlŷ ρ V (k) measured close to 40 • can be used in the least-squares algorithm without introducing too much error. It should also be noted that the noise figure is now already calculated based onâ H (40) so the least-squares algorithm can be directly applied to estimateâ V (40) without engaging the iterative process of Algorithm 1. Although this would again introduce some error, simulations suggest that the error can be considered negligible. Multiple simulations suggest that the retrieved field is erroneous due to underdetermination in the least-squares algorithm when less than eight SNR measurements from 40 • (ρ V (40) toρ V (51.2)) are used. Figure 3c shows the retrieval ofâ V (40) for the same 10 voxels as in Figure 3b. Here, eight SNR measurements are used but the relative errors for individual voxels are much higher than in the retrieval of a H . For instance, for voxel 5, the relative error of estimatedâ V (40) is 39%. This indicates that using estimatedâ V (40)/â H (40) of individual voxels to calculate the slope parameter Λ will lead to very large errors. As a result, we propose a new model for retrievingâ V (40) and thusâ V (40)/â H (40), taking advantage of the assumption that Λ does not change within a layer.  For vertically polarized links, the change in with the elevation angle will be too large to ignore. Table 4 shows the value of ( )/ (40) for different values of . It can be seen that (54.34) already exceeds (40) by over 4% when is 2.8 mm −1 . As a result, only ( ) measured close to 40 ° can be used in the least-squares algorithm without introducing too much error. It should also be noted that the noise figure is now already calculated based on (40) so the least-squares algorithm can be directly applied to estimate (40) without engaging the iterative process of Algorithm 1. Although this would again introduce some error, simulations suggest that the error can be considered negligible. Multiple simulations suggest that the retrieved field is erroneous due to underdetermination in the least-squares algorithm when less than eight SNR measurements from 40° ( (40) to (51.2)) are used. Figure 3c shows the retrieval of (40) for the same 10 voxels as in Figure 3b. Here, eight SNR measurements are used but the relative errors for individual voxels are much higher than in the retrieval of . For instance, for voxel 5, the relative error of estimated (40) is 39%. This indicates that using estimated (40)/ (40) of individual voxels to calculate the slope parameter Λ will lead to very Before retrivingâ V (40), the entire field is divided into three layers according to Λ and then nine areas as shown in Figure 3d, with each of area being 2 km in width. The same least-squares algorithm is applied to retrieve the specific attenuation for vertical polarization for each area. Figure 3d shows the retrievedâ V (40) field in one simulation, in which SNR measurements at eight different elevation angles (ρ V (40) toρ V (51.2)) are used. For each area,â H (40) is then calculated by averaging the previously retrievedâ H (40) over all of the voxels within the area. For example,â H (40) for area I is the average of eight voxels in the top left in Figure 3a and for area V it is the average of 16 voxels in the center. The calculatedâ H (40) for each of the areas is also shown in Figure 3d. Note that the relative errors inâ H (40) for the center 16 voxels are very low soâ H (40) for area V is expected to be very accurate. Hence, usingâ V (40)/â H (40) of area V is the best choice for estimating the slope parameter Λ 2 .
In order to achieve the best result of solvingâ V (40), we need to use only SNR measurements close to a 40 • elevation angle but also to ensure there are enough data points for the least-squares algorithm to work. Nine groups of simulations were carried out to find the optimal balance between the two requirements. In each group, a certain number of SNR measurements from 40 • onwards are used to retrieveâ V (40), i.e., two measurements (ρ V (40) andρ V (41.4)) are used for the first group and sequentially ten measurements (ρ V (40) toρ V (54.3)) for the ninth group. Each group contains 100 independent simulations. Figure 4 shows the mean and standard deviation ofâ V (40) for area V in the 100 simulations for each group. It can be seen that the standard deviation is very large when less than five SNR measurements are used, indicating that the retrievedâ V (40) could be erroneous. These simulations suggest that using 7 to 10 SNR measurements will lead to the best overall results for retrievedâ V (40).

Retrieval of DSD Parameters
As the ratio of a V (40)/a H (40) holds the key to retrieving the slope parameter Λ, we examine the retrieved ratio for different slope parameters. Because the retrieval of specific attenuations is relatively more accurate for areas with higher values, we will focus on area V in Figure 3d and evaluate the performance of retrieving Λ 2 . In total 15 simulation groups are designed, in each of which parameter Λ 2 takes values from 2.8 mm −1 to 4.2 mm −1 in 0.1 mm −1 increments. In each of the simulation groups, 100 independent simulations are carried out. Figure 5a shows the retrievedâ V (40)/â H (40) for area V, in which the blue circles mark the average retrieved ratio across 100 simulations and the error bars show the standard deviation. It is found that the relationship between true Λ 2 and the mean of the retrievedâ V (40)/â H (40) is almost linear. Using a linear least-squares regression fit of the blue circles, this linear relationship is generated and shown by the black line and the equation in the figure.
Our following simulations utilize the linear relationship above to infer Λ 2 from the retrievedâ V (40)/â H (40). In the new 141 groups of 100 independent simulations, the true parameter Λ 2 takes values from 2.8 mm −1 to 4.2 mm −1 in 0.01 mm −1 increments, and the retrieved Λ 2 is calculated from the retrievedâ V (40)/â H (40) using the linear equation in Figure 5a. Figure 5b shows the retrieval outcome, in which the blue line is the average retrieved Λ 2 and the black line is a reference line with a gradient of 1. The red bars show the standard deviation across 100 simulations when the true Λ 2 is at some selected values.
Furthermore, the intercept parameter N 0 for each of the 16 voxels in area V of Figure 3d is calculated from the retrieved a H (40) and Λ 2 , and then compared with its true value. In Figure 5c, the retrieved N 0 is compared with its true value for voxel A in Figure 2a in 15 different Λ 2 (true Λ 2 being from 2.8 mm −1 to 4.2 mm −1 in 0.1 mm −1 increments) with each one containing 20 independent simulations. The black line is also a reference line with a gradient of 1. The same was carried out for voxel B in Figure 2a where the true N 0 is less than that for voxel A, and the results are shown in Figure 5d. It can be concluded from the two figures that the retrieved N 0 has a close agreement with the true N 0 for the 15 different Λ 2 . More simulations confirm that the retrieval outcomes are very similar for all of the 16 voxels in area V. As the retrievedâ H (40) is quite accurate, the level of error in retrieved N 0 is mainly due to the fluctuation in retrieved Λ 2 , which can be seen in Figure 5b. the optimal balance between the two requirements. In each group, a certain number of SNR measurements from 40° onwards are used to retrieve (40), i.e., two measurements ( (40) and (41.4)) are used for the first group and sequentially ten measurements ( (40) to (54.3)) for the ninth group. Each group contains 100 independent simulations. Figure 4 shows the mean and standard deviation of (40) for area V in the 100 simulations for each group. It can be seen that the standard deviation is very large when less than five SNR measurements are used, indicating that the retrieved (40) could be erroneous. These simulations suggest that using 7 to 10 SNR measurements will lead to the best overall results for retrieved (40).

Retrieval of DSD Parameters
As the ratio of (40)/ (40) holds the key to retrieving the slope parameter , we examine the retrieved ratio for different slope parameters. Because the retrieval of specific attenuations is relatively more accurate for areas with higher values, we will focus on area V in Figure 3d and evaluate the performance of retrieving . In total 15 simulation groups are designed, in each of which parameter takes values from 2.8 mm −1 to 4.2 mm −1 in 0.1 mm −1 increments. In each of the simulation groups, 100 independent simulations are carried out. Figure 5a shows the retrieved

Simulations of Another Rain Event
A second rain event is utilized to further verify the performance of the proposed method. Shown in Figure 6a, it is an idealized convective storm generated by WRF using the same resolution settings as before. The vertical attenuation field measures 8 km (width) by 5 km (height). The highest specific attenuation (for horizontal polarization at 40 • ) in all voxels is larger than the frontal rain field in Figure 2a. The DSD profile for each voxel is generated following the exact same process. The entire field is redivided into 4 × 3 (marked by white solid lines in Figure 6a) areas, two of which (marked by white dashed lines in Figure 6a) are used to retrieve the DSD parameters. Figure 6b shows the retrieved a V (40)/â H (40) when Λ 2 takes different values, in which the blue circles mark the average retrieved ratio across 100 simulations and the error bars show the standard deviation. Again, using a linear least-squares regression fit of the blue circles, a linear relationship is identified and shown by the black line and the equation in Figure 6b. Figure 6c shows the retrieved Λ 2 using the identified linear relationship, in which the blue line is the average retrieved Λ 2 and the black line is a reference line with a gradient of 1. The red bars show the standard deviation across 100 simulations when the true Λ 2 is at the same selected values as in Figure 5b. Figure 6d is the scatter plot of the retrieved N 0 for voxel A in Figure 6a.
The results above suggest that the performance of retrieval is consistent for both the synthetic convective storm and the simulated frontal rainfall. The proposed methods display great potential in retrieving the slope and intercept parameters of DSD. Although only the retrieval of vertical DSD fields is discussed in this paper, it can be inferred that 3-D retrieval can be achieved through the aggregation of a series of retrieved vertical fields.

Conclusions
In this paper, we use computer simulations to explore the feasibility of using dualpolarized microwave signals from LEO satellites to measure the raindrop size distribution. Our analysis suggests that measurements of attenuations on both horizontally and vertically polarized microwave links offer the opportunity to estimate two independent DSD parameters. As the difference between the two specific attenuations is crucial for retrieving the slope parameter and becomes less distinct as the elevation angle of the link goes higher, its retrieval needs to be performed at low elevation angles, preferably close to the minimum elevation angle of the receivers. Simulation results show that using SNR measurements from all elevation angles to retrieve the specific attenuation for horizon-

Conclusions
In this paper, we use computer simulations to explore the feasibility of using dualpolarized microwave signals from LEO satellites to measure the raindrop size distribution. Our analysis suggests that measurements of attenuations on both horizontally and vertically polarized microwave links offer the opportunity to estimate two independent DSD parameters. As the difference between the two specific attenuations is crucial for retrieving the slope parameter and becomes less distinct as the elevation angle of the link goes higher, its retrieval needs to be performed at low elevation angles, preferably close to the minimum elevation angle of the receivers. Simulation results show that using SNR measurements from all elevation angles to retrieve the specific attenuation for horizontally polarized signals through a least-squares algorithm is feasible. Particularly for voxels with heavy rain, the retrieval error is less than 2%. It is also confirmed that the specific attenuation for vertically polarized signals can be retrieved using a few SNR measurements taken close to the minimum elevation angle. In order to achieve the precision needed for estimating the DSD parameters, the retrieval algorithm is performed on a large area that covers several voxels. It is suggested that the specific attenuation ratio of vertically to horizontally polarized signals can be used through a linear relationship to estimate the slope parameter of the DSD, and that the intercept parameter can also be retrieved using the estimated slope parameter. The retrieved values for both parameters closely agree with their true values.
As the estimation of DSD parameters demands high accuracy in the specific attenuation measurements, the following sources of error in the system need to be further examined in future work. The first source of error comes from the definition of the voxels, which inevitably introduces a difference from the true field. This difference can be reduced by using a finer resolution grid, but it also means that the retrieval task will become more difficult. The second source of error is the unknown sky noise. The attenuation field is initially retrieved by assuming the sky noise is zero. Although the iterative process can update the sky noise and fine-tune the retrieval result, the final error is still significant. The third is the SNR estimation error, which can be reduced by designing good estimators but cannot be completely eliminated. We expect that the proposed approach can provide even better DSD estimations if these sources of error can be further mitigated.