1. Introduction

On March 1st, 2002 the European Space Agency (ESA) launched the ENVISAT satellite on a polar orbit. The satellite payload includes the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS, [1]). MIPAS sounds the mid-infrared atmospheric limb-emission spectrum. This region contains the signatures of many atmospheric constituents. ESA set up a retrieval algorithm [2, 3] to determine pressure at tangent points and geolocated profiles of temperature and of Volume Mixing Ratios (VMR) of some atmospheric constituents, among which H₂O, O₃, HNO₃, CH₄, N₂O and NO₂.

The full resolution MIPAS measurements acquired until March 2004 were extensively validated by several research groups (see Atmos. Chem. Phys. 2006 special issue on MIPAS). Oscillations beyond the error bars were observed in several MIPAS profiles, particularly in CH₄ and N₂O VMR [4]. The optimized spectral resolution adopted by MIPAS since January 2005 (see for instance [5]) is characterized by a finer vertical sampling step of 1.5 km in the upper troposphere - lower stratosphere region. The field of view of the instrument is approximately 3 km in the vertical, so the atmosphere turns out to be oversampled. The retrieval grid of the ESA processor coincides with the tangent points of the limb measurements. Therefore the oversampling of the optimized resolution measurements amplifies the unphysical oscillations already present in the full resolution measurements.

The EC (error consistency) [6] Tikhonov regularization scheme has been implemented in the ESA retrieval algorithm to damp the profile oscillations. The choice of the scalar, altitude-independent Tikhonov parameter, which determines the trade-off between the smoothing of the oscillations and the preservation of small-scale features, is rather conservative. This choice guarantees that small-scale profile features in the altitude domain are preserved. However, the performances of this method (as for any scalar regularization scheme) are degraded when the profile changes by several orders of magnitude in the retrieval altitude range. For this reason the regularization of water vapor profiles has been turned off in the operational ESA retrievals.

In [7] we proposed the VS (variable strength) regularization scheme, a self-adapting and altitude-dependent approach that detects whether the actual observations contain information about small-scale profile features, and determines the strength of the regularization accordingly. While representing an optimal solution, this method has two drawbacks. First, the method relies on an external minimization routine that needs to be carefully tuned to obtain optimal performances. Second, the computational requirements of the method, and specifically of the optimization routine, lead to an increase of the retrieval time of about 20%.

In this paper we propose the IVS (iterative variable strength) method, an alternative to the VS scheme. The IVS, while based on the same rationale of the VS, does not use any minimization routine. As a consequence the method is more robust, and the implementation is easier. The additional computational effort required amounts only to about 1.5% of the total retrieval time. This feature makes the IVS method suitable also for implementation in the near real-time processor of MIPAS data.

In Sections 2 and 3 we outline the theoretical background of the developed regularization scheme. In Section 4 we highlight the aspects specific to our implementation of the proposed method. In Sections 5 and 6 we present some tests illustrating the self-consistency of the method when synthetic measurements are used. In Section 7 we show results from real measurements. Finally in Section 8 we draw the conclusions.

For the readers’ convenience, we report here a list of acronyms used in the text to specify various techniques. This list supplements the definitions given in the text.

2. Theory

Tikhonov regularization is often used to improve the conditioning of atmospheric profile inversion. Smoother profiles are obtained by penalizing the oscillating solutions in the inversion formula. Let y = f(x) be the forward problem, where y is the m–dimensional vector of the observations with error covariance matrix S_y, f is the forward model, function of the n–dimensional atmospheric state vector x, whose components represent the unknown profile at altitudes z = z_j, j = 1,...,n. The Tikhonov solution is the state vector x_t minimizing the following cost function:

The first term of the right side of Eq. (1) is referred as χ² and represents the cost function minimized in the least–squares (LS) approach. The vector x_s is an a-priori estimate of the solution, L is a h×n matrix operator, usually approximating a linear combination of the i – th order vertical derivatives (i = 0,1,2). The h×h matrix Λ is diagonal, positive semi-definite and drives the strength of the regularization. Note that normally h ≤ n. Assuming that ${\left(\mathbf{L}\mathbf{x}\right)}_j\sim {\sum _{i=0}^2{c}_i\frac{d^{i}x}{d{z}^i}|}_{z={\tilde{z}}_j}$, we may think of Λ_jj as the regularization strength at z = z̃_j, j = 1,...,h. The altitudes z̃_j are determined by the choice of L, and may not necessarily coincide with the retrieval grid z_j, j = 1,...,n. Thus we may speak of a vertical profile of Λ.

