Next Article in Journal
Learning Hierarchical Representations with Spike-and-Slab Inception Network
Next Article in Special Issue
Rainfall Intensity and Quantity Estimation Method Based on Gamma-Dose Rate Monitoring
Previous Article in Journal
SiamMFC: Visual Object Tracking Based on Mainfold Full Convolution Siamese Network
Previous Article in Special Issue
On the Use of Dynamic Calibration to Correct Drop Counter Rain Gauge Measurements
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Communication

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

Department of Electrical, Electronic and Computer Engineering, The University of Western Australia, Perth 6009, Australia
*
Author to whom correspondence should be addressed.
Sensors 2021, 21(19), 6389; https://doi.org/10.3390/s21196389
Submission received: 12 August 2021 / Revised: 21 September 2021 / Accepted: 22 September 2021 / Published: 24 September 2021
(This article belongs to the Special Issue Rain Sensors)

Abstract

:
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.

1. 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-to-noise 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:
N ( D ) = N 0 D μ e Λ D ,
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 path-averaged 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 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.

2. Model and Methods

2.1. 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 satellite-to-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]
ρ ( k ) = C ( k ) F n ( k ) A I ( k ) A F ( k ) ,
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.

2.2. Specific Attenuation and DSD

The relationship between the specific attenuation ( a in dB/km) and the DSD of (1) is given by [35]
a = 4.343 × 0 N ( D ) Q ( D ) d D ,
where Q ( D ) is the extinction cross section of the raindrop (in m2), 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.
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
y x = 1.03 0.062 D ,
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.

2.3. 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
a V a H = 0 D μ e Λ D Q V ( D ) d D 0 D μ e Λ D Q H ( D ) d D   ,
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).

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.

3.1. 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 Dmin = 0.01 mm to Dmax = 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 ) .

3.2. 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 a ^ H ( 40 ) and a ^ 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 a ^ H , which will be regarded as a ^ 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 a ^ H ( 40 ) . Figure 3a shows the estimated a ^ 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 a ^ 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, only ρ ^ 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 a ^ H ( 40 ) so the least-squares algorithm can be directly applied to estimate a ^ 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 a ^ 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 a ^ V ( 40 )   is 39%. This indicates that using estimated a ^ V ( 40 ) / a ^ 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 a ^ V ( 40 ) and thus a ^ V ( 40 ) / a ^ H ( 40 ) , taking advantage of the assumption that Λ does not change within a layer.
Before retriving a ^ 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 a ^ 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, a ^ H ( 40 ) is then calculated by averaging the previously retrieved a ^ H ( 40 ) over all of the voxels within the area. For example, a ^ 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 a ^ H ( 40 ) for each of the areas is also shown in Figure 3d. Note that the relative errors in a ^ H ( 40 ) for the center 16 voxels are very low so a ^ H ( 40 ) for area V is expected to be very accurate. Hence, using a ^ V ( 40 ) / a ^ H ( 40 ) of area V is the best choice for estimating the slope parameter Λ 2 .
In order to achieve the best result of solving a ^ 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 a ^ 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 a ^ 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 a ^ V ( 40 ) could be erroneous. These simulations suggest that using 7 to 10 SNR measurements will lead to the best overall results for retrieved a ^ V ( 40 ) .

3.3. 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 a ^ V ( 40 ) / a ^ 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 a ^ V ( 40 ) / a ^ 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 a ^ V ( 40 ) / a ^ 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 a ^ V ( 40 ) / a ^ 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 a ^ 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.

3.4. 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 ) / a ^ 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.

4. Conclusions

In this paper, we use computer simulations to explore the feasibility of using dual-polarized 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.

Author Contributions

Conceptualization, D.D.H. and X.S.; methodology and simulation, X.S.; writing—original draft preparation, X.S.; writing—review and editing, D.D.H. and X.S.; supervision, D.D.H.; project administration, D.D.H. Both authors have read and agreed to the published version of the manuscript.

Funding

