## Abstract

We present a general analysis of the error budget in the spectral inversion of atmospheric radiometric measurements. By focussing on the case of an occultation experiment, we simplify the problem through a reduced number of absorbers in a linearized formalism. However, our analysis is quite general and applies to many other situations. For a spectrometer having an infinite spectral resolution, we discuss the origin of systematic and random errors. In particular, the difficult case of aerosols is investigated and several inversion techniques are compared. We underline the importance of carefully simulating the spectral inversion as a function of the target constituent to be retrieved, and the required accuracy level.

© Optical Society of America

## 1 Introduction

In remote sensing experiments, atmospheric constituents are retrieved from measurements of a spatially integrated optical property. Remote probing can be carried out either in a passive or in an active way and using many spectral regions. From the near-ultraviolet (UV) to the near-infrared (IR) domain of the spectrum, the total column amount of a specific absorbing gas can be inferred from ground-based measurements combined with the Langley method. On the other hand, space-borne instruments using the occultation technique offer three major advantages: they possess a natural self-calibration property, they allow an accurate retrieval of vertical profiles of number density, with a possible quasi-global coverage.

However, all these techniques face a common problem: measured vertical or slant-path optical thicknesses generally result, at any wavelength, from additive light extinctions (absorption or scattering) by different atmospheric constituents. The optical thickness separation into its different contributions is often called the spectral inversion as opposed to the vertical inversion where the optical thickness is decomposed into the contribution of individual atmospheric layers.

In the literature, many methods have been proposed [1],[2],[3],[4] to solve this problem in a pragmatic way. Although these retrieval schemes turned out to be reasonably successful, they are often obscure in their genesis and do not say very much about the fundamental limitations of the measurement. Even more, there are few assertions about the possible superiority of a particular method and the error budget is mostly discussed only with respect to random errors.

A major reason for this situation (at least in the stratosphere), arises from the presence of aerosols which can contribute to scattering in the near-UV to near-IR range according to the Mie theory. A first difficulty comes from their variability; an increase factor of 100 may be observed between background aerosols and the huge loading produced by major volcanic eruptions such as Mount Pinatubo in 1991. Another difficulty lies in the a priori unknown wavelength dependence of the aerosol optical thickness. This is mainly due to the fact that the aerosol extinction is a macroscopic quantity related to an (unknown) underlying particle size distribution. In contrast with the other absorbing gases whose optical thicknesses at different wavelengths are related by their absorption cross sections, the aerosol wavelength dependence is a product of the spectral inversion and is not known in advance.

The purpose of the present paper is to discuss the problems associated with the spectral inversion of an occultation experiment. For instance, the GOMOS experiment [5](on board Envisat) is a clear example that illustrates the difficulties encountered although the analysis can also be easily extended to other cases, including ground-based experiments or even reflected solar radiation in the Earth’s atmosphere. Therefore we will first reduce the problem to a simplified formalism by means of reasonable approximations. All essential features that will be explained hereafter are necessarily present in more sophisticated and (possibly) exact inversion schemes. In section 2, we will discuss the issue of bias and random errors and we will suggest a criterion for the optimal description of the aerosol optical thickness. Finally, in section 3, we will extend the analysis to the problem of microwindows and filtering. No definitive optimal solution will be proposed to the reader. Instead, we will underline the importance of carefully simulating the spectral inversion as a function of the target constituent to be retrieved, and the required accuracy level.

## 2 Spectral inversion for an ideal occultation radiometer

According to Beer-Lambert law, the atmospheric transmittance measured by an ideal monochromatic occultation radiometer working at wavelength λ can be written:

where *h* is the tangent altitude of the measured ray of light and *τ*(λ) is the slant path optical thickness due to absorption or scattering. It is well known that refractive effects of the Earth’s atmosphere have to be taken into account below altitudes of about 40 km. In particular, the diverging atmospheric lens induces a dilution of the incoming flux (for light sources smaller than the instrument field of view) and a shift in the tangent altitude. The optical path length increase due to ray bending is however very small and, as most of the extinction occurs in the neighbourhood of the tangent point, we can approximate the slant path optical thickness by taking the logarithm of the dilution corrected transmittance at the physical (refracted) tangent altitude.

If we restrict our analysis to the near UV - near IR range (λ ∈ [0.2,1.0] *μ*m), the true optical thickness (unknown to the experimenter) may be described by:

