Dynamic rainfall monitoring using microwave links

In this work, we propose a sparsity-exploiting dynamic rainfall monitoring methodology using rain-induced attenuation measurements from microwave links. To estimate rainfall field intensity dynamically from a limited number of non-linear measurements, we exploit physical properties of the rainfall such as spatial sparsity and non-negativity along with the dynamics of rainfall intensity. We develop a dynamic state estimation algorithm, where the aforementioned spatial properties are utilized as prior information. To exploit spatial sparsity, we use a basis function to tailor the sparse representation of the rainfall intensity. The basis is selected based on some criteria for sparse reconstruction such as orthonormality and mutual coherence. The tuning parameter that controls the sparsity in the spatial rainfall distribution is dynamically updated at every correction step. The developed methodology is applied to dynamically monitor the rainfall field intensity in an area with a specified spatial resolution using less number of simulated non-linear measurements than pixels. The proposed methodology can be generalized for any dynamic field reconstruction, where the limited number of non-linear measurements are field intensities integrated over a linear path.


Introduction
Spatial rainfall mapping from the measurements of raininduced attenuations collected from microwave links (used by cellular telecommunication networks) is an emerging technology which can serve as an alternative to traditional approaches like rain gauges and weather radar [1]. The practicability of the method is illustrated in [2] by comparing its performance with rain gauges and radar. The motivation behind this methodology is to utilize existing systems such as cellular networks to improve the quality of rainfall estimates using rain gauges and radar, as well as to use it as an independent rainfall measuring unit in areas, where traditional measuring modalities are scarce. The attenuation measurements from microwave links can also be used for monitoring snowfall, fog, and humidity [3]. Seminal works in this domain include tomographic rainfall mapping [4] and a stochastic implementation of the microwave tomographic inversion technique (MTIT) [5]. Recently, it has been observed that signal processing algorithms like a modified *Correspondence: v.roy@tudelft.nl 1 Faculty of Electrical Engineering, Mathematics, and Computer Science, Delft University of Technology, Mekelweg 4, 2628 CD Delft, the Netherlands Full list of author information is available at the end of the article weighted least squares method can be implemented to spatially map the rainfall intensity on a regular grid, using microwave link attenuation measurements [6]. Also, a direct spatial reconstruction from non-linear measurements using a variable grid size is exhibited in [7]. The robustness of a practical application of such techniques is illustrated in [8], where a country-wide (the Netherlands) rainfall mapping is shown to be possible using link attenuation measurements using a data set of 12 days (with a temporal resolution of 15 min). However, in order to achieve some desired spatial resolution of the rainfall field estimate (in terms of number of pixels), the number of microwave links, i.e., the number of attenuation measurements, is always much smaller than the number of pixels in a given service area. In this case, to dynamically monitor the rainfall intensity, physical properties of rainfall like spatial sparsity and non-negativity can be exploited as extra information. In [9], a sparse reconstruction of the rainfall field from a limited number of non-linear measurements is presented. In [10], a sparsity as well as a ridge-penalized, non-negativity constrained, ordinary least squares method is used to estimate the spatial rainfall map from linear path-averaged rainfall intensities, albeit for a single snapshot. Furthermore, incorporating the non-linearity of the measurements as well as a state-space model, a spatio-temporal rainfall monitoring method using an extended Kalman filter (EKF) is described in [11]. Recently, a linear Kalman filter is used for the reconstruction of rainfall maps inspired by object tracking algorithms [12]. However, none of the above dynamic rainfall monitoring methods exploits structural properties of the rainfall field like sparsity or non-negativity.
Commingling the concepts of the aforementioned literature, estimating a spatio-temporally evolving rainfall field can be viewed as a dynamic sparse field estimation problem, where the spatial sparsity of the rainfall field can be tailored by representing it as a sparse signal in a suitable "sparsifying" basis [13]. Such a dynamic estimation of sparse signals, also known as sparsity-aware Kalman filtering, is a well-studied problem in the field of signal processing with quite a number of applications like target tracking and video coding. Next to the spatial sparsity, also, the temporal sparsity can be exploited in the state estimation [14]. Sparsity penalties lead to a faster convergence than a clairvoyant Kalman filter, as illustrated in [14]. Also, a non-negativity constrained sparsity-aware Kalman filter is applied to the target tracking problem in [15]. In [16], the "dynamic filtering" is implemented by introducing an iterative re-weighted 1norm penalty. In that work, a Bayesian hierarchical model is used for the dynamically varying sparse coefficients of the signal. Also, in [17], the convergence of the aforementioned approach has been illustrated, formulating it as a basis pursuit denoising (BPDN) problem. Another notable approach of tracking a sparse signal in an underdetermined measurement scenario is viewing sparsity as a pseudo-measurement and implementing a parallel state and covariance update scheme for this extra measurement [18]. In the Bayesian paradigm, a sparsity-aware state estimation can be formulated as a constrained maximum a posteriori estimator (MAP) [19].
In this work, we assume that the spatial rainfall intensity can be represented as a sparse environmental signal. We assume two scenarios for the spatio-temporal evolution of the rainfall field. In the first case, we assume that the dynamics of the rainfall field are perfectly known. In this case, we use a linear but non-stationary dynamical model for the space-time evolution of the rainfall event, which incorporates physical phenomena like advection, diffusion, and convection [20,21]. In the second case, we assume that the information regarding the dynamics is not perfectly known. In this case, we approximate the spatiotemporal evolution by a simple Gaussian random walk model.
We develop a complete structured framework to dynamically monitor the rainfall intensity exploiting a priori knowledge, i.e., the spatial sparsity and the non-negativity of the rainfall field. The overall dynamic rainfall monitoring setup is pictorially represented in Fig. 1. The proposed setup accepts attenuation measurements, in a given service area at any given snapshot from the operating links, whose geometry and operating frequencies are known. Accumulating these non-linear measurements, the spatial rainfall intensity in the given service area is computed in a centralized approach with a specified resolution. The developed dynamic sparsity-aware rainfall monitoring algorithm has the following salient features: • The used measurement model is non-linear, underdetermined, and also time-varying. We perform a dynamic linearization, followed by a state estimation where sparsity and non-negativity is utilized, in order to achieve a stable solution from the underdetermined measurement setup. • The selection of the tuning parameter controlling the sparsity in the solution can be dynamically updated in every correction step. • The algorithm is also generalized to dynamically select the representation basis that minimizes the mutual coherence between the basis matrix and the measurement matrix at a particular time instance, which represents the geometry of the available link measurements at that time instance.
The rest of the paper is organized as follows. In Section 2, the measurement model, state model, and the spatial covariance structure of the rainfall field are presented. In Section 3, the dynamic rainfall monitoring algorithm is discussed. Issues such as the selection of an appropriate sparsifying basis and the selection of the tuning parameter controlling the sparsity is discussed in Section 4. In Section 5, simulation results are presented. Here, we mention that the true value of the rainfall, i.e., the gauge adjusted radar images and the location of the microwave links are known to us. But we do not have real attenuation measurements from the microwave links. We simulate the measurements using the ground truth, available link locations using the non-linear measurement model mentioned in Section 2 and add additive white Gaussian noise (AWGN) of known variance. Section 5 summarizes this paper and looks at future directions.
Notations: Matrices are in upper case bold while column vectors are in lower case bold. The symbol [X] ij is the (i, j)-th entry of the matrix X and [ x] i is the i-th entry of the vector x. The identity matrix of size N × N is denoted by I N . The transpose operator is denoted by (·) T , x is the estimate of x, defines an entity, and x p =

