Interpolating Sparsely Corrupted Signals in Micrometeorology

In real applications where data acquisition is carried out under extreme conditions, post-processing techniques for systematic corrections are of critical importance. In micrometeorological studies, it is often the case that acquired data contains both missing information and impulse noise due to instrumentation failure, data transmission and data rejection for quality assurance. In this work, we propose a simple algorithm based on an 1 − 1 variational formulation that simultaneously suppresses impulse noise and interpolates missing information. Our approach consists of relaxing the objective function in the variational formulation with a strictly convex and continuously differentiable function that depends on a regularization parameter. We solve a sequence of strictly convex optimization subproblems as the regularization parameter goes to zero, converging to the solution of the original problem. Numerical experiments on real micrometeorological datasets are conducted showing the effectiveness of the proposed approach. Furthermore, a convergence analysis is presented providing theoretical guaranties of our method.


Introduction
In Micrometeorology, the tower-based Eddy Covariance (EC) method is widely used to measure the exchange of carbon dioxide, water vapor, and energy between terrestrial ecosystems and the atmosphere.From EC measurements, patterns and controls of ecosystem exchange are typically deduced and the potential for environmental impacts such as desertification of arid rangelands, deforestation of tropical forest, and warming of tundra landscapes are assessed.
Importantly, the EC method can be limited by environmental variability (hilly areas, heterogeneous canopy height, downwind turbulence), extreme environmental events (freezes, snow, dust storms, hail, lightening) and general instrumentation requirements (insufficient or none power supply, incorrect calibration, transmission), thereby producing incomplete and sparsely corrupted (spiky) datasets (Moncrieff et al., 1996).
In order to compute accurate representations of land-atmosphere trace gas fluxes and energy exchange so that reliable conclusions can be drawn, raw EC data needs to be corrected and processed.As a consequence, gap filling and de-spiking procedures need to be established (Falge et al., 2001).
Several gap filling methods have been developed and tested by the micrometeorological community.These methods include non-linear regression, Kalman filtering, artificial neural networks, mean diurnal variation, and the recently proposed wavelet sparse representation (Ramirez et al., 2012).A comprehensive comparison of such methods can be found in (Moffat et al., 2007).These methods, however, only address the interpolation process, and consider the removal of spikes as a separate task.In this work, we propose a simple algorithm based on an an 1 − 1 variational formulation that simultaneously suppresses the impulse noise (i.e., sparse corruption) and interpolates the missing information.
Although the 1 − 1 variational model has been a subject of interest for the signal processing community in recent years, most of this interest has been focused on image processing applications such as denoising and deblurring in the presence of impulse noise (Bar & Kiryati, 2006;Nikolova, 2004).To the best of our knowledge, research aimed at interpolating signals in the presence of impulse noise has been very limited, perhaps with the exception of Wang et al. (2011) andYan (2011) which, nonetheless, fall in the image processing setting.This paper is organized as follows.In Section 2, we present the mathematical formulation of the problem.In Section 3, we describe our algorithm and give theoretical guaranties that motivate the numerical implementation.In Section 4, we present numerical results involving real micrometeorological data supporting the proposed method.In Section 5, we present the relevance of this study and how it relates to previous work.Finally, Section 6 presents some concluding remarks.

Problem Formulation
An observed signal b ∈ R m is obtained after the original data set u ∈ R n is subject to a data loss process followed by a selective corruption in some elements of b.This degradation process can be modeled as an underdetermined linear system b = Hu + ν, where H ∈ R m×n (m < n) is constructed by removing from the n × n identity matrix, the n − m rows associated with the missing data.The vector ν ∈ R m is additive impulse or sparse noise in which the nonzero elements are drawn from a normal distribution with zero mean and standard deviation σ 0 .Thus, the difference Hu − b must be sparse at the solution, so that the fidelity term Hu − b 1 becomes a suitable choice for measuring the discrepancy between b and Hu.Furthermore, if u is assumed to be sparse in a given basis Ψ, that is u = Ψx with x sparse, this prior knowledge can be imposed at the solution with an 1 penalty term.Therefore, in order to recover or estimate the original signal (data), we formulate the following unconstrained optimization problem where A = HΨ and λ ∈ 0, A T sig(b) ∞ is a penalization parameter that balances the fidelity of the solution in the 1 norm sense, while promoting sparsity in x.The operator sig(•) is applied component wise, and returns the sign of the argument.

