## Abstract

It is well known that optical satellite remote sensing is mostly based on measurements of the radiance exiting the top of the earth's atmosphere and is interpreted using a comparison of the data with precalculated tables for some “pure” aerosol models. However, such an approach seems to meet some difficulties connected to both the necessity to search in the multidimensional parameter space and to increase the volume of the precalculated database, if the number of basic aerosol components increases. This problem becomes more real now because modern satellite sensors [multiangle imaging spectroradiometer (MISR), Polarization and Directionality of the Earth's Reflectances (POLDER)] provide more information by performing multiangle radiance measurements over the same pixel and more detailed characteristics of the atmospheric aerosol profile (not only the aerosol optical depth and Ångström coefficient) are expected to be retrieved.

Recently, it has been shown that the perturbation technique was an ideal tool to solve such atmospheric physics problems as an investigation of the effect of a variation in aerosol profile on fluxes. Moreover, it has shown its efficiency in analysis of the optical remote sensing in the simplest case of the sounding of a homogeneous medium.

In the present paper the perturbation technique is applied to a realistic atmosphere–ocean model. Simulations were performed for a stratified slab medium with an underlying surface. Perturbations considered include variations of the profiles of extinction coefficient and single scattering albedo, and also the adding of a small cloud contamination. It was shown that the perturbation technique allows one to predict the effect of such variations of the atmospheric profile to the exiting radiance with high accuracy. The nature of the errors is analyzed and discussed.

## 1. Introduction

It is well known that the atmospheric aerosol exerts a significant impact on the earth's climate, both directly (d'Almeida et al. 1991), by absorbing and scattering radiation, and indirectly (Charlson et al. 1992; Curry 1995) by their effect on cloud optical and microphysical properties. To better understand these effects it is necessary to know not only the aerosol optical depth, but also the aerosol scattering and absorption properties, and their vertical profiles. In other words, an adequate model of global aerosol distribution, taking into account their seasonal and spatial variability, has to be developed.

The integrated field experiments being conducted [e.g., ASTEX/MAGE in June 1992 (Huebert et al. 1996), SCAR-B in August–September 1995 (Kaufman et al. 1998), ACE-1 in November–December 1995 (Bates et al. 1988), TARFOX in July 1996 (Russell et al. 1999)] can specify many features of local aerosol atmosphere models and their results are the best touchstone for any kind of retrieval algorithms, especially for satellite remote sensing.

The only approach that could provide the necessary spatial and temporal coverage is earth observation from space. During the last 5 years several sensors to fulfill this task have been launched. Moreover, some of these modern sensors provide much more information than those of previous generations, by viewing the same surface area under different angles [POLDER on board the *Advanced Earth Observing Satellite* (*ADEOS*, 1996) has up to 14 angles (Deschamps et al., 1994), MISR on the *Earth Observing System* (*EOS*, 1999) has 9 angles (Diner et al. 1991)]. Substantial increasing of the available information makes it possible to retrieve not only the aerosol optical depth, but also aerosol stratification, detection of absorbing aerosol (over the oceans), or estimation of the bidirectional reflection function of the surface (over the land) (Kaufman et al. 1997). However, this leads to the retrieval algorithms becoming more sophisticated and causes the number of base aerosol models to increase. Nevertheless, the actual number of aerosol models that are really taken into account is usually limited to constrain the size of the precalculated database (Martonchik et al. 1998). This could be a possible source of retrieval error, if nonstandard aerosol contributions occur, or in the case of an undetected small cloud contamination of a pixel, but no thorough sensitivity study has yet been performed because of the multiple-factors nature of the problem. However, the papers by Mishchenko and Travis (Mishchenko and Travis 1997a,b) should be noted in which, despite of the simplest aerosol–ocean model being considered, the important conclusion that the algorithm based on multiple viewing angle measurements performs far better than that based on single viewing angle radiance measurements is made.

Usually, to simulate satellite remote sensing the approach based on developing of precalculated database (Rao et al. 1989; Mishchenko and Travis 1997a; Martonchik et al. 1998) is used. In this paper a fundamentally different technique is applied to the problem of multiangle satellite observation of a stratified atmosphere. This approach, known as the perturbation technique, has been adapted from nuclear reactor theory (Marchuk 1964; Gerstl 1980), and developed comprehensively for the purpose of the radiative transfer in atmosphere in the papers of Box and his coauthors (Box et al. 1988a,b, 1989a,b; Trautmann et al. 1992). This approach shows its efficiency to analyze the possibility of the retrieval of the optical thickness, single scattering albedo, and phase function (Sendra and Box 2000) in the case of the sounding of a homogeneous medium. In this paper it will be applied to the sounding of a stratified medium. The perturbation approach also seems to be a perfect base to develop the retrieval algorithm, because it provides the necessary derivatives and hence the standard techniques of the optimization of multivariable nonlinear functions (Fletcher 1987) may be successfully implemented.

## 2. Perturbation technique

The perturbation technique of radiative transfer theory is based on the joint solution of both the direct and adjoint equations, and has been used for more than 50 years. Initially it was successfully applied to nuclear reactor problems (Lewins 1965; Bell and Glasstone 1970). Marchuk (1964) was the first who has introduced this approach as a general technique to interpret optical measurements into the field of the atmospheric optics. However, the most successful applications of the perturbation technique started more then a decade later. This approach was applied to compute atmospheric radiative effects (Gerstl 1980; Box et al. 1988a, 1989a,b; Trautmann et al. 1992), to consider the problems of thermal sounding (Ustinov 1990) and the photometric observation of solar radiation reflected from the optically thick vertically inhomogeneous planet atmosphere (Ustinov 1991a,b, 1992), and to investigate the sensitivity of exiting radiances to optical characteristics in the simplest case of a homogeneous atmosphere (Box and Sendra 1995).

Let us introduce briefly the basic results. We consider only the case of radiation transfer in a plane parallel horizontally homogeneous atmosphere with solar illumination. In the most common cases, the radiative effects (e.g., fluxes) or optical measurements are given by

where *I*(*z,* Ω) is the radiance at altitude *z* and in the direction Ω, *R*(*z,* Ω) is the response function that to simulate satellite measurements has the form

where *μ* is the cosine of the zenith angle, *δ*(*z*) is the Dirac *δ* function, *z*_{t} is the altitude of the atmosphere “top,” *μ*_{R} and *ϕ*_{R} are the zenith angle cosine and azimuth angle of the satellite receiver.

The *I*(*z,* Ω) satisfies the radiative transfer equation (RTE)

where ɛ(*z*) and *σ*(*z*) are the extinction and scattering coefficients, respectively, and *Q*(*z,* Ω) is the source function, which has the form

where *μ*_{0} and *ϕ*_{0} are the zenith angle cosine and azimuth angle of the sun, and *χ*(Ω, Ω′) is the phase function, normalized by

For the sake of simplicity it is convenient to rewrite Eq. (3) in the operator form

where

The notation ⊗ is used to indicate that the final term is an integral operator, not a simple definite integral.

Let us introduce an adjoint operator *L̂*^{+} , which is defined by requiring that

where 〈…〉 = ∫∫ (…) *dz **d*Ω denotes a scalar product and *I*^{+} is the solution of the adjoint equation

Taking into account the common boundary conditions for the radiance (Lenoble 1985, p. 179), that is,

where *ρ*(Ω, Ω′) is the bidirectional reflection function of the underlying surface and the sunlight is taken into account in the radiative transfer equation itself, we can obtain from (6) the explicit form for *L̂*^{+} (Marchuk 1964; Box et al. 1988a; Ustinov 1991a),

and the boundary condition for *I*^{+}

Let us consider two atmospheres whose optical properties differ slightly (the first one we denote as “base” and the other as “perturbed”) and let

It can be shown that if only the most significant contributions are taken into account the effect perturbation *δE,* which corresponds to the perturbation *δL̂* of the atmosphere, is given by (Marchuk 1965; Box et al. 1988a; Ustinov 1991a)

More detailed results and discussion concerning higher-order terms can be found in (Box et al. 1988b). Unfortunately, it is difficult to provide simple criteria to determine whether a given atmosphere perturbation is small enough and hence to estimate the accuracy of (12), but some simulation results provided below will serve to give a rough idea.

## 3. Simulation technique

We introduce a Cartesian coordinate system with the *z* axis directed to atmosphere “top” along the normal to its boundary and the *x* axis in the plane of sunlight incidence. The direction of the light propagation is specified by Ω = (*ϑ, **ϕ*), where *ϑ* is the zenith angle measured from the positive *z* axis and *ϕ* is the azimuth angle measured clockwise from the positive *x* axis.

To perform the perturbation calculation we must evaluate integral (12). To simplify the integration, let both *I*(*z,* Ω) and *I*^{+}(*z,* Ω) be represented as a series of associated Legendre functions:

Here, *A*_{nm} = (2*n* + 1)(2 − *δ*_{0m}), *μ* = cos(*ϑ*), and *χ*(*z, **β*) [*β* = (Ω, Ω′) is the scattering angle] may also be expanded in a series of Legendre polynomials:

Note that in the most common case the atmosphere perturbation can be represented in the form

where the first member is due to the extinction coefficient variation *δσ*_{e}, the second correspondents to the single scattering albedo variation *δω*_{0}, and the third describes the influence of the phase function variation *δχ.* Therefore, the corresponding expression for the effect perturbation can be obtained by substitution of (13)–(15) into (12), and takes the form

It is interesting to point out that if we need to estimate the variation of the exiting radiances due to adding some aerosol characterized by the parameters *σ*^{P}_{e}, *σ*^{p}_{s}, *χ*^{P}, which are independent of altitude, Eq. (16) is simplified to

and the integral over *z* should be calculated only once, even if the influence of different aerosols is to be investigated.

The natural way to calculate *ψ*_{nm} is the spherical harmonics approximation (SHA) (Dave 1975). Moreover, it allows one to get both *ψ*_{nm} and *ψ*^{+}_{nm} during the same calculation without considerable expense of computer time. It is well known that SHA allows all integrals from the radiance over angles to be estimated with high accuracy, but unavoidable oscillations in the radiance angle distribution appear near the boundary. To improve the accuracy the iteration procedure (Lenoble 1985, p. 35) is used. This idea is based on using the formal solution of the radiative transfer equation (4), which has the form

where the optical thickness *τ*(*z*) is defined in the form

and the SHA solution *I*_{SHA}(*z,* Ω) is being used to get

There are two important reasons to improve the SHA solution. First, it is necessary to simulate satellite measurements. Second, to estimate some effects of the perturbation, the integral