Signal model
In this section, we first describe the non-linear measurement model which can be used to dynamically estimate the rainfall intensity with a prescribed spatial resolution by the time-varying attenuation measurements. Next, we describe the structure of the spatial covariance matrix of the rainfall field. After that, we describe the dynamic state model, based on the physics behind the movement of a rainstorm both spatially and temporally.

Non-linear measurement model
The geometry of the microwave links deployed by any telecommunication service provider in any area is fixed. These links can be viewed as a fixed network of sensors to monitor rainfall since the received signal level (RSL) measurements related to these links depend on the rainfall. Note that the signal attenuation on a microwave link is not only due to rainfall but also depends on other atmospheric effects like humidity, wet antenna attenuation, and propagation loss [6]. For simplicity, we assume that the attenuation caused by these other effects (except precipitation) can be precomputed, e.g., during "dry periods", and subtracted from the recorded RSL measurements. So, the effective measurements only include the rain-induced attenuation. The conventional empirical relationship between the rain-induced specific attenuation and the path-averaged rainfall rate is given by y s = ar b , where y s is the specific attenuation of the link (dB/km), and r is the path-averaged rainfall rate over the link (mm/h) [22]. If L is the length (km) of the microwave link, then the total rain-induced attenuation over the link is y = ar b L dB. The parameters a and b are related to the drop size distribution (DSD) of rain, the polarization and frequency of the transmitted electromagnetic wave, the length of the link, the ambient temperature, etc. It has been extensively studied and shown in several works that variations of the aforementioned environmental and non-environmental parameters affect the estimate of the path-averaged rainfall rate. A quantitative analysis of DSD related errors in estimating the path-averaged rainfall from direct rain-induced attenuation measurements is illustrated in [23,24]. It can be observed that the attenuation for links operating in frequencies around 35 GHz can be treated as a linear measurement of the path-averaged rainfall rate [23]. A detailed analysis of the effects of the frequency, DSD, link length, and temporal sampling in estimating the path-averaged rainfall rate has been presented in [25,26]. Also, in a wide coverage area, the link (measurement) availability in different hours of the day may significantly vary. All of these aforementioned studies advocate a dynamic tuning of the a and b coefficients in order to better monitor the rainfall from link attenuations.
The non-linear attenuation measurements from the microwave links in any given service area for a fixed time can be used to estimate the spatial rainfall intensity over the same area. Let us consider a uniform discretization of the specified service area A (square) into N pixels where we would like to estimate the rainfall intensity. Here, we make the assumption that the rainfall intensity is constant within any pixel. This assumption is flexible as any resolution can be attained by tailoring the size of the square pixels. Let us assume that there are M links in the given service area. The length of the i-th link can be written as L i = N j=1 l ij , where l ij is the length of the i-th link passing through the j-th pixel, where i = 1, . . . , M. If the i-the link does not pass through the j-th pixel then l ij = 0; otherwise, it is computed by the link and the pixel coordinates. The total attenuation over a link can be modelled as the sum of the attenuations over the link segments [6]. Using this, the attenuation over the i-th link at time t can be expressed as y i,t ≈ N j=1 y ij,t , where y ij,t is the attenuation over the linksegment of length l ij . Using the power-law relationship for the attenuations over the link segments, the measurement model can be constructed in the following way, where y i,t is the attenuation measurement of the i-th link and u j,t is the intensity of the rainfall field in the j-th pixel at time t. The power-law coefficients of the i-th link at time t are given by a i,t and b i,t . The measurement model in (1) is a generalized time-varying non-linear tomographic measurement model. In this work, we consider the fact that all the M links are operated in the same frequency and that the other environmental conditions (e.g., DSD, temperature) are fixed for all t. Based on these assumptions, the aforementioned measurement model can be simplified as The measurement noise incurred at the i-th link measurement at time t is given by e i,t . The measurements are corrupted by errors which are mainly due to quantization but also other sources of noise exist. A more detailed description of the statistical nature of the measurement noise can be found in [6]. For the sake of simplicity, let us assume that e i,t is zero-mean spatio-temporally white Gaussian noise with variance σ 2 e . Further, we assume that e i,t is uncorrelated with u j,t .
Combining all the measurements from the M links at time t, we can construct the following non-linear measurement model at time t: where y t ∈ R M stacks the measurements from the M links at time t, whereas e t ∈ R M does the same for the noise. The vector u t ∈ R N gathers the rainfall intensities for all of the N pixels at time t, i.e., it is the parameter to be estimated dynamically. The non-linear mapping between the rainfall intensities and the attenuation measurements is given by : R N → R M and it is assumed to be perfectly known. The elements of u t are given by [ where u t (x) represents the continuous rainfall field at any arbitrary location x ∈ R 2 and x j =[ x j , y j ] T is the centroid of the j-th pixel of the service area. The measurement noise components associated with the M link measurements are characterized by

Spatial variability of u t
At any snapshot t, the spatial rainfall intensity u t (x j ), for j = 1, . . . , N can be viewed as a wide-sense stationary (WSS) random process. In spatial statistics, . . , N in the service area) and if the spatial covariance between any two points is dependent only on the distance between them (i.e., isotropic) [27]. The parameter μ t is the mean/trend of the rainfall field. A variogram model can be used to represent the spatial variations. The variogram (i.e., 2ξ(h)) or semivariogram (i.e., ξ(h)), as a function of h [27]. Generally, several variogram models are used as it is computationally hard to calculate the spatial dependency for every lag distance h. Some statistical functions like a Gaussian, exponential, or empirically fitted models like spherical functions are often used as variogram models [28]. From the analysis of [29], the spherical variogram model is seen to be an appropriate model to describe the spatial variability of rainfall. This is given as, The parameters that characterize a variogram model are the sill N 0 + S 0 of the variogram (ξ(h) for h → ∞) with S 0 as the partial sill, the nugget N 0 (non-zero value of ξ(h) for h → 0), and the range d (value of h for which the variogram reaches the sill). The advantage of the spherical variogram model is that the parameters S 0 , N 0 , and d can be well estimated in hourly scales for a specific day of the year [29]. Now, the spatial covariance function Using the second-order stationarity of the random process u t (x j ), the semivariogram can be related to the spatial covariance function [27]. Now, the elements of the spatial covariance matrix u can be computed as The covariance matrix, constructed in this way, is symmetric and positive definite.