Methodology
The major challenge of solving problem (1) is dealing with the non differentiability of the 1 norm.While several numerical methods have proliferated for solving (1) in which the fidelity term is measured with the 2 norm, e.g.(Chen et al., 2001;Kim et al., 2007;Figueiredo et al., 2007;Wright et al., 2008;Argaez et al., 2011;Becker et al., 2010), research efforts aimed at solving (1) as it is presented here have been more limited.Such works include Fu et al. (2006), where problem (1) is posed as a linear programming problem and solved using interior point methods; Nikolova (2004), who considers the denoising case, i.e.A = I, and proposes a regularization term based upon local regularity replacing the second term in (1) by where N i is a neighborhood of x i .Finally, Bar and Kiryati (2006) consider the image deconvolution problem in the presence of impulse noise using an 1 fidelity term, and a Mumford-Shah functional to promote piece-wise regularity as well as edge preserving features.
In order to solve (1), we extend our ideas reported in (Ramirez & Argaez, 2013) where the non differentiability of the 1 norm is overcome by approximating the absolute value with a continuously differentiable and strictly convex function that depends on a regularization parameter μ k .More specifically, |x| ≈ x 2 + μ k where μ k → 0 as k → ∞.Therefore, our strategy consists of solving a sequence of strictly convex optimization subproblems of the form as k → ∞.
The main difference of this formulation and the one presented in (Ramirez & Argaez, 2013) is that the latter assumes normality in the distribution of residuals r i = a T i x − b i whereas in this work the residual vector r is modeled as sparse vector due to the presence of spikes in the measurements b.Now, notice that for any μ k > 0, Problem (2) is strictly convex since the Hessian matrix In fact, for any non-zero d ∈ R n , we have that Therefore for μ k > 0, the optimal set of ( 2) is nonempty and has at most one element.Put differently, the global minimizer x μ k of (2) exists and is unique.

Optimality Conditions
The partial derivative of f μ k (x) with respect to x i is Using matrix notation, the gradient of f μ k (x) is written more compactly as Consequently, the unique minimizer of (2) is characterized by the following system of nonlinear equations (3) Noticing that λD x is a positive definite matrix, we can write (3) as a fixed point equation of the form where Therefore, given an initial point x , we solve for x the following linear system where D r = diag((r ) 2 i + μ k ) with (r ) i = a T i x − b i , and D x = diag((x ) 2 i + μ k ).We continue this process iteratively until two consecutive points x and x are sufficiently close.This fixed point scheme motivates the algorithmic framework presented below.

The Algorithm
We propose a modification of the Fixed Point Least Squares-Preconditioned Conjugate Gradient (FPLS PCG) algorithm presented in (Ramirez & Argaez, 2013) for solving Problem (1).Step 5.
Using a Preconditioned Conjugate Gradient (PCG) solve for x : Step 8.
Improve approximate solution: x = x Step 10.
Remain inner loop: go to Step 6 Step 11.

end-if
Step 12. end-for Before concluding this section, we want to comment on the following two aspects.Firstly, it is important to realize that the main difference between the FPLS PCG algorithm presented here and the one in (Ramirez & Argaez, 2013) is in the Steps 6 and 7.In these Steps there is an extra re-weighting matrix D r associated with the residual vector r = Ax − b which accounts for the non-Gaussian noise present in the measurement vector b.Secondly, although any monotonically decreasing sequence {μ k } k is theoretically a guarantee for convergence to the optimal set of (1) as proved in next section, numerical experiments have demonstrated that the choice μ k = 1 2 k+15 performs well in practice for a broad set of problems.This is, however, not a general rule and the choice of {μ k } k can be considered a non-sensitive parameter.