The standard scalar Tikhonov regularization is obtained when Λ = λI. There are many different methods to determine λ, a good review may be found in [8]. The variable regularization of atmospheric profiles has been proposed in [9], where combinations of 0th, 1st, and 2nd order Tikhonov constraints were considered. The method was applied to nadir retrievals from the Tropospheric Emission Spectrometer (TES) in [10]. In [11], Λ = λS_h and S_h is a diagonal matrix containing the reciprocal of the a-priori estimation of the profile. The VS method was introduced and tested on MIPAS measurements in [7], where a comprehensive comparison of scalar and altitude-dependent regularization schemes may also be found.

In this paper we test an altitude-dependent regularization method determining a profile of Λ as the result of an iterative process. We choose to apply the regularization a-posteriori to obtain a better convergence rate, as reported in [12] and [7]. In other words, we first find the minimum of the χ² via a Gauss-Newton iterative method:

The Gauss-Newton scheme of Eq. (2) may fail if for some p ≤ k the matrix ${\mathbf{K}}_p^T{\mathbf{S}}_{\text{y}}^{-1}{\mathbf{K}}_p$ is singular or too much ill–conditioned to be inverted with satisfactory accuracy. In this case, the optimal estimation (OE) and/or Levenberg–Marquardt (LM) modifications (see e.g. [13]) can be used. For example, if OE is used, Eq. (2) becomes:

Let k be the iteration count at convergence, thus x_OE ≡ x_k+1, the unregularized solution. Let A_OE be the averaging kernel of x_OE and S_OE its measurement error covariance matrix. The ESA retrieval algorithm uses the LM solution, and A_OE and S_OE can be calculated with alternative algorithms of different sophistication as explained in [12] and [14]. The theory explained hereafter does not depend on the method used to determine A_OE and S_OE.

We can compute the regularized solution x_Λ as the Gauss-Newton iterate for the minimization of:

3. Altitude-dependent regularization

In the VS method [7] the Λ-profile is determined as the minimizer of the following target function:

The same objectives of the VS method can be pursued with an iterative approach, avoiding the expensive minimization of ψ_VS. In the IVS approach we define a Λ-profile λ(z) on a vertical grid so fine that we can consider it a continuous function. We start with a large λ⁽⁰⁾(z) = λ_max constant profile and decrease it iteratively until the following requirements are fulfilled:

Fix a threshold λ_min. Let J ⊂ {1,...,n} the set of indices of the altitudes z_j for which λ^(l)(z_j) > λ_min and:

4. Implementation details

We implemented the IVS regularization method in the Optimized Retrieval Model (ORM, see [2, 3]), the scientific prototype of the retrieval algorithm used by ESA for near real–time inversion of MIPAS data. In this code the LM modification is applied to the Gauss-Newton minimization with ${\mathbf{S}}_{\text{a}}^{-1}=\alpha \text{diag}\left({\mathbf{K}}_p^T{\mathbf{S}}_{\text{y}}^{-1}{\mathbf{K}}_p\right)$, and the EC and VS regularization methods are also implemented. The regularization method is selectable via a switch. In the IVS implementation we used h = n – 2, with the operator L defined so that:

We chose the thresholds λ_min = 10⁻², and λ_max = 10. A notable exception is the regularization of H₂O profiles where we set λ_max = 10³, and implemented a damping scheme which rapidly reduces λ⁽⁰⁾(z) from λ_max to λ_min below the tropopause. The aim of this modification is to avoid introducing a systematic bias in the H₂O profile in this altitude range. Any regularization in this region, even if within the error bars, would introduce a bias around or below the knee of the profile. This effect is due either to the large second derivative values necessary to model the profile knee around the tropopause (when L = L₂) or to the large values of the profile and its first derivative below the tropopause (when L = L₀ or L = L₁).

Unlike the IVS case, the VS computation load depends heavily on the number of grid points used to represent the Λ-profile. For this reason, the Λ-profile of the VS method was defined in our tests on a coarse grid of 9 points and then interpolated to the grid z̃_j, j = 1,..., n. Finally we set (w_e, w_r) = (1, 5) for all retrieval targets, according to the results obtained when testing the VS method on real observations [7]. The same values were used for both VS and IVS, since the parameters have the same physical meaning in both methods. Finally note that the EC method is implemented in the ORM with L = L₁.

5. Self-consistency test with synthetic observations from a single limb scan

In this section we test the self–consistency of the IVS method and its capability to preserve sharp profile features measured by the instrument. We repeated the test O₃ retrieval based on synthetic observations already used to assess the VS method [7]. To avoid the forward model error, these observations were generated by the forward model included in the ORM, using the reference atmosphere model of [15]. The O₃ profile was modified with a sharp bump in the 18–24 km altitude range, reflecting the double–peak sometimes observed in the real O₃ profiles in pre ozone hole conditions (see [16]).