State model 2.3.1 A Kernel-based state model
One standard approach of modelling the spatio-temporal evolution of any environmental field is based on the integro-difference equation (IDE) [28]. Following this approach, the dynamics of the rainfall field for any specific temporal sampling interval δ t can be modelled as the following discrete time IDE Here, g(x, x ; θ ) is the space-time interaction function parameterized by θ , which can be deterministic or random and dependent on the temporal sampling interval δ t . The quantity q t (x) is the process noise which is generally modelled as independent in time but correlated in space.
The space-time interaction function g(·) can be modelled as a parameterized Gaussian dispersal kernel which captures the underlying physical processes behind the spatio-temporal evolution of rainfall, i.e., diffusion, advection, and convection [20,21]. In this case, the spacetime interaction function is given as , i.e., a Gaussian kernel. The translation parameter of the kernel, i.e., w t ∈ R 2 models the time-varying advective displacement, i.e., the spatial drift of the rain storm, and the dilation parameter of the kernel, i.e., D ∈ S 2 ++ models the diffusion. Note that w t can also vary with space but we assume that it is averaged over the entire area and fixed. The diffusion coefficient D can be used to model isotropic as well as anisotropic diffusion. The amount and the directions of the spatial anisotropy can be introduced by D. The parameter D can also vary with time but this is not considered here. The scalar scaling parameter α ∈ R ++ is used to control the stability (i.e., to avoid the explosive growth) of the dynamic process.
Here, the entire service area is uniformly discritized into N pixels. We assume a state transition matrix H t ∈ R N×N whose elements are modelled by the aforementioned simple 2D Gaussian kernel. After proper vectorization of the field intensities and state noise for N pixels, we obtain where the elements of the state transition matrix H t are given by and q t is the spatially colored yet temporally white Gaussian state noise vector. The quantity w t is the advective displacement during the temporal sampling interval δ t , which can be represented more precisely as where v t is the advection velocity. Note that the aforementioned model is non-stationary when the advection vector w t changes with time, which happens in many real scenarios [20]. If there is no advection, i.e., w t = 0 and D = I, the model is stationary and isotropic. We assume that the dynamic model, i.e., the state transition matrix H t is perfectly known through the parameters w t , D, and α which are considered to be deterministic and known. Without loss of generality, we follow the assumptions of [20] and [30] that the distribution of q t is given by But this assumption is not true in practical scenarios because the rainfall process cannot be negative. In the simulation section, after generating u t using the sate model of (6), we set the negative elements of u t to 0. This is a modelling approximation.
One notable advantage of the model in [20] is the linear relation of the rainfall intensities in one snapshot with the ones in the previous snapshot.

Gaussian random walk model
In the last section, we assume that the parameters of the state model are perfectly known. But in many practical scenarios for a large N, it can be computationally intractable to estimate the N 2 elements of the state transition matrix H t using the available data. In this case, without any prior knowledge regarding the parameterization of H t , one way to approximate the dynamics is by assuming that the process follows a Gaussian random walk model [31]. In this case, we assume that H t = H = I and the process model is given by The benefit of a Gaussian random walk model is that it has very few model parameters rather than a parameterized process model as mentioned in Section 2.3.1.
Note that the parameterized state model of (6) can be viewed as a random walk model by incorporating negligible diffusion, i.e., D = I, where 1 and no advection, i.e., w t = 0. In this case, we have H t ≈ I assuming α = 1.

