A self-adapting and altitude-dependent regularization method for atmospheric profile retrievals

MIPAS is a Fourier transform spectrometer, operating onboard of the ENVISAT satellite since July 2002. The online retrieval algorithm produces geolocated profiles of temperature and of volume mixing ratios of six key atmospheric constituents: H 2O, O3, HNO3, CH4, N2O and NO2. In the validation phase, oscillations beyond the error bars were observed in several profiles, particularly in CH 4 and N2O. To tackle this problem, a Tikhonov regularization scheme has been implemented in the retrieval algorithm. The applied regularization is however rather weak in order to preserve the vertical resolution of the profiles. In this paper we present a self-adapting and altitudedependent regularization approach that detects whether the analyzed observations contain information about small-scale profile features, and determines the strength of the regularization accordingly. The objective of the method is to smooth out artificial oscillations as much as possible, while preserving the fine detail features of the profile when related information is detected in the observations. The proposed method is checked for self consistency, its performance is tested on MIPAS observations and compared with that of some other regularization schemes available in the literature. In all the considered cases the proposed scheme achieves a good performance, thanks to its altitude dependence and to the constraints employed, which are specific of the inversion problem under consideration. The proposed method is generally applicable to iterative GaussNewton algorithms for the retrieval of vertical distribution profiles from atmospheric remote sounding measurements. Correspondence to: M.Ridolfi (Marco.Ridolfi@unibo.it)


MIPAS (Michelson Interferometer for Passive Atmospheric
Sounding, Fischer et al. 2008) is a Fourier transform spectrometer operating onboard of ENVISAT, a satellite launched by the European Space Agency (ESA) on 1 March 2002 in a polar orbit.MIPAS measures the atmospheric limb-emission spectrum in the middle infrared (from 685 to 2410 cm −1 ), a spectral region containing the signatures of the vibrational transitions of many atmospheric constituents.Beyond pressure at the tangent points, the ESA online retrieval algorithm (Ridolfi et al., 2000;Raspollini et al., 2006) produces geolocated profiles of temperature (T ) and of Volume Mixing Ratios (VMR) of six key atmospheric constituents: H 2 O, O 3 , HNO 3 , CH 4 , N 2 O and NO 2 .
The MIPAS measurements from July 2002 to March 2004, consisting of scan patterns of 17 sweeps in the 6-68 km altitude range with 3 km steps in the lower atmosphere, 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 4 and N 2 O VMR (Payan et al., 2009).Starting from January 2005 MIPAS is operated at a reduced spectral resolution with a nominal scan pattern consisting of 27 sweeps in the 6-68 km altitude range with 1.5 km steps in the lower atmosphere.The field of view of the instrument is approximately 3 km in the vertical, so the atmosphere turns out to be oversampled.Since the ESA retrieval grid coincides with the tangent altitudes of the measurements, the finer sampling of the vertical profiles is expected to amplify the unphysical oscillations already present in the measurements acquired until March 2004.
To tackle this problem, a Tikhonov regularization scheme has been implemented in the ESA retrieval algorithm.The choice of the Tikhonov parameter determines the trade-off M. Ridolfi and L. Sgheri: A self-adapting regularization method between the smoothing of the oscillations and the preservation of small-scale features.In the ESA retrieval the adopted choice for the strength of the regularization is rather conservative, to guarantee that small-scale profile features in the altitude domain are preserved (Ceccherini, 2005).
In this paper we present a self-adapting and altitudedependent regularization approach that detects whether the actual observations contain information about small-scale profile features, and determines the strength of the regularization accordingly.The objective of the method is to smooth out artificial oscillations as much as possible, while preserving the fine detail features of the profile when related information is detected in the observations.
In Sect. 2 we outline the theoretical background of the developed regularization scheme.In Sect. 3 we show a series of tests focusing on a single MIPAS limb scan.In Sect. 4 we analyze the performance of the method on the basis of a full MIPAS orbit, first using simulated data and then real measurements.Finally in Sect. 5 we draw conclusions and outline the future developments.