Instrument configuration, including field of view, vertical scan pattern and spectral line–shape, was adjusted to the MIPAS optimized resolution scenario [17]. Spectral measurement noise was added to synthetic observations. For altitudes ≤ 40 km the noise was chosen consistent with MIPAS specifications; for altitudes > 40 km the noise was amplified by a factor 20 in order to obtain amplified oscillations in the unregularized retrieved profile. The initial guess profile was obtained by multiplying the climatological profile by a factor of 1.3 and with no bump modification. In this test the LM term was kept as small as possible (α = 10⁻³) in order to limit its regularization effect. Because of the artificially amplified errors, the LM term is however necessary to guarantee the convergence of the minimization sequence.

Figure 1 illustrates the results of the test. Panel (a) shows the reference profile (solid gray), the initial guess profile (dashed black), and the LM solution (solid blue). The same panel shows also the solutions regularized with the VS (solid green), EC (solid orange) and IVS (solid purple) methods. Panel (b) shows the LM retrieval errors calculated as the square roots of the diagonal elements of the covariance matrix S_OE (dashed blue). The same panel shows also the differences between the LM (solid blue), VS (solid green), EC (solid orange), IVS (solid purple) solutions and the reference profile, i.e. the actual errors.

 

Fig. 1 Simulated O₃ retrieval with amplified noise above 40 km and artificial bump added from 18 to 24 km: (a) Reference, initial guess and retrieved profiles; (b) estimated LM retrieval error and actual differences between retrieved and reference profiles.

Download Full Size | PDF

From Fig. 1 we see that all the considered regularization methods are able to preserve the sharp bump feature of the reference profile in the 18–24 km altitude range. Above 40 km however, the artificial oscillations induced by the large noise are better smoothed out by the VS and IVS methods. In fact, while all the presented methods are self-adaptive (i.e. they adjust the regularization strength on the basis of the error bars of the unregularized profile), the altitude-dependence of the VS and IVS methods permits to achieve a stronger regularization localized in the altitude range where the error bars of the unregularized profile are particularly large (above 40 km). From panel (b) we see that on average the VS and IVS profiles are consistent with the LM profile within a fraction w_e = 1 of the LM error bars as required by the algorithms. Note that the relatively large errors obtained in this test retrieval are mainly due to the artificial amplification of the measurement noise that we applied above 40 km. Therefore the results of this test, while useful to assess the consistency of the IVS method, should not be considered as representative of the real MIPAS performance.

6. Results of retrievals from a simulated full MIPAS orbit

In this section we show the results of the IVS method applied to a full orbit of synthetic MIPAS measurements. Since inspection of a large number of profiles is unpractical, we introduce a few scalar quantifiers to characterize the performance of the retrieval.

The consistency of the final simulated spectra with the synthetic measurements is measured by ${\overline{\chi }}_R^2$. This is the arithmetic mean (on the orbit) of the normalized chi-square ${\chi }_R^2$ (see [18]) related to individual profile retrievals.

We introduce an oscillation quantifier Ω₂ to evaluate the profile smoothness. For any profile x_i = x(z_i), i = 1,..., n we define

Any regularization scheme with strength Λ aims at a reduction of Ω̄₂, at the expenses of an increase of ${\overline{\chi }}_R^2$. To measure the efficiency of a regularization scheme, we introduce the quantifier

Finally we consider the number of degrees of freedom (DoF) of the retrieval, calculated as the trace of the averaging kernel matrix [13]. This quantifier represents the overall reduction of the vertical resolution as a consequence of the applied constraints. We divide this number by the number n of points of the retrieved profile, since this latter can vary from scan to scan due to cloud contamination. We then take the arithmetic mean $\overline{\text{DoF}/n}$ on the orbit. Note that $1-\overline{\text{DoF}/n}$ measures the fraction of the vertical resolution lost (with respect to the LS solution) because of external constraints. Therefore $1-\overline{\text{DoF}/n}$ is also an estimator of the overall strength of the combined effect of the LM technique and the a-posteriori regularization.

Since in this test the true atmosphere is known, we can also characterize the error of the retrieved profiles with the mean ( $\overline{\Delta x}$) and the standard deviation (σ) of the difference between the retrieved and the true profile at each retrieval grid point along the orbit.