Structure of the state error covariance matrix
It is assumed that the state error, i.e., q t , is a spatially colored but temporally white Gaussian process. Assuming spatial isotropy and stationarity of the state error q t , the elements of the covariance matrix Q t = Q can be represented using the Matern covariance function as where (·) is the Gamma function, K p (·) is the modified Bessel function of the second kind, and γ is a positive shaping parameter [28]. With p → ∞ and p = 1/2, (8) becomes the squared exponential and the exponential covariance functions, i.e., , and

Dynamic rainfall mapping
We dynamically estimate the rainfall intensities at the N pixels, i.e., u t at t = 1, . . . , T snapshots from the attenuation measurements y t at t = 1, . . . , T. The measurement and state models can be represented in the following forms A standard practice to estimate the rainfall intensity u t at every time t = 1, . . . , T from the measurement and state equations of (9) and (10) is the non-linear semblance of the standard Kalman filter, i.e., the extended EKF [32]. Note that we have non-linearity only in the measurements.
As one of the criteria for the optimal behavior of the Kalman fiter, we assume that the measurement and the state noise statistics are completely known. The measurement and the state noises are characterized by e t ∼ N (0 M , R), and q t ∼ N (0 N , Q), respectively. The dimension of the measurement noise covariance matrix depends on the number of the available measurements at time t. As the state model is a linear function of u t , the standard Kalman fiter prediction steps are given bŷ where the prediction of u t from the last t − 1 observations is given byû t|t−1 with the error covariance matrix The prediction based on the state model is corrected by the measurements. But here, we have a non-linear measurement model. To linearize that model, let us introduce the M × N Jacobian matrix computed at u t =û t|t−1 as The elements of the Jacobian matrix are given by . . , M, and j = 1, . . . , N. A first order Taylor series expansion of the non-linear measurement function aroundû t|t−1 is then given as Substituting this in (9), we obtain the following linearized measurement equation: Note that here, we have less observations than unknowns, i.e., the number of links (M) is much smaller than the number of pixels (N), i.e., the dimension of u t . Hence, in the correction step, to utilize the measurements along with the state model, we need to solve the underdetermined system (13) in order to updateû t|t−1 leading toû t|t . After the dynamic linearization, the state estimates can be obtained using a standard Kalman filter. In this case, both the expressions for the state estimateû t|t and its state error covariance M t|t can be obtained in closed form [32].

Limitations of a standard EKF
The estimation of u t from only M measurements using an ordinary EKF has the following uncertainties.
• First of all, the quality of the estimate strongly depends on the degree of non-linearity and the accuracy of the linearization [32]. Also, for a highly underdetermined (M N) and unpredictable measurement matrix (many rows of J t can be zero for anyû t|t−1 ), the solution can be highly inaccurate and dependent mainly on the predictions using the state model and the initialization. • In the above case, if the available information regarding the dynamics are incomplete or imperfectly known, then the prediction using the state model will be inaccurate. In this case, an ordinary EKF may produce unrealistic estimates in the presence of high measurement noise. • Also, there is no guarantee that an ordinary EKF will always produce non-negative estimate ofû t|t . For instance, let us assume that an element of the predicted value, i.e., [û t|t−1 ] j (predicted using (11)) is less than 0 at any t. In that case, if l ij = 0, we may have an imaginary is a fractional quantity. As mentioned in [22], the standard values for b mainly lie in the interval of 0 < b < 2.
In these circumstances, any further prior information about u t (beyond the dynamics) is desirable to achieve a stable and more accurate solution.

Available prior knowledge regarding rainfall field
Prior information about u t can be acquired from the physical properties of rainfall such as sparsity and nonnegativity. In a given area, the rainfall intensity itself can be assumed to be a sparsely distributed environmental field over the entire service area [9,33]. But sparsity can also be introduced by representing u t in an orthonormal basis t , which can in principle be time-varying. When rainfall itself is sparse, we simply have t = I. Denoting u t = t z t , sparsity is measured by the number of non-zero entries in z t , i.e., z t 0 .
As the rainfall intensity cannot be negative, another prior knowledge about u t is the non-negativity of the rainfall field. For N pixels, this is represented as u t ≥ 0 N .
Comment: Here, we mention that the prior information regarding sparsity and non-negativity along with the measurements can be efficiently utilized to monitor the rainfall over multiple snapshots. For this, we do not need any information regarding the dynamics. This can be implemented for both linear [10] as well as non-linear [9] measurement models. However, one limitation of this dynamics-agnostic method is that the rainfall events should occur in areas where microwave links are present for accurate estimation. Otherwise, the effect of the measurement noise can be dominant. In this case, we need other spatial/temporal information (e.g., covariance structure, dynamics) to interpolate the rainfall field over the entire service area.
In the next section, we illustrate iterative approaches to dynamically estimate the state of u t for t = 1, . . . , T, exploiting both sparsity and non-negativity.

Estimation of u t
A simple Kalman estimation step without the sparsity and the non-negativity constraint can be formulated as the following weighted least squares optimization problem [34]: This estimation step is not aware of sparsity or nonnegativity. The sparsity information can be incorporated in the optimization problem of (14), by adding an 1penalty that enforces sparsity. Note that here, we use the 1 norm as a convex relaxation of the non-convex 0 norm. Using the sparse representation of u t , i.e., z t , the optimization problem of (14) can be formulated as a sparsity and non-negativity constrained optimization problem. This can be given aŝ where λ t is the tuning parameter that controls sparsity. The standard error covariance update ofû t|t for the estimation step (14) is given by [32]. This expression of M t|t can be used to update the covariance of the estimate of (15) but is an approximation as it is not aware of the sparsity and the non-negativity constraint. If we do not consider to propagate the second-order statistics of the estimate, like in the traditional Kalman filter, the state noise minimization term in (15) can also be regularized by Q instead of M t|t−1 . This can be viewed as a weighted least squares problem to estimate u t using the measurement (13) and the state Eq. (10) constrained by sparsity and non-negativity. In this case, the simple iterative state estimates are given bŷ Note that this is a suboptimal approach to dynamically estimate the states u t avoiding the computation of M t|t . As mentioned in [14], different penalties (like 2 or 1 ) can be applied to the state error minimization term in (18) depending on the nature of the sparse state (u t ) and/or the state noise (q t ).

Formulation of a constrained MAP estimator for u t
Using the representation of u t in the t = domain, the measurement and state equation of (13) and (10) can be written asỹ N (0, Q). From the above measurement and state equations, we can derive the conditional probability density functions p(ỹ t |z t ) ∼ N (J t z t , R) and p(z t |z t−1 ) ∼ N (H t z t−1 , T Q ). Using Bayes' rule, the posterior pdf p(z t |ỹ t ) can be given as p(z t |ỹ t ) ∝ p(ỹ t |z t )p(z t |z t−1 ). So, a MAP estimator for z t can be formulated as, where z t−1 is computed from the previous time step. However, there is no guarantee that the estimator in (22) will produce a sparse estimate of z t . On the other hand, the representation of u t in the domain is targeted to exploit sparsity. So, the estimator of (22) can be formulated as a constrained MAP estimator by adding the sparsity and non-negativity constraint in the optimization problem of (22). After substituting the pdfs, following the same approach as used in (15), the sparsity and non-negativity constrained MAP estimator can be given aŝ whereλ t controls the sparsity in the estimateẑ t . From the solution of (23), the state estimate is given asû t|t = ẑ t . It is seen that resorting to a Bayesian paradigm, the developed constrained MAP estimator of (23) has a structure similar to the optimization problem of (18).

Selection of the basis ( t )
If the spatial rainfall distribution is physically sparse, then we simply solve the optimization problem of (18) for t = I. But in the absence of physical sparsity, which is a more general case in many practical scenarios, we use a basis t . Revisiting the celebrated theory of compressive sampling, we know that some of the properties of both t and J t , (or t = J t t ) like mutual coherence and restricted isometry property (RIP) are important in the framework of sparse reconstruction [13]. Let us denote the quantity μ( t ) as the mutual coherence of the matrix t , which is the maximum absolute inner product of different columns of t [35]. Without the state error minimization term and the non-negativity constraint in (18), the problem is a simple BPDN problem. As derived in [36], if a suitable sparse representation of z t is possible, which is given by z t 0 < 1 2 (1 + 1 μ( t ) ), then a "stable" solution with the standard BPDN algorithm can be obtained with a bound on the estimation error. Along the same lines, the cost function of (18) can be viewed as a BPDN problem by augmenting the measurement and the state noise minimization terms into a single least squares term [17]. In [17], the convergence guarantees of the aforementioned BPDN problem are also derived based on some assumptions on the dynamics and the measurement matrix (here J t ).
In our application, the design of the measurement matrix J t , in every snapshot, is dictated by the link locations and the predicted state estimates (û t|t−1 ). So, to maximally exploit the sparsity information of the rainfall field, we focus mainly on a suitable sparse representation of the state u t . Standard orthonormal bases such as a discrete cosine transform (DCT) basis C or wavelet basis W are quite popular in sparse signal representation for communications as well as image processing. Also, a Gaussian basis function can be used to sparsely represent environmental signals [37]. However, an orthonormal basis can also be constructed using the spatial covariance matrix of the rainfall field. An orthonormal basis can be constructed by the spatial covariance matrix u described in Section 2.2, by simply choosing t = U, where u = U U T is the eigenvalue decomposition of u with U T U = I and a diagonal matrix. In this case, z t = U T u t is similar to applying a Karhunen-Loeve transform (KLT), which is also advocated as a sparse representation technique [38].
We choose a basis t that has a minimum mutual coherence with the measurement matrix J t . Mutual coherence can be measured for the overall dictionary t = J t t . In this case, μ( t ) can be quantified as the maximum magnitude off-diagonal element of D t =˘ T t˘ t , wherȇ t is obtained by normalizing the columns of t . In this case, the mutual coherence can be defined as μ( 39]. So, given a set U of U sparsifying basis matrices, the minimal coherence basis matrix at time t can be selected by solving the following optimization problem, Note that to optimally select the basis, we need to solve this optimization problem on every snapshot as the matrix J t is recomputed at every time step t. In our simulations, we specifically use U = {C, U}.

Selection of the tuning parameter (λ t )
The sparsity regulating parameter λ t in the optimization problem of (18) can be adapted dynamically. It can also be kept fixed for multiple snapshots for monitoring a short period of rainfall, within which the sparsity pattern can be assumed to be fixed. An upper bound on λ t is given by λ t = λ max t , which gives the sparsest solution, i.e.,ẑ t = 0 N orû t = tẑt = 0 N . Note that the cost function of (18) is non-differentiable but convex for λ t > 0. So, following the methodologies of [15] and [40], we use a subdifferential based approach to compute λ max t . The subgradient of the non-differentiable cost of (18) with respect to z t can be written as where∇ z t is the subgradient operator towards z t . Using the first-order optimality condition, we have where j = 1, . . . , N. Now, we consider the case z t = 0. Substituting this in the above equation, the optimal value of λ max t can be selected as λ max Traditional approaches to select the tuning parameter for an 1 -penalized regression problem are crossvalidation and generalized cross-validation (GCV) [41]. Recent methods suggest information theoretic approaches like Mallow's Cp type criterion [42], Akaike information criterion (AIC) [43], and Bayesian information criterion (BIC) [44] to find an optimal λ t . In all of these approaches, the optimal tuning parameter is selected that minimizes a cost function which depends upon the estimate of u t using a set of {λ k t } K k=1 , where the length of the search grid for λ t is K. In this case, we need to solve the optimization problem of (18) K times in every iteration and select the optimal λ k t that minimizes any of these aforementioned model selection criteria. After that, (18) needs to be solved again with the selected λ k t , in order to estimate u t . This seems to be computationally unrealistic for an online application of the dynamic rainfall monitoring over a large service area (large N). To circumvent this problem, the K optimization problems can be solved once and the selected λ t can be used for multiple snapshots for a short term monitoring application.
It is clear that increasing the value of λ t , i.e., the term z t 1 becomes smaller and vice versa. So, if an approximation of z t , i.e., z approx t , is available, it can be related to the tuning parameter by λ t ∝ 1/ z approx t 1 . Following this, a coarse but relatively fast approach to dynamically tune λ t could be selecting λ t = ν( T tû t|t−1 1 ) −1 , whereû t|t−1 can be regarded as an approximation ofû t|t and ν > 0 is a proportionality constant. In our simulations, we use ν = 1.
For the sake of completeness, we summarize the steps of the two proposed dynamic rainfall monitoring algorithms. In algorithm 1, we follow the standard steps of dynamic state estimation, but we do not update the second-order statistics of the estimate. In algorithm 2, we use the approximate approach, where we use the standard Kalman covariance update (unaware of sparsity and nonnegativity).
The performance of both these algorithms strongly depends upon the initializationû 0|0 . One should avoid initializations like an all zero vector or anû 0|0 that consists of negative elements. If we consider the initialization u 0|0 = 0 N , it will produce J t = 0 M×N , asû 1|0 = H tû0|0 = 0 N . It is mentioned in [22], that the standard values for b mainly lie in the interval of 0 < b < 2. It should also be noted that, for b < 1 (for frequencies in the range 1 − 3 GHz, or frequencies above 40 GHz [22, Table II This problem can be circumvented by replacing the 0 rainfall scenario by a very small value (close to zero) like in the order of 10 −4 mm denoting a no rainfall event and the non-negativity constraint can be replaced by u t ≥ 10 −4 1 N in the optimization problems. However, in our simulations, we use b > 1.

Simulation results
In this section, we present some simulation results to test the developed methodologies to dynamically monitor the rainfall in a given area. Here, we perform numerical experiments for three scenarios. In the first case, we assume that the dynamics/state model, i.e., H t is perfectly known through the parameters α, D, and w t . In the second case, we consider that the dynamics are not perfectly known and we assume that the state model is a Gaussian random walk. In the third case, we consider the scenario where we do not have any information regarding the state model/dynamics. The simulations for these three scenarios are presented below.

Ground truth
We consider two sets of ground truth/true values of the rainfall vector {u t } 8 t=1 . In the first case, we assume that the state model is perfectly known. In the next case, we assume that the state model is not perfectly known. given a, b, l ij , (i = 1, . . . , M; j = 1, . . . , N), y t , R t , H t , Q

Ground truth with known dynamics
The ground truth is used from a practical rainfall event in an area of 25 × 25 km 2 in Amsterdam, the Netherlands 1 . We take one spatial map of 15 min gauge adjusted radar rainfall depth (mm) of the same area of the day June 11, 2011, which is shown as the first state, i.e., u 1 in Fig. 2.
We assume that the state transition matrix, i.e., H t and the process noise covariance matrix, i.e., Q t = Q is perfectly known in this case. The parameters of the state transition matrix H t = H are given as w t = w =[ 1, 0] T , α = 0.33, and D = I 2 (isotropic diffusion) for all t = 2, . . . , T snapshots, where T = 8. We assume the temporal sampling interval, i.e., δ t = 15 min. The parameter w represents a constant advective displacement in 15 min.
The covariance matrix of the state noise, i.e., Q is assumed to have an exponential structure given as , with σ 2 s = 10 −3 and γ = 3.33. The state noise vector q t at every snapshot is generated from the distribution q t ∼ N (0 N , Q). After generating the states of u t using the state model of (6), we set the negative elements of u t to 0. This is a modelling approximation adopted to avoid the generation of the negative rainfall values for very low rainfall intensities. The total number of pixels is given as N = 25 × 25 = 625, each of size 1 km 2 . Using this, we generate the states u 2 . . . u 8 using the state model mentioned in (6). Based on these parameters, the space-time evolution of rainfall over t = 1, . . . , 8 snapshots (each of 15 min, i.e., in total 120 min) is shown in Fig. 2. The unit of the rainfall field is millimeter (mm).

Ground truth with approximate dynamics
In this section, we consider eight consecutive snapshots of 15-min radar rainfall depths of the same day and area as mentioned in the previous section. The eight snapshots of the true gauge adjusted radar rainfall depths are shown in Fig. 3. Here, we assume that the state model is a Gaussian random walk, i.e., H t = I. The process noise covariance matrix Q is assumed to be the same as before.

Measurements
In this work, we simulate the measurements using the locations of the microwave links from a network of 151 microwave links. The total number of measurements for all t = 1, . . . , 8 is M = 151. The locations of the microwave links in the service area along with the (1 km 2 ) pixels, where we would like to estimate the rainfall intensities, are shown in Fig. 4.
We would like to mention that most of the microwave links in the Netherlands (specially in the urban areas) are operated at 38 GHz. In that case, b ≈ 1, i.e., the measurement model becomes linear. To check the estimation performance of the developed algorithms in a non-linear measurement framework, we intentionally choose b = 1. In this case, we assume the rain temperature to be 20 • , which corresponds to [22, Table II]. We select the operating frequency to be 15 GHz, with a = 3.28 × 10 −2 and b = 1.173. Using these, we simulate the measurements at eight snapshots, i.e., {y t } 8 t=1 , using the non-linear measurement model of (2), where u j,t 's are the true values for j = 1, . . . , 625, and t = 1, . . . 8 as mentioned in the previous section.
Here, we generate two different sets of measurements. The first set of measurements is for the seven snapshots, i.e., {y t } 8 t=2 . These measurements are computed using the ground truth where the dynamics, i.e., H t = H for t = 2, . . . , 8 are perfectly known (Fig. 2). The second set of measurements are computed using the exact radar rainfall maps (Fig. 3) for eight snapshots, i.e., {y t } 8 t=1 whose dynamics are unknown.
The parameters L i and l ij are known from the geometry of the links as shown in Fig. 4. It is assumed that the M measurements are collected in every 15-min interval which are corrupted by additive white Gaussian noise characterized by e t ∼ N (0 M , σ 2 e I M ).

Dynamic rainfall monitoring
The noisy sets of measurements are used to estimate the rainfall depths at N = 625 pixels over T = 8 snapshots. In this section, we perform simulations for three Fig. 3 Spatio-temporal evolution of the rainfall (mm) (unknown dynamics) different scenarios which are perfectly known dynamics, a Gaussian random walk dynamics, and no dynamics.

Perfectly known dynamics
The measurements {y t } 8 t=2 are computed using the true values shown in Fig. 2, i.e., {u t } 8 t=2 . Here, we use σ 2 e = 10 −3 . These measurements are used to estimate the states {û t } 8 t=2 . The parameters of the spherical variogram model are computed for the particular day of the year, i.e., June 11, 2011 [29], to compute u . The sill (S 0 ) and the range (d) parameters are computed using (11) of [29], whose parameters are taken from [29, Table 5]. The 15-min time interval is rescaled in hourly scales, i.e., 0.25 h. We assume that the nugget is N 0 = 0. The value of the range (d) is 17.4675 km, and the sill (S 0 ) is 5.3328 mm 2 .
Based on the predictions and the available link locations, it is seen that the DCT matrix exhibits minimal coherence with J t , in every iteration. We initializê u 1|1 =μ1 N , whereμ is computed by empirically averaging the ground truth of the first state u 1 over N pixels for both algorithms. In a real application, an appropriate initialization can be computed using the trend of the rainfall field, which is generally available from the climatological information of the area. For algorithm 2, we initialize M 1|1 = I N . We use the software CVX [45] (parser CVX, solver SeDuMi [46]) to solve the convex optimization problems (i.e., (18) for algorithm 1 and (15) for algorithm 2).
In Fig. 5, we show the reconstructed spatial rainfall map for the statesû 2 andû 8 using the algorithm 1. The same estimates are shown in Fig. 6 using the algorithm 2.
We plot the pixel-wise comparisons of the estimates with the true values for all the seven snapshots, i.e., a total of 625 × 7 = 4375 pixels for the algorithms 1 and 2 in Fig. 9a and Fig. 9b respectively. The dark black lines in Fig. 9a and Fig. 9b, represent the y = x line. It is observed that algorithm 2, i.e., using the standard Kalman covariance update exhibits better estimation performance than algorithm 1, where the second order statistics are not updated.

Gaussian random walk
In this section, the state model is considered to be a Gaussian random walk, i.e., H t = I for t = 1, . . . , 8. The process noise statistics are considered to be the same as before. The measurements for eight snapshots, i.e., {y t } 8 t=1 are generated using the true radar rainfall depths shown in Fig. 3, using the measurement model of (2) with the same a, b coefficients as the previous case. The measurements are different from the previous known dynamics case as the true values of {u t } 8 t=1 , i.e., the true radar rainfall depths (as shown in Fig. 3) are different. In this case, the measurement noise variance is reduced to σ 2 e = 10 −5 . Due to better estimation performance (as seen in Section 5.3.1), we select algorithm 2 to estimate the states {û t } 8 t=1 using the measurements generated by the true radar rainfall depths. Now, as the predictions using the sate model are not accurate in this case, we do not perform the tuning of λ t based on the predictions. However, the tuning of λ t , in this case can be performed using the standard methods mentioned in Section 4.2. In the current setup, to exploit the sparsity prior on every snapshot, we fix λ t = λ = 2 for the sake of simplicity. The initializations are given bŷ u 0|0 =μ1 N and M 0|0 = 1 N .
In Fig. 7, we show the estimated spatial rainfall maps of the statesû 1 andû 8 assuming that the state model is a Gaussian random walk. In Fig. 10a, we compare the estimation performance of the estimates of the 625 pixels over eight snapshots, i.e., a total of 625 × 8 = 5000 pixels with the true gauge adjusted radar rainfall depths.

No dynamics
In the last case, we assume that we do not use any prior information regarding the dynamics. On every snapshot, we estimateû t using the measurements and exploiting the prior information regarding the sparsity and the nonnegativity. In this case, for the sake of simplicity, we assume that the measurement model is linear, i.e., the case when the links are operated around 38 GHz. Here, we use a, b = 1.
The linear measurement model is given by On every snapshot, we solve the sparsity-aware nonnegativity constrained optimization problem given aŝ Fig. 7 Estimate of the spatial rainfall (mm) map with unknown dynamics (Fig. 3). (a)û 1 , (b)û 8 ; (algorithm 2 assuming Gaussian random walk dynamics) However, this can be easily extended to a non-linear measurement model by adopting an iterative linearization with respect to a suitable initial guess. Like in the previous case, here, we also fix λ = 2 and the used basis is DCT matrix on every t. However, an upper bound on λ can be easily computed in this case by using the same methodology discussed in Section 4.2.
In Fig. 8, we show the estimated statesû 1 andû 8 by solving the optimization problem of Eq. (26). The measurement noise variance is set as σ 2 e = 10 −5 . In Fig. 10b, we compare the estimation performance of the estimates of total 625 × 8 = 5000 pixels with the true gauge adjusted radar rainfall depths.

Analysis of the results
The following inferences can be drawn from the aforementioned simulation studies.
• The estimation performances in the first two cases are highly dependent on the accuracy of the state model, availability of the measurements in any region, and the initialization of the algorithm. In Fig. 6b, the estimate of the stateû 8 is much better than the same estimate in Fig. 7b, as the dynamics are perfectly known in the first case. Also, the estimation performance is improved with time as the state error is minimized with temporal iterations (Fig. 6b). • As seen in Fig. 4, there are many regions without any microwave links/measurements but where a rainfall field is present. In these regions, the estimates are mainly dependent on the predictions. On the other hand, if there is no rainfall over any link, the rainfall can be overestimated or underestimated in those regions, due to the effects of the measurement noise. This effect severely impairs the estimation performance in the case when we do not have an accurate prediction (or no prediction). • There is always a trade-off between the estimation performance and the "availability" of the measurements and/or the "accuracy" of the predictions. • The reasons behind the scatter plots being not very symmetric are due to the biased estimates in the measurement-void regions and the rainfall-void links. • In the case of an imperfectly known or unknown state model, the prior knowledge of sparsity and nonnegativity help to achieve a stable solution (Figs. 7, 8,  10a, b).

Performance metrics
To compare the estimation performances of the developed methods, we use some performance metrics, which are described in the following part of this section. The performance of a rainfall monitoring method can be quantified by root mean square error (rmse in mm), the mean bias (mb in mm), and the correlation coefficient (ρ) [7]. We quantify the overall estimation performances of all the above scenarios for N pixels over T snapshots using the aforementioned metrics. If the true value and the estimate of the rainfall field at any t are given by u t andû t , respectively, then the performance metrics can be defined in the following ways [ u t ] j are the sample means of the estimated and the true values of rainfall for N pixels over T snapshots. The performance metrics are computed for the estimates using algorithm 2 for the scenarios of perfectly known dynamics and Gaussian random walk dynamics. For both of these scenarios, we also estimate the rainfall depths using a simple EKF without any sparsity and nonnegativity constraint. To avoid the negative estimates produced by the EKF, we set the negative estimates to 0. While computing the performance metrics, we fix the process and the measurement noise variances to σ 2 s = 10 −4 and σ 2 e = 10 −3 , respectively, in all the cases. When the dynamics are perfectly known, then the performance metrics are computed for the estimates of 625 pixels for seven snapshots and averaged over 20 different measurement noise realizations. In Table 1, we present the performance metrics computed for algorithm 2 and a simple EKF (with the thresholding) for the perfectly known dynamics case.
In Table 2, we present the aforementioned performance metrics computed for the estimates of 625 pixels for eight snapshots using algorithm 2 (with fixed λ t = λ = 2) and an EKF (with the thresholding), where the state model is assumed to be a Gaussian random walk. Here, we also average the performance metrics for 20 different measurement noise realizations. In all of these realizations, the measurements are generated using the true radar rainfall depths as shown in Fig. 3.

Analysis of the performance metrics
• For a perfectly known state model, extra information like sparsity and non-negativity does not play any significant role in terms of estimation accuracy. This is clear from Table 1 where it is shown that the performance improvement over a simple EKF (with setting the negative estimates to 0) is negligible. • When the information regarding the state model is unknown and approximated as a Gaussian random walk model, then the performance of a simple EKF is very poor. In this case, the sparsity and non-negativity information along with the measurements improve the estimation performance (Table 2). • The last mentioned observation is quite useful in practical cases, where the availability of an accurate state model is scarce.
The computation times for both algorithm 1 and algorithm 2 including the basis selection part is less than a minute for N = 625 pixels using the aforementioned off the shelf solvers. The computation time is increased if N is higher than 625. Algorithm 1 is computationally simpler than algorithm 2 because there is no covariance update state. But the price we pay is in terms of the estimation performance. However, the speed of the developed algorithms can be increased by using a projected subgradient method [47].

Conclusions
We have developed a generalized dynamic rainfall monitoring algorithm from limited non-linear attenuation measurements by utilizing the spatial sparsity and nonnegativity of the rainfall field. We have formulated the dynamic rainfall monitoring algorithm as a constrained convex optimization problem. The performance of the developed algorithm is compared with the standard approaches like an EKF for the scenarios, where we have both perfect knowledge about the state model and an approximate state model. Numerical experiments show that the developed approach outperforms a simple EKF in scenarios, where the state model is not perfectly known. The proposed methodology can be equivalently implemented for dynamic field tracking in tomographic applications like MRI and microwave tomography, where we have path-integrated Gaussian measurements. However, tackling more complicated dynamics of the rain (possibly non-linear and highly time-varying), and non-Gaussian measurement noise could be possible future extensions of this work. In that case, both the state and the measurement models are non-linear. This triggers one to use an unscented Kalman filter (UKF), particle filtering based algorithms, or other heuristic approaches. Estimation of the underlying dynamics of rainfall from the available ground truth and using it for real-time dynamic monitoring is also a part of the future research. A realtime selection of the most informative attenuation measurements from the available links could be interesting in order to reduce the processing time and computational complexity.

Competing interests
The authors declare that they have no competing interests.