Theoretical basis
Ill-conditioning is a common feature of many inverse atmospheric problems.In the case of the retrieval of vertical atmospheric profiles from spectroscopic limb measurements, ill-conditioning produces oscillations in the retrieved profiles beyond the error margins defined by the mapping of the measurement noise into the solution.Altitude dependent systematic errors such as the forward/reverse differences analyzed in Kleinert et al. (2007) may also trigger oscillations.Tikhonov regularization is often used to improve the conditioning of the 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 mdimensional vector of the observations with error covariance matrix S y , f is the forward model, function of the ndimensional atmospheric state vector x.The Tikhonov solution is the state vector x t minimizing the following cost function: where • is the 2 norm, x s is an a-priori estimate of the solution, L is a l×n matrix operator, usually approximating the i−th order vertical derivative (i=0, 1, 2).Note that normally l=n−i.Finally λ is the non-negative scalar Tikhonov parameter driving the strength of the regularization.The first term of the right side of Eq. ( 1) is referred as χ 2 and represents the cost function minimized in the least-squares approach.
The choice of λ is a crucial step.If the selected λ is too small, a poor regularization will be achieved, whilst if λ is too large, Lx t will be strongly biased toward Lx s .Many general methods have been proposed for the selection of λ, such as cross validation (Allen, 1974), generalized cross validation (GCV) (Wahba, 1977 andGolub et al., 1979), the Lcurve method (LC) (Lawson andHanson, 1974 andHansen, 1992) and the discrepancy principle (Morozov, 1993).See e.g.Choi et al. (2007) for a recent paper on the comparison of the various techniques.See also the monographic issue of June 2008 of the Inverse Problems journal.
On the other hand better results may be expected if the operator L and the value of λ are adapted to the problem under investigation.The following references deal with the inversion of atmospheric state parameters.The LC method has been adopted in Schimpf and Schreier (1997) and more recently in Doicu et al. (2004).A-priori estimates of the degrees of freedom or of the retrieval error have been used by Steck (2002) to get λ.Alternatively Sofieva et al. ( 2004) tested both the discrepancy principle and vertical resolution requirements for the determination of λ.The error consistency (EC) method proposed by Ceccherini (2005) determines λ analytically by imposing the consistency of the difference between the regularized and the unregularized profiles with the error bars of the regularized profile.
In this paper we propose an altitude-dependent regularization scheme.Though there are more general mathematical formulations we only treat the case of a diagonal l×l positive semi-definite matrix .Assuming that d i x dz i z=z j ∼ (Lx) j , we may think of jj as the regularization strength at altitude z j .Thus we may speak of a vertical profile of .Then Eq. (1) becomes: Historically, the first idea of a variable regularization, the so-called localized Tikhonov regularization has been successfully used in the case of Volterra integral equations (see Lamm (1999) for a good survey).To the best of our knowledge, only few papers deal with variable regularization in other fields.In a recent paper Modarresi and Golub (2007) show that a vectorial version of the GCV achieves better results than the ordinary scalar version for an image reconstruction problem.
A method for calculating altitude dependent Tikhonov constraints for atmospheric retrievals was previously introduced in Kulawik et al. (2006), where combinations of 0th, 1st, and 2nd order Tikhonov constraints were considered.The polynomial weights of the constraints were selected by optimizing a function of the mean a posteriori error covariance and degrees of freedom.Pre-determined altitudedependent 0th and 1st order Tikhonov constraints are currently used for temperature, water, ozone, methane, and carbon monoxide nadir TES retrievals, as described in Bowman et al. (2006).In Steinwagner and Schwarz (2006), =λS h and S h is a diagonal matrix containing the reciprocal of the a-priori estimation of the profile.
Methods like LC, EC and the altitude-dependent methods proposed in this paper derive taking into account the achieved χ 2 .The χ 2 quantifies the compliance of the regularized profile with the available observations through the matrix S y .Oscillations produced by systematic errors not accounted for in S y may be nevertheless smoothed out, even if is determined regardless of them.
In this paper we test altitude-dependent regularization methods determining a profile of as the result of the minimization of a target function.After the -profile is obtained, Eq. ( 2) is solved via an iterative Gauss-Newton scheme.
Let k be the iteration count, K be the m×n Jacobian matrix of f in x k , then the Gauss-Newton iteration for the minimization of Eq. ( 2) is: (3) For the determination of , we will use the unregularized iterate solution x LS ≡ x k+1 ( =0).To get x LS , the inversion of K T S −1 y K is required.If this matrix is singular or too much ill-conditioned optimal estimation (OE) and/or Levenberg-Marquardt (LM) solutions (see e.g.Rodgers 2000) can be used instead of x LS .For example, if OE is used, Eq. ( 2) becomes: where x a is an a-priori estimate of the profile x with covariance matrix S a .Both x a and x s are estimates of the solution, however usually x a constrains the values of the profile, while x s constrains its derivatives.Therefore two different symbols are used.The iterative solution of Eq. ( 4) is: (5) When Eq. ( 5) is used, set x OE ≡ x k+1 ( =0).The LM solution with damping factor α can also be represented with Eq. ( 5), by setting S −1 a =αI and x a =x k .The OE/LM modifications permit to apply the method also in the case of singular or severely ill-conditioned K T S −1 y K matrices by adding the term S −1 a .In fact, with the OE/LM modifications the matrix to be inverted is K T S −1 y K+S −1 a .Since we want to use Tikhonov regularization, the only purpose of S −1 a is to permit the inversion of K T S −1 y K+S −1 a with reasonably small numerical errors.Therefore this term should be chosen positive definite, diagonal or diagonally dominant and kept as small as possible.
Fix any and let x be the profile minimizing Eq. ( 4).For moderately non-linear problems and a suitable initial guess, x k converges to x .The covariance matrix S mapping the measurement error S y into the solution x is given by: In the linear approximation, the spatial response function of x is represented (Rodgers, 2000) by the averaging kernel matrix (AK) A given by: Vertical resolution is a measure of the dispersion of the signal, usually calculated via the averaging kernel A .Still following Rodgers (2000), there are many practical ways of measuring the vertical resolution, such as the full width at half height (FWHH) of the AK rows: where z j , j =1, . . ., n are the altitudes, and . Note that AK rows not peaking at the diagonal element are penalized by Eq. ( 8), which in this case provides an overestimate of the FWHH.Throughout this paper, we use a modified version of Eq. ( 8), with |A | ij in place of (A ) ij in order to penalize also the negative lobes of the averaging kernel.When A =I, Eq. ( 8) provides the vertical step z i =(z i−1 −z i+1 )/2 of the retrieval grid.The Backus-Gilbert spread (Rodgers, 2000) is an alternative measure of the vertical resolution.However, in our tests it provided similar results while being more demanding from the computational point of view.

