Ozone profile smoothness as a priori information in the inversion of limb measurements

In this work we discuss inclusion of a priori information about the smoothness of atmospheric profiles in inversion algorithms. The smoothness requirement can be formulated in the form of Tikhonov-type regularization, where the smoothness of atmospheric profiles is considered as a constraint or in the form of Bayesian optimal estimation (maximum a posteriori method, MAP), where the smoothness of profiles can be included as a priori information. We develop further two recently proposed retrieval methods. One of them – Tikhonov-type regularization according to the target resolution – develops the classical Tikhonov regularization. The second method – maximum a posteriori method with smoothness a priori – effectively combines the ideas of the classical MAP method and Tikhonov-type regularization. We discuss a grid-independent formulation for the proposed inversion methods, thus isolating the choice of calculation grid from the question of how strong the smoothing should be. The discussed approaches are applied to the problem of ozone profile retrieval from stellar occultation measurements by the GOMOS instrument on board the Envisat satellite. Realistic simulations for the typical measurement conditions with smoothness a priori information created from 10-years analysis of ozone sounding at Sodankyl ä and analysis of the total retrieval error illustrate the advantages of the proposed methods. The proposed methods are equally applicable to other profile retrieval problems from remote sensing measurements.


Introduction
The problem of profile retrieval from remote sensing measurements is always under-determined: a continuous Correspondence to: V. F. Sofieva (viktoria.sofieva@fmi.fi) profile of an atmospheric constituent is reconstructed from a finite number of measurements. The problems of this kind belong to the class of ill-posed problems. Measurement data themselves do not uniquely determine the solution; therefore, some kind of prior constraint must be used to make the problem well-posed. Following Rodgers (2000), it can be done via: (a) A discrete representation. Dividing the atmosphere into layers according to a measurement structure and making certain assumptions about the profile within the layers (constant, linear, polynomial altitude dependence) transform the problem to a well-determined one.
(b) An ad hoc constraint such a smoothness of atmospheric profiles. It can be applied in the form of the Tikhonovtype regularization.
(c) Prior information about atmosphere, which includes a mean state and its covariance, or more generally, a prior probability density function of the state.
It is also possible to use a representation via orthogonal polynomials or wavelets, but these approaches are not discussed here.
A discrete representation appears in almost all practical solutions of inverse problems, and often the case (a) is referred (also in this paper) to as having no a priori. Although the discretization makes the problem well-determined, it can remain ill-posed or noise-amplifying. It means that noisy data lead to instability of inversion: the reconstructed profile has non-physical oscillations. Application of Tikhonov-type regularization (also referred to as the Twomey method) (Tikhonov and Arsenin, 1977;Twomey, 1977;Rodgers, 2000) helps to recover the stability of inversion. This method assumes a smoothness of an atmospheric constituent profile and therefore we can consider regularization as a kind of prior information.
If a priori information about the state is available (case c), the maximum a posteriori method (MAP) usually allows a significant improvement of retrievals. In the case of Gaussian noise and a priori, it gives the minimum variance solution.
The Gaussian distribution is often assumed in practice, thus introducing elements of "ad hoc guess" into the case (c).
In the current paper, we develop further two recently proposed methods which can be applied if the information needed for the MAP method is unavailable. The first method considers the smoothness of atmospheric profiles not as an ad hoc constraint, as it is done in Tikhonov regularization, but as a priori information. This method (Haario et al., 2004) "combines" cases (b) and (c) and can therefore be called the MAP method with smoothness a priori. The second method  chooses the regularization parameter in a Tikhonov-type scheme according to a target resolution, which is determined from requirements to resolve fine structures of the profiles.
For both methods we propose a grid-independent formulation, which is very important for measurements unevenly distributed in altitude. For example, the sampling vertical resolution of stellar occultation measurements depends on altitude and on the obliqueness of occultations (i.e. angle between the orbital plane and the direction to the star).
In Sect. 2 we describe the stellar occultation measurements by the GOMOS instrument and discuss the need for advanced data analysis. The advanced inverse methods that use a priori information of a different kind, as well as the availability of this information are discussed in Sect. 3. Comparison of the proposed inversion schemes based on realistic simulation is presented in Sect. 4. The a priori information used in simulation is created from analysis of 10 years of ozone sonde measurements at Sodankylä. The analysis of the total retrieval error, which includes components due to measurement noise and due to smoothing properties of inversion, illustrates advantages of the proposed methods.

