Rainfall Map from Attenuation Data Fusion of Satellite Broadcast and Commercial Microwave Links

The demand for accurate rainfall rate maps is growing ever more. This paper proposes a novel algorithm to estimate the rainfall rate map from the attenuation measurements coming from both broadcast satellite links (BSLs) and commercial microwave links (CMLs). The approach we pursue is based on an iterative procedure which extends the well-known GMZ algorithm to fuse the attenuation data coming from different links in a three-dimensional scenario, while also accounting for the virga phenomenon as a rain vertical attenuation model. We experimentally prove the convergence of the procedures, showing how the estimation error decreases for every iteration. The numerical results show that adding the BSL links to a pre-existent CML network boosts the accuracy performance of the estimated rainfall map, improving up to 50% the correlation metrics. Moreover, our algorithm is shown to be robust to errors concerning the virga parametrization, proving the possibility of obtaining good estimation performance without the need for precise and real-time estimation of the virga parameters.


Introduction
In recent years, the higher and higher occurrence of extreme phenomena related to climate change has considerably spurred the demand of accurate and real-time rainfall maps. To this end, classical methods for rainfall measurement, i.e., surface sensors, weather radars and satellite systems have been intensively used. These solutions, however, require non-negligible installation and operating costs, while the measurements they offer have a temporal and spatial resolution that is often not sufficient for the tasks of interest.
In the last decades, the use of opportunistic sensors for rainfall estimation has opened the way to various methods based on signal processing techniques applied to attenuation measurements of existing commercial microwave links (CMLs). Several different strategies have been published on this topic such as [1,2], and in the sequel, [3][4][5]. As an alternative option, it has also been proposed to obtain the rain attenuation contribution from the overall attenuation measured at the CML receiver [6][7][8]. In addition, exploiting the measurements of the attenuation data available at the CML receiver, the generation of rainfall maps through the inverse distance weighting (IDW) algorithm or tomographic estimation has been addressed in [9][10][11] and in [1,12,13], respectively.
Recently, a significant research effort has been devoted also to the estimation of the rainfall rate from the received signal level at the ground station of direct-to-home (DTH) broadcast satellite links (BSLs) [14][15][16][17][18][19][20]. Hence, taking advantage of the limited cost and ease of installation of the commercial-grade BSL receivers for DTH broadcast, the rationale of our current work consists of effectively merging together the attenuation data coming from both the CML and BSL receivers. Differently from our preliminary approach [21], (to the best of our knowledge) the first rainfall estimation system based on fusing together such measurements, here, we propose a modified GMZ algorithm [9] which properly handles the BSL and CML attenuation data over a three-dimensional (3-D) scenario of interest, where the vertical variations of rain intensity affecting the satellite links are considered.
The main contributions are as follows.
• We consider the virga phenomenon, i.e., a variation of rainfall rate with respect to the height due to a gradient of environmental parameters such as humidity, which may cause evaporation or sublimation of rain [22]; • Based on the virga model, we propose a hybrid 3-D modified version of the GMZ algorithm [9,11] which provides the rainfall rate estimation by properly merging together the attenuation measurements collected at the BSL and CML receivers; • We numerically show that the root mean square estimation error steadily decreases at each step of the iterative procedure, thus leading to a stable solution for the rainfall intensity although the underlying optimization problem is not convex; • We prove that a few BSLs placed in locations scarcely covered by CMLs can act as gap-fillers, with the result of notably improving the rainfall estimation performance with respect to conventional schemes based on either CMLs or BSLs only; • The robustness of the proposed approach is achieved even in case of non-ideal knowledge of the parameters describing the virga phenomenon, which is significant indication that the hard task of its real-time characterization [22,23] can be avoided.
The organization of the paper is as follows. In Section 2, the environmental scenario is described, i.e., we define the geometrical position of CMLs and BSLs in Section 2.1, the rainfall model for both horizontal and vertical planes in Section 2.2, and the quantization procedure to obtain the rainfall map in Section 2.3. In Section 3, we outline the hybrid iterative optimization algorithm, in Section 4, we discuss the numerical results, and finally, some conclusions are drawn in Section 5.