This research was funded in part by Australian Research Council (ARC), under Grant DP190100786.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Atlas, D.; Ulbrich, C.W. Path- and Area-Integrated Rainfall Measurement by Microwave Attenuation in the 1–3 cm Band. J. Appl. Meteorol. 1977, 16, 1322–1331. [Google Scholar] [CrossRef] [Green Version]
  2. Ippolito, L.J. Propagation Effects Handbook for Satellite Systems; NASA: Washington, DC, USA, 1989; Volume 1082.
  3. Giuli, D.; Toccafondi, G.; Gentili, B.; Freni, A. Tomographic Reconstruction of Rainfall Fields through Microwave Attenuation Measurements. J. Appl. Meteorol. 1991, 30, 1323–1340. [Google Scholar] [CrossRef] [Green Version]
  4. Messer, H.; Zinevich, A.; Alpert, P. Environmental Monitoring by Wireless Communication Networks. Science 2006, 312, 713. [Google Scholar] [CrossRef] [Green Version]
  5. Leijnse, H.; Uijlenhoet, R.; Stricker, J.N.M. Rainfall measurement using radio links from cellular communication networks. Water Resour. Res. 2007, 43, W03201. [Google Scholar] [CrossRef]
  6. Goldshtein, O.; Messer, H.; Zinevich, A. Rain Rate Estimation Using Measurements from Commercial Telecommunications Links. IEEE Trans. Signal Process. 2009, 57, 1616–1625. [Google Scholar] [CrossRef]
  7. Zinevich, A.; Pinhas, A.; Messer, H. Estimation of rainfall fields using commercial microwave communication networks of variable density. Adv. Water Resour. 2008, 31, 1470–1480. [Google Scholar] [CrossRef]
  8. Messer, H.; Zinevich, A.; Alpert, P. Environmental sensor networks using existing wireless communication systems for rainfall and wind velocity measurements. IEEE Instrum. Meas. Mag. 2012, 15, 32–38. [Google Scholar] [CrossRef]
  9. Messer, H.; Sendik, O. A New Approach to Precipitation Monitoring: A critical survey of existing technologies and challenges. IEEE Signal Process. Mag. 2015, 32, 110–122. [Google Scholar] [CrossRef]
  10. Barthès, L.; Mallet, C. Rainfall measurement from the opportunistic use of an Earth–space link in the Ku band. Atmos. Meas. Tech. 2013, 6, 2181–2193. [Google Scholar] [CrossRef] [Green Version]
  11. Giannetti, F.; Reggiannini, R.; Moretti, M.; Adirosi, E.; Baldini, L.; Facheris, L.; Antonini, A.; Melani, S.; Bacci, G.; Petrolino, A.; et al. Real-Time Rain Rate Evaluation via Satellite Downlink Signal Attenuation Measurement. Sensors 2017, 17, 1864. [Google Scholar] [CrossRef]
  12. Gharanjik, A.; Bhavan Shankar, M.R.; Zimmer, F.; Ottersten, F. Centralized Rainfall Estimation Using Carrier to Noise of Satellite Communication Links. IEEE J. Sel. Areas Commun. 2018, 36, 1065–1073. [Google Scholar] [CrossRef] [Green Version]
  13. Arslan, C.H.; Kuletgin, A.; Urbino, J.V.; Dyrud, L. Satellite-Link Attenuation Measurement Technique for Estimating Rainfall Accumulation. IEEE Trans. Geosci. Remote Sens. 2018, 56, 681–693. [Google Scholar] [CrossRef]
  14. Su, Y.; Liu, Y.; Zhou, Y.; Yuan, J.; Cao, H.; Shi, C. Broadband LEO Satellite Communications: Architectures and Key Technologies. IEEE Wirel. Commun. 2019, 26, 55–61. [Google Scholar] [CrossRef]
  15. Pachler, N.; Del Portillo, I.; Crawley, E.F.; Cameron, B.G. An Updated Comparison of Four Low Earth Orbit Satellite Constellation Systems to Provide Global Broadband. In Proceedings of the 2021 IEEE International Conference on Communications Workshops (ICC Workshops), Montreal, QC, Canada, 14–23 June 2021. [Google Scholar]
  16. Lemorton, J.; Castanet, L.; Lacoste, F.; Riva, C.; Matricciani, E.; Fiebig, U.-C.; Van de Kamp, M.; Martellucci, A. Development and validation of time-series synthesizers of rain attenuation for Ka-band and Q/V-band satellite communication systems. Int. J. Satell. Commun. Netw. 2007, 25, 575–601. [Google Scholar] [CrossRef] [Green Version]
  17. Huang, D.; Lu, X.; Feng, X.; Wang, W. A hypothesis of 3D rainfall tomography using satellite signals. J. Commun. Inf. Netw. 2016, 1, 134–142. [Google Scholar]
  18. Kak, A.C.; Slaney, M. Principles of Computerized Tomographic Imaging; SIAM: Philadelphia, PA, USA, 2001. [Google Scholar]
  19. Shen, X.; Huang, D.D.; Song, B.; Vincent, C.; Togneri, R. 3-D Tomographic Reconstruction of Rain Field Using Microwave Signals from LEO Satellites: Principle and Simulation Results. IEEE Trans. Geosci. Remote Sens. 2019, 57, 5434–5446. [Google Scholar] [CrossRef]
  20. Shen, X.; Huang, D.; Vincent, C.; Wang, W.; Tigneri, R. A Differential Approach for Rain Field Tomographic Reconstruction Using Microwave Signals from Leo Satellites. In Proceedings of the ICASSP 2020—2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Barcelona, Spain, 4–8 May 2020. [Google Scholar]
  21. Song, B.; Huang, D.; Xi, S.; Togneri, R. Performance Analysis for Path Attenuation Estimation of Microwave Signals Due to Rainfall and Beyond. In Proceedings of the ICASSP 2020—2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Barcelona, Spain, 4–8 May 2020. [Google Scholar]
  22. Xu, L.; Huang, D.; Feng, X.; Wang, W. Tomographic reconstruction of rainfall fields using satellite communication links. In Proceedings of the 2017 23rd Asia-Pacific Conference on Communications (APCC), Perth, Australia, 11–13 December 2017. [Google Scholar]
  23. Ulbrich, C.W. Natural Variations in the Analytical Form of the Raindrop Size Distribution. J. Appl. Meteorol. Climatol. 1983, 22, 1764–1775. [Google Scholar] [CrossRef] [Green Version]
  24. Holt, A.R.; Upton, J.W.F.; Willis, G.J.G.; Rahimi, A.R.; Baxter, P.D.; Collier, C.G. Measurement of rainfall by dual-wavelength microwave attenuation. Electron. Lett. 2000, 36, 2099–2101. [Google Scholar] [CrossRef]
  25. Rincon, R.F.; Lang, R.H. Microwave link dual-wavelength measurements of path-average attenuation for the estimation of drop size distributions and rainfall. IEEE Trans. Geosci. Remote Sens. 2002, 40, 760–770. [Google Scholar] [CrossRef]
  26. Ruf, C.S.; Kultegin, A.; Mathur, S.; Bobak, J.P. 35-GHz Dual-Polarization Propagation Link for Rain-Rate Estimation. J. Atmos. Ocean. Technol. 1996, 13, 419–425. [Google Scholar] [CrossRef] [Green Version]
  27. Aydin, K.; Daisley, S.E.A. Relationships between rainfall rate and 35-GHz attenuation and differential attenuation: Modeling the effects of raindrop size distribution, canting, and oscillation. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2343–2352. [Google Scholar] [CrossRef]
  28. Song, K.; Liu, X.; Gao, T.; He, B. Raindrop Size Distribution Retrieval Using Joint Dual-Frequency and Dual-Polarization Microwave Links. Adv. Meteorol. 2019, 2019, 7251870. [Google Scholar] [CrossRef]
  29. Kozu, T.; Nakamura, K. Rainfall Parameter Estimation from Dual-Radar Measurements Combining Reflectivity Profile and Path-integrated Attenuation. J. Atmos. Ocean. Technol. 1991, 8, 259–270. [Google Scholar] [CrossRef] [Green Version]
  30. Kumar, L.S.; Lee, Y.H.; Ong, J.T. Two-Parameter Gamma Drop Size Distribution Models for Singapore. IEEE Trans. Geosci. Remote Sens. 2011, 49, 3371–3380. [Google Scholar] [CrossRef]
  31. Munchak, S.J.; Tokay, A. Retrieval of Raindrop Size Distribution from Simulated Dual-Frequency Radar Measurements. J. Appl. Meteorol. Climatol. 2008, 47, 223–239. [Google Scholar] [CrossRef]
  32. Pruppacher, H.R.; Beard, K.V. A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air. Q. J. R. Meteorol. Soc. 1970, 96, 247–256. [Google Scholar] [CrossRef]
  33. Mishchenko, M.I.; Travis, L.D.; Mackowski, D.W. T-matrix computations of light scattering by nonspherical particles: A review. J. Quant. Spectrosc. Radiat. Transf. 1996, 55, 535–575. [Google Scholar] [CrossRef]
  34. Shen, X.; Defeng, D.D.; Wang, W.; Prein, A.F.; Togneri, R. Retrieval of Cloud Liquid Water Using Microwave Signals from LEO Satellites: A Feasibility Study through Simulations. Atmosphere 2020, 11, 460. [Google Scholar] [CrossRef]
  35. Hogg, D.C.; Ta-Shing, C. The role of rain in satellite communications. Proc. IEEE 1975, 63, 1308–1331. [Google Scholar] [CrossRef]
  36. Beard, K.V.; Bringi, V.N.; Thurai, M. A new understanding of raindrop shape. Atmos. Res. 2010, 97, 396–415. [Google Scholar] [CrossRef]
  37. Mishchenko, M.I.; Travis, L.D. Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers. J. Quant. Spectrosc. Radiat. Transf. 1998, 60, 309–324. [Google Scholar] [CrossRef]
  38. Waterman, P.C. Matrix formulation of electromagnetic scattering. Proc. IEEE 1965, 53, 805–812. [Google Scholar] [CrossRef]
  39. Shen, X.; Huang, D.D.; Xu, L.; Tang, X. Reconstruction of vertical rainfall fields using satellite communication links. In Proceedings of the 2017 23rd Asia-Pacific Conference on Communications (APCC), Perth, WA, Australia, 11–13 December 2017; pp. 1–6. [Google Scholar]
  40. Tang, Q.; Xiao, H.; Guo, C.; Feng, L. Characteristics of the raindrop size distributions and their retrieved polarimetric radar parameters in northern and southern China. Atmos. Res. 2014, 135, 59–75. [Google Scholar] [CrossRef]
  41. Yang, Q.; Dai, Q.; Han, D.; Chen, Y.; Zhang, S. Sensitivity analysis of raindrop size distribution parameterizations in WRF rainfall simulation. Atmos. Res. 2019, 228, 1–13. [Google Scholar] [CrossRef]
  42. Van Leth, T.C.; Leijnse, H.; Overeem, A.; Uijlenhoet, R. Estimating raindrop size distributions using microwave link measurements: Potential and limitations. Atmos. Meas. Tech. 2020, 13, 1797–1815. [Google Scholar] [CrossRef] [Green Version]
  43. Ryzhkov, A.; Pinsky, M.; Pokrovsky, A.; Khain, A. Polarimetric Radar Observation Operator for a Cloud Model with Spectral Microphysics. J. Appl. Meteorol. Climatol. 2011, 50, 873–894. [Google Scholar] [CrossRef]
  44. Somerville, W.R.C.; Auguié, B.; Le Ru, E.C. Smarties: User-friendly codes for fast and accurate calculations of light scattering by spheroids. J. Quant. Spectrosc. Radiat. Transf. 2016, 174, 39–55. [Google Scholar] [CrossRef] [Green Version]
  45. Kay, S.M. Fundamentals of Statistical Signal Processing: Estimation Theory; Prentice-Hall: Hoboken, NJ, USA, 1993. [Google Scholar]