where *τ*
_{N2}, *τ*
_{O3}, *τ*
_{NO2} respectively stand for slant path optical thicknesses of air, ozone and nitrogen dioxide. For a standard atmosphere, *a*
_{N2}, *a*
_{O3} and *a*
_{NO2} are set to 1. The individual optical thicknesses will not be translated into trace gas columns and it should be kept in mind that this step could introduce supplementary uncertainties if the absorption or scattering cross sections vary along the line-of-sight. *τ*
_{A}(λ) is the aerosol optical thickness of which the wavelength dependence is always smooth, due to integration of the Mie cross section over the particle size distribution. In order to simplify the discussion hereafter, we have deliberately removed contributions from water vapour, A and B oxygen bands as well as other minor constituents. In Fig. 1, we have represented the partial slant path optical thicknesses for a tangent altitude of 20 km and from climatological values (*a*
_{N2} = *a*
_{O3} = *a*
_{NO2} = 1). The λ^{-1} aerosol wavelength dependence of *τ*
_{A}(λ) in a moderate volcanic situation [6] is considered as our standard case. It is clear that ozone is characterized by two strong spectral signatures in the Huggins (λ < 340 nm) and Chappuis (λ ≃ 600 nm) bands, and this constituent should normally be the easiest to retrieve. Nitrogen dioxide is not very absorbing but possesses a useful differential structure superimposed on a broad continuum centered at about λ ≃ 400 nm. Finally, aerosol and Rayleigh scattering (∝ λ^{-4}) are clearly mutual competitors both in magnitude and in spectral shape. Considering that the natural smoothness of aerosol optical thickness suggests its representation by a low order (*n*) polynomial, the spectral inversion consists of retrieving the partial optical thicknesses by using a model *R*(λ) defined by (we remove the variable h for the sake of readiness):

where *τ*_{i}
(λ) ∝ (λ - λ_{0})^{i} and λ_{0} is an arbitrary wavelength reference (the *x*’s are the retrieved parameters). The experimental data will always be measured at a defined noise level depending on the light source (star, planet or Sun) and the detector intrinsic sensitivity. For the sake of simplicity we will only consider here shot noise, which is proportional to the square root of the signal for a Poisson distribution. If we call *S* the detector sensitivity (the dynamic range expressed in photons or electrons per wavelength unit), we may write for the error Δ*T* on transmittance (Eq. 1):

which leads, by first order propagation of errors, to

Here we fully agree with [1] that the true a priori distribution of the model parameters is only gaussian in a first approximation, which may cause difficulties in low absorption (when noise makes *T* occasionally higher than 1) and saturation regimes (with possible negative *T* and troubles in computing *τ*(λ)). These are regions where the information content of the measurement is questionable and problems can be partially cured by appropriate screening and filtering. We will assume that Eq. 5 is valid during the complete occultation even if more accurate approximations could be used.

It is useful to consider the problem for an infinite resolution spectrometer because it is the limiting case of modern instruments where spectral resolutions of 1 nanometer can be achieved nowadays. Therefore the matching between the modelled and measured optical thicknesses will be evaluated by means of the following merit function *M*:

with λ_{1}=0.2 *μ*m and λ_{2} =1 *μ*m as default values. The minimization of *M* is obtained by imposing

which leads to the linear system

where *x*⃗is the column vector of unknowns and *C* is the design or covariance matrix:

whose elements are the inner products:

The column vector *y*⃗ contains the respective projections of the experimental *τ*(λ) onto the basis functions *τ*_{i}
(λ) as

*C* being a real symmetric matrix, it can be diagonalized by using an orthogonal transformation *U*:

where Λ is the diagonal matrix of eigenvalues and the solution of the spectral inversion reads

whereas the covariance matrix of the solution is simply *C*
^{-1}. Therefore, the expected random error *e*_{r}
on N_{2},O_{3},NO_{2} is:

The expected random error error on the aerosol optical thickness is

The relative random errors are defined by:

and

It is worth discussing the previous formulas. It is clear that the random error of the solution is driven by the *C*
^{-1} matrix and hence by the smallest eigenvalues of the covariance matrix *C*. A first consequence is that any aerosol that can be expressed as a quasi-linear combination of the spectral dependence of the other constituents will makes *C* almost singular and will induce huge random errors in the solution. A trivial example could be a λ^{-4} dependence that would make aerosol and Rayleigh scattering physically undistinguishable by the optical experiment. A second consequence is that increasing the polynomial degree *n* necessarily implies more and more small eigenvalues in the spectrum of *C* resulting in enhanced random error. On the other hand, increasing the spectral range [λ_{1}, λ_{2}] or the detector sensitivity *S* will increase the elements of *C* and ameliorate the random error.

Apart from random errors, the spectral inversion may also suffer from systematic or bias errors. Indeed, if the true unknown aerosol is not completely described by a polynomial expansion, there will also subsist an incompleteness error which cannot be cured by a statistics of a large number of inversions as it is possible for random errors. If we call Δ*τ*
_{A} the difference between the true aerosol and a pure polynomial expansion, this difference will project onto all the basis functions, producing a *δy*⃗

which, in turn, will induce a positive or negative bias *δx*⃗

It is clear that all elements have, a priori, non-zero values. Furthermore, depending on the ill-conditioning of *C* (small eigenvalues), these elements will be amplified and mixed. This means that the aerosol error will cause a bias error *e*_{b}
for all retrieved constituents and this kind of bias has been often disregarded in the literature. For N_{2}, O_{3} and NO_{2}, the relative bias errors are

and for aerosols, we have

and

In particular, this can lead to the naive idea that increasing the spectral range may surely ameliorate the spectral inversion. Mostly, the contrary will occur because the elements of *δy*⃗ will generally increase. In Fig. 2, we have computed the systematic error produced by different ”true” aerosols generated from

which reduces to the default aerosol for *γ* = 0. Clearly, when *γ* is varied, within the range of realistic values corresponding to already observed volcanic modes, the systematic error of all constituents evolves considerably. Notice that the ozone errors are always quite small due to the pronounced spectral signature of the ozone cross section.

A possible way to decrease Δ*τ*
_{A} is to increase the polynomial degree *n* but this will increase the random error. Hence, we are led to the conclusion that there must exist an optimal *n* that realizes a trade-off between the random and bias errors. This is illustrated in Fig. 3 where this optimal value is found to be different from constituent to constituent. An important consequence of this analysis is that a global inversion of all constituents together cannot be optimal and that the trade-off should depend on the geophysical objectives (e. g. an optimal accuracy for independent ozone retrievals). This has not been noticed by [1] because they performed their spectral inversion by assuming that the aerosol wavelength dependence is known a priori, which is not the case in reality.

## 3 Exploring alternative solutions to global spectral inversion

#### 3.1 Do microwindows help?

Looking at Eqn. 6, it is worth investigating the role played by Δ*τ*(λ) which acts as a spectral range cut-off when transmittance gets too weak, avoiding oversensitivity to the UV part. In view of the above-mentioned influence of the constituent non-orthogonality, it is interesting to consider the possibility of selecting only a subrange of pixels to perform the optimal inversion of a particular constituent. A simple way to investigate this possibility is the inclusion of a filter function *F*(λ) into the merit function as

and, for instance, *F* can be given a simple gaussian dependence in order to simulate a single spectral window:

where *c*
_{1} and *c*
_{2} respectively play the role of center and width of the window. In Fig. 4, we present the isopleths of bias and random errors for the four retrieved constituents when both parameters are varied. As expected, the random error (right column) increases when a smaller window is used (low *c*
_{2} values) but not uniformly: air and nitrogen dioxide are better retrieved toward the UV range while ozone and aerosol are favored by the center of the spectral range. Interesting patterns emerge from the bias error in the left column. Air retrieval is more effective in the UV-part in opposition to aerosol for which no clear improvement can be obtained with respect to full scale window (*c*
_{2}=∞). Ozone and nitrogen dioxide exhibit very structured domains. In particular, isolines of zero error bias exist that can remove possible systematic errors caused by the aerosol incomplete description. Indeed, if we respectively call *e*_{r}
(*c*
_{1},*C*
_{2}) and *e*_{b}
(*c*
_{1},*C*
_{2}) the random and bias error as function of the spectral window, a natural strategy could be to minimize *∊*_{r}
under the constraint *∊*_{b}
= 0. This can be easily accomplished by constructing the Lagrangian function *L* as

and numerically solving the nonlinear system of equations

For example, one gets for NO_{2} a solution at the location (${c}_{1}^{*}$, ${c}_{2}^{*}$) = (0.450,0.148) where (${\mathit{\u220a}}_{b}^{*}$,${\mathit{\u220a}}_{r}^{*}$) = (0,49.2) [%] instead of (12.5,30.2) [%] when *c*
_{2}=∞. Although the random error is somewhat larger, the benefit of having a zero bias is evident when data will be averaged or binned in a set of different measurements because the random error will also decrease by a factor equal to the square root of the number of measurements. For only one measurement however, the total error ${\u220a}_{t}=\sqrt{{\u220a}_{b}^{2}+{\u220a}_{r}^{2}}$ may be the most important quantity to consider. It is worth mentioning that the parameters of the optimal spectral windows were found to be quasi-independent of the aerosol spectral shape (the value of *γ* in Eqn. 23) and of the polynomial order *n*. The same procedure can be applied for an optimal ozone retrieval.

In order to ameliorate the aerosol retrieval, we have also implemented a multi-window retrieval, where *F*(λ) is no longer a simple Gaussian but is allowed to take multiple maxima over the whole spectral range and fully optimized by simulated annealing. For *n* = 2, a weak amelioration is observed with respect to the Gaussian window, and the optimal windows are mainly centered at locations where some constituents have the same order of magnitude (for instance, there is a maximum around λ = 0.5*μm* where Rayleigh, aerosol and ozone extinctions have similar contributions to the optical thickness). When using a *n* = 3 polynomial, the optimal window is the full window (*F*(λ) = 1) and no more progress can be expected from this approach.

#### 3.2 Regularization

A possible solution to increase the accuracy of the inversion, with simultaneous control of the sensitivity to random errors, could be achieved by a constrained linear inversion [7]. In this regularization technique, we construct a more elaborate merit function *M*
_{1} from Eqn. 24 as:

where the second term of the left-hand side constrains the aerosol optical thickness to be smooth whatever the polynomial order *n* is. The latter equation consists of the sum of two quadratic forms and leads to a modified linear system similar to Eqn. 8. The relative contribution of the constraint with respect to *M* is controlled by *ρ*, whose optimal value can be determined by generalized cross validation [8]. However, it is interesting to investigate the outcome of the constrained inversion for a large range of *ρ* values. The results are presented in Fig. 5 for aerosol bias, random and total errors when *n* is varied from 2 to 7. Clearly, when *n* becomes larger, there is a competition between smoothing (a factor of random error reduction) that increases with *ρ* and accuracy (that reduces bias error) when *ρ* is very small. The total error shows a clear minimum when *n* ≥ 4 but the important point is that the value of the total error at this minimum is not better than for *n* = 2 or *n* = 3. We deduce that constrained linear inversion does not improve the aerosol retrieval.

An other approach is to consider the possibility of exploiting the low and high frequency components of the measured optical thickness, i.e. we construct the merit function *M*
_{2} as:

where the symbol prime refers to the first derivative (in practice, the differentiation is performed numerically) with respect to wavelength and *χ* is an adjustable parameter in the range [0,1]. Inversion results are presented in Fig. 6. As expected, the use of the higher frequencies (*χ* ≃ 1) increases accuracy and noise sensitivity. Notice, however, the drastic reduction of the bias error on the NO_{2} retrieval which is the basis of the differential absorption spectroscopy (DOAS) method[9].

#### 3.3 Filtering

On the other hand, we know that the total atmospheric optical thickness is a relatively smooth function of wavelength on which measurement noise is superimposed.

Although it is generally not recommended to pre-filter before data fitting , we have to consider that inverse problems are extremely sensitive to high frequency noise which suggests the use of low-pass filters capable of preserving aerosol information. Amongst innumerable methods, we have applied the Savitzky-Golay (hereafter referred to as SG) filtering technique ([10]) which is known to preserve moments of the underlying signal ([11]). Briefly, the SG filter is characterized by *n*_{sg}
, the number of neighbouring points that determine the filtered value (discretized with Δλ = 0.001 *μ*m), and *m*_{sg}
, the associated polynomial order that locally fits the data in the sliding window. It is important to apply the filter on the total optical thickness already weighted by *T*(λ)*S* as in Eqn. 10. In Fig. 7, we present the isopleths of the total error *∊*_{t}
(*n*_{sg}
,*m*_{sg}
) when using the *n* = 2 and *n* = 3 aerosol polynomial order representations. A clear minimum appears in both cases, for *m*_{sg}
≃ 2 – 3 and *n*_{sg}
≃ 51. The amelioration is not negligible. For *n* = 2, the total error is decreased from *e*_{t}
≃ 13.5% to *e*_{t}
≃ 8.5%, and, surprisingly, is smaller than for *n* = 3. This is due to the unavoidable action of the filter on the true signal itself and not only on the noise.

However, it is important to realize that, until now, we have investigated the inversion error budget by using the true aerosol optical thickness and comparing it to the retrieved value. In the real world, it would then be necessary to have a good approximation of the solution to estimate, via a table look-up procedure, the optimal filter parameters.

## 4 Conclusions

This work has proven that the optimal retrieval of a particular atmospheric constituent may require a specific strategy instead of a global inversion.

For aerosols, we conclude the above-mentioned approaches by the following consideration: the retrieval faces an apparently untractable problem: it is necessary to increase *n* in order to diminish the bias error but in doing so we always increase the noise sensitivity. In the near-UV to near-IR wavelength range that we have considered here, the number of aerosol parameters that can be retrieved from radiometric measurements can hardly exceed 4 or 5, depending on the signal to noise ratio of the spectrometer. This should be kept in mind when trying to retrieve bimodal particle size distributions (requiring at least 6 parameters) from the wavelength dependence of the aerosol optical thickness.

## Acknowledgements

This work was partly performed within projects ”SADE (Prodex 6)” and “Measurement, understanding and climatology of stratospheric aerosols (grant MO/35/004)” funded by the SSTC/DWTC service of the Belgian Government.

## References and links

**1. **E. Kyrola, E. Shivola, Y. Kotivuori, M. Tikka, and T. Tuomi, “Inverse Theory for OccultationMeasurements: 1. Spectral Inversion,” J. Geophys. Res. **98**, 7367–7381 (1993). [CrossRef]

**2. **D. E. Flittner, B. M. Herman, K. J. Thome, J. M. Simpson, and J. A. Reagan, “Total Ozone and Aerosol Optical Depths Inferred from Radiometric Measurements in the Chappuis Absorption Band,” J. Atmos. Sci. **50**, 1113–1121 (1993). [CrossRef]

**3. **M. King, “Sensitivity of Constrained Linear Inversions to the Selection of the Lagrange Multiplier,” J. Atmos. Sci. **39**, 1356–1369 (1982). [CrossRef]

**4. **W. P. Chu, M. P. McCormick, J. Lenoble, C. Brogniez, and P. Pruvost, “SAGE II InversionAlgorithm,” J. Geophys. Res. **94**, 8839–8351 (1989). [CrossRef]

**5. **J. L. Bertaux, G. Megie, T. Widemann, E. Chassefiere, R. Pellinen, E. Kyrola, S. Korpela, and P. Simon, “Monitoring of ozone trend by stellar occultations: the GOMOS instrument,” Advances in Space Research **11**, 237–242 (1991). [CrossRef]

**6. **D. Fussen, F. Vanhellemont, and C. Bingen, “Evolution of stratospheric aerosols in the post-Pinatubo period measured by the occultation radiometer experiment ORA,” Atmos. Env. **35**, 5067–5078 (2001). [CrossRef]

**7. **S. Twomey, “Comparison of Constrained Linear Inversion and an Iterative Nonlinear Algorithm Applied to the Indirect Estimation of Particule Size Distributions,” J. Comput. Phys. **18**, 188–200(1975). [CrossRef]

**8. **G. H. Golub and C. F. Van Loan, “Matrix Computations,” (The Johns Hopkins University Press1996).

**9. **U. Platt, “Air monitoring by Spectroscopic Techniques, Chapter 2,” (John Wiley and Sons 1994).

**10. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Numerical Recipes in FORTRAN, Second Edition,” (Cambridge University Press, Cambridge 1992).

**11. **M. U. Bromba and H. Ziegler, “Applications Hints for Savitzky-Golay Smoothing Filters,” Analytical Chemistry **53**, 1583–1586 (1981). [CrossRef]