Altitude-dependent regularization methods
In this paper we compare a new altitude-dependent approach for the determination of that we call variable strength (VS) with two other altitude-dependent methods.A) In the VS method is determined as the minimiser of the following target function: where the bar over a vector stands for the average of the vector elements, and a superscript + stands for the positive part of a function.Finally, w e and w r are tunable parameters.Formula (9) selects a -profile minimizing the error of the regularized profile (first term of ψ VS ), with penalization terms that are effective when: (i) the χ 2 of the regularized solution (χ 2 ( )) increases beyond a threshold nw 2 e with respect to the χ 2 of the unregularized solution (second term of ψ VS ), and/or (ii) at some altitude z j the vertical resolution of the regularized solution is degraded beyond a factor w r with respect to z j (third term of ψ VS ).
In other words, we aim at the strongest regularization which leads to a profile (i) compatible with the available observations and (ii) with a vertical resolution not degraded beyond a pre-defined margin.
The calculation of ψ VS requires the evaluation of χ 2 ( ).This quantity is known only after the forward model f (x ) is calculated.Since this is a very time consuming operation, we use an approximation of χ 2 ( ) in Eq. ( 9).We have: It is possible to linearise f about x k , obtaining: Inserting Eq. ( 11) in Eq. ( 10), after some algebraic manipulations we obtain: where S −1 x ≡ K T S −1 y K.When no LM or OE modifications are employed, x OE =x LS and Eq. ( 12) may be further simplified.If we extract the term K T S −1 y (y−f (x k )) from Eq. (3) with =0 and plug it in Eq. ( 12) we obtain: Expression ( 13) shows the meaning of the factor w e in Eq. ( 9): on average the regularized and the unregularized profiles should differ by less than a fraction w e of the error bar of the unregularized profile.The averaging of residuals at different altitudes involved in the total χ 2 may in principle cause over-regularization if an isolated profile bump is encountered.Therefore we also tested some more restrictive versions of the second term of ψ VS , in which the χ 2 increase is penalized at each individual altitude, similarly to the vertical resolution.The results however did not change significantly.Therefore we preferred to stay with the formulation of Eq. ( 9) which checks the overall increase of the χ 2 , consistently with the actual implementation of the convergence criteria of the retrieval algorithm.Note that in the first term of Eq. ( 9) the errors (S ) jj are not individually normalized with (x ) j , to avoid singularities when the profiles approach zero.This happens frequently for example in the case of HNO 3 profiles above 30 km.Should this choice cause problems in case of a retrieved profile with a pronounced vertical variability (e.g.H 2 O), one of the following techniques may be used.The sum may be performed over (S ) jj /(x ) j (as suggested in the discussion phase of the present paper) and/or coefficients (either in the S or in the matrices) may be damped in the altitude regions where x is much larger than x .
The parameters w e and w r drive the strength of the regularization.As outlined above, these parameters reflect general requirements on the retrieval and therefore they do not depend on the shape of the actual profile.
B) In the vectorial version of the GCV approach the optimal value of is obtained as the minimiser of the target function ψ GCV defined as in Modarresi and Golub (2007).Within our framework ψ GCV becomes: This expression shows that the GCV method selects a -profile with the smallest possible number of degrees of freedom for the retrieval (given by trace(A )) compatibly with a small χ 2 ( ).In our implementation we calculate χ 2 ( )=χ 2 (0)+ χ 2 , with χ 2 given by Eq. ( 12) as in the VS method.The vertical resolution of the regularized profile is factored in the GCV method only through the χ 2 ( ).However, the χ 2 ( ) may be not sensitive to vertical resolution, e.g. when attempting the retrieval of a constant vertical profile.In this case the GCV approach produces profiles with dramatically degraded vertical resolution.On the other hand, when m n as in our case, the variation of the denominator in Eq. ( 14) may be marginal compared with that of the numerator.Therefore, even with a mild dependence of χ 2 ( ) on , the regularization produced by the GCV method may be very weak.
C) To overcome the drawbacks of the GCV approach, we also tested a scaled GCV method (SGCV).In this method we first find a -profile 0 with the GCV approach.Secondly we determine a scalar factor s 0 by minimizing ψ VS (s 0 ) with respect to s.Then we take s 0 0 as the final -profile.With this strategy we select the shape of the -profile determined with the GCV method, subsequently we tune the overall strength of the regularization according to the constraints of the VS method that are more specific to the inversion problem under consideration.