needs to be calculated. If we take into account that both *I*(*z,* Ω) and *I*^{+}(*z,* Ω) are contributed significantly by the singular component *S*^{(+)} exp[−*τ*(*z*)/*μ*_{0}]*δ*(Ω − Ω_{0}) at small *z* (“0” denotes either the sun illumination or receiver angles), we meet again the necessity of the exact radiance calculation.

## 4. Numerical results

We consider the widely used model of the atmosphere as a stratified slab with an underlying surface. Although there is no restriction on the reflection properties of the underlying surface, for the sake of simplicity we consider the Lambertian case. Moreover, it is a reasonable approximation to model the real atmosphere–ocean system. In these simulations three profiles (Lenoble 1985) are chosen as appropriate for our purpose. The characteristics of the layers are given in Table 1. Note that model III is very similar to model I except for the large cloud. We introduce this model to check the applicability of the perturbation technique in the case of an optically thick atmosphere. The *U.S. Standard Atmosphere* is chosen as a model of the molecular atmosphere. We restrict our simulation to a single wavelength of 0.55 *μ*m, at which the contribution of the molecular atmosphere cannot be neglected, but it is still less than the aerosol component.

Let us consider a typical problem of the aerosol satellite remote sensing. We would like to investigate how the exiting radiances change if the basic characteristics of the lowermost layer, such as the optical thickness or the single scattering albedo, are changed. Also we are interested to estimate the influence of a small cloud contamination to it. Let us solve this problem using two approaches: the perturbation technique and the direct solution of the RTE, which is the common approach to do it. Taking into account that the perturbation technique assumes a linear dependence of the effect on the perturbation (12), it is convenient to evaluate the corresponding derivatives instead of the absolute values. Moreover, they describe the sensitivities of the exiting radiance to variation of these parameters.

Using the perturbation technique we calculate the following characteristics:

where *τ*_{C} is the cloud optical thickness, *ω*^{C}_{0} is the cloud single scattering albedo and *χ*^{C}_{n} are the expansion coefficients of the cloud phase function in Legendre polynomials. The direct solution of the RTE provides us with the corresponding estimates

We performed the calculation assuming the sensor azimuth angle is 90° and the surface albedo is 0.05 for several values of Δ*σ*_{e}, Δ*ω*_{0}, and Δ*τ*_{C}. This allows us to estimate the applicability range of the perturbation technique over the parameters variation area.

The results of our calculation for the atmosphere models are given in Tables 2–9. For convenience they also contain the relative error, which, for example, for the extinction coefficient perturbation has the form

Tables 2–4 show how the perturbation technique can predict the effect of the extinction coefficient variation. The range was chosen to be rather wide from −0.5*σ*_{e} to *σ*_{e} (i.e., from half the original value to double it). Tables 5–7 contain results for the single scattering albedo variations. Tables 8–9 show the sensitivity of the exiting radiance to a relatively small cloud contamination.

In all the tables we can see that the linear approximation (12) of the dependence of the effect on the parameter perturbation is reasonable in most of the case considered, especially for atmosphere models I and II, which are not optically thick. The tables give us the area of the perturbation technique applicability. For distinctness let us set 10% as a permissible limit of the relative error. We can see that for atmosphere model I, which describes clear atmosphere, the perturbation technique and direct simulation results coincide very well in practically all considered cases, except for a very few. This can be explained by the small multiple scattering contribution to the exiting radiances. For model II, which corresponds to a rather polluted atmosphere, the accuracy of the perturbation technique estimation becomes unacceptable at extinction coefficient variation |Δ*σ*_{e}| ≫ 0.1*σ*_{e}, but is still tolerable for all considered variations of the single scattering albedo, except Δ*ω*_{0} = −0.5.

Consideration of the results for model III shows that we meet an interesting situation, where noting that the optical thickness of the lowermost layer is not small (*τ* = 10), and its single scattering albedo is close to 1 (*ω*_{0} = 0.998), it is clear that multiple scattering contributes significantly to the exiting radiance, but the nonlinearity of its dependence on the extinction coefficient variation decreases slightly in comparison to model II. In contrast, Table 7 shows us that to predict the effect of the single scattering albedo variation on the exiting radiance we can use the perturbation technique only at |Δ*ω*_{0}| < 0.01 to meet the above criteria. Are there some contradictions? Let us recall some analytical results of the asymptotic theory of an optically thick weakly absorbing medium. For this case the exiting radiance for a semi-infinite layer is given by (Zege et al. 1991):

Here *I*_{0}(*z*_{t}, Ω, Ω_{0}) is the exiting radiance in the case of *ω*_{0} = 1, *γ* = 3(1 − *g*)(1 − *ω*_{0}), *q* = 1/[3(1 − *g*)], and *u*_{0}(*μ*_{0}) = [1 + (3/2)*μ*_{0}]/2. If we substitute into (23) the approximate values of the parameters (*g* = 0.86, *u*_{0}(*μ*_{0}) ≈ 1, *μ*_{0} = 1, *I*_{0}(*z*_{t}, Ω, Ω_{0}) ≈ 0.5), we find that *I* depends on *ω*_{0} as exp(−10(1 − *ω*_{0})). Therefore, the dependence of Δ*I* on Δ*ω*_{0} could be linear for *τ* ≫ 1 only when Δ*ω*_{0} ≪ 1 − *ω*_{0}, and that is why (1 − *ω*_{0} = 0.002) the relative error at |Δ*ω*_{0}| > 0.01 in Table 7 is so high. The perturbation technique points to it indirectly by giving the very high value for *dI*/*dω*_{0}, which is approximately 10 times greater than the exiting radiance.

Let us now look closely at Tables 8 and 9. They contain the results for the important case of small cloud contamination, which is difficult to detect with a high level of confidence. The tables show that the perturbation technique allows the effect to be predicted with the above-mentioned accuracy, except in the case of very large angles of illumination and observation. It is not surprising that the error of the prediction is less for model II than for model I at the same optical thickness of the cloud component because of the optical thickness of the lowermost layer of the model II is 5 times greater than for model I. However, despite such a difference, the sensitivity of the exiting radiance to a cloud contamination is substantially higher for model II than for model I.

## 5. Conclusions

The above comparison shows that the theoretically established relationship (12) between a small variation of the atmospheric parameters and the corresponding changes in the exiting radiance can be used for both clear and polluted aerosol atmospheres, even with the presence of small cloud contamination. Moreover, it can be the basis of a retrieval algorithm. Let us consider the basic idea of the MISR retrieval algorithm. The aerosol atmosphere is considered as a stratified medium, whose each layer can be contributed to by 11 possible basic aerosols (each model has its own range of altitudes). Starting from some initial model and following a sophisticated technique, based on using a precalculated radiance database, the aerosol optical thickness is retrieved. The perturbation technique can be incorporated into this algorithm in a natural way. Let us introduce a vector **T** = (*τ*_{1}, … , *τ*_{11}) (*τ*_{k} is the optical thickness of the *k*th pure component), which characterizes our aerosol atmosphere model. We also have the measured radiances, *I*_{M}, by satellite sensors, and the estimated ones, *I*_{E}, for the same geometry, but for some model **T**_{0}. The perturbation technique provides the necessary correction, Δ**T**, to the model **T**_{0}, on the basis of the comparison between *I*_{M} and *I*_{E}, which has the form

Note that the matrix ∂*I*/∂**T**_{0} is calculated together with *I*_{E} without noticeable addition to the computer time. It is difficult to estimate for the very first pixel how many iterations it takes to get the final solution, because of it strongly depends on the quality of the initial model. However, when the retrieval is started for the neighboring pixels, the situation is changed. Taking into account that the aerosol profiles at the neighboring pixels differ slightly as a rule, a few iterations should be enough to get the final solution. Additionally, Eq. (24) also allows us to obtain natural estimation of the retrieval errors for the each pure component if the calibration error is known.

However, the main goal of this paper is to show the applicability of the perturbation technique to describe how the variation of the aerosol atmosphere characteristics effects the exiting radiance for realistic atmosphere models. The comparison of the perturbation technique predictions with the results of the direct solution of the RTE shows its high accuracy. We emphasize that to use the perturbation technique we have to solve the RTE only once, plus some extra calculation. Note that if we need to analyze the dependence on a single parameter the time costs of both the perturbation technique and direct simulation are more or less similar, but if we need to investigate multiparameter dependences, the advantages of the perturbation technique are incontestable. That is why the perturbation technique can be successfully used to make sensitivity analyses of a remote sensing experiment on the basis of consideration of very complex models and, as mentioned above, as a core of the retrieval algorithm. These topics have not been covered in this paper, but will be considered in future work.

## Acknowledgments

This work is supported by the Australian Research Council.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Michael A. Box, School of Physics, University of New South Wales, Sydney 2052, Australia. Email: m.box@unsw.edu.au