{'fraction_non_alphanumeric': 0.035247518133502694, 'fraction_numerical': 0.012183708259518, 'mean_word_length': 4.7335185894978915, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational Downscaling, Fusion and Assimilation of Hydrometeorological States via Regularized Estimation

Improved estimation of hydrometeorological states from down-sampled observations and background model forecasts in a noisy environment, has been a subject of growing research in the past decades. Here, we introduce a unified framework that ties together the problems of downscaling, data fusion and data assimilation as ill-posed inverse problems. This framework seeks solutions beyond the classic least squares estimation paradigms by imposing proper regularization, which are constraints consistent with the degree of smoothness and probabilistic structure of the underlying state. We review relevant regularization methods in derivative space and extend classic formulations of the aforementioned problems with particular emphasis on hydrologic and atmospheric applications. Informed by the statistical characteristics of the state variable of interest, the central results of the paper suggest that proper regularization can lead to a more accurate and stable recovery of the true state and hence more skillful forecasts. In particular, using the Tikhonov and Huber regularization in the derivative space, the promise of the proposed framework is demonstrated in static downscaling and fusion of synthetic multi-sensor precipitation data, while a data assimilation numerical experiment is presented using the heat equation in a variational setting.


Introduction
In parallel to the growing technologies for earth remote sensing, we have witnessed an increasing interest to improve the accuracy of observations and integrate them with physical models for more accurate estimates and forecasts of environmental states. Remote sensing observations are typically noisy and coarse-scale 1 arXiv:1210.0854v1 [physics.data-an] 2 Oct 2012 representations of a true state variable of interest, lacking sufficient accuracy and resolution for finescale environmental modeling. In addition, environmental predictions are not perfect as models often suffer either from inadequate characterization of the underlying physics or their initial conditions. Given these limitations, several classes of estimation problems present themselves as continuous challenges for the atmospheric, hydrologic and oceanic science communities. These include: (1) Downscaling which refers to the class of problems to enhance the resolution of a measured field or produce a fine-scale representation of a coarse-scale modeled field; (2) Data fusion, to produce an improved estimate from a suite of noisy observations at different scales; and (3) Data assimilation which deals with estimating initial conditions in a predictive model consistent with the available noisy observations and the model dynamics.
In this paper, we revisit the problems of downscaling (DS), data fusion (DF), and data assimilation (DA) focusing on a common thread between them as inverse estimation problems. Proper regularizations and solution methods are proposed to efficiently handle large-scale data sets while preserving key statistical and geometrical properties of the underlying fields of interest. Here, we only examine hydrometeorological problems with particular emphasis on land-surface applications.
In land surface hydrologic studies, downscaling of precipitation and soil moisture signals has received considerable attention, using a relatively wide range of methodologies. DS methods in hydrometeorology and climate studies generally fall into three main categories namely, dynamic downscaling, statistical downscaling, and variational downscaling. Dynamic downscaling often uses a regional physical model to reproduce fine-scale details of the state of interest consistent with the large-scale scale observations or outputs of a global circulation model [e.g., Reichle et al., 2001a;Castro et al., 2005;Zupanski et al., 2010]. Statistical downscaling methods encompass a large group of methods that typically use empirical multiscale statistical relationships, parameterized by observations or other environmental predictors, to reproduce realizations of fine-scale fields. Precipitation and soil moisture statistical downscaling has been mainly approached via spectral and (multi)fractal interpolation methods, capitalizing on the presence of a power law spectrum and a statistical self-similarity/self-affinity in precipitation and soil moisture fields [among others, Lovejoy and Mandelbrot, 1985;Lovejoy and Schertzer , 1990;Gupta and Waymire, 1993;Kumar and Foufoula-Georgiou, 1993a,b;Perica and Foufoula-Georgiou, 1996;Veneziano et al., 1996;Wilby et al., 1998b,a;Deidda, 2000;Kim and Barros, 2002;Rebora et al., 2005;Badas et al., 2006;Merlin et al., 2006]. In variational approaches, a direct cost function is defined whose optimal point is the desired fine-scale field which can be obtained via using an optimization routine. Recently along this direction, Ebtehaj et al. [2012] cast the rainfall DS problem as an inverse problem using sparse regularization to address the intrinsic rainfall singularities and non-Gaussian statistics. This variational approach belongs to the class of methodologies presented and extended in this paper.
The DF problem has also been a subject of continuous interest in the precipitation science community mainly due to the availability of rainfall measurements from multiple spaceborne (e.g., TRMM and GOES satellites) and ground-based sensors (e.g., the NEXRAD network). The accuracy and space-time coverage of remotely sensed rainfall are typically conjugate variables. In other words, more accurate observations are often available with lower space-time coverage and vice verse. For instance, low-orbit microwave sensors provide more accurate with less space-time coverage compared to the high-orbit geo-stationary infrared (GOES-IR) sensors. Moreover, there are often, multiple instruments on a single satellite platform (e.g., precipitation radar and microwave imager on TRMM), each of which measures rainfall with different footprints and resolutions. A wide range of methodologies including weighted averaging, regression, filtering, and neural networks has been applied to combine microwave and Geo-IR rainfall signals [e.g., Adler et al., 2003;Huffman et al., 1995;Sorooshian et al., 2000;Huffman et al., 2001;Hong et al., 2004;Huffman et al., 2007]. Furthermore, a few studies have addressed methodologies to optimally combine the products of the TRMM precipitation radar (PR) with the TRMM microwave imager (TMI) using Bayesian inversion and weighted least squares approaches [e.g., Masunaga and Kummerow , 2005;Kummerow et al., 2010]. From another direction, Gaussian filtering methods on Markovian tree-like structures, the so-called scale-recursive-estimation (SRE), have been proposed to merge spaceborne and ground-based rainfall observations at multiple scales [e.g., Gorenburg et al., 2001;Tustison et al., 2003;Bocchiola, 2007; Van de Vyver and Roulin, 2009;Wang et al., 2011], see also [Kumar , 1999] for soil moisture application.
Recently, using the Gaussian scale mixture probability model and an adaptive filtering approach Ebtehaj and Foufoula-Georgiou [2011b] proposed a fusion methodology in the wavelet domain to merge TRMM-PR and ground-based NEXRAD measurements, aiming to preserve the non-Gaussian structure and local extremes of precipitation fields.
Data assimilation has played an important role in improving the skill of environmental forecasts and has become by now a necessary step in operational prediction models [see Daley, 1993]. Data assimilation amounts to integrating the underlying knowledge from the observations into the first guess or the background state, typically provided by a physical model. The goal is then to obtain the best estimate of the current state of the system, the so called analysis. The analysis with reduced error is then used to forecast the state at the next time step and so on [see Daley, 1993;Kalnay, 2003, for a comprehensive review].
Note that, in practice, the present background state, used in each analysis cycle, is commonly a forecast of the state by the underlying model, initialized by the analysis state in the previous time step. One of the most common approaches to the data assimilation problem relies on variational techniques [e.g., Sasaki, 1958Sasaki, , 1970Lorenc, 1986;Talagrand and Courtier , 1987;Courtier and Talagrand , 1990;Parrish and Derber , 1992;Zupanski, 1993;Courtier et al., 1994;Reichle et al., 2001b, among many others]. In these methods, one explicitly defines a cost function, typically quadratic, whose potentially unique minimizer is the analysis state. Very recently Freitag et al. [2012] proposed a regularized data assimilation scheme to improve assimilation results in advective sharp fluid fronts.
What is common in the DS, DF, and DA problems is that, in all of them, we seek an improved estimate of the true state given a suite of noisy and down-sampled observations and/or uncertain model-predicted states. Specifically, let us suppose that the unknown true state in continuous space is denoted by x(t) and its indirect observation (or model-output), by y(r). Let us also assume that x(t) and y(r) are related via a linear integral equation, called the Fredholm integral equation of the the first kind, as: where H(r, t) is the kernel relating x(t) and y(r). Recovery of x(t) knowing y(r) and H(r, t) is a classic linear inverse problem. Clearly, the deconvolution problem is a very special case with the kernel of the form H(r − t), which in its discrete form, plays a central role in this paper. Linear inverse problems are by nature ill-posed, in the sense that they do not satisfy at least one of the following three conditions: (1) Existence; (2) Uniqueness; and (3) Stability of the solution. For instance, when due to the kernel architecture, the dimension of the observation is smaller than that of the true signal, infinite choices of x(t) lead to the same y(r) and there is no unique solution for the problem. For the case when y(r) is noisy and has a larger dimension than the true state, the solution is typically very unstable, because, the high frequency components in y(r) are typically amplified and spoil the solution. In fact, the higher the frequency, the larger the amplification in the solution, which is often called the inverted noise; see, e.g., Hansen [2010] for a comprehensive account on linear inverse problems. A common approach to make an inverse problem well-posed is via the so-called regularization methods. The goal of regularization is to reformulate the inverse problem aiming to obtain a unique and sufficiently stable solution. The choice of regularization typically depends on the continuity and degree of smoothness of the state variable of interest, often called the regularity condition. For instance, some state variables are very regular with high degree of smoothness (e.g., pressure) while others might be more irregular and suffer from frequent and different sorts of discontinuities (e.g., rainfall). In fact, the choice of regularization not only yields unique and stable solutions but also reinforces the underlying regularity of the true state in the solution. It is important to note that, different regularity conditions are theoretically consistent with different statistical signatures in the true state, a fact that may guide proper design of the regularization, as explored in this study.
The central goal of this paper is to propose a unified framework for the class of DS, DF, and DA problems by recasting them as discrete linear inverse problems using relevant regularizations in the derivative space aiming to solve them more accurately compared to the classic formulations. Examples on rainfall DS and DF are presented to quantitatively elaborate on the effectiveness of the presented methodologies.
In these examples, regularization is performed consistent with the non-Gaussian statistics of rainfall in the derivative space, which might be generalized to other land-surface signals such as soil moisture. For the DA family of problems, the promise of the presented framework, is demonstrated via an elementary example using the heat equation. It turns out that the accuracy of the analysis and forecast cycles can be improved compared to the classic DA methods, especially when the initial state is discontinuous. This simple example provides insight into how the new formulations can outperform the classic methods and be of special interest in hydrometeorological applications.
This paper is structured as follows. Section 2 provides conceptual insight into the discrete inverse problems.
Section 3 describes the DS problem in detail, as a primitive building block for the other studied problems. Figure 1: Schematics of under-determined (a) and over-determined (b) linear discrete observation model in (2). (a) The observation matrix is fat with more columns than rows and x has a greater dimension than y.
(b) The observation matrix is skinny with more rows than columns and size of y exceeds the dimension of x.
Important classes of regularization methods are explained and their statistical interpretation is discussed from a Bayesian point of view. Examples on rainfall downscaling are presented in this section, taking into account the specific regularity and statistical distribution of the rainfall fields in the derivative space.
Section 3 is devoted to the regularized DF class of problems with examples and results again on rainfall data. The regularized DA problem is discussed in Section 4 with an illustrative example using the heat equation. Concluding remarks and future research perspectives are discussed in Section 5.

Discrete Inverse problems: Conceptual Framework
In this section, we briefly explain the conceptual key elements of discrete linear inverse estimation relevant to the problems at hand and leave further details for the next sections. Analogous to equation (1), linear discrete inverse problems typically amount to estimating the true m-element state vector x ∈ R m from the following linear observation model: where y ∈ R n denotes the measurement, e.g., output of a sensor, H ∈ R n×m is an n-by-m observation matrix and v ∼ N (0, R) is the Gaussian error in R n . Note that, the observation operator, which is a discrete representation of the kernel in (1), and the noise covariance are supposed to be known or properly calibrated. Depending on the relative dimension of y and x, this linear system can be under-determined (m n) or over-determined (m n) (see, Figure 1). Provided that H is full rank, in the first case, there are infinite different x's that satisfy (2) while for the second case a single exact solution does not exist.
As is evident, the DS class of problems belongs to the under-determined class because the sensor output is a coarse-scale and noisy representation of the true state. However, the class of DF and DA problems fall into the category of over-determined problems, as the total size of the observations and background state exceeds the dimension of the true state ( Figure 1).
In each of the above cases, we may naturally try to obtain a solution with minimum error variance by solving a linear least squares (LS) problem. However, for the under-determined case the solution still does not exist, while for the over-determined case it is commonly ill-conditioned and sensitive to the observation 5 noise (see, Section 4). Therefore, the minimum variance LS treatment can not properly make the above inverse problems well-posed. To obtain a unique and stable solution, beyond the LS minimization of the error, the basic idea of regularization is to further constrain the solution. For instance, among many solutions that fit the observation model in (2), we can obtain the one with minimum energy, mean-squared curvature or total variation. The choice of this constrain or regularization, highly depends on a priori knowledge about the underlying regularity of x. For sufficiently smooth x we naturally may promote a solution with minimum mean-squared curvature to impose smoothness on the solution. However, if a state is non-smooth and contains frequent jump discontinuities, a solution with minimum total variation might be a better choice. In subsequent sections, we explain these concepts in more detail for the DS, DF, and DA problems with hydrometeorological examples.

Problem Formulation
To put the DS problem in a linear inverse estimation framework, we recognize that in the observation model of equaiton (2), the true state x ∈ R m has a larger dimension than the observation vector y ∈ R n , m n. Note that throughout this work, a notation is adopted in which the vector x ∈ R m may also represent, for example a 2D field X ∈ R √ m× √ m which is vectorized in a fixed order (e.g., lexicographical).
As explained in the previous section, the DS problem naturally amounts to obtaining the best weighted least squares (WLS) estimatex of the high-resolution or fine-scale true state aŝ where, x 2 A = x T Ax denotes the quadratic-norm, while A is a positive definite matrix. In the above notation, 1 2 y − Hx 2 R −1 is called the cost function whose minimum is attained atx. Due to the ill-posed nature of the problem, this optimization does not have a unique solution. Specifically, taking the derivative and setting it to zero we get in which H T R −1 H is definitely singular. Indeed, many choices ofx lead to the same right-hand side.
To narrow down all possible solutions to a stable and unique one, a common choice is to regularize the problem by constraining the squared Euclidean norm of the solution to be less than a certain constant, that is Lx 2 2 ≤ const., where L is an appropriately chosen transformation and x 2 2 = i x 2 i denotes the Euclidean l 2 -norm. Note that, by putting a constraint on the Euclidean norm of the state, we not only narrow down the solutions but also implicitly suppress the large components of the inverted noise and reduce their spoiling effect on the solution.

6
Using the theory of Lagrangian multipliers the dual form of the constrained version of the optimization in where λ > 0 is the Lagrangian multiplier or the so-called regularizer.
The first term in (5) measures how well the solution approximates the given (noisy) data, while the second term imposes a specific regularity on the solution. In effect, the regularizer plays a trade off role between making the fidelity to the observations sufficiently large, while not imposing too much regularity (in this case, smoothness) on the solution. The smaller the value of λ, the more weight is given to fitting the (noisy) observations resulting in solutions that are less regular and prone to overfitting. On the other hand, the larger the value of λ the more weight is given to the regularity of the solution. The goal is to find a regularized solution, by finding a good balance between the two terms such that the it is sufficiently close to the observations while obeying the underlying property of the true state. As is evident, the Ltransformation of the state also plays a key role in the solution of equation (5). For instance, choosing an identity matrix implies that we are looking for a solution with a small Euclidean norm (energy), the so called least-norm solution, while if L represents an approximation of the second order derivative (i.e., Laplacian transform), the Lx 2 2 is a measure of the mean-square curvature of the state which imposes extra smoothness on the solution.
The problem in equation (5) is a smooth convex quadratic programming problem since the cost function is differentiable and its Hessian H T R −1 H + 2λL T L is always positive definite for any λ > 0, provided that L T L is positive definite. This problem is known as the Tikhonov regularization with the following analytical solutionx [Tikhonov et al., 1977]. It is easy to show that the covariance of the estimate is the inverse of the Hessian of the cost function in (5) which is H T R −1 H + 2λ L T L −1 . Notice that, for large scale regularization problems, this closed form solution can not be directly computed and efficient iterative methods (e.g., Newton's method) need to be employed.
Depending on the selected transformation and the intrinsic regularity (degree of smoothness) of the underlying state, other choices of the regularization term are also common. For example, in case the Ltransformation projects the major body of the state vector onto (near) zero values, a common choice is the l 1 -norm instead of the l 2 -norm in (5) [e.g., Chen et al., 1998]. Such a property is often referred to as sparse representation in the L-transformation space and gives rise to the following formulation of the regularized DS problem:x = argmin where, the l 1 -norm is x 1 = i |x i |. Note that, the problem in (7) is a non-smooth convex optimization as the regularization term is non-differentiable and the conventional iterative gradient descent methods 7 are no longer applicable in their standard forms. Several optimization techniques are currently under development to directly solve the non-smooth inverse problem in (7) that may include: (1) The iterative shrinkage thresholding algorithms pioneered by Tibshirani [1996]; (2) The basis pursuit by Chen et al. [1998]; (3) Constrained quadratic programing [e.g., Figueiredo et al., 2007]; (4) Proximal gradient based methods [e.g., Beck and Teboulle, 2009]; and (5) Interior point methods [e.g., Kim et al., 2007].
One of the common approaches to treat the non-differentiability in (7) is to replace the l 1 -norm with a smooth approximation, the so called Huber norm, and T denotes a non-negative threshold ( Figure 2). The Huber norm is a hybrid one that behaves similar to the l 1 -norm for values greater than the threshold T while for smaller values it is identical to the l 2norm. From the statistical regression point of view, the sensitivity of a norm as a penalty function to the outliers depends on the (relative) values of the norm for large residuals. If we restrict ourselves to convex norms, the least sensitive ones to the large residuals or say the outliers are those with linear behavior for large input arguments (i.e., l 1 and Huber). Because of this property, these l 1 -like norms are often called robust norms, [Huber , 1964[Huber , , 1981Boyd and Vandenberghe, 2004]. Throughout this paper, for solving l 1 -regularized inverse problems, we use the Huber relaxation due to its simplicity, efficiency and adaptivity to all of the concerning classes of DS, DF, and DA problems.  Figure 2: The l 1 , l 2 and Huber norms. The Huber penalty is a smooth relaxation of the l 1 -norm which acts quadratically for input values smaller than the threshold while it behaves linearly for larger inputs. For inputs with heavy tail behavior, linear penalization in the regularization term is advantageous compared to the quadratic penalty functions in which the cost function becomes dominated by a few large values in the tail of the distribution.
In downscaling of state variables containing frequent jumps and isolated singularities (e.g., local rainfall 8 extremes) using regularization in the derivative space (i.e., Lx is an approximate derivative measure), the l 1 -like norms typically preserve rapid variations and lead to improved and sharper solutions, compared to the overly smooth results of the l 2 -regularization; see, examples presented by Ebtehaj et al. [2012] for rainfall downscaling and other examples in the next section.

Statistical Interpretation
From the frequentist statistical point of view, it is easy to show that the WLS solution of (3) is equivalent to the maximum likelihood estimator (ML) given that the likelihood density is Gaussian, taking − log(·), one can find the minimizer of the negative log-likelihood function − log {p(y|x)} as, which is identical to the WLS solution of (3).
It is important to note that in the ML estimator, x is considered to be a deterministic variable (fixed) while y has a random nature.
On the other hand, in the Bayesian perspective, the regularized solution of equation (5) is equivalent to the maximum a posteriori (MAP) estimator where both x and y are considered of random nature. Specifically, using Bayes theorem, ignoring the constant terms in x and applying − log(·) on the posterior density p(x|y), we get The first term, − log p (y|x), is just the negative log-likelihood as appeared in the ML estimator and the second term is called the prior which accounts for the a priori knowledge about the density of the state vector x. Accordingly, the proposed Tikhonov regularization in (5) is equivalent to the MAP estimator assuming that the state, or the linear transformed state Lx, is a multivariate Gaussian of the form where the covariance is Q = L T L [e.g., Tikhonov et al., 1977;Elad and Feuer , 1997;Levy, 2008].
Clearly, the choice of the l 1 -norm in equation (7), implies that log p(x) ∝ Lx 1 or say the transformed state can be well explained by a multivariate Laplace density with heavier tail than the Gaussian case [e.g., Tibshirani, 1996;Lewicki and Sejnowski, 2000]. Similarly, selecting the Huber norm can also be interpreted which is equivalent to assuming the Gibbs density function as the prior probability model of the state [Geman and Geman, 1984;Schultz and Stevenson, 1994]. The equivalence between the regularized estimation approach, which imposes constraints on the regularity of the solution, and its Bayesian counterpart, which imposes constraints on the prior probability of the state is very insightful.
This relationship establishes an important duality which can guide the selection of regularization depending on the statistical properties of the state in the real or derivative space.

Problem Formulation and Settings
To downscale a remotely sensed hydrometeorological signal using the explained discrete regularization methods, we need to have proper mathematical models for the downgrading operator and also a priori knowledge about the form of the regularization term. Clearly, the downgrading operator in the presented framework needs to be a linear approximation of the sampling property of the sensor. If a sensor directly measures the state of interest, while its maximum frequency channel is smaller than the maximum fre- In the filtering operation, we basically slide the kernel over the field and sum the element-wise multiplication of the kernel elements with those of the field and then put the results at the center of the current position of the kernel. However, in the convolution operation we first rotate the kernel by 180 • and follow the same procedure.
As is evident, the smoothing kernel needs to be estimated for each sensor, possibly by learning from a library of coincidental high and low-resolution observations (e.g., Ebtehaj et al., 2012), or through a direct minimization of an associated cost. In the absence of prior knowledge, one possible choice is to assume that the sensor observes a coarse grained (i.e., non-overlapping box averaging) and noisy version of the true state. In other words, to produce a field at the grid-scale of s c × s c from a 1 × 1, this assumption is equivalent to selecting a uniform smoothing kernel of size s c × s c followed by a downsampling ratio of s c ( Figure 4a).
The discrete (high-pass) generalized Laplacian filter of size 3 × 3, where κ is a parameter ranging between 0 to 1. The Laplacian coefficients, obtained by filtering the 2D state with the Laplacian kernel, are approximate measures of the second order derivative. Throughout this paper, we choose κ = 0.5 which corresponds to the 2D standard second order differencing operation.
The choice of the regularization term also plays a very important role on the accuracy of the DS solution. This well behaved non-Gaussian structure in the derivative space mainly arises due to the presence of spatial coherence (correlation) in the rainfall fields and abrupt occurrences of piece-wise discontinuities (large gradients). In effect, over the large areas of uniform rainfall intensity, a measure of derivative translates rainfall values into a large number of (near) zero values; however, over the less frequent jumps and isolated high-intensity rain-cells, values of the derivative measure are markedly larger than zero and form the tails. Note that this non-Gaussianity is due to the intrinsic spatial structure of rainfall reflectivity fields and can not be resolved by a logarithmic or power-law transformation (e.g., Z-R relationship). It is easy to see that after applying the Z-R relationship on the reflectivity fields, the shape of the rainfall histogram remains non-Gaussian and still can be explained by the Laplace density (not shown here). In practice, the histogram of the derivatives may exhibit a thicker tail than the Laplace density, requiring a heavier tail model such as the Generalized Gaussian Density (GGD) of the form p(x) ∝ exp (−λ |x| p ), where p < 1 [see, Ebtehaj and Foufoula-Georgiou, 2011a]. However, using such a prior model gives rise to a non-convex optimization problem in which convergence to the global minimum can not be guaranteed.
Hence, the choice of the l 1 -norm (the Laplace prior) is indeed the closest convex regularization that can partially fulfill the strict statistical interpretation, while the actual prior might still be better explained by a heavier tail model than the Laplace density.
Following our observations related to the distribution of the rainfall derivatives, here we direct our attention to the Huber penalty function as a smooth approximation of the l 1 -regularization, Minimization of the above cost function can be easily achieved using first order efficient gradient descent based methods. However, as the rainfall fields are non-negative, we used the gradient projection (GP) method [Bertsekas, 1999, pp. 228], to solve the above problem in a constrained mode such that x 0 (see Appendix A).

Results
The same rainfall snapshot shown in Figure 5 has been used to examine the performance of the proposed regularized DS methodologies. Throughout the paper, to make the reported parameters independent of the range of intensity values, the rainfall reflectivity fields are first scaled into the range between 0 and 1; however, the results and Figures are presented in the true range.
To demonstrate the performance of the proposed regularized DS methodology, the NEXRAD high-resolution observation x was assumed as the true state while the low-resolution observations y were obtained by smoothing x with an average filter of size s c × s c followed by a downsampling ratio s c . Given the true state and constructed observations, we can quantitatively examine the effectiveness of the presented DS methods. In fact, selecting some common quality metrics, we demonstrate the effective improvement of those measures using the presented regularized DS framework.
Both the Huber and Tikhonov regularizations were examined to downscale the observations from scales 4 × 4 and 8 × 8 km down to 1 × 1 km ( Figure 6). A very small amount of white noise (i.e., standard deviation of 1e-3) was added to the low-resolution observations (equation 2), giving rise to a diagonal error covariance matrix. In both of the regularization methods, the regularization parameter λ was set to 5e-3 and 1e-2 for downscaling from 4-to-1 and 8-to-1 km in grid spacing, respectively. These values are selected through trial and error; however, there are some formal methods for automatic estimation of this parameter which are left for future work [e.g., the L-curve; see Hansen, 2010]. In our experiments, it turned out that small values of the Huber threshold T , typically less than 10% of the field maximum range of variability, lead to a successful recovery of isolated singularities and local extreme rainfall cells ( Figure 7).
In the studied snapshot, coarse graining of the rainfall reflectivity fields to the scales of 4 × 4 and 8 × 8 kilometers was equivalent to loosing almost 20 and 30 percent of the rainfall energy in terms of the relative Root Mean Squared Error (RMSE), RMSE r = x−x 2/ x 2 (see, Table 1). Note that, to compute the RMSE of the low-resolution observations, the size of those fields was extended to the size of the true field using the nearest neighborhood interpolation, that is, each low-resolution pixel was replaced with s c × s c pixels with the same intensity value. In addition to the relative RMSE measure, we also used three other metrics: (1) Relative Mean Absolute Error (MAE), MAE r = x−x 1/ x 1 ; (2) A logarithmic measure often called the peak signal-to-noise ratio (PSNR), PSNR = 20 log 10 ( max(x) /std(x−x)) where std (·) denotes the standard deviation and; (3) The Structural Similarity Index (SSIM) by Wang et al. [2004]. The PSNR in decibel (dB), represents a measure that not only contains RMSE information but also encodes the recovered range.
The latter measure varies between -1 and 1 and the upper bound refers to the case where the estimated and reference fields are perfectly matched. The SSIM metric is popular in the image processing community as it takes into account not only the marginal statistics such as the RMSE but also the correlation structure between the estimated and reference (true) image. This metric seems very promising for analyzing the forecast mismatch with observations in hydro-meteorological studies, especially when the systematic errors in the large scale features of the predicted state (e.g., displacement error) might be more dominant than the random errors; see Ebtehaj et al. [2012] for applications of SSIM in rainfall downscaling.
On average, it was seen that one third of the lost relative energy of the rainfall reflectivity fields can be restored via the regularized DS (Table 1). The l 2 -regularization led to smoother results and as the scaling ratio grows, this regularization was almost incapable to recover the peaks and the correct variability range of the rainfall field ( Figure 7). Typically, as expected, the Huber regularization results were slightly better than the Tikhonov ones. For large scaling ratios (i.e., > 4 × 4 km) the results of those methods tended to coincide in terms of the global quality metrics such the RMSE. However, using the Huber prior, the recovered range was markedly better than that by the Tikhonov regularization as reflected in the PSNR metric.    Table 1: Results showing the effectiveness of the proposed regularized DS by reducing the estimation error and increasing the accuracy of the rainfall fields. The first two columns refer to the values of the quality metrics obtained by comparing the constructed low-resolution observations with true 1 × 1 km reflectivity field. The performance of the Huber prior is slightly better than the l 2 -regularization, especially for the small scaling ratios (i.e., ≤ 4 × 4 km). † RMSE r : relative root mean squared error; MAE r : relative maximum absolute error; SSIM: structural similarity; and PSNR: peak signal to noise ratio.

Problem Formulation
Analogous to the DS problem in the previous section, here we focus on the formulation of the DF problem.
In the DF class of problems, typically, an improved estimate of the true state is sought from a series of lowresolution and noisy observations. Let x ∈ R m be the true state of interest while a set of N downgraded measurements y i ∈ R n i , i = 1, . . . , N , are available through the following linear observation model where n i m, H i ∈ R n i ×m and v i ∼ N 0, R i denotes uncorrelated Gaussian error in Compared to the DS family of problems, DF is more constrained in the sense that usually there are more equations than the number of unknowns, Σ N i n i m, giving rise to an overdetermined linear system.
As previously explained, naturally the linear weighted least squares (WLS) estimate of the true state, given the series of N observations, amounts to solving the following optimization problem: Notice that the solution of the above problem not only contains information about all of the available observations (Fusion) but also, with proper design of the observation operators, allows us to obtain an estimate with higher resolution than any of the available observations (Downscaling). Clearly, the inverse of each covariance matrix in (16) Figure 6g. The results indicate a better performance of the Huber regularization, especially for smaller scaling ratios. The Huber regularization yields sharper results and is more capable to recover high-intensity rainfall cells and the correct range of variability; see Table 1 for quantitative comparison using a suit of metrics.
For notational convenience, the above system of equations can be augmented as follows: where, the concatenated error vector v has the following block diagonal covariance matrix, Therefore, the DF problem can be recast as the classic problem of estimating the true state from the augmented observation model of y = Hx + v. Setting the gradient of the cost function in equation (16) to zero, yields the following linear system: In this case, the left hand side H T R −1 H is likely to be very ill-conditioned giving rise to an unstable LS solution, highly sensitive to any perturbation in the observation vector in the right-hand side (Figure 8c) [see, e.g., Elad and Feuer , 1997;Hansen, 2010].
One possible remedy for stabilizing the solution is regularization. Recalling the formulation discussed in the previous section, a general regularized form of the DF problem can be written aŝ where the convex function ψ L (x) can take different penalty norms such as: the smooth Tikhonov Lx 2 2 ; the non-smooth l 1 -norm Lx 1 ; and the smooth Huber-norm Lx Hub .
For the case of the linear Tikhonov regularization, computation of the solution requires to invert the Hessian (H T R −1 H + 2λ L T L) of the objective function similar to equation (6). Analogous to the DS problem, the covariance of the posterior distribution is H T R −1 H + 2λ L T L −1 .
Based on the selected type of regularization, statistical interpretation of the DF regularized class of problems is also similar to what was explained in Section 3.2. In other words, given the augmented classic observation model in (17), it is easy to see that the solution of (16) is the ML estimator while (20) can be interpreted as the MAP estimator with a prior density depending on the form of the regularization term.

Example on Rainfall DF
To quantitatively analyze the effectiveness of the regularized DF for rainfall data, we reconstructed two synthetic low-resolution and noisy observations from the original high-resolution NEXRAD reflectivity snapshot. To resemble different sensor constraints we chose different smoothing and downsampling operations for each of the reconstructed field. The first observation field y 1 was produced, at resolution 6 × 6 km, using a simple averaging filter of size 6 × 6 followed by a downsampling ratio of s c = 6. Analogously, the second field y 2 was generated at scale 12 × 12 km using a Gaussian smoothing kernel of size the 12 × 12      with a standard deviation of 4. A white Gaussian noise, with standard deviation of 1e-2 and 2e-2 was also added, respectively, to resemble the measurement random error. Roughly speaking, this selection of the error magnitudes implies that the degree of confidence (relative weight) on the observation at 6 × 6 km is twice as large as for the one at 12 × 12 km scale. According to the selected smoothing and downsampling operations, the downgrading operators H 1 and H 2 are designed to produce a high-resolution field at the scale of 1 × 1 km. To solve the DF problem, we have used the same settings for the Gradient Projection (GP) method as explained in Appendix A.
The solution of the ill-conditioned WLS formulation or the ML estimator in (16) is blocky, out of range and severely affected by the amplified inverted noise (Figure 8c). On the other hand, the Huber regularization can properly restore a fine-scale and coherent estimate of the rainfall field. The results show that almost 30% of the uncaptured subgrid energy of the rainfall field can be restored through solving the regularized DF problem ( Table 2). Improvements of the selected fidelity measures in the DF problem is more pronounced than the results of the DS experiment. This naturally arises, because more observations are available in the DF problem than the DS one and thus, the results are likely to be more accurate.  Table 2: Values of the selected fidelity metrics in the rainfall DF experiment using the Huber regularization, see the text for the definitions. Here, metrics refer to comparison of the low-resolution (6 × 6 and 12 × 12 km) observations the DF results with the true field (1 × 1 km).

Problem Formulation
Compared to the previously explained problems of downscaling and data fusion, the data assimilation (DA) problem is more involved in the sense that we need to address the evolution of a dynamical system and the available (low-resolution) observations at the same time in the estimation process. Despite the increased complexity, DA shares the same principles with the explained formulations of the DS and DF problems, from the estimation point of view. Here, we briefly explain the linear 3D and 4D-VAR data assimilation schemes and extend their formulations to a regularized format. Sample results of the regularized variational data assimilation problem are illustrated on the estimation of the initial conditions of the linear heat equation in a 3D-VAR setting.
The 3D-VAR is a memoryless assimilation method. In other words, at each time step, the best estimate of the true state or say analysis is obtained based only on the present-time noisy observations and background information of the dynamical system. The analysis is then being used for forecasting the state at the next time step an so on. Suppose that the true state of interest at discrete time t k is denoted by x k ∈ R m , a single noisy observation is y k ∈ R n , and x b k ∈ R m represents the background state. In the linear 3D-VAR data assimilation problem, obtaining the analysis state x a k ∈ R m amounts to finding the minimum point of the following cost function: In equation (21), B ∈ R m×m is the background error covariance matrix, H denotes the observation operator, while the analysis is the optimal solution: x a k = argmin Here we assume that the observation matrix is translation invariant. Clearly, this 3D-VAR problem is the WLS problem which has the following analytic solution Because the error covariance matrices are positive definite, the matrix B −1 + H T R −1 H is always positive definite and hence invertible. Thus, solution of the 3D-VAR requires no rank or dimension assumption on H. However, this problem might be very ill-conditioned depending on the architecture of the covariance matrices and the measurement operator.
The classic 4D-VAR is, indeed, an extension to the explained 3D-VAR which exploits the temporal memory of the system by constraining the solution (analysis) to the underlying dynamics. Assuming that we are interested in estimating the initial condition of a dynamical model at previous discrete time t k , to be used for an improved forecast in future time. The 4D-VAR assimilation method is formulated in such a way that uses the initial background state x b k and also the observations within a finite discrete time interval [t k , . . . , t k+T ]. Accordingly, this assimilation method amounts to obtaining the minimum point of the following cost function: while constraining the solution to the underlying dynamics by assuming that The model operator M k,i is a discrete linear representation of the system dynamics that evolves the state from t k to t i [e.g., Courtier and Talagrand , 1990;Daley, 1993;Zupanski, 1993;Kalnay, 2003]. Note that, in the above formulation it is implicitly assumed that the background and all measurement errors are mutually uncorrelated. Clearly, the linear 4D-VAR optimization contains a background cost which measures the weighted Euclidean distance of the true initial condition to the background state at the beginning of the interval and an accumulated cost, associated with all of the noisy observations within the selected time interval. Taking into account the imposed constraint by the model operator, the 4D-VAR cost function can be recast as, where its optimal solution, x a k = argmin is the analysis state at k th time-step. Thus, it is easy to see that, the linear 3D and 4D-VAR problems are in the category of the classic WLS problems which might be very ill-conditioned.
Analogous to the previous discussions, the generic regularized form of the linear 3D-VAR under the predetermined L-transformation might be considered as follows: x a k = argmin where ψ L (x) can take any of the explained regularization penalty functions including: the smooth Tikhonov Lx 2 2 ; the non-smooth l 1 -norm Lx 1 ; and the smooth Huber-norm Lx Hub . It is easy 21 to realize that the regularized 4D-VAR estimation of the the initial condition can also take the following form: x a k = argmin In the above regularized formulations, the solution not only becomes close to the background and observations, in the weighted Euclidean sense, but it is also enforced to follow a regularity imposed by the ψ L (x). Here, we emphasize that the regularization typically yields a stable and improved solution with less uncertainty; however, this gain comes at the price of introducing bias in the solution whose magnitude can be kept small, by proper selection of the regularizer [Hansen, 2010].

Statistical Interpretation
Statistical interpretation of the classic variational DA problems is a bit tricky compared to the DS and DF class of problems, mainly because of the involvement of the background information in the cost function.
Lorenc [1986] derived the 3D-VAR cost function using Bayes theorem and called it the ML estimator [see, e.g., Lorenc, 1988;Bouttier and Courtier , 2002]. More recently, it has been argued that the 4D-VAR, and thus as a special case the 3D-VAR cost function, can be interpreted via the Bayesian MAP estimator [Johnson et al., 2005;Freitag et al., 2010;Nichols, 2010]. For notational convenience, here we only explain the statistical interpretation of the 3D-VAR and its regularized version which can be easily generalized for the case of the 4D-VAR problem.
As discussed earlier, the ML estimator is basically a frequentist view to estimate the most likely value of an unknown deterministic variable from an (indirect) observation with random nature. The ML estimator intuitively requires to find the state that maximizes the likelihood function asx M L = argmax x p (y|x). Let us assume that, at time step t k , the background x b k is just a (random) realization of the true deterministic state x k . In other words, we consider x b k = x k +w, where the error w can be well explained by a zero mean Gaussian density N (0, B), uncorrelated with the observation error, E wv T = 0. Here, the background state is treated similar to an observation with random nature. Thus, let us recast the problem of obtaining the analysis as a classic linear inverse problem by augmenting the available information in the from of Notice that R is block diagonal because the background and observation errors are uncorrelated. Following the augmented representation and applying − log (·), we have − log p(y|x k ) ∝ 1 /2(y − Hx k ) T R −1 (y − Hx k ) and thus it is easy to see that the ML estimator in terms of the augmented observations, x a k = argmax x k p y|x k , is equivalent to minimizing the 3D-VAR cost function in (21). Therefore, following this statistical interpretation, the classic 3D-VAR, can be derived via the frequentist ML estimator.
On the other hand, from the Bayesian perspective, the state of interest and the available observations are considered to be random and the MAP estimator is the optimal point which maximizes the posterior density asx M AP = argmax x p (x|y). Let us assume a priori that the (random) state of interest has a Gaussian density with the mean x b and covariance B, that is p(x k ) ∼ N (x b k , B). More formally, this assumption implies that the deterministic background is the central (mean) forecast and is related to the random true state via N (0, B). Therefore, using Bayes theorem; see equation (12), it immediately follows that the 3D-VAR is the MAP estimator with the assumed Gaussian prior for the true state, x a k = argmax x k p (x k |y).
In conclusion, the regularized 3D-VAR in (25) might be interpreted as the MAP estimator, x a k = argmax with the prior density, p(x k ) ∝ λ ψ L (x k ), when we follow the frequentist approach and use the augmented notation to interpret the classic 3D-VAR as the ML estimator. On the other hand, taking the MAP interpretation for the classic 3D-VAR, the regularized version might be understood as the MAP estimator which also accounts for an extra and independent prior on the distribution of the state under the L transformation.

Heat Equation Example
The promise of the proposed regularized 3D-VAR data assimilation methodology, is shown via assimilating For a space-time representation of a 1D scalar quantity x (s, t), the well-known heat equation is x(s, 0) = x 0 (s), where ∞ < s < ∞, 0 < t < ∞. For mathematical treatment, let us assume that γ = 1 [ L 2 /T]. It is well understood that the general solution of the heat equation at time t is then given by the convolution of the 23 initial condition with the fundamental solution (kernel) as where We can see that x (s, t) is obtained via convolution of the initial condition by a Gaussian kernel with the standard deviation of σ = √ 2t. Clearly, estimation of the initial condition x 0 (s) only from the diffused and possibly noisy observations x(s, t) is an ill-posed deconvolution problem (see equation 1).
To reconstruct a 3D-VAR assimilation experiment, we assume that the true initial condition in discrete space is a vector with 256-elements (x ∈ R m where m = 256) as follows: the so-called top-hat initial condition. We added a white Gaussian noise with σ w = 0.05 to the true initial condition for reconstructing the background state x b 0 for the assimilation experiment.
We assumed that the observation vector is a downgraded version of the true state, while the sensor can only capture the mean of every four neighbor elements of the true state. In other words, the observation is a noisy and low-resolution version of the true state with one quarter of its size (Figure 9). To this end, using the linear model in (2), we employed the following architecture for the observation operator: and an n-element white Gaussian observation error with σ v = 0.03, where n = 64.
The top-hat initial condition is selected to emphasize the role of regularization, especially those with linear penalization (i.e., the Huber or l 1 -norm ). Clearly, the first order derivative of the above initial condition is very sparse. In other words, the first order derivative is zero everywhere on its domain except at the location of the two jumps, resembling a heavy-tailed and sparse statistical distribution. This the L-transform in (25), as follows: Note that, instead of incorporating a measure of curvature in the regularization term to impose smoothness on the solution, here we seek a solution with minimize total variation. Figure 9 shows the inputs of the assimilation experiment (left panels) and the results of the analysis cycle (right panels) using the classic versus the regularized 3D-VAR estimators. In this example, it is clear that the classic solution overfits, while slightly damps the noise. Indeed, the 3D-VAR is unable to effectively damp the high-frequency error components and impose the underlying regularity of the true state, because, its cost function is only about minimizing the weighted squared error. On the other hand, in the regularized assimilation methods not only the error term, but also a cost associated with the state regularity is also minimized. The Tikhonov regularization (T3D-VAR), i.e.,ψ L (x) = Lx 2 2 led to a smoother result compared to the classic one with slightly better error statistics. However, result of the Huber regularization (H3D-VAR), i.e., ψ L (x) = Lx Hub , is the best. The rapidly varying noisy components are effectively damped in this regularization, while the sharp jump discontinuities have been preserved better than the T3D-VAR. The quantitative metrics in Table 3, indicate that in the analysis cycle, the RMSE and MAE metrics are improved dramatically, up to 85% in the H3D-VAR scheme.
As previously explained, there is no unique and universally accepted methodology for automated selection of the regularization parameters, namely λ and T . Here, to select the best parameters in the above assimilation examples, we performed a few trial and error experiments. In other words, over a feasible range of parameter values, we computed the solutions of the regularized assimilation methods and obtained the RMSE measure by comparing those solutions with the true initial condition x 0 . Figure 10 shows the RMSE for different choices of regularization parameters in both T3D-VAR and H3D-VAR. Note that, the true initial condition is definitely not available in practice; however, here we picked the optimal values of the regularization parameters in the RMSE sense for comparison purposes and demonstrating the importance of a proper regularization. In the T3D-VAR, as expected, larger values of λ typically damp rapidly varying error components of the noisy background and observation; however, they may give rise to an overly smooth solution with larger bias and RMSE (Figure 10a). Here, for the T3D-VAR experiment, we picked the value λ T = 0.05 associated with the minimum RMSE (Figure 9). In the H3D-VAR, in addition to the regularizer λ H , we also need to choose the optimal thresholding value T of the Hubernorm. A contour plot of the RMSE values for different choices of λ H and T is shown in Figure 10b. By inspection, we roughly picked λ H =35 and T = 1.5e-3 for the H3D-VAR assimilation experiment presented in Figure 9.
The main purpose of the DA process is, indeed, to increase the quality of the forecast. Given the initialtime analysis state, we can forecast the profile of the scalar quantity, x (s, t), at any future time step through the heat equation. One important property of the heat equation is its diffusivity. In other words, naturally, noisy components and rapidly varying perturbations are damped but become more correlated as the profile of the initial condition evolves. Thus, rapidly-varying uncorrelated error components become low-varying and correlated features which their detection and removal are naturally more difficult than  Table 3: The root mean squared error (RMSE) and the mean absolute error (MAE) for the studied classic and regularized 3D-VAR in the analysis cycle (A) and forecast step (F) . The forecasts based on the classic 3D-VAR and the T3D-VAR almost coincide while the regularized result is marginally better. This behavior arises because neither of those methods could properly eliminate the noisy features in the analysis cycle and hence low-varying error components are appeared in the forecast profile. However, the quality metrics in Table 3 indicate that using the H3D-VAR, the RMSE and MAE of the forecast are almost improved more than 50% compared to the other methods.

Conclusion
In this paper, we presented a new direction in approaching hydrometeorological estimation problems by taking into account important continuity and statistical properties of the underlying states such as the presence of sharp jumps, isolated singularities (i.e., local extremes) and statistical sparsity in the derivative space. We started by explaining the concept of regularization and discussed about the common points of the hydrometeorological problems of downscaling (DS), data fusion (DF), and data assimilation (DA) as discrete linear inverse problems. We argued about the importance of proper regularization which not only makes an inverse problem well-posed but also imposes the desired regularity and statistical property on the solution. Regularization methods were theoretically linked to the underlying statistical structure of the states and it was shown how the statistical information about the density of the state, or its derivative, can be used for proper selection of the regularization method. Specifically, we emphasized  , Table 3). three types of regularization, namely, the Tikhonov, l 1 -norm and Huber regularizations. We theoretically showed that these methods are statistically equivalent to the maximum a posteriori (MAP) estimator while respectively assuming the Gaussian, Laplace and Gibbs prior density for the state. It was argued that piece-wise continuity of the state and the presence of frequent jumps are often translated into heavy-tailed distributions in the derivative space that favor the use of l 1 -regularization.
Informed by the non-Gaussian distribution of rainfall derivatives, the effectiveness of the regularized DS and DF problems was tested via analysis of remotely sensed precipitation fields. Then the performance of the regularized DA was studied via assimilating noisy observations into evolution of the heat equation.
We showed that regularized assimilation methods outperform the classic 3D-VAR method, especially for the case where the initial condition exhibits a sparse distribution in the derivative space (e.g., first order derivative of the top-hat initial condition).
The presented frameworks can be potentially applied to other hydrometeorological problems, such as soil moisture downscaling and fusion. Clearly, proper selection of the regularization method requires careful statistical analysis of the underlying state of interest. Moreover, the problem of rainfall or soil moisture retrieval from satellite microwave radiance, can be considered as a non-linear inverse problem.
This nonlinear inversion may be cast in the presented context, provided that the nonlinear kernel can be (locally) linearized with sufficient accuracy. Application of regularization in the DA problems is in its infancy [e.g.; see, Freitag et al., 2012, for a recent study] and is expected to play significant role over the next decades, especially for highly non-linear chaotic dynamical systems with frequent jumps. problem needs to be solved on the non-negative orthant {x| x i ≥ 0 ∀ i = 1, . . . , m}, We have used one of the primitive gradient projection (GP) methods to solve the above constrained DS problem [see, Bertsekas, 1999, pp. 228]. Accordingly, to obtain the solution of (A.4), it amounts to obtaining the fixed point of where α is a stepsize and denotes the Euclidean projection operator onto the non-negative orthant. As is evident, the fixed point can be obtained iteratively as Thus, if the descent at step k is feasible (i.e., x k − α k ∇J (x k ) 0) the GP iteration becomes an ordinary unconstrained steepest descent method, otherwise the result is mapped back onto the feasible set by the projection operator in (A.6). In effect, the GP method finds iteratively the closest feasible point, in the Euclidean distance, to the solution of the original unconstrained minimization.
In our study, the stepsize (α k ) was selected using the Armijo rule, or the so called backtracking line search, that is a convergent and very effective stepsize rule and depends on two constants 0 < ξ < 0.5 , 0 < ς < 1.
In this method, the step size is assumed α k = ς m k , where m k is the smallest non-negative integer for which J (x k − α k ∇J (x k )) ≤ J (x k ) − ξα k ∇J (x k ) T ∇J (x k ).