Geometrical Model
The area of interest has a parallelepiped shape P with a square base of area B and height equal to the rain height h 0 , defined as in [24] and assumed to be constant over B. The coordinates (measured in km) of each point inside P are referenced as p = [x, y, z] T ∈ P, [·] T denoting the transpose of a vector, with The above-defined geographical scenario to be monitored includes N CML horizontal CMLs along with N BSL slanted BSLs, whose available data measurements are properly gathered up and fused together; consider Figures 1 and 2 as examples. After indicating with n the link index, either terrestrial or satellite, and denoting the overall number of links as N = N CML + N BSL , the n-th link is geometrically delimited by p n,1 = [x n,1 , y n,1 , z n,1 ] T ∈ P and p n,2 = [x n,2 , y n,2 , z n,2 ] T ∈ P, 1 ≤ n ≤ N.
Exploiting the former equations, it can be further remarked that: (iv) θ n = 0 for all the CMLs; (v) The BSLs are characterized by a slanted path between the satellite transmitter and the ground receiver; (vi) Due to the distance of the geostationary satellite from the scenario and assuming that √ B is around few kilometers, θ n = θ, for all the BSLs, i.e., the elevation angle can be assumed constant for all the links pointing to the same satellite. (The case of using BSLs from different satellites visible at different azimuth and elevation angles is under study.)