Convergence Analysis
We show that any limit point of the sequence of minimizers {x μ k } μ k →0 of Problem (2) is contained in the optimal set of Problem (1).
Notice that the family of global minimizers {x μ } as μ → 0 is bounded.This is a consequence of the bound where the last inequality follows because f μ (x μ ) ≤ f μ (x) for all x ∈ R n , in particular for x = 0. Therefore, and the family {x μ } μ→0 is bounded.Hence, there exists a convergent subsequence of {x μ k } k which for simplicity will be denoted by {x μ k } k , and whose limit is written as where x * = arg min f (x) is a solution of (1).
Taking the limit as μ k → 0 in the above expression and by continuity it yields Consequently, f ( x) ≤ f (x * ) and therefore x is a minimizer of f (x).In other words, any limit point of the sequence of minimizers {x μ k } μ k →0 of Problem (2), is contained in the optimal set of (1).

Numerical Results
In  Our aim is to simultaneously estimate the missing information while suppressing the impulse or sparse noise.We further assume that the variable of interest can be sparsified by the Discrete Cosine Transform (DCT).This is a reasonable assumption in micrometeorology, since the variable of interest derives most of its features from seasonal and diurnal changes.Therefore, we solve Problem (1) where Ψ is the DCT transform, H is the subsampled identity matrix associated with the missing information, and b is the observed data.Figure 1 (left) shows the incomplete and corrupted datasets corresponding to CO 2 , H 2 O, AT and PRESS respectively.Figure 1 (right) shows a successfully recovered dataset after running our FPLS PCG algorithm for all four variables.
In order to assess the error associated with the increase of missing and corrupted elements, we take a 512-length fraction of uncorrupted data denoted by u * , and solve the recovery problem after artificial gaps and impulse noise are randomly added.The added artificial gaps range from 5% up to 20%, and the impulse noise from 1% up to 3%.In each case, we compute the 2-norm error defined by where û is the recovered signal.This experiment set up is run one hundred times, and the averaged results are reported in Table 1, which shows that even with 20% of missing pixels and 3% of impulse noise, a successful recovery is achieved.Figure 2 illustrates the recovery for the AT data set where 10% of the data is missing, and 2% of the elements are affected by impulse noise.

Input
Matrix A, and vector b.Output Return x = arg min { Ax − b 1 + λ x 1 } Parameters Penalization parameter λ, and maximum number of outer iterations N. Step 1.Initial point: x = A T (AA T ) −1 b Step 2. Outer loop Step 3. for k = 1 : N this section we use datasets from the Eddy Covariance tower established by UTEP's Systems Ecology Lab (SEL) and Cyber-ShARE center at the Jornada Basin Experimental Range, Las Cruces, New Mexico (USA).We select a 31-month data set of four 30-min averaged micrometeorological variables (Carbon Dioxide-CO 2 , Water Vapor-H 2 O, Air Temperature-AT and Atmospheric Pressure-PRESS).The variables H 2 O and AT have 44.953measurements where a total of 5.303 are missing or rejected, and nearly 1% are corrupted by impulse noise.On the other hand, the variables CO 2 and PRESS whose sensor platform was deployed at a later date, have 42.363 measurements where a total of 4.229 are missing or rejected, with the same percentage of impulse noise.

Figure 1 .
Figure 1.(Left) Corrupted and incomplete datasets of CO 2 , H 2 O, AT and PRESS respectively.(Right) Recovered datasets after applying the FPLS PCG algorithm

Figure 2
Figure 2. (Top) Original Signal (AT).(Middle) Incomplete and corrupted signal.10% of the data has been removed, and a 2% of impulse noise has been added.(Bottom) Recovered signal after applying the FPLS PCG algorithm with a relative 2-norm error of 0.075

Table 1 .
Scalability 2-norm error for different percentages of artificial gaps and impulse noise in a 512-length signal