GOMOS (Global Ozone Monitoring by Occultation of Stars)
is a stellar occultation instrument operating on board the Envisat satellite (http://envisat.esa.int/instruments/gomos) launched 1 March 2002. The GOMOS spectrometers measure the stellar spectrum from ∼140 km down to ∼15 km. The products retrieved from the UV-Visible spectrometer measurements are ozone, NO 2 , NO 3 , aerosol and air density vertical profiles. The basis for the geophysical data retrieval from GOMOS measurements is the transmission function (Kyrölä et al., 1993). In the GOMOS data processing, the inversion is split into two parts: the spectral inversion part and the vertical inversion part (Kyrölä et al., 1993). In the spectral inversion, horizontal column densities are retrieved from the atmospheric transmission data from which refractive effects and scintillation modifications have been removed. In the vertical inversion, vertical profiles are reconstructed from the horizontal column densities. In this paper, we will deal only with the vertical inversion. The vertical inversion of GOMOS is linear. By discretization of the atmosphere into layers according to the measurement structure (the number of unknowns is equal to the number of measurements), the forward model connecting the measurements (horizontal column densities N ) and vector of unknowns (profile ρ), can be written as where ε presents a vector of measurement noise. The discretized GOMOS vertical inversion problem (1) is well-conditioned (the condition number of the forward model matrix K for a typical occultation is ∼25), so it can be solved with the usual matrix inversion. The inversion is slightly noise-amplifying by a factor of ∼2 . However, in the case of dim stars, the reconstructed profiles are significantly contaminated with noise.

Characterization of retrieved profiles. Vertical resolution
The error estimate and the vertical resolution are two quantities which completely characterize retrieved profiles. The resolution of the retrieved profiles depends, in addition to the sampling resolution, on smoothing properties of a retrieval algorithm. In profile reconstruction, the resolution can be studied by computing the averaging kernel (Backus and Gilbert, 1970;Rodgers, 1990): where ρ is the retrieved profile and ρ is the true profile. In linear problems the averaging kernel matrix can be computed as where G is the inversion matrix and K ∞ is the forward model matrix in (infinitely) dense grid. One commonly used measure of the resolution of the retrieved profile (or a measure of the width of the averaging kernels) is the Backus-Gilbert spread s(z) (e.g. Rodgers, 2000): In the following, we use the Backus-Gilbert spread (4) as a measure of resolution. The GOMOS instrument is capable of retrieving the atmospheric profiles with a very good resolution. The sampling resolution for different line-of-sight azimuth angles is shown in Fig. 1. The measurement grid becomes denser in the lower part of the atmosphere due to refraction. The sampling resolution depends on obliqueness (an angle between orbital plane and direction to the star) of the occultation. For very oblique occultations it can be nearly twice better than for vertical ones. The averaging kernels of the pre-launch GOMOS inversion are sharply peaked (Fig. 2). For a typical occultation, the resolution (Backus-Gilbert spread) is ∼2 km in the mesosphere and upper stratosphere, and is less than 1 km in the lower stratosphere and troposphere.

Need for advanced data analysis
The signal-to-noise ratio in stellar occultation measurements strongly depends on stellar parameters (visual magnitude and effective temperature), so does the error of the vertical profile reconstruction.
In pre-launch simulations we observed significant nonphysical oscillations of reconstructed ozone profiles for dim stars (visual magnitude >2.5) but only for altitudes above 80 km and below 18 km (Fig. 3). Therefore, no a priori information or regularization were explicitly used in the GOMOS baseline inversion, as the lowermost altitude of the GOMOS measurements was expected to be ∼18 km.
The first year of GOMOS validation has shown that the noise level is higher than expected. Additional errors come from incomplete scintillation correction.
Analysing one reference data set, which includes more than 600 occultations in dark limb, we have found that a significant share of occultations (∼10%) is terminated at altitudes below 10 km. This enables GOMOS to probe the troposphere, but advanced inversion methods are required because of a low signal-to-noise ratio.
The vertical resolution achieved in the GOMOS measurements is better than the characteristic vertical scale of the ozone fine structures (1.0-1.4 km in the troposphere and in the lower stratosphere, according to the recent study of .
These features prompt us to apply advanced methods in the GOMOS inversion. Unfortunately, application of the MAP method or Tikhonov-type regularization degrades the vertical resolution. The question therefore is: how to achieve the stability of the inversion without critical degradation of resolution?

Inversion methods
Let us consider the linear forward model (1). All the methods considered below can be applied to any linearized inverse problem. They are also valid for general nonlinear problems, but nonlinear techniques like iterative minimization or the MCMC (Markov chain Monte Carlo ) method (Tamminen and Kyrölä, 2001), are required instead of matrix multiplication.
The essence of inclusion of a priori information is expressed by Bayes' formula describing the posterior probability density function (pdf) P (ρ|N) via a likelihood function P (N |ρ) and a prior pdf P prior (ρ): Provided that the prior pdf is known and assuming it is Gaussian, with a mean ρ a and a covariance C a :ρ a ∼N (ρ 0 , C a ), the maximum a posteriori (MAP) estimation of the retrieved profile ρ MAP can be presented in the following form (e.g. Rodgers, 2000): Here a measurement noise is also assumed to be Gaussian: ε∼N (0, C ε ). The validity of this assumption is verified in Tamminen (2004) with the MCMC method: the intermediate products of the GOMOS inversion -horizontal column densities -have a distribution close to Gaussian, provided that the instrumental noise is Gaussian. The matrix C ε is diagonal due to the independence of instrumental noise at each altitude.
A priori information on the profile smoothness can be included in the MAP method via correlation of prior uncertainties. Then the prior covariance C a can be presented, for example, in the form (Rodgers, 2000)  where r i denotes an altitude grid point, σ i is the standard deviation for a priori profile at the point r i and L is the characteristic scale. Alternatively, a Gaussian or a triangular shape of the correlation function can be used. This method includes the prior information in the most complete form. It is also grid-independent, i.e. the profile smoothness is defined independently of the discretization grid. The actual resolution of retrieved profiles depends only on the noise level: it is worse for larger measurement noise.
Unfortunately, at the moment this a priori information is often not available, even for ozone. The climatological data (Fortuin and Kelder, 1998) can, in principle, be used as an a priori estimate of the mean of ozone profile, but the inter-annual variability of the climatology does not reflect that of the individual profiles. It is illustrated in Fig. 4, where the mean and the standard deviation of ozone profile obtained from analysis of 10 years (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) of ozone sonde data at Sodankylä are compared with those of the monthly mean ozone profile from the Fortuin-Kelder climatology. The variability of individual ozone profiles is 1.5-2 times larger than the variability of monthly mean profiles.
Nevertheless, in a few cases when an occultation is located near an ozone observation station, useful a priori information can be obtained. Recently,  analysed 11-year ozone sonde data at Sodankylä and showed that the characteristic scale of the ozone fine structure is ∼1.0 km in the troposphere and ∼1.4 km in the lower stratosphere (up to 25 km). It was also found that the characteristic scale of ozone fine structure is a relatively stable atmospheric characteristic without significant seasonal or interannual variations at this location.

MAP method with smoothness a priori
In most cases a priori information includes only a general measure of smoothness of atmospheric profiles. For such situations, we develop a method recently proposed by Haario et al. (2004). It is based on the assumption that the neighboring discretized values of a retrieved profile cannot be too different. This a priori can be presented as are mutually independent Gaussian random variables with the zero mean and the diagonal covariance matrix C reg and h i is the discretization grid. Alternatively, first order or higher order differences can be used in the left-hand side of Eq. (8). The regularization Equation (8) can be expressed in the matrix form as where tri-diagonal matrix H approximates second derivatives: It corresponds to the prior distribution (smoothness a priori) (P prior ∝ exp(− 1 2 ρ T H T C −1 reg Hρ). If the matrix H T C −1 reg H is invertible, the prior distribution is Gaussian: ρ a ∼N (0, (H T C −1 reg H) −1 ). The MAP estimation ρ sm can be written as The only information needed for application of this method are the uncertainties of the second differences C reg . It can be also created from analysis of high-resolution profile measurements, such as ozone sonde data. Analogously to the classical MAP estimation, this method efficiently combines the measurements and a priori information, applying additional smoothing only when it is required by a low signal-to-noise ratio. It is also grid-independent, i.e. the chosen amount of smoothing is defined by the ratio of a measurement noise and a priori uncertainty, but not by the discretization grid.
3.2 Tikhonov-type regularization. Choosing regularization parameter according to the target resolution The classical Tikhonov regularized solution of the problem (1) was originally derived as a minimizer of the functional Here λ is the regularization parameter and H is the matrix representing first, second (10) or higher order differences (which are assumed to be bounded, thus characterizing the smoothness of the solution), and · is 2 -norm. The Tikhonov regularized solution of (12) exists, and it is unique. It is given by the formulâ It is equivalent to the MAP solution (6), provided that the prior distribution is P prior ∝ exp(− 1 2 ρ T H T Hρ) and the measurement error ε is Gaussian ε∼N (0, σ 2 I), where I is the unit matrix.
The optimal choice of the regularization parameter λ is a central issue in the literature discussing the Tikhonov regularization. It can be chosen, for example, according to the Morozov's discrepancy principle (e.g. Morozov, 1993;Hansen et al., 2000), which states that λ should be selected so that the residual (difference between model and measurement data) should be of the same value as the measurement noise: Application of the Tikhonov regularization, as well as the MAP methods discussed above, leads to a certain degradation of resolution. The regularization parameter can also be chosen according to some target resolution, if the optimal value of the regularization parameter does not meet the resolution requirements. Then the actual vertical resolution does not depend on the instrumental noise and the discretization grid. This simplicity of profile characterization makes this method attractive, despite some disadvantages discussed further in Sect. 4.

GOMOS measurements: comparison of regularization schemes
In this section we apply methods discussed above to the problem of reconstruction of ozone profiles from the GOMOS measurements. The following methods are used: 1. Onion peeling: this is the standard GOMOS inversion without any a priori. In discretization, a constant density inside each layer is assumed. Onion peeling can be formulated as a matrix inversion or, equivalently, the problem can be solved sequentially, starting from the uppermost layer. We used the matrix inversion in our simulation.
2. Tikhonov regularization in grid-dependent formulation with regularization parameter λ=10 15 , which is chosen according to the Morozov's discrepancy principle for the typical GOMOS star of visual magnitude m=2 (The horizontal column densities are assumed to be in 1/cm 2 and distances needed for computing K and H are in cm.).
3. Tikhonov-type regularization in grid-independent formulation and altitude-dependent regularization parameter λ. The regularization parameter is chosen according to the target resolution (Table 1), that corresponds to the characteristic scale of ozone fine structure determined from Sodankylä ozone sonde data below 30 km. The values of the target resolution above 30 km are taken from climatology.
4. MAP method with the smoothness a priori described in Sect. 3.1. The covariance of second order finite differences was created from the Sodankylä ozone sonde data.
5. MAP method (classical). The Fortuin-Kelder climatology (Fortuin and Kelder, 1998) was taken as a priori mean ozone profile below 53 km. The upper part of the prior mean ozone profile was simulated with the LIMBO simulator (Kyrölä et al., 1999). A priori standard deviation was calculated from 10 years of ozone sonde data (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) at Sodankylä. The correlation scale is taken in accordance with   (Table 1). In the upper atmosphere the prior standard deviation is assumed to be a linear growing function from 20% to 200% for altitudes 30-100 km with the correlation scale 3 km.

Simulation setup
The realistic simulation was performed for the comparison of the methods. The measurements grid was taken coincident with a real GOMOS measurement grid, with a typical angle between the orbital plane and the star direction (∼18 • ). The simulations were performed for stars of different magnitudes and of the typical effective temperature 10 000 K. Both instrumental noise and modelling errors (Kyrölä et al., 1993;Sofieva and Kyrölä, 2003) were taken into account in computing the measurement error (i.e. error of horizontal column densities). The "true" profiles used in the simulation coincide in the lower atmosphere (below 30 km) with ozone profiles measured at Sodankylä (http://fmiarc.fmi.fi). The upper part of ozone profile was simulated with the LIMBO simulator and perturbed in accordance with the "error pattern method" (Rodgers, 2000). Figure 5 shows an example of reconstructed ozone profile in the middle stratosphere and troposphere for the star of magnitude m=2 and effective temperature T =10 000 K. We will not show, discuss or compare the inversion at high altitudes, because no a priori information is available there. All the advanced methods give almost the same results in the mesosphere, as they use the same ad hoc assumptions on profile smoothness. The lower part of the "true" profile used in the simulation coincides with the ozone profile measured at Sodankylä on 23 November 2002.

An exemplary ozone profile retrieval
The ozone profile reconstructed with the onion peeling method, which gives the results almost identical to the standard GOMOS vertical inversion, looks very noisy below 15 km. The Tikhonov-type regularization, with a parameter chosen according to the discrepancy principle, smoothes out fine structures of ozone profile. If the regularization parameter is chosen according to the target resolution, the fine structure of the ozone profile is well detected in the altitude range 15-30 km. For altitudes below 10 km, both Tikhonov-type regularizations give the same results and they somewhat improve the reconstruction, but the retrieved profiles are still far away from the real profile. The ozone profile is reconstructed most accurately with maximum a posteriori methods (classical and with smoothness a priori). However, the fine structure of the ozone profile may stay unresolved with these methods. In the considered example, the lamina around 13 km is not detected with the MAP methods either.

Error analysis
The efficiencies of the discussed methods characterized by estimated error of retrieval and resolution are compared in Fig. 6. The error estimates in Fig. 6 are square roots of diagonal elements of the posterior covariance matrix (Given the retrieved profileρ in the form ofρ=GN , where G is the inversion matrix, the posterior covariance matrix is C post =GC ε G T .). The prior uncertainty of the ozone profile is also shown in this figure. The original vertical resolution of the GOMOS measurements (it corresponds to onion peeling method) is very good: ∼2 km for altitudes above 40 km and better than 1 km below 20 km. The accuracy of the reconstruction can be significantly improved by advanced methods at a price of degraded vertical resolution.
However, the error estimate, as presented in Fig. 6, is not the total error of ozone profile reconstruction. The total error (i.e. a deviation of a retrieved profileρ from the true one ρ) iŝ Here A is the averaging kernel matrix, I is the unit matrix and G is the inversion matrix. The latter term in Eq. (15) is the error due to noise in measurements, while the former one describes the smoothing error caused by deviation of averaging kernels from delta-functions. The covariance of the total error Eq. (15) is (Rodgers, 2000) where C e is the covariance of an ensemble of real profiles about the mean profile. For the occultation located near Sodankylä C e can be obtained from the ozone sonde measurements, thus enabling one to estimate the total error of the profile reconstruction and compare the real efficiency of different inversion algorithms. In our work, we used 10 years of sonde data (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999), in order to estimate the ozone variability (i.e. to obtain the ensemble covariance C e ). Figure 7 shows the total rms error of reconstruction by different inversion algorithms for stars of various magnitudes.
The following conclusions can be drawn from this analysis: -The standard Tikhonov regularization with one constant smoothing parameter chosen according to the discrepancy principle has the largest error for altitudes above 15 km. This method seems not to be applicable to ozone profile reconstruction, because the chosen smoothing violates the resolution requirements, thus leading to the largest total error.
-For the other four methods, the error of reconstruction is almost the same for occultations of bright stars at altitudes >15 km. This confirms that it is not necessary to apply advanced methods in these conditions. -For altitudes below 15 km the accuracy can be significantly improved by advanced inverse methods. The classical MAP method has the smallest total error, as theoretically predicted. The MAP with smoothness a priori is almost as good as the classical MAP method. The target resolution method also improves the inversion at the lower atmosphere, especially for dim stars.
In the Tikhonov regularization methods, the vertical resolution is predetermined, while in the MAP methods, it depends on the noise level (Fig. 8). The application of the MAP methods regularizes the inversion only when data are seriously contaminated with the noise, while using the target resolution method, some smoothing is always applied, even if it is not necessary. We have to recall that the classical MAP method cannot be applied at the moment, because the variability of ozone profile is not known, except for locations coincident with ozone sounding stations. This method was considered in order to illustrate its potential applicability in the future. The features of the inversion methods, their advantages and drawbacks are collected in Table 2. -A priori is not available globally.
MAP method (classical) -mean profile -The inversion uses a -Actual resolution depends on -standard deviation priori information most stellar spectrum. -characteristic completely.
-A priori is not available globally. scale -Solution is statistically optimal.

Summary and conclusion
In this paper we have discussed how the prior information about smoothness of the atmospheric profiles can be used in the profile retrievals. It can be applied either in the form of the smoothness constraint (Tikhonov-type regularization) or in the form of a priori information (maximum a posteriori methods). The methods can be used when the classical Bayesian MAP approach is not applicable due to insufficient information about prior distribution. The grid-independent formulations were proposed for all the considered methods, thus providing results independent of a discretization grid.
In the GOMOS measurements, the low signal-to-noise ratio for dim stars prompted the adoption of the Bayesian approach, which allows one to efficiently combine measured data with a priori information about the smoothness of the atmospheric profiles.
These methods were applied to reconstruction of the ozone profile from GOMOS measurements. The realistic simulations for typical measurement conditions and analysis of the total error of reconstruction revealed advantages and drawbacks of each proposed procedure.
The following main conclusions can be drawn: 1. The standard Tikhonov regularization with smoothing parameter chosen according to discrepancy principle is not recommented for the ozone profile retrieval from the GOMOS measurements, because the "optimal" smoothing violates the resolution requirements: almost all fine structures of ozone profiles are smoothed by this method.
2. The regularization with the choice of the smoothing parameter according to the resolution requirements is the most attractive method because of the predetermined vertical resolution, independence from the measurement grid and from stellar properties. It is a kind of "minimal guaranteed strategy", but it is not optimal in the statistical sense.
3. MAP with a smoothness a priori allows one to include the information about smoothness of the atmospheric profile in the form of Bayesian optimal estimator. It efficiently combines the measurement data and information on smoothness of ozone profiles. It is a good alternative to the classical MAP method, when the mean and/or standard deviation of retrieved quantities are not known.
4. The successful application of the classical MAP method was demonstrated. The a priori information about the mean and the standard deviation of the ozone profile was created from the Sodankylä ozone sonde data. This method includes the a priori information in the most complete form. Probably, this method will be used widely in the future algorithms of ozone profile reconstruction, but now it can be applied only in rare cases.
In summary, the proposed inversion methods, which use a priori information on smoothness of atmospheric profiles, improve significantly the quality of retrievals. These methods can be used in reconstruction of atmospheric profiles not only from stellar occultation but also from other remote sensing measurements.