Rainfall Model
The specific rain attenuation, in dB/km, experienced at a location p lying along the nth wireless communication link with carrier frequency f n in the 10-30 GHz frequency range depends (with a good approximation) on the local rainfall rate r(p), in mm/h, according to the power law formula [25] α n (p) = a n [r(p) where a n and b n are empirical coefficients (assumed to be known throughout the paper) relying on f n and on the polarization, i.e., horizontal, vertical or circular, [14,26,27]. It is also worth emphasizing that the rain attenuation represents only one of the contributions to the total attenuation affecting the n-th link. State-of-the-art works, however, have shown that the rainfall contribution can be reliably extracted from the measurements of the total attenuation [7,14,15]. Furthermore, for the sake of simplicity, we assume the ITU model based on stratiform rain with two layers only, i.e, solid and liquid [24] (A more accurate model considering a third layer, named melting layer between the solid and the liquid layers, was proposed and investigated in [14,16].), as shown in Figure 1, where, for instance, the link 1 is a CML having length L 1 , attenuation A 1 and elevation angle θ 1 = 0, and the link 2 is a BSL with wet path length, i.e., the portion of the path inside the liquid layer, L 2 , attenuation A 2 and elevation angle θ 2 . It can be remarked that the solid layer marginally influences the attenuation, thus explaining the reason why the BSL link length is defined as the wet-path only.
We model the rainfall rate as the spatially-continuous random process, expressed in mm/h, r(p) ∈ P, i.e., with values depending on the position within the 3-D scenario. Hence, according to (4), the total rain attenuation, in dB, experienced by the n-th link reads as As detailed hereafter, the rainfall will be modeled differently in a generic x − y horizontal (H) plane at height z, denoted as π H (z), and in a generic x − z vertical (V) plane at coordinate y, denoted as π V (y).

The π H (z) Plane
Let us consider K points lying on the same H-plane at height z, collect them into the set K, i.e., p k = [x k , y k , z] T ∈ K, and assume that the relevant rainfall rates r(p k ), 1 ≤ k ≤ K, are known. Then, let us consider p u = [x u , y u , z] T , which lies on the same plane π H (z), but with an unknown rainfall rate. The rainfall rate at p u can be estimated according to the Shepard's inverse distance weighting (IDW) method [10,28,29] The adimensional weights W u,k ∈ [0, 1] in (6) are expressed as where [·] + = max{0, ·}, the constant Γ is the radius of influence, i.e., the radius of a circumference centered on p u and lying on π H (z), which is suitably set (In [9], the radius of influence depends on the density of the data points and is adaptively chosen so as to include at least five data points.) to encompass those points on the H plane whose rainfall rates are assumed to appreciably contribute to the evaluation of the rainfall rate at p u , and denotes the distance between the point with unknown rainfall rate and the k-th point of the set K. Consequently, the rainfall rate in (6) is evaluated considering only those points for which d u,k < Γ.

The π V (z) Plane
Let us take into consideration the following two points, both lying on the same vertical line, i.e., p k = [x, y, z k ] T , where the relevant rainfall rate r(p k ) is known, and p u = [x, y, z u ] T , where the relevant rainfall rate is instead unknown. Additionally, we consider the presence of a vertical gradient of the rainfall rate, the so-called virga phenomenon. An accurate yet analytically involved expression of the upwards vertical variation of the rainfall rate is provided by [30]. Such analytical model, however, involves many parameters, which vary according to the type of the precipitation, and so, in general, are hard to estimate [22]. Nevertheless, experimental results for temperate climates presented in [23] show that the dependence of the rainfall rate with the height can be accurately modeled by a simple linear law, i.e., where ν(p u , p k ) is the rainfall rate evaluated at p u by applying the linear model from the knowledge of the rainfall rate at p k . By convention, the gradient of the rainfall rate g(x, y), expressed in mm/h/km, is assumed to be positive if the rainfall rate increases with the altitude, so that the highest value is attained at height h 0 . Further, if the model (9) yields a null rain rate at a given height, the rain rates will be zero for all the points below, down to the ground. For the sake of simplicity, we assume a constant gradient all over the scenario, i.e., g(x, y) = g, ∀x, y such that p = [x, y, z] T ∈ P. (The proposed procedure can be easily generalized also to the case of a non-uniform gradient over the area of interest).

Overall Model
Merging together the assumptions on the planes π H (z) and π V (y), along with assuming that the rainfall rates r(p k ) of the set of K points p k are known, the rainfall rate at the generic point p u ∈ P can be estimated as

Quantization
The distribution of the rainfall rate at the base of the volume P, i.e., at z = 0, is obtained by building a spatially quantized two-dimensional (2-D) map of J × J pixels, each pixel with an area of ∆ = B/J 2 . The center of each pixel is denoted with c j , 1 ≤ j ≤ J 2 , while the grid points consists of the overall set C = {c j } J 2 j=1 . Hence, assuming that the rainfall rate does not significantly change over ∆, r(c j ) stands for the rainfall rate for each pixel. Therefore, the overall rainfall map the algorithm yields is given by r(C) = [r(c 1 ), . . . , r(c J 2 )] T .
Following the approach of [11], the quantization process is applied to the links as well. We subdivide each link into segments with length D where the rainfall rate can be assumed to be nearly constant, thus obtaining for the n-th link a number of Q n intervals equal to where L n cos(θ n ) is the projection of the BSL or CML link onto the base, and · takes the nearest lower integer of the argument. The center of the q-th segment, 1 ≤ q ≤ Q n , for the n-th link, 1 ≤ n ≤ N, is called data point, with coordinates d n,q = [x n,q , y n,q , z n,q ] T . All the data points of the n-th link are then collected in the set Q n = {d n,j } Q n j=1 , being the corresponding rainfall rates denoted as r(Q n ) = [r(d n,1 ), . . . , r(d n,Q n )] T , [mm/h].

Estimation Algorithm
An iterative estimation algorithm of the rainfall rate for the data points r(Q n ) and grid points r(C) within the scenario of interest is outlined by combining the IDW estimation algorithm in (10) with a proper constrained optimization problem (OP).
As a first step, let us consider the estimation of the rainfall rate of the data points of the n-th link for the i-th iteration. We employ the rainfall rate corresponding to the data points not belonging to the n-th link obtained at iteration i − 1, i.e., r(Q ) (i−1) , 1 ≤ ≤ N, = n, to estimate the rainfall rate of the n-th link through (10), denoted asr(Q n ), 1 ≤ n ≤ N. Specifically, the distance between the desired r(Q n ) and the one estimated from the other linksr(Q n ) is minimized, under the constraint that the overall rainfall attenuation A n measured over the n-th link be having approximated the integral (5) as the summation over all the intervals the n-th link has been subdivided into, whose centers are located at the corresponding data points. Hence, the OP can be formalized as arg min A n Q n a n L n cos θ n − as a 3-D generalization of the OP described in [9]. Since the constraint is not affine, the OP (14) is not convex, and thus, we argue that the optimal solution is not unique. Our approach is to apply the gradient descent-based method to converge to (at least) a local optimal solution [31], i.e., r(Q n ) (i) . As a second step, the whole optimal set of data points is obtained, running the OP (14) for each link n, until convergence is reached, or equivalently, when e i < , with an arbitrary > 0, being e i the error of the optimization procedure at iteration i-th, i = 1, 2, . . ., defined as As a third and final step, using the rainfall rate at the data points as input, the rainfall rate estimation on the grid points is performed using (10), thus obtaining the map r(C).
The proposed iterative algorithm is outlined in Algorithm 1. The rainfall rates are initialized with the values r(Q n ) (0) = R n , 1 ≤ n ≤ N, where R n = A n a n L n 1 bn is obtained by inverting (4) and assuming as an initial guess of the specific attenuation α n (p) for all the data points of the n-th link.

Numerical Results
To corroborate the effectiveness of the proposed approach, a set of simulations were run for the scenario of interest. The simulation makes use of a synthesized rain profile, and the network topology shown in Figure 2, whose parameters are presented in Table 1.
(The evaluation of the performance using experimental data is under investigation.) We simulate a square-shaped scenario with side length √ B, which contains a given number of CMLs and BSLs at known locations. For the sake of simplicity, we assume that all the links are operating with the same carrier frequency f n , and the same polarization. (We assume frequency values typically employed by CMLs [9] to be used also for the BSLs, for simplicity. While BSLs' carriers are usually in the Ku band (i.e., 10-13 Hz) [16], we remark that the effect of the operational frequency is completely described by the parameters a n and b n [27]; therefore, we may use the same frequency for both the set of links without loss of generalization.) All the (horizontal) CMLs have 0°elevation angle. Moreover, we also assume that all the BSLs are pointed toward the same the satellite, laying on the local meridian, so that, in a city located within 43°and 44°parallels (e.g., Pisa, Italy), the corresponding elevation angle result is approx. θ n = 39.5°; this value is set for any BSL terminal. Based on (4), the coefficients a n , b n , computed according to [27], are used to model the rain attenuation in the whole scenario, as in [14,26]. The value of the radius of influence Γ is computed in order to have non-zero weight for at least five other data points in the evaluation of Equation (7), according to [9]. The simulated rainfall intensity in mm/h is synthesized as a 2-D Gaussian-shaped spatial distribution on the x − y plane with standard deviation σ G , and peak value R mm/h located at p G = [x G , y G , 0] T , as in [10]. Moreover, the simulated precipitation is assumed to experience a fixed vertical gradient g G , in mm/h/km, on the x − z plane, to take into account for the virga phenomenon. The resulted precipitation rainfall rates are generated at the same position of the grid points of the output estimated map C (see Section 2.3), and it is denoted asr(C) = [r(c 1 ), . . . , r(c J 2 )] T . To obtain the quantized data points, the length D where the rainfall rate can be assumed constant (see Section 2.3) is set lower than 1/10 of the diameter of the rainfall phenomenon, as suggested in [9]. The overall accuracy performance of the estimation algorithm is quantified by both the root mean square error (RMSE) ε RMS and the (adimensional) correlation coefficient ρ ρ = cov{r(C), r(C)} std{r(C)}std{r(C)} (19) where cov{·} and std{·} denote the covariance and the standard deviation operators, respectively. The accuracy contribution provided by the slanted satellite links in the estimate of the rain intensity map is affected by how the degree of accuracy of the gradient model. It is worth recalling that accurate real-time measurements of the vertical rain gradient are difficult to achieve as they would require costly and complex instrumentation; see, e.g., [22]. To evaluate the performance of the algorithm under realistic conditions, the gradient g alg is adopted, which is not necessarily equal to the true value g G used in the generation of the synthetic rain. Numerical results are presented, indeed, under both the assumptions g alg = g G , i.e., ideal knowledge of the current virga phenomenon, and g alg = g G , i.e., imperfect knowledge. Convergence behavior of the algorithm. To quantify the convergence behavior through an analytical approach was found to be too complex, and therefore, we resort to simulations. The error e i (15) is averaged for 1000 simulation runs. The resulting average error e i is then plotted in Figure 3 for the first 100 iterations. As becomes apparent, the average error steadily decreases, thus providing the experimental evidence about the convergence of the proposed procedure. Rainfall rate estimation for a single realization. A single realization of the rainfall rate is estimated for the topology depicted in Figure 2, in the presence of a given 2-D Gaussian rain intensity profile with peak value R = 15 mm/h. Figure 4a shows the π H (0) horizontal plane of the scenario at ground level, including the following links: (i) A mesh of 21 CMLs (as white lines); (ii) The ground projections of 8 BSLs (as black lines). The rain cell is also visible at the top left of the map. The black circumference, which is centered on the peak of the precipitation and has radius σ G , is shown for reference as a core area of the precipitation. Figure 4b shows the π V (y G ) vertical plane of the scenario, where the CMLs appear as near-ground horizontal white lines, while the BSLs are the slanted black lines reaching h 0 . Also depicted is the vertical profile of the synthetic precipitation, characterized by a virga effect with gradient g R = 5 mm/h/km. Figure 4c,d offers the precipitation maps generated by the proposed algorithm, assuming a perfectly estimated vertical gradient, i.e., g alg = g G . In both cases, the rainfall intensity distribution and the number (N CML = 21) along with the geometry of the CMLs are the same. In case (c), there are no BSLs; thus, the estimation procedure relies on the pure GMZ algorithm [9,10], whose performance results in ε RMS = 5.715 and ρ = 0.470. In case (d), there are eight BSLs supplementing the measurements of the CMLs. Employing the proposed algorithm, the performance improves to ε RMS = 1.981, i.e., −34%, and ρ = 0.934, i.e., +50%. Hence, fusing BSL data through our approach leads to better estimation performance, proving the possibility of using, either already installed or purposely installed, satellite receivers as gap-fillers in CML networks. Rainfall rate estimation average performance. To quantify the average performance, 1000 different scenarios are considered, each of them with randomly positioned BSLs, while the CMLs are randomly selected from the network topology shown in Figure 2. For every scenario, 50 different random positions of the rain cell p G are generated. Figure 5 shows the near-ground overall performance in terms of ε RMS (Figure 5a) and ρ (Figure 5b) as a function of N BSL . The results are for different values of N CML , while keeping fixed R = 15 mm/h and g alg = g G . Again, when N BSL = 0, the algorithm coincides with the pure GMZ algorithm. On the other hand, we present also the results for N CML = 0 to show the performance obtained by collecting measures from pure BSL approaches, such as [14][15][16]. For both RMSE and ρ metrics, the boost in performance given by adding BSLs is greater than the one obtained by adding CML, as long as a minimum number of CML is present in the scenario. For example, for N CML = 7, we can halve the RMSE by increasing the number of BSLs from 8 to 13 (i.e., inserting five BSL terminals in the scenario). To obtain a similar boost in performance with the CML, we need to triple the number of microwave links in the network. On the contrary, when no CMLs are present in the network, the estimation procedure is not able to reach good estimation performance. This phenomenon is partially due to the limited length of the projection of the BSL on the x − y plane. In fact, to cover the whole scenario, a large number of satellite terminals are needed. Another reason for this drop in performance may be due to the (approximate) parallel links obtained by pointing all terminals toward the same satellite, which could generate errors in the estimation of the data points' rain rate. Nevertheless, the quantitative impact of this phenomenon on estimation performance is still under study. The results further support the choice of CML and BSL data fusion for improving the estimation performance.   Figure 6 illustrates ρ as a function of g alg for N CML = 21, a different number of N BSL , and with: (a) g G = 0, (b) g G = 3, (c) g G = 6 mm/h/km. The plots emphasize that the performance is slightly influenced by the error of the gradient, which gives minimal effects in the case g G = 3 mm/h/km. This proves that the proposed rainfall estimation technique does not require a precise real-time assessment of the rain gradient. In fact, just the use of a suitably selected constant value of the rain gradient (e.g., the average or a typical value) is enough to guarantee a small estimation error of the rainfall intensity.

Conclusions
In this paper, a novel hybrid procedure was illustrated to obtain the rainfall rate map estimation from fusing together the attenuation data collected at the receivers of both commercial microwave and broadcast satellite links. The proposed algorithm consists of the modified version of the GMZ one which has been properly extended to: (i) Apply the estimation procedure to a three-dimensional scenario; (ii) The virga phenomenon experienced by the BSLs; (iii) Merge the attenuation data from the mixed available CMLs and BSLs; (iv) Take into account the experimentally proven convergence of the iterative estimation algorithm. The numerical results show the boost in the accuracy performance provided by employing BSL terminals together with CMLs. Moreover, we show that the proposed algorithm is robust to possible errors in the parameters used to model the virga effect, reducing the need for precise real-time estimation of those parameters. Combined with the low installation cost of BSL terminals, the proposed algorithm ensures that the BSLs can act as as gap-fillers of an pre-existing CML network.