Minimization of the target function ψ
While the diagonal matrix has always dimension l×l, it is possible to represent the vertical -profile with fewer base points.The required (z j ) strengths are then calculated via linear interpolation in altitude between base points.This approach has the advantage of reducing the number of unknowns in the minimization of the ψ functions (ψ VS and ψ GCV defined in Sect.2), thus shortening the calculation time.The number of points used to represent the -profile however should be sufficient to allow an adequate altitude variability of the regularization strength.On the other hand it is useless to employ more than l points for the representation of the -profile.Note that the piecewise linear representation has the advantage that a change in a single coefficient produces only a localized change in the -profile.In our implementation this feature helps stabilizing the minimization process, in the sense that the minimal -profile does not depend on the initial guess.
Due to the large amount of local minima, analytical methods like conjugate gradients do not perform well when applied to the minimization of ψ.We found better results using the simulated annealing method (see e.g.Press et al., 1992, Sect. 10.9).
For the efficiency of the algorithm, we allow negative elements of the -profile and take the absolute values, instead of bounding them to be positive.Fine tuning of the -profile is not rewarded by the inversion procedure, therefore the minimization can be stopped as soon as the location of the minimum is approached.In this way it is possible to limit the computational overhead required by the minimization.The -profile corresponding to the minimum of ψ depends on the vertical shape of the actual profile x OE which, in turn, exhibits a large variability in the atmosphere.The lack of a preferred shape for the -profile makes it impossible to predict an a-priori annealing temperature for which the process should be stopped.To avoid useless calculations, the process should also be stopped when repeatedly failing to reduce significantly the target function.
We tried several implementations of the simulated annealing method, and we found the best results with the routine SA of Goffe (1994).The settings of this routine were optimized according to the guidelines mentioned above.In this way we achieved a much faster convergence compared to the standard settings suggested by the authors.