In this section we still use synthetic MIPAS measurements generated assuming the climatological atmosphere of [15]. We however apply an interpolation within the tabulated latitude bands to avoid latitude intervals with constant atmosphere and large jumps at the edges. In this test retrieval errors due to interfering species and pressure and temperature error propagation in VMR are avoided by assuming the true atmospheric state except for the retrieval target. For the retrieval target the initial guess is set equal to the true profile multiplied by 1.3. We used the same convergence criteria as the ESA Level 2 Processor version 6.0. These consist in checking the inter-iteration variation of the state-vector and of χ². The thresholds were adjusted to keep the convergence error smaller than 10% of the error due to the measurement noise. A more comprehensive description of the convergence criteria may be found in the technical note [19].

We compare the IVS method with the LM, EC and VS methods. As mentioned in Sect. 1, the EC regularization is not used in the H₂O retrievals, hence the related results are missing.

In Table 1 we report $\overline{\Delta x}$ and σ for each target parameter and retrieval method. For reference, we also report the standard deviation of the LS profiles. Since the pure LS method does not always converge, we calculate the σ of the LS profiles (σ_LS) as the average of $\sqrt{{\left({\mathbf{S}}_{\text{LS}}\right)}_{jj}}$ over the altitudes z_j and the scans of the orbit.

Table 1. Summary of retrieval errors from a full orbit of synthetic MIPAS measurements. Average ( $\overline{\Delta x}$) and standard deviation (σ) of the differences between retrieved and true profile (K for temperature and ppmv for VMR). The standard deviation of the LS profile is estimated from the diagonal of the S_LS matrix, and is reported for reference.

View Table | View all tables in this article

To check the altitude behavior of the errors, we also broke down the differences between retrieved and true profiles into altitude bins centered around nominal MIPAS tangent altitudes. For each bin we calculated the mean and the standard deviation of the sample. The plot of Fig. 2 shows the altitude behavior of the standard deviations for H₂O retrieval. The σ of the LS method (σ_LS) is calculated by binning $\sqrt{{\left({\mathbf{S}}_{\text{LS}}\right)}_{jj}}$. Figure 3 shows the retrieved H₂O profiles for a sample scan.

 

Fig. 2 Water vapor retrieval with LM, VS and IVS methods. Standard deviations of the difference between the retrieved and true profiles. The differences are binned around nominal MIPAS tangent altitudes.

Download Full Size | PDF

 

Fig. 3 Water vapor retrieval with LM, VS and IVS methods. Initial guess, reference and regularized profiles for a typical scan of simulated orbit 15451.

Download Full Size | PDF

In Table 2 we report, for each target species and retrieval method, the values of ${\overline{\chi }}_R^2$, Ω̄₂, E and $\overline{\text{DoF}/n}$, except for the reference (REF) profile, for which only the Ω̄₂ is defined. The last row of the table contains the efficiency E averaged over the retrieval targets.

Table 2. Summary of retrieval performances for simulated MIPAS measurements of orbit 15451.

View Table | View all tables in this article

From Table 1 we can see that, for all retrieved profiles $\overline{\Delta x}$ is smaller than the related σ, thus indicating that the bias of the retrieved profiles is not significant. We also note that in general the tested regularization methods achieve a reduction of the σ with respect to the LM technique. The estimated σ contains the contributions of both the smoothing error and the retrieval noise [13]. Hence the smoothing error possibly introduced by the regularization is generally more than compensated by the achieved reduction of the noise error.

The standard deviation of the retrieved profiles is much smaller than σ_LS except for the H₂O case, where σ_LM is already greater than σ_LS. In the H₂O case the differences between retrieved and reference profiles, while mostly constant if normalized to the profile value, become very large below the tropopause. As a consequence the global σ_LM remains quite large, because its value is driven by the squares of the differences. The same effect does not apply to σ_LS, which is simply the arithmetic average of the $\sqrt{{\left({\mathbf{S}}_{\text{LS}}\right)}_{jj}}$. This can be seen from Fig. 2 where, conversely, each binned σ_LM is smaller than the related σ_LS. The effect of the damping of λ⁽⁰⁾(z) in the water retrieval can also be seen in Fig. 2, where σ_LM ∼ σ_IVS below approximately 12 km.

From Table 2 we see that a good reduction of the Ω̄₂ is achieved at the expenses of a marginal increase of the ${\overline{\chi }}_R^2$. Because of their altitude-dependent strength, the VS and IVS methods have a better efficiency in reducing the oscillations of the profiles. Despite the theoretical suboptimality, the IVS achieves comparable performances with respect to the VS method. The only marginal increase of the efficiency in the water vapor case is motivated by the fact that here, the Ω̄₂ quantifier is dominated by the contributions of the steep part of the profile below the tropopause, where the profile is only marginally regularized because of the damping we applied to λ⁽⁰⁾(z). Nevertheless the profile oscillations in the stratosphere are well smoothed out by the VS and IVS methods as illustrated in Fig. 3, even if this effect is not properly factored by Ω̄₂.

