H. Berger , T. Hay and E. Margerum
USA Topographic Engineering Laboratory, Alexandria, VA 22315
ABSTRACT
Most atmospheric correction
codes are based in part or in whole on the radiative transfer equation (RTE),
which is an integrodifferential equation.
It is well known to be an ad hoc equation which can and has produced
incorrect answers. This paper initiates
a program for a new way of exploring where the RTE can produce unphysical
answers in parameter-ratio ranges of genuine concern. The program begins with formulating a new technique for
rigorously transforming the scalar RTE, without approximation, into a “pure” partial
differential equation (PDE), i.e. one involving only partial derivatives of
finite and relative small order. The
virtue of this approach is that there are only a small number of analytical and
numerical techniques for dealing with integrodifferential equations compared to
the vast array of techniques for PDEs.
We develop a variety of tools more powerful than needed for the
particular physical problem to demonstrate the robustness of the
technique. We then consider an
atmosphere where Rayleigh scattering is dominant and derive its PDE (We are unaware of anyone having done this
before). We also transform a class of nonlinear integrodifferential
equations into linear PDEs and solve for a multiplicity of solutions.
Keywords: radiative transfer equation, atmospheric
correction, Rayleigh scattering atmosphere, linear integrodifferential
equations, linear partial differential equations, nonlinear integrodifferential
equations, annihilation operators, Whitham’s method
1.
INTRODUCTION
The motivation for this work
is the concern about the very real possibility that most major atmospheric
correction codes for optical wavelengths may have parameter-ratio ranges where
they give unphysical results. This is
because most are based on the scalar radiative transport equation (RTE) (or one
of its generalizations) which has long been known to be an ad hoc equation.1
The RTE is an
integrodifferential equation of such a character that the analytical and
numerical techniques available to investigate its solutions nature are
relatively few compared with the immense array of techniques developed over the
years for partial differential equations (PDE). For this reason this paper demonstrates that the RTE can be
rigorously reformulated, without approximation, as a “pure” linear PDE, i.e.
one involving only partial derivatives.
There is, however, a different PDE corresponding to each phase function
(defined later).
The paper applies a new
technique developed to derive the RTE-based PDE for a Rayleigh scattering
atmosphere, a fourth order linear PDE.
To obtain a clear sense of what this has accomplished, the reader is
invited to inspect the traditional set of
coupled integral equations with variables that are integrals of the ones
of interest, for Rayleigh scattering, that are derived from the RTE, given in
the Appendix. Note that Rayleigh
scattering occurs even in “clear” air.
The ultimate purpose of this
program is 1) to broaden the search for such unphysical parameter-ratio ranges
in optical-wavelength atmospheric analyses for parameter values of any real
concern and 2) to explore the power of the PDE approach to atmospheric
problems. The basic program is to
develop solutions of the RTE for different atmospheric phenomena of interest
using the new technique and compare them to rigorous solutions, where
available, from electromagnetic theory.
This should expose any such
unphysical parameter ranges. Rayleigh
scattering is a good candidate for such an approach because of the availability
of electromagnetic solutions and because a comparison with the traditional
coupled integral equations formulation for Rayleigh scattering gives some
measure of the simplicity introduced by the PDE approach.
There have been a small
number of cases in the literature considering special, specific circumstances
under which the RTE has been demonstrated to give the correct answer or a good
approximation. This has strengthened
confidence that it can often be relied on to produce accurate results. The
literature, however, observes that incoherent fields are assumed so that power
additivity takes place rather than field additivity or some other
relationship. It is known from the Van
Cittert-Zernike theorem of optical coherence that even spatially incoherent
sources can generate highly coherent fields over large regions of space. The literature also relates the solution of
the RTE as a statistical average of a randomly fluctuating Poynting vector
whose Fourier transform is approximately the mutual coherence function. 1
Ishimiru also observes that although the description
of interference and diffraction effects are included in the scattering and
absorption characteristics of a single particle, transport theory itself does
not include diffraction effects.1 As a specific example of the
problem with using the RTE, we cite a recent work in which a careful
field-based analysis of scattering in a dense medium with particle correlations
considered produced a governing equation in the same general form as the RTE
but differing in some details.2
The RTE has been recently demonstrated to have parameter-ratio ranges
with the potential for large errors in microwave frequency non-atmospheric
investigations for a few special cases.
In those it was possible to use numerical solutions of the rigorous
electromagnetic theory for comparison.
That work was reported in a number of journal articles.
To illustrate how the transformation technique works, Section 2 discusses in considerable detail an example of an integrodifferential equation is solved so that the new method will be transparent to almost any reader. This example can be followed easily, but nevertheless contains the essence of the difficulties to be surmounted.
The use of the technique is
illustrated again in Section 5 to solve an example of a nonlinear
integrodifferential equation, which satisfies the defining equation and the
boundary conditions. If [f(y)] appears
in the equation raised to the Nth power (N = integer), then N boundary
conditions will be needed. Two boundary
conditions are required here.
Withham discussed a technique
for converting some integrodifferential equations into PDEs, but works for more
limited circumstances than are of interest here.3 Ochs and Kristensson considered a problem
which involved the propagation of electromagnetic waves in two types of dispersive
dielectrics.4 The annihilation operators required were second order
constant coefficients involving ordinary derivatives. The functions to be annihilated in the integrands were
separable. There was one single
integral in each partial differential equation.
Those authors expressed an
awareness that the annihilation operator method they used, which is a valuable
one and clearly useful here, could be generalized but gave no sense this could
include nonlinear equations.4 Unlike the technique used in reference 4, the method developed
here uses the defining equation directly to determine the additional
relationships needed for a unique solution.
Furthermore, it can be used whether ordinary or partial derivatives are
present.
Previously one of the authors
(H. Berger) presented, in an unpublished 30 March 1998 memorandum, a method for
converting the scalar RTE into a second order linear partial differential
equation, free of any integrals, for the special case of isotropic, and possibly lossy, scatterers.
The method discussed here
represent a method capable of attacking much more general scattering
problems and can be seen to include the technique this author used earlier.
There are few, if any,
universities in the USA which offer full length courses on integrodifferential
equations. There is only a tiny
community of people skilled in dealing with them. Yet almost all universities offer multiple full-length courses on
PDEs. Most physics and mathematics
majors are required to take multiple full length courses on them in their
education and most electrical, aeronautical and mechanical engineering students
are required to take at least one such course.
Thus transforming the RTE into PDEs opens the area to exploration by
technical professionals from many backgrounds including some we have not
listed.
Further, the solution of PDEs
is currently a very active area with a great many practitioners, and has been
for decades, in many engineering disciplines.
The three IEEE Transactions, Microwave Theory and Techniques, Antennas
and Propagation and Electromagnetic
Compatibility issue research and exploratory reports that are dominated by
articles authored by people that apply PDEs in new and innovative ways to solve
engineering problems. Of course the
various mathematics societies offer many journals devoted exclusively to
PDEs. This brief description in no way
comes close to exhausting the list of journals that include many articles in
this area. Finally the modern monograph
and textbook literature on PDEs so completely dwarfs what is available on
integrodifferential equations that comparison is difficult.
While we are interested in
solving the radiative transfer equation, to illustrate how the basic technique
works, a much simpler example is first offered here involving a linear
integrodifferential equation with an ordinary derivative. The method has been used when partial
derivatives are present and the more complex RTE is considered later. We will consider the integrodifferential
equation for f(x) given in eq.(1). The
basic idea is to try to convert this equation into an ordinary linear
differential equation and to use the defining equation to constrain the
relation between what appears to be an excess of independent constants to
insure that the correct and unique solution is obtained. The equation and boundary condition are
2
df/dx = ò g(x) f(y) dy, g(x) = x, f(0) = A = constant (1)
0
For this specific function
g(x) = x we can find a particular
differential operator (D2) that does the following.
(D2)
g(x) = ( d2 / dx2 ) g(x) =
0
(2)
This equation also defines (D2)
that is said to annihilate g(x).
If we apply the operator (D2)
to both sides of eq.(1) we quickly find
2
d [(D2)
f] / dx = ò (D2) g(x)
f(y) dy
(3)
0
Since the operator, (D2),
operates only on functions of x and (D2) g(x) =
0. then the integral in eq.(3)
becomes zero, and the equation reduces to the linear ordinary differential equation
d [(D2) f] / dx
= d3f / dx3 =
0. (4)
At first glance, this is a simpler equation to solve than eq.(1). Of course more boundary conditions are needed with eq.(4) than with eq.(1). The general solution to eq.(4) is obviously
f(x) = B1x2
+ B2x + B3 (5)
B1, B2 = undetermined constants, B3 = A
= boundary condition constant
This process ends,
incomplete, unless we can generate more boundary conditions. To do this eq.(1) is rewritten as
df / dx = C1 x
(6)
where C1 is the
integral in eq.(1) after x is factored outside of it. Eq.(6) integrates to
f(x) = (C1/2)
x2 + C2, C2
= constant
(7)
Comparing eqs.(5) and (7) it
is clear that B2 = 0 and A =
B3 = C2.
Now we are not free to generate a new boundary condition because the
relation between C1 and A are constrained by eq.(1). Substituting Eq.(7) into Eq.(1)
produces
2
C1 = ò [ (C1/2) y2 + C2
] dy =
( 8/6) C1 + 2 C2 (8)
0
Carrying out the remaining
simple algebra produces the solution
f(x) = - 6A [ (x2/2)
- (1/6) ]
(9)
This is a valid solution
because it is seen to satisfy the original equation, eq. (1), and the boundary
condition. The RTE certainly appears
more challenging than this example but the same basic approach will allow its
transformation into a PDE whose form will depend on the form of the phase
function (to be defined later).
Integrodifferential equations involving partial derivatives have been
solved by this method as well.
The annihilation
operator technique employed in Reference 4 uses linear constant-coefficient
differential operators of first and second order and the integrands can be
expressed as a sum of separable polynomial functions. We will generalize that here in a fairly obvious way, and to an
extent beyond what is needed to consider the Rayleigh scattering problem that
we focus on in later sections. The
reason we do this is to make clear that this is a robust method capable of
ultimately handling much more general scattering problems. But first we try to make clear what the
terms used in the first sentence of this paragraph mean.
The second order linear
constant-coefficient differential operator has the form
DA2 = A0 + A1
(d /dx) + A2 (d2 /dx2)
(10)
where the A coefficients are constants. The integrands in Reference 3 have the form
I22(x,y) = S g2(x) f2(y)
(11)
where
S is the summation sign and g2(x) and f2(y)
are both quadratic (second-order) polynomials.
It
is clear from the relatively simple example in the Introduction, eqs. (1) –
(9), that the same kind of results are obtainable for an arbitrary function
g(x), so long as there exists an N th order constant-coefficient
operator
N
DAN = S An
(dn / dxn) (12)
n=0
such
that (DAN) g(x) = 0. If we
could find such an operator, i.e. an annihilation operator for g(x), then its
use, when we have an arbitrary g(x) in eq.(1) instead of just x, would, in
place of eq.(4), lead to
dNf / dxN
= 0.
(13)
It
is probably also reasonably clear that if instead of having g(x) in the
integrand there appeared some g(x,y) the method would still work so long as (DAN)
g(x,y) = 0. This is in fact the case
for Rayleigh scatters as we show and make use of later. The next step in generalizing the
annihilation operator method is to
consider the operator in eq.(10) when the An are not constant but
rather functions of x, i.e. An = An(x). This
can be demonstrated in a straightforward manner.
More
generally if the operator in eq. (12) has coefficients An which are
functions of x then the arguments in the paragraph immediately above,
supporting eqs. (7–8), still apply. We
then see that the annihilation operator technique applies for all g(x,y) and An(x)
for which
(DAN) g(x,y) = 0
(14)
and for which integrals like those in the equations exist.
Now
consider an integrodifferential equation involving two integrals. For clarity we use the relatively following
simple case as an example.
b d
df/dx = ò g(x,y) f(y) dy
+ ò h(x,z) f(z) dz (15)
a c
If
(DAN)
g(x,y) = 0 and (DBM) h(x,z) =
0 (16)
(DAN)
(DBM) = (DBM) (DAN)
then applying the two operators to eq. (15) we
get the ordinary differential equation
d[(DAN)
(DBM) f(x)] / dx = 0 (17)
which is free of integrals,
the desired result.
This idea is easily generalized in a trivial, but important, way to the more complex circumstance where the integrand of any one integral in an integrodifferential equation involves a kernel function that is the sum of terms. Each of which are the products of functions and for one of each of these the corresponding annihilator operators are known. This frequently comes up when the phase function in the RTE is expanded in a Legendre polynomial series.
Finally these ideas are applied to integrodifferential equations
involving one or more multiple integrals.
We illustrate the idea through another relatively simple example. Consider
¶f(x,y)/¶x = òò g(x,y,w,z) f(w,z,y) dw dz + òò h(x,y,w,z) f(w,z,y) dw dz, (18)
where the limits of the
integrations are fixed numbers. Then
if (DAN) g(x,y,z) = 0,
and (DBM) h(x,y,z) = 0 and are communtaive operators and both
of these operators are applied to eq. (18)
we now get a partial differential equation, free of integrals, that appears as
¶ [ (DAN) (DBM) f(x,y) ] / ¶x = 0. (19)
Most of the generalizations of the annihilation operator technique in this section can be further generalized and are more advanced than what is needed to investigate the Rayleigh scattering problem. However the intent was to demonstrate that the method could be generalized to the point where a wide variety of atmospheric scattering problems, possibly some much more complex than Rayleigh scattering, could be investigated with this technique.
4. THE PDE FOR
THE RTE AND A RAYLEIGH SCATTERING ATMOSPHERE
In this section we briefly
consider the application of the annihilation operator technique to the scalar
radiative transfer equation. That
equation, for azimuth-independent plane-parallel media, is the
integrodifferential equation (Reference 5, pg. 15)
1
m ¶I(t, m)
/ ¶t = I(t,m) + (1/2) ò p(m,m’)
I(t,m’)
dm’ (20)
-1
where I = specific intensity,
m = cosJ , J = zenith angle and t = optical thickness. If p(m,m’) is annihilated by some operator, (Dm), then it follows from the previous discussion that
¶ [(Dm) mI] / ¶t = (Dm) I
(21)
In this section some of the
potential virtues of using the annihilation operator technique may become
clearer. The atmospheric phase function
for Rayleigh scattering is well known and is given on page 121 in reference 6
(after integration over the azimuth angles j and j’) as
p(m,m’) = (3/4) a [ 1 + m2m’2 + (1/2) (1 - m 2)(1 - m’2) ]
(22)
Then straightforward
differentiation yields
¶3p/¶m3 =
0 [23]
From this we see that the
annihilation operator for the Rayleigh phase function is ¶3 /¶m3. Applying this annihilation operator to the
RTE in eq. (20) with p( m,m’ ) given in eq. (22), we find that
3¶3I / ¶m 2¶t + (m) ¶4I / ¶m3¶t = ¶3I / ¶m3
(24)
This fourth order, linear,
variable coefficient PDE is then
equivalent to the radiative transfer equation for Rayleigh scattering. Generating an additional two boundary
conditions here may not be straightforward but this may not be necessary when
the method in the elementary example is used..
5. SOLVING NONLINEAR INTEGRODIFFERENTAL EQUATIONS
As mentioned in the Introduction, the scalar radiative transport equation has been derived by some authors as an approximation, indirectly from some forms of the Boltzmann transport theory and directly from some random media theories. Some of these theories take nonlinear forms that have remained largely intractable. Thus it is of some interest that we have unexpectedly found that the generalized method of annihilation operators can be used to solve certain kinds of nonlinear integrodifferential equations that it converts to linear ordinary or linear partial differential equations.
We discuss the basic idea by
outlining the solution method of a particular example. The accuracy of the solution is
verified. It also seems as though an
additional boundary condition is needed for uniqueness that is not required in
the linear case. Additionally the
method seems like it may be applicable to a nonlinear form of the radiative
transport equation.
The integrodifferential
equation and boundary condition we consider
is
2
df/dx = ò x [ f(y) ]2 dy.
f(0) = A [25)
0
The method of solution runs
parallel to that given for eq.(1) but with one key aspect different. Following the kind of analysis given for
eq.(1) we find with some messier algebra that eq. (25) yields
f(x) = - A { [x2 /
K1,2 ] - 1 }
(26)
where
K1,2 = -
{ (4/3) ± Ö(26/90) }
Thus we obtain two solutions from eq. (30), corresponding to the ± signs that satisfy the equation and the boundary condition.
This type of result is not
seen for linear problems and strongly suggests that an additional boundary
condition is needed for a unique solution.
This needs to be explored further..
The preceding approach seems as though it should also be successful for
integrands involving f(x) raised to the Nth power, N = any integer, which
should produce N distinct solutions satisfying the defining equation and
the single boundary condition
used. Here it would seem that N
boundary conditions are needed to obtain a unique solution.
6: WHITHAM’S
METHOD FOR OBTAINING
A PDE FROM AN
INTEGRODIFFERENTIAL EQUATION
There exists in the
literature a special case where an equation that is similar to the RTE is
converted into a partial differential equation (PDE). We will briefly comment on this and observe that the approach
used is far more limited than the annihilation operator approach. Whitham considers the equation
µ
¶y(x,t)/¶t + ò K(x-z) {¶y(z,t)/¶z}dz = 0 (27)
-µ
where µ represents infinity.2
Whitham observes that because
this is a linear equation then one would expect elementary solutions of the
form
y = exp (kx – wt)
(28)
He indicates that this is
possible if the phase velocity, given
by the ratio
c(k) = w/k,
(29)
is the Fourier transform of
the kernal function, K(x). He further
indicates that if c(k) can be represented by a finite length polynomial
involving k raised only to even integer powers then one could derive a partial
differential equation, which he displays, that is equivalent to the original
integrodifferential equation.
One could use the same kind
of reasoning for the RTE and arrive at the same kind of conclusion as well as
an equivalent partial differential equation under the same very special
conditions. This approach, which is
very interesting, is not equivalent in any way to the annihilation operator
approach. It is very limited in its
applicability while the annihilator approach will yield a partial differential
equation under far more general circumstances, specifically whenever one can
find the appropriate annihilation operator.
We infer from the ease of doing so for Rayleigh scattering, above, that
it may not be too difficult for many other practical problems.
7. CONCLUSIONS
This paper has reminded the
atmospheric correction community of the potential for serious problems in the
use of current major atmospheric correction codes and discussed an approach to
investigating that problem. The problem
arises from the use of the radiative transfer equation (RTE), which is widely
regarded as an ad hoc usage in atmospheric problems. Recently published work, carried out for microwave frequency
non-atmospheric problems, have shown in special cases where rigorous results
from electromagnetic theory could be obtained for comparison, parameter ranges
in which the use of the RTE produced serious errors.
The approach to uncovering similar parameter ranges for optical
wavelength atmospheric problems we are following is the same: find instances where solutions to the RTE
could be found and compared with rigorous solutions from electromagnetic
theory. A method has been developed
which seems to make it far more likely that this approach could actually be
persued on a phenomenon-by-phenomenon basis starting with the Rayleigh
scattering problem.
.
This paper has considerably extended and augmented the annihilation operator technique demonstrated in Reference 4 for the equations of electromagnetic wave propagation in dispersive dielectrics. It makes use of the original defining equation as a constraint and can be applied to both integrodifferential equations with ordinary and partial derivatives present. It has applied this to the scalar radiative transfer equation usually used in atmospheric correction codes. Although the techniques developed here should be useful for many different kinds of scatterers in the atmosphere we have applied it here to Rayleigh scattering which occurrs even in clear air.
A pure partial differential equation (PDE) has been rigorously derived without approximation, i.e. one involving only partial derivatives. It is the first such RTE-based equation for Rayleigh scattering to our knowledge. The virtue of having transformed the problem into this form is that there is an immense solution technique literature for PDEs and a very small one for the RTE beyond numerical methods. Although the PDE for Rayleigh scattering is a fourth-order PDE with variable coefficients some progress toward solving it has already been made and will be reported elsewhere. Other RTEs with different phase functions will transform into different PDEs.
The work reported here
includes the extension of the above methods to certain types of non-linear
integrodifferential equations. This has
been done by way of a relatively elementary example. This example suggests, if
generalizations are valid, it may be useful in dealing with certain more
complex theories regarded by some as more accurate expressions of transport
theory than the RTE for atmospheric problems.
REFERENCES
1. A. Ishimaru, “Wave Propagation and Scattering in Random Media,”
IEEE Press, N.Y. and Oxford Univ. Press, Oxford, Chap. 7, 1997.
2. B. Wen, L Tsang, D.P. Winebrenner and A. Ishimaru, “Dense Medium
Radiative Tranfer Theory: Comparison
with Experiment and Application to Microwave Remote Sensing and Polarimetry,
IEEE Trans. Geosci, Remote Sensing , 28.
pp. 46-59, Jan, 1990.
3. G. B. Whitham, “Linear and Nonlinear Waves”, John Wiley, New
York, p. 368, 1974.
4. R. L. Ochs, Jr., and G. Kristensson, “Using Local Differential
Operators to Model Dispersion in Dielectric Media,” J. Opt. Soc. Am., Part
A, 15,
pp. 2208-2215, Aug. !998.
5. S. Chandrasekhar,
“Radiative Transfer”, Dover Pubs., New York, 1960.
6. M. E. Thomas and D. D. Duncan, “Atmospheric Transmission”, Chap.
1 in “Atmospheric Propagation of Radiation”, F. G. Smith, Ed., Vol. 2 of “The
Infrared and Electro-Optical Systems Handbook,
J. S. Aceetta, et. al., Ed., SPIE Optical Engineering Press, 1993.
APPENDIX: THE TRADITIONAL COUPLED INTEGRAL
EQUATIONS FOR
RAYLEIGH SCATTERING
Chandrasekhar considers
Rayleigh scattering in a plane-parallel azimuth-averaged homogeneous atmosphere
and derives the integral equation for it as5
1 1
m ¶I(t, m) / ¶t = I(t,m) - (3/16)[ (3- u2)
ò I(t,m’) dm’ +
(3u2 – 1) ò (m’)2
I(t,m’)
dm’] (30)
-1 -1
He then observes that this
integral equation with multiple integrals can be reformulated as a fairly
complex set of coupled linear integral equations which we give below.5 But first we must define the key symbols and
how they relate to I. These are
J(t, m
) =
1/(4Õ) ò I(t,m ) dW (31)
where dW is a differential element of solid angle and the
axial-symmetric K-integral is
1
K(t) = (1/2) ò I(t , m ,
j) [m]2 dm (32)
-1
It is also necessary to
define a type of exponential integral with an integer index, n, which takes on
the values 1,3 and 5 in the coupled integral equations. It is
µ
En(y) = ò {[ exp(-xy) ] / xn } dx , y = ½t -