Figure 1. An oblate raindrop has different extinction cross sections and causes different attenuations for vertically and horizontally polarized signals.
Figure 1. An oblate raindrop has different extinction cross sections and causes different attenuations for vertically and horizontally polarized signals.
Sensors 21 06389 g001
Figure 2. (a) The true specific attenuation field for horizontally polarized signals with 40° elevation angle. Warmer colors indicate higher specific attenuations. (b) The true specific attenuation field for vertically polarized signals with 40° elevation angle.
Figure 2. (a) The true specific attenuation field for horizontally polarized signals with 40° elevation angle. Warmer colors indicate higher specific attenuations. (b) The true specific attenuation field for vertically polarized signals with 40° elevation angle.
Sensors 21 06389 g002
Figure 3. (a) The retrieved specific attenuation field for horizontally polarized signals at 40 degrees. (b) Retrieval results of horizontal specific attenuations for the ten voxels in the 4th row from the bottom in 100 simulations: the true specific attenuations are marked by red crosses and the average retrieved specific attenuations are marked by red circles. The blue triangles show the relative error. (c) Retrieval results of vertical-specific attenuation for the ten voxels in the 4th row from the bottom in 100 simulations: the true specific attenuations are marked by red crosses and the average retrieved specific attenuations are marked by red circles. The blue triangles show the relative error. (d) The newly divided nine areas and the retrieved specific attenuations for both horizontal and vertical polarizations.
Figure 3. (a) The retrieved specific attenuation field for horizontally polarized signals at 40 degrees. (b) Retrieval results of horizontal specific attenuations for the ten voxels in the 4th row from the bottom in 100 simulations: the true specific attenuations are marked by red crosses and the average retrieved specific attenuations are marked by red circles. The blue triangles show the relative error. (c) Retrieval results of vertical-specific attenuation for the ten voxels in the 4th row from the bottom in 100 simulations: the true specific attenuations are marked by red crosses and the average retrieved specific attenuations are marked by red circles. The blue triangles show the relative error. (d) The newly divided nine areas and the retrieved specific attenuations for both horizontal and vertical polarizations.
Sensors 21 06389 g003
Figure 4. Using different numbers of SNR measurements to retrieve the specific attenuation for vertically polarized link at 40 degrees. The circles mark the average retrieved specific attenuations for area V across 100 independent simulations and the error bars show the standard deviation.
Figure 4. Using different numbers of SNR measurements to retrieve the specific attenuation for vertically polarized link at 40 degrees. The circles mark the average retrieved specific attenuations for area V across 100 independent simulations and the error bars show the standard deviation.
Sensors 21 06389 g004
Figure 5. (a) The retrieved specific attenuation ratio when the slope parameter for area V takes different values. Each blue circle represents the mean for 100 simulations, and each error bar shows the standard deviation. The black line is the linear regression result over the blue circles, which is also given by the equation. (b) The retrieved slope parameter using the linear relationship. The blue line shows the average for 100 simulations and the red bars show the standard deviation at certain values. (c) Scatter plot of the retrieved intercept parameter for voxel A. (d) Scatter plot of the retrieved intercept parameter for voxel B.
Figure 5. (a) The retrieved specific attenuation ratio when the slope parameter for area V takes different values. Each blue circle represents the mean for 100 simulations, and each error bar shows the standard deviation. The black line is the linear regression result over the blue circles, which is also given by the equation. (b) The retrieved slope parameter using the linear relationship. The blue line shows the average for 100 simulations and the red bars show the standard deviation at certain values. (c) Scatter plot of the retrieved intercept parameter for voxel A. (d) Scatter plot of the retrieved intercept parameter for voxel B.
Sensors 21 06389 g005
Figure 6. (a) Second rain event: the true specific attenuation field for horizontally polarized signals with 40° elevation angle. (b) The retrieved specific attenuation ratio for different slope parameter values. Each blue circle represents the mean for 100 simulations, and each error bar shows the standard deviation. The black line is the linear regression result over the blue circles, which is also given by the equation. (c) The retrieved slope parameter using the linear relationship. The blue line shows the average for 100 simulations and the red bars show the standard deviation at certain values. (d) Scatter plot of the retrieved intercept parameter for voxel A.
Figure 6. (a) Second rain event: the true specific attenuation field for horizontally polarized signals with 40° elevation angle. (b) The retrieved specific attenuation ratio for different slope parameter values. Each blue circle represents the mean for 100 simulations, and each error bar shows the standard deviation. The black line is the linear regression result over the blue circles, which is also given by the equation. (c) The retrieved slope parameter using the linear relationship. The blue line shows the average for 100 simulations and the red bars show the standard deviation at certain values. (d) Scatter plot of the retrieved intercept parameter for voxel A.
Sensors 21 06389 g006
Table 1. The extinction cross section for horizontal polarization calculated by the T-Matrix method for various raindrop sizes and elevation angles.
Table 1. The extinction cross section for horizontal polarization calculated by the T-Matrix method for various raindrop sizes and elevation angles.
Q H ( D )
( × 10 8   m 2 )
Elevation Angle (°)
40557090
D
(mm)
0.51.13371.13361.13351.1335
1.023.63423.49923.38923.327
2.0509.52507.46505.76504.82
4.04019.04139.84243.74302.9
6.09064.29318.29487.69562.0
Table 2. The extinction cross section ratio of vertical to horizontal polarizations calculated by the T-Matrix method for various raindrop sizes and elevation angles.
Table 2. The extinction cross section ratio of vertical to horizontal polarizations calculated by the T-Matrix method for various raindrop sizes and elevation angles.
Q V ( D ) / Q H ( D ) Elevation Angle (°)
40557090
D
(mm)
0.50.99890.99940.99981
1.00.97000.98310.99391
2.00.90450.94610.98071
4.00.85850.92200.97271
6.00.83660.91910.97461
Table 3. The specific attenuation for a horizontally polarized link changes while the elevation angle increases from 40 degrees. Nonetheless, the maximum change is only about 1%.
Table 3. The specific attenuation for a horizontally polarized link changes while the elevation angle increases from 40 degrees. Nonetheless, the maximum change is only about 1%.
a H ( k ) /
a H ( 40 )
Λ ( mm 1 )
2.83.33.84.2
k
(°)
42.801.00000.99960.99920.9991
54.341.00000.99770.99610.9955
66.681.00010.99590.99330.9921
88.721.00030.99430.99090.9893
Table 4. The specific attenuation for a vertically polarized link changes while the elevation angle increases from 40 degrees. The change is much more significant compared to that for a horizontally polarized link.
Table 4. The specific attenuation for a vertically polarized link changes while the elevation angle increases from 40 degrees. The change is much more significant compared to that for a horizontally polarized link.
a V ( k ) /
a V ( 40 )
Λ ( mm 1 )
2.83.33.84.2
k
(°)
42.801.00901.00641.00521.0044
54.341.04361.03441.02671.0216
66.681.07731.06001.04671.0381
88.721.10571.08191.06391.0523
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Shen, X.; Huang, D.D. Retrieval of Raindrop Size Distribution Using Dual-Polarized Microwave Signals from LEO Satellites: A Feasibility Study through Simulations. Sensors 2021, 21, 6389. https://doi.org/10.3390/s21196389

AMA Style

Shen X, Huang DD. Retrieval of Raindrop Size Distribution Using Dual-Polarized Microwave Signals from LEO Satellites: A Feasibility Study through Simulations. Sensors. 2021; 21(19):6389. https://doi.org/10.3390/s21196389

Chicago/Turabian Style

Shen, Xi, and Defeng David Huang. 2021. "Retrieval of Raindrop Size Distribution Using Dual-Polarized Microwave Signals from LEO Satellites: A Feasibility Study through Simulations" Sensors 21, no. 19: 6389. https://doi.org/10.3390/s21196389

APA Style

Shen, X., & Huang, D. D. (2021). Retrieval of Raindrop Size Distribution Using Dual-Polarized Microwave Signals from LEO Satellites: A Feasibility Study through Simulations. Sensors, 21(19), 6389. https://doi.org/10.3390/s21196389

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop