# A Direct Numerical Solution to theInverse 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.0τ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).(eki - ek+1i ) =
k=1
N
:= Σ S(τk).Ai,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 Ai,k) as given).

It is well known however that this system of equations is ill-conditioned due to the nature of the coefficients Ai,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).Ai,k ]/Ai,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 Ai,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 ai to the diagonal elements Ai,i, and, in order to keep the equation algebraically correct, a term ai.S(τi) to the angular bracket. Basically, all sufficiently large values ai 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 ai just equal to the sum over the absolute values of the off-diagonal row elements i.e.

N
(5)       ai = Σ|Ai,k|     ,
k=1(≠i)

and thus

N
(6)       S(τi) = [ I(μi) +ai.S(τi) -Σ S(τk).Ai,k ]/(Ai,i+ai)     .
k=1(≠i)

Note that Eq.(6) is formally exactly identical to Eq.(4) (as evident by multiplying the equation through by Ai,i+ai ), 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 S0i) and additionally adding a random noise of a given relative magnitude to it. These intensities I0i) 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 S0i)= τi+1 ), the other two tables assume an exponentially increasing source function (S0i) = 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.

Table 1: linear source function; noise = 10-2
``` i      S0(τi)      S(τi)            μi          I0(μi)       I(μi)
1    0.50000   0.52859      0.52632   0.81878   0.82140
2    1.00000   1.05124      0.55556   0.84951   0.84592
3    1.50000   1.29770      0.58824   0.87492   0.87354
4    2.00000   1.93520      0.62500   0.91341   0.90483
5    2.50000   2.32458      0.66667   0.94817   0.94043
6    3.00000   2.91554      0.71429   0.99188   0.98115
7    3.50000   3.45683      0.76923   1.03834   1.02787
8    4.00000   4.19160      0.83333   1.09479   1.08164
9    4.50000   5.20870      0.90909   1.16424   1.14354
10  5.00000   6.20795      1.00000   1.23911   1.21459
```
after 36 iterations in 0.025 sec

Table 2: linear source function; noise = 10-3
``` i      S0(τi)      S(τi)            μi          I0(μi)       I(μi)
1    0.50000   0.51371      0.52632   0.81511   0.81512
2    1.00000   0.95334      0.55556   0.84193   0.84201
3    1.50000   1.50472      0.58824   0.87240   0.87221
4    2.00000   2.06561      0.62500   0.90663   0.90626
5    2.50000   2.47072      0.66667   0.94452   0.94479
6    3.00000   3.03972      0.71429   0.98820   0.98850
7    3.50000   3.63486      0.76923   1.03790   1.03821
8    4.00000   4.11129      0.83333   1.09392   1.09478
9    4.50000   4.41105      0.90909   1.15696   1.15907
10  5.00000   4.89185      1.00000   1.22935   1.23180
```
after 233 iterations in 0.15 sec

Table 3: exponential source function; noise = 10-2
``` i      S0(τi)      S(τi)            μi          I0(μi)       I(μi)
1    1.00000   0.94527      0.52632   1.22080   1.22183
2    1.28403   1.36377      0.55556   1.24668   1.24598
3    1.64872   1.88247      0.58824   1.27682   1.27322
4    2.11700   2.31424      0.62500   1.30797   1.30412
5    2.71828   2.41485      0.66667   1.33586   1.33936
6    3.49034   2.95783      0.71429   1.37623   1.37978
7    4.48169   4.02274      0.76923   1.43132   1.42636
8    5.75460   4.90534      0.83333   1.48807   1.48022
9    7.38906   6.38219      0.90909   1.56242   1.54257
10  9.48774   8.03012      1.00000   1.64666   1.61458
```
after 50 iterations in 0.074 sec

Table 4: exponential source function; noise = 10-3
``` i      S0(τi)      S(τi)            μi          I0(μi)       I(μi)
1    1.00000   0.99339      0.52632   1.21714   1.21712
2    1.28403   1.27555      0.55556   1.23997   1.24006
3    1.64872   1.70467      0.58824   1.26653   1.26646
4    2.11700   2.23120      0.62500   1.29731   1.29702
5    2.71828   2.59056      0.66667   1.33215   1.33262
6    3.49034   3.31145      0.71429   1.37386   1.37432
7    4.48169   4.29783      0.76923   1.42322   1.42339
8    5.75460   5.72540      0.83333   1.48214   1.48132
9    7.38906   7.33265      0.90909   1.55112   1.54976
10  9.48774   9.73347      1.00000   1.63368   1.63042
```
after 287 iterations in 0.31 sec

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 Ai,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 Ai,k as input).