7. Retrieval from real MIPAS data

In order to verify that the good performance of the IVS method persists also when real MIPAS measurements are processed, in this section we report the results of a test based on real observations. We selected 12 orbits of MIPAS measurements from the year 2007. The measurements cover the four seasons to ensure that the test includes sufficient atmospheric variability. Table 3 shows the results of the test, with the same format of Table 2. In this case the reported quantities are averages on all the scans (1128 in total) of the considered orbits. Note that the H₂O retrieval in the EC case has been carried out with the LM method. The difference in the H₂O results with respect to the full LM case is mainly due to the method used to retrieve the temperature.

Table 3. Summary of retrieval performances for real MIPAS measurements of a sample of 12 orbits.

View Table | View all tables in this article

We may compare the results of the full orbit retrieval from synthetic (Table 2) and real (Table 3) measurements. The numbers reported in Table 3 show clearly that, globally, the performance of the regularization is preserved also when real measurements are processed. We note that the Ω̄₂ of the LM retrieved profiles is generally smaller in the synthetic case. This is due to the combination of two causes. First, synthetic measurements do not include systematic model errors which are present in the real observations. Second, the reference model atmosphere of the synthetic test retrieval is probably smoother than the actual atmosphere sounded by MIPAS in the selected orbits. The only exception is H₂O. In this case the presence of clouds in the sample of real measurements makes the average bottom altitude of the profiles higher than in case of orbit 15451, used to model the synthetic measurements. As already mentioned, in the case of H₂O the Ω̄₂ quantifier is dominated by the contributions of the profile below the tropopause, hence the smaller value obtained with real measurements.

8. Conclusions

In this paper we present a new approach to the variable strength (VS) Tikhonov regularization that we already proposed in a earlier paper [7]. The new approach avoids the two drawbacks of the VS method: the dependence on an external specifically tuned minimization routine and the increase of computing time by about 20%.

The new method we propose, the iterative variable strength (IVS), has the same rationale of the VS method. The aim is to apply the strongest possible regularization that, within pre-defined margins, preserves the compatibility of the regularized profile with the observations (only a small increase of chi-square is permitted) and does not destroy the profile vertical resolution. Instead of minimizing a target function, the IVS method starts with a large (strong) Tikhonov constraint at all altitudes, then decreases it locally and iteratively, until the user requirements on chi-square and vertical resolution are met.

We prove the self-consistency of the proposed algorithm on the basis of synthetic limb-scanning observations. In particular we show that the IVS method is able to preserve sharp profile features, such as ozone double peaks, that might be unexpected and therefore not modeled in the initial guess of the retrieval. Based on the analysis of a statistically significant set of synthetic observations we also quantify the bias and the smoothing error introduced by the regularization. It turns out that the bias is one order of magnitude smaller than the noise error. Furthermore, in total, the smoothing and the noise error components of the regularized profile never exceed the noise error of the unregularized solution.

We measure the performance of the regularization on the basis of its capability to damp the profile oscillations with a marginal chi-square increase. Based on both synthetic and real MI-PAS measurements we evaluate the performance of the proposed IVS method in comparison to the original VS algorithm and to the error consistency (EC) method currently implemented in the on-line ESA Level 2 processor for MIPAS. In all the tested cases the IVS method, while sub-optimal with respect to the VS, achieves similar good performance. Compared to scalar regularizations, the altitude-dependence of the VS and the IVS methods permits to better smooth-out profile oscillations when these are associated with large retrieval errors localized in specific altitude ranges. The self-adaptability of the VS and IVS methods permits to obtain a sufficiently strong regularization and, at the same time, the risk of over-smoothing sharp profile features is avoided when related information is present in the analyzed observations.

Compared to the VS, the IVS method does not rely on external minimization routines, therefore the method is more robust, its implementation is easier and the required computational effort is much smaller. For the tests presented in this paper the additional computing time required by the IVS method amounts to only 1.5% of the total retrieval time. The proposed method can be implemented in any Gauss-Newton-type algorithm for the retrieval of vertical distribution profiles. Currently the IVS algorithm is coded in a standard FORTRAN routine both in a stand-alone version and in a version integrated within the ORM code. The routine can be easily interfaced with any existing inversion software. The authors will freely supply the IVS routine to scientists who would like to test the algorithm in their inversion codes, for no-profit purposes.