A Direct Numerical Solution to the
Inverse Radiative Transfer Problem

In many applications, it is desirable to directly retrieve physical parameters from measured data, rather than indirectly obtaining the former by trying to find best-fit solutions of model calculations. For instance, in the theory of radiative transfer in an emitting/absorbing atmosphere, the intensity observed outside the atmosphere (t=0) for a given line of sight i is given by the equation (assuming a plane-parallel geometry):

(1)       I(μ_i) = 1/μ_i^.₀∫^τ_vdτ' S(τ')^.e^-τ'/μ_i     ,

where

(2)       μ_i = cos(η_i)     ,

with η_i the local zenith angle (i.e. μ_i =1 for a vertical line of sight), τ_v the total vertical optical depth, and S(τ') the source function (emissivity) of the atmosphere at level τ' (note that the exponential function under the integral implies a continuum (i.e. frequency independent) absorption; in case of radiative transfer in spectral lines for instance, this factor would take on a different form).

If one subdivides the atmosphere into N layers and assumes S(τ') as constant within each layer, one can write Eq.(1) in the discrete form (after performing the integration over the exponential)

                       _N
(3)       I(μ_i) = Σ S(τ_k)^.(e^-τ_k/μ_i - e^-τ_k+1/μ_i ) =
                      ^k=1
                       _N
                  := Σ S(τ_k)^.A_i,k     .
                      ^k=1

If one has now N measurements of the intensity I(μ_i) for different values of μ_i , Eq.(3) constitutes a system of N equations for the N unknowns S(τ_k) (assuming the total vertical optical depth τ_v (and thus the coefficient matrix A_i,k) as given).

It is well known however that this system of equations is ill-conditioned due to the nature of the coefficients A_i,k, and a direct inversion of Eq.(3) proves in general impossible due to numerical instabilities. Therefore, the solution is usually obtained by an implicit technique where forward model calculations (i.e. an evaluation of Eq.(1) for an a priori assumed source function) are compared to the measured intensities and successively improved by means of certain algorithms until the agreement is satisfactory.

By means of casting Eq.(3) into a different form, one can improve the conditions of the equation however, and the source function can be directly obtained by using a straightforward iteration scheme. First of all, in order to bring the system of equations into a form suitable for an iterative solution, one solves it for S(τ_i), i.e.

                                  _N
(4)       S(τ_i) = [ I(μ_i) -Σ S(τ_k)^.A_i,k ]/A_i,i     .
                                 ^k=1(≠i)

This alone is not sufficient yet however, as the equation is still ill-conditioned and the iteration will indeed not converge. It is evident from the equation that a potential problem can arise if the diagonal element A_i,i is too small, as this will lead to large values of the source function S and then possibly on the subsequent iteration to negative values for S. This suggests to add some constant a_i to the diagonal elements A_i,i, and, in order to keep the equation algebraically correct, a term a_i^.S(τ_i) to the angular bracket. Basically, all sufficiently large values a_i turn out to work in this case (i.e. enable convergence of the iteration), but larger values also mean more iterations, and the most efficient choice here is to make a_i just equal to the sum over the absolute values of the off-diagonal row elements i.e.

                   _N
(5)       a_i = Σ|A_i,k|     ,
                   ^k=1(≠i)

and thus

                                                 _N
(6)       S(τ_i) = [ I(μ_i) +a_i^.S(τ_i) -Σ S(τ_k)^.A_i,k ]/(A_i,i+a_i)     .
                                                ^k=1(≠i)

Note that Eq.(6) is formally exactly identical to Eq.(4) (as evident by multiplying the equation through by A_i,i+a_i ), but the point is that S(τ_i) on the right hand side is now the source function value of the previous iteration. So in practice the iteration equation is different, and only upon convergence approaches the original equation.

The following tables give some numerical results I obtained. In all cases a total vertical optical depth τ_v=5 was chosen, subdivided into 10 layers (i.e. the optical depth of each layer i was thus 0.5). By means of Eq.(3), synthetic intensities were then first produced (for 10 viewing directions μ_i) by assuming given source functions S₀(τ_i) and additionally adding a random noise of a given relative magnitude to it. These intensities I₀(μ_i) were then used for the inverted version of Eq.(3) (which solves for the retrieved source functions S(τ_i) and this system of equations was then iterated until a set convergence criterion was achieved. In each case, the convergence value was twice the noise level (it does not make much sense to iterate further as either convergence problems could arise, or the solution would become random and reflect the noise rather than the systematic data).

The first two tables were obtained by assuming the source function as linearly increasing with optical depth (in fact S₀(τ_i)= τ_i+1 ), the other two tables assume an exponentially increasing source function (S₀(τ_i) = e^τ_i/2). In each case, noise levels of 10^-2 and 10^-3 were considered. The retrieved values are in bold blue print here.

Note that in all cases the source function values agree less well than the intensities, especially as the bottom of the atmosphere is approached. This is an expected behaviour caused by the exponential factor in Eq.(1), which strongly suppresses any contributions of layers at large optical depths to the measured intensities. Relatively large source function variations at deeper regions will therefore have little or no impact on the resulting intensities measured outside the top of the atmosphere (given the limited accuracy of the latter associated with the data noise).

Also note that in the above cases the optical depth assumed for the inversion was exactly identical to that used in producing the synthetic data. In reality this will of course not be exactly known but only approximately (as determined by independent measurements). But as long as the errors associated with this are not too drastic, the results should be reasonably close to the correct case (an artificial error of 10% in the assumed vertical optical proved to be insignificant for the accuracy of the retrieved source function throughout the whole range in case of the 10^-2 noise level results above).

The number of layers/data-pairs N can of course be arbitrarily increased in principle, but this will obviously be associated with a corresponding increase in computation time (and might in many cases anyway not make sense due to the inherent inaccuracies associated with the data noise).

The procedure proved to work as well if the coefficient matrix A_i,k is not given by the exponentials as shown above (which would be appropriate for a continuum (i.e. frequency-independent) absorption coefficient), but by the corresponding kernel function appropriate for the radiative transfer in spectral lines (in fact, the convergence is even better as the function varies much less strongly with τ).
This suggests that the particular mathematical method used here to enable a direct inversion of the problem is not restricted to the inversion of the radiative transfer equation, but should be applicable to other remote sensing applications and in general to inverse problems considered to be 'ill-conditioned'.

I have listed the C/C++ program code used here separately on the page Program Code for Solution (Inversion) of Linear System of Equations (this is not restricted to this particular problem, so it requires the coefficient matrix A_i,k as input).
I also made an Online-Demo available where you can put in your own data (this is however limited to the problem described on this page) (note that this page is password-protected; use username:demo, password:demo).

Print Version

Home

Search Site