Retrievals from a single limb scan
We implemented the VS regularization method in the Optimized Retrieval Model (ORM, see Ridolfi et al. 2000;Raspollini et al. 2006) that is used by ESA for near real-time inversion of MIPAS data (Fischer et al., 2008).For comparison purposes, in the same code we also implemented the GCV and SGCV methods with a selectable switch.All of these methods can be applied either after each Gauss-Newton iteration or as a final step after the convergence of the inversion.However, in general we found that applying regularization after each iteration leads to heavier calculations and a slower convergence rate, with no benefits on the results (in agreement with findings reported in Ceccherini et al. 2007).As a consequence in all the test cases presented in this paper we applied the regularization only after reaching the convergence of the inversion.In all the tests presented we selected the regularization operator L=L 2 , the second derivative operator, with an exception for the EC method which is im-plemented in the ORM with L=L 1 .The choice L=L 2 produces slightly better results when the profile varies almost linearly with altitude.The regularization schemes take into account the LM approach employed by the ORM, as outlined in Sect. 2. Since in practical cases it is always difficult to find reliable a-priori profile estimates, in this paper we always select x s =0.The joint choice of L=L 2 and x s =0 implies that retrieved profiles with shape deviating from a straight line will be penalized by the Tikhonov regularization term.
3.1 Self-consistency test with synthetic observations from a single limb scan In this subsection we test the self-consistency of the VS method and its capability to detect possible sharp profile features measured by the instrument.For this purpose we carried out a test O 3 retrieval starting from synthetic observations.These observations were generated by the forward model included in the ORM (thus avoiding forward model errors), using the reference atmosphere model of Remedios et al. (2007) with the O 3 profile modified with a sharp bump in the 18-24 km altitude range.This modification reflects the double-peak feature sometimes observed in the real O 3 profiles for instance in pre ozone hole conditions (see Nemuc and Dezafra, 2005).Instrument features such as field of view, vertical scan pattern and spectral line-shape were adjusted to the MIPAS configuration adopted for the nominal reduced resolution measurements acquired from January 2005 onward (Dudhia, 2008).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 VS regularization was applied after the convergence of the unregularized (LS) solution, using parameters (w e , w r )=(1, 5).The altitude grid of the retrieved profiles consisted of 27 points, coinciding with the tangent points of the limb measurements, while the -profile consisted of 27−2=25 points.In this particular test case we disabled the LM modification in the ORM.
Figure 1 shows the results of the test.In panel (a) we show the reference profile (solid grey), the initial guess profile (dashed black), the unregularized LS solution (dotted red) and the regularized VS solution (solid blue).The initial guess profile was obtained by multiplying the climatological profile by a factor of 1.3 and with no bump modification.The VS method was able to distinguish between the oscillations of the LS solution due to lack of stability (mainly in the 40-70 km height range, where the error has been artificially amplified) and the real bump present in the reference profile.In the 40-70 km range the oscillations were smoothed out thanks to the large error bars of the LS solution.On the other hand the real bump was retained since the relatively small error bars in this altitude region prevented a strong smoothing.As required by the VS method, on average the VS profile is www.atmos-chem-phys.net/9/1883/2009/consistent with the LS profile within a fraction w e =1 of the LS error bars.This result is illustrated in panel (b) of Fig. 1 which shows the percentage retrieval errors of the LS (dotted red) and VS (dashed blue) solutions (obtained from Eq. 6) and the actual percentage difference between the VS and the reference profiles (solid blue), i.e. the actual error.We note that this difference is mostly consistent with the error of the VS solution.Only below 18 km the regularization introduces a noticeable smoothing error, which is not included in Eq. ( 6).This error is however consistent with the LS error bounds and is quite small in absolute value (<0.1 ppmv), the profile itself being very close to zero in this altitude range.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 VS method, should not be considered as representative of the real MIPAS performance.
Panel (c) shows the LS (solid red) and VS (solid blue) vertical resolutions.The LS vertical resolution, as mentioned in Sect.2, coincides with the vertical limb scanning step of the measurements.The dashed red line shows the maximum allowed vertical resolution for the VS solution, i.e. w r =5 times the LS vertical resolution.While this upper bound is never violated, we note that in the 20-37 km range this bound is not even approached.This is due to the simultaneous occurrence of small error bars of the LS solution and the changing slope of the reference profile.This combination prevents a stronger VS regularization by triggering the χ 2 penalization term in the ψ VS target function.Panel (d) shows the obtained -profile for the VS solution.Note that small values of are obtained whenever the above combination occurs.
Figure 2 shows the rows of the AK of the regularized profile.The number of degrees of freedom obtained for the VS profile (the trace of the averaging kernel) was 14.7.
The AK plotted in Fig. 2 is calculated with Eq. ( 7), thus assuming that the -profile derived with the VS method does not depend on the actual VMR profile encountered in the atmosphere.Therefore, the AK of Fig. 2 represents only locally (i.e. for the current ) the spatial response function of the measuring system.We point out that the large width of the AK rows for altitudes >40 km is due to the strong regularization triggered by the artificially amplified noise in the synthetic observations.

Tests with real observations from a single limb scan
In this subsection we present the results of retrievals based on real MIPAS measurements related to a single limb scan.For this analysis we selected scan number 060 of ENVISAT orbit 15451 from 12 February 2005.The approximate average latitude of the tangent points is 82 • South.This scan shows low stratospheric temperatures, hence a reduced S/N ratio that triggers oscillations in the unregularized retrieval.Moreover, it includes limb views with tangent altitudes penetrating the cloud-free upper troposphere.

Selection of VS parameters
The choice of the parameters driving the strength of the regularization is often a critical step when they have to be determined by the user on the basis of a tuning procedure.In fact the tuning is necessarily based on some assumed profile and therefore the results may not be optimal when there is a substantial difference between the actual and the assumed profile.
In the case of the VS method, the strength of the regularization (i.e. the magnitude of the -profile) is indirectly driven by the w e and w r coefficients.In principle w e and w r may be chosen arbitrarily and independently from each other.However, there is a positive correlation between allowed χ 2 increase and vertical resolution degradation, therefore not all the couples (w e , w r ) are equally meaningful.For instance, if a small w r is imposed, the difference between the regularized and unregularized profiles will be small, and therefore the increase of χ 2 will also be small, so that a large value of w e would make ineffective the related constraint.
This concept is illustrated in Fig. 3, which is a color map of the logarithm of the minimum of the target function ψ VS of Eq. ( 9) as a function of w e and w r for CH 4 retrieval.From this map we see that well chosen couples (w e , w r ) are those for which variations of the target function minimum occur for small variations of any of the two parameters.This situation occurs in Fig. 3 around the diagonal from bottom-left to top-right.Analogous maps for the other MIPAS retrieval targets show the same behavior for roughly the same values of (w e , w r ).Therefore we find that well chosen couples do not depend on the actual atmosphere and target profile.
From the previous considerations one may argue that a single parameter w e or w r could be sufficient to control both the vertical resolution degradation and the χ 2 increase.While this is true for most MIPAS retrieval targets, the double con- straint in ψ VS ensures, with very little overhead added, a good behavior even in some pathological conditions.These include, for instance, the case of an almost linear profile versus altitude, or the case of very large relative error bars such as in the case of NO 2 retrievals above 60 km.In both cases a single limitation on the χ 2 increase would lead to profiles with a dramatically degraded vertical resolution.In the NO 2 case, the regularized profiles become straight lines above 40 km, a physically unacceptable shape.On the other hand a single constraint on the vertical resolution leads to the loss of detailed features of the profile also in the case of relatively small error bars, such as in the double-peaked O 3 profile retrieval considered in Sect.3. Figure 4 reports the obtained CH 4 profiles for some well chosen (w e , w r ) couples.The related χ 2 increase with respect to the LM method is reported in the legend.We see that the couple (w e , w r )=(1, 5) permits to achieve a strong regularization with a marginal increase (0.56%) of the χ 2 .Considering the large oscillations detected in MIPAS profiles during the validation phase (see Sect. 1), we would suggest this choice for routine MIPAS retrievals and we will use this couple for the tests reported hereafter in this section.
Note however that a softer regularization (such as (w e , w r )=(0.6,3)) may be advisable for specific applications of the retrieved profiles.In fact strong regularization degrades significantly the vertical resolution, generally making more difficult the comparison to independent observations (Rodgers and Connor, 2003;Ridolfi et al., 2006Ridolfi et al., , 2007) ) or model results (Lahoz et al., 2007).

Comparison of altitude-dependent regularizations
In this subsection we briefly compare the VS method with the other altitude-dependent techniques (GCV and SGCV) introduced in Sect. 2. The purpose of this comparison is twofold.On one side we show that the -profiles obtained with the VS method have some correlation with those obtained with other more general methods such as the vectorial version of GCV.On the other hand we also show that the VS method achieves better results by implementing constraints specific to the inversion problem under consideration.
Figure 5 illustrates the results of the comparison for the retrieval of CH 4 .The obtained -profiles, reported in panel (d), show similar shapes as a function of altitude.As shown in panel (c) the GCV method produces a dramatic degradation of the vertical resolution in the 25-40 km altitude range.To restore the vertical resolution constraint of the VS method, the scaling factor of SGCV is less than 0.001.As a consequence the regularization achieved by the SGCV method is very weak, as confirmed by panel (a) and (b), reporting profiles and errors, respectively.Despite the generally large degradation in vertical resolution, the GCV method is not able to smooth out the feature of the LM profile in the 10-15 km range.On the other hand this objective is achieved by the VS method with only the marginal χ 2 increase mentioned above (0.56%).

Comparison of VS method with self adapting scalar regularizations
In this subsection we compare the VS method with the LC and EC scalar regularization methods already introduced in Sect.2, using the retrievals of CH 4 , O 3 and H 2 O as test cases.
The results are illustrated in Figs. 6, 7 and 8, respectively.We see that with the rather strong choice of (w e , w r )=(1, 5) the VS method is able to smooth out quite large oscillations, such as those in the H 2 O profile above the tropopause.Due to the large variability of the water profile across the retrieval altitude range, these oscillations could not be smoothed by any of the scalar methods considered.
On the other hand, in the ozone retrieval small error bars suggest that the feature in the 20-26 km range may be real.
In this case, both the EC and VS methods are able to preserve this feature, while the LC method smooths it out badly.
These results indicate that the VS method, due to its adaptive capability, is able to achieve a strong regularization while preserving small-scale profile features when the LM profile errors are small compared with the amplitude of the feature itself.On the other hand the structures at 15-20 km in the H 2 O LM profile and at 25-30 km in the CH 4 LM profile are smoothed out by the VS method.We do not know if these structures are real, however the point is: can we believe an oscillation or feature of the LM profile if its amplitude is comparable with the error bars?The answer depends on the specific application for which the profile is used.A smaller (w e , w r ) couple would maintain these structures in the regularized profiles, as explained in Sect.3.2.1.We believe that ultimately, in order to resolve the mentioned structures we would need a thinner instrument field of view and/or better signal to noise ratio in the measurements.

Results of retrievals from a full MIPAS orbit
In this section we analyze the performance of the VS method based on measurements from a full MIPAS orbit.Visual inspection of individual profiles from a large sample of scans is unpractical, so we introduce some quantifiers to characterize the average performance of the retrieval.
The first quantifier we consider is χ2 R , which is the arithmetic mean (on the orbit) of the normalized chi-square χ 2 R (see Bevington and Robinson, 2003) related to individual profiles.
To measure the smoothness of a profile we introduce an oscillation quantifier 2 that, for a single profile x i =x(z i ), i=1, . . ., n is defined as The quantity 2 represents the root mean square distance between each profile point x i and the linear interpolation at z i from the two adjacent points x i−1 and x i+1 .The factor 100 is introduced for better readability of the actual numbers.Note that 2 =0 if and only if the profile is a line.Moreover, when the z i are equispaced, 2 is proportional to the 2 norm of the discrete second derivative of the profile.We then take the arithmetic mean ¯ 2 (on the orbit) of the 2 related to individual profiles.
The last quantifier we consider is the number of degrees of freedom (DoF) of the retrieval.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 DoF/n on the orbit.
We compare the VS method with different (w e , w r ) couples with the LM (no regularization, only LM modification) and the EC methods.For each of the VS tests, -profiles with 9 base points have been used.We found that the LC  method poses some problems when there is no user supervision of the individual retrievals.In fact, the L-curve is a log-log plot between the squared norm of the constraint (second term of Eq. 1) and the χ 2 (first term of Eq. 1) for a range of values of the regularization parameter λ.Generally the shape of the L-curve is similar to that of the "L" letter, (hence the name) showing a single corner with maximum curvature.The LC method selects the value of λ corresponding to the corner of the L-curve.However, for the problem under investigation, we found that the L-curve is not always really L-shaped.In these cases the values of the λ parameter obtained for the maximum of curvature may be meaningless.
In this section we will use both synthetic and real MI-PAS measurements.Real measurements refer to the full EN-VISAT orbit 15451, already considered in the single scan tests.The orbit consists of 79 nominal scans (Dudhia, 2008).
Table 1.Summary of retrieval errors from a full orbit of synthetic MIPAS measurements.Average ( 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 x matrix, and is reported for reference.Several measurements related to scan 4 are however corrupted, therefore the retrieval is performed only on 78 scans.The synthetic measurements, generated from a known atmosphere, emulate the same acquisition scenario.

LS
The computational overhead introduced by the VS method depends on the number of base points used for the -profile, and on how often the -profile is updated (i.e.how often the minimization of ψ VS is carried out).Within our setup (9 base points and -profile updated every scan) the overall runtime increase is less than 20% with respect to the LM method.This is a quite encouraging result, considering that in operational retrievals fewer base points might be sufficient to achieve a good regularization and also that the -profile could be updated only when strictly necessary.Update of is in fact necessary only when the actual unregularized VMR profile encountered in the atmosphere is significantly different (i.e.beyond a few error bars) from the VMR profile used for the last calculation of .

Retrieval from synthetic MIPAS data
In this subsection we use synthetic measurements that assume the observation geometry of the full MIPAS orbit 15451.The synthetic measurements were generated using the same forward model setup as in Sect.3.1, but without any O 3 bump or noise amplification.Since the true atmosphere is known, we can characterize the error of the retrieved profiles with the mean ( x) and the standard deviation (σ ) of the difference between the retrieved and the true profile at each retrieval grid point along the orbit.This data is reported in Table 1 for each target species and retrieval method, the parameters (w e , w r ) used in the VS method are shown in parenthesis.Since the pure LS method did not always converge, we calculated the σ of the LS method (σ LS ) as the average of (S x ) jj over the altitudes z j and the scans of the orbit.
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.We report in Figs. 9, 10 and 11 the plots of the standard deviations for CH 4 , O 3 and H 2 O retrievals, respectively.The σ of the LS method (σ LS ) is calculated by binning (S x ) jj .
We report in Table 2, for each target species and each profile type, the values of χ2 R , ¯ 2 and DoF/n, except for the reference (REF) profile, for which only the ¯ 2 is defined.The last row of the table contains the percentage variation of χ2 , ¯ 2 and DoF/n with respect to the LM method, averaged over the retrieval targets.
From Table 1 we can see that, for all retrieved profiles x is smaller than the related σ , thus indicating that the bias of the retrieved profiles is not significant.Furthermore, the standard deviation of all the retrieved profiles is smaller than σ LS .The estimated σ contains the contributions of both the smoothing error and the retrieval noise (Rodgers, 2000).Hence the smoothing error possibly introduced by the regularization is more than compensated by the achieved reduction of the noise error.We also note that both the EC and VS methods achieve a significant reduction of the σ with respect to the LM technique.Further comments arise from the inspection of the altitude behaviors of σ reported in Figs. 9, 10 and 11.As expected, the VS(1,5) achieves a stronger regularization than the VS(0.6,3).Therefore it obtains a smaller error in the altitude regions where the reference profile is close to a line (see e.g.CH 4 profile or the H 2 O profile above the tropopause).We note however that the σ of the profiles retrieved with VS(1,5) is smaller than σ LS in all altitude ranges, consistently with the choice of w e =1.
From Table 2 we see that the weakest VS regularization considered (w e , w r )=(0.6, 3) already provides on average both a smaller χ2 R increase and a larger ¯ 2 reduction with respect to the EC scalar method.A further reduction of the ¯ 2 is achieved by the VS method with (w e , w r )=(1, 5) at the expenses of a quite large (0.97%, i.e. double of that of the EC method) χ2 R increase.Because of this large χ 2 R increase, stronger regularizations (such as VS with (w e , w r )=(2, 8)) were not attempted with synthetic measurements.

Retrieval from real MIPAS data
In this subsection we use real MIPAS measurements related to the full ENVISAT orbit 15451.Table 3 shows the results of the test, with the same format of Table 2.
As in the case of synthetic measurements, we note that the weakest VS regularization considered (w e , w r )=(0.6, 3) already provides on average both a smaller χ2 R increase and a larger ¯ 2 reduction with respect to the EC scalar method.A further reduction of the ¯ 2 is achieved by the VS method  with (w e , w r )=(1, 5) with χ 2 R values close to those of EC.The VS method with (w e , w r )=(2, 8) achieves a further reduction of the ¯ 2 at the expenses of a quite large χ2 R increase.The advantages of the VS method are particularly noticeable in the case of the H 2 O, CH 4 and N 2 O target species.The H 2 O profile probably gets a particular benefit from different strengths of regularization that are applied above and below the tropopause.Above the tropopause a strong regularization can be applied since the profile is almost linear with altitude.Below the tropopause only a weak regularization can be applied since the profile deviates significantly from linearity.In the case of CH 4 and N 2 O, there are quite large altitude inter- vals where the profiles behave almost linearly so that the VS method can apply a strong regularization without significant χ2 R increase.We note that these are the two MIPAS species for which unphysical oscillations were reported in the validation phase (see Payan et al., 2009).
We may compare the results of the full orbit retrieval from synthetic (Table 2) and real (Table 3) measurements.We note that the ¯ 2 of the LM retrieved profiles is 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 orbit 15451.There is however a reassuring similarity in the behavior of the regularization methods.This can be seen from the averages reported in the last row of the tables, i.e. the same regularization method achieves similar variations of ¯ 2 and χ2 R with respect to the LM method.

Conclusions
In this work we introduce a new self-adapting method (VS) for determination of the altitude dependent strength of Tikhonov regularization.The method can be applied to the retrieval of vertical distribution profiles from observations sounding the atmosphere either at the limb or vertically.
We first prove the self-consistency of the implemented algorithm on the basis of synthetic limb-scanning observations.Secondly we test the method using both synthetic and real MIPAS observations.We compare the performance of the method with that of some scalar (LC and EC) and altitudedependent (GCV, SGCV) regularization schemes available in the literature.In all the tested cases the VS method achieves a better performance than the other methods, thanks to its altitude dependence and to the constraints employed, which are specific of the inversion problem under consideration.The self-adaptability of the VS method 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.Future work will include a further assessment of the performance of the VS method on the basis of difficult but realistic situations, such as polar winter CH 4 retrieval in the case of vortex air masses sounded in a limited vertical region or polar winter NO 2 retrieval in the presence of descended NO x produced by particle precipitation.
An additional task will be the optimization of the algorithm for operational MIPAS data analysis and its extension to 2-D retrieval schemes.
The proposed method can be implemented in any Gauss-Newton-type algorithm for the retrieval of vertical distribution profiles.Currently the VS algorithm is coded in a standard FORTRAN routine both in a stand-alone version and in a version interfaced with the ORM code.The routine can be easily interfaced with any existing inversion software.The authors will be happy to freely supply the VS routine to scientists that would like to test the algorithm in their inversion codes, for no-profit purposes.

Fig. 3 .
Fig. 3. Map of the logarithm of the minimum of ψ VS ( ) as a function of w e and w r for CH 4 retrieval.

Fig. 6 .
Fig. 6.CH 4 retrieval with LM (reference), VS, EC and LC methods.Profiles (left), errors (center) and vertical resolution (right).All VMR profiles but the leftmost are horizontally shifted for a clearer representation.

Fig. 7 .
Fig. 7. O 3 retrieval with LM (reference), VS, EC and LC methods.Profiles (left), errors (center) and vertical resolution (right).All VMR profiles but the leftmost are horizontally shifted for a clearer representation.

Fig. 8 .
Fig. 8. H 2 O retrieval with LM (reference), VS, EC and LC methods.Log-scale plot of profiles (left), errors (center) and vertical resolution (right).All VMR profiles but the leftmost are horizontally scaled for a clearer representation.

Fig. 10 .
Fig. 10.O 3 retrieval with LM, EC and VS regularization methods: Standard deviations of the difference between the retrieved and true profiles.The differences were binned around nominal MIPAS tangent altitudes.

Fig. 11 .
Fig. 11.H 2 O retrieval with LM, EC and VS regularization methods: standard deviations of the difference between the retrieved and true profiles.The differences were binned around nominal MIPAS tangent altitudes.
Comparison of CH 4 profiles retrieved with the VS method for various (w e , w r ) couples.The retrieved LM profile is shown for reference.All profiles but the leftmost are horizontally shifted by 0.4 ppmv each for a clearer representation.Legend includes percentage χ 2 increase with respect to the LM method. www.atmos-chem-phys.net/9/1883/2009/

Table 2 .
Summary of retrieval performances for a full orbit of synthetic MIPAS measurements.

Table 3 .
Summary of retrieval performances for real MIPAS measurements of orbit 15451.