MASTER'S THESIS
Semi-Infinite Solutions to the Radiative Transfer Equation Applied on Snow
Benjamin Friberg 2014
Master of Science in Engineering Technology Engineering Physics and Electrical Engineering
Luleå University of Technology
Institutionen för teknikvetenskap och matematik
Semi-Infinite Solutions to the Radiative Transfer Equation Applied on Snow
Benjamin Friberg
Lule˚ a University of Technology
Department of Engineering Sciences and Mathematics Division of Fluid and Experimental Mechanics
March 3, 2014
Abstract
The radiative transfer equation has been solved numerically for semi-infinite plane parallel cases using Chandrasekhar’s H-functions, as well as by implementing a more general matrix formulation.
The solutions have then been compared to a parametric model, derived by Rosendahl.
The comparisons showed that the solution to the radiative transfer equation is non-linear in its behaviour with respect to the single scattering albedo and that a more general method for com- paring the two models is needed. A more general method is therefore suggested and implemented as an inverse problem formulation.
The results from the numerical simulations using the inverse problem formulation, showed that neither the use of the Henyey Greenstein phase function, nor a more general low order phase function exhibits similar behaviour as the parametric model. This thesis shows that it is possible to modify the parametric model to better suit radiative transfer approximations.
Finally, frequency resolved measurements has been done on white, fine grained snow. The inverse problem formulation is then used to obtain the corresponding solution to the radiative transfer equation. The results from the measurements show that it is possible to model the reflec- tion distribution from a snow surface using the Henyey Greenstein phase function. Better results can be obtained, however, by using a more general low order phase function. Comparisons with the modified parametric model showed that the modified model has difficulties in resolving the direct reflection that arise from the snow surface. When measurements were done on an ice covered snow surface, the results showed that more complex modelling is necessary.
preface
The making of this thesis concludes my final year on my MSc in Applied Physics and Electrical Engineering. Now - at the very end of my studies, it feels like this thesis project has given me the opportunity to test my problem solving skills to the fullest.
I would like to take the opportunity to thank Prof. Mikael Sj¨odahl and Dr. Johan Casselgren for their support and encouragement. I also would like to thank Dr. Per Gren and Dr. Henrik Lycksam for their invaluable help during the experimental part and Dr. Sara Rosendahl for her help with questions concerning the parametric model.
Contents
1 Introduction 3
1.1 Outline . . . 3
1.2 Frequently used variables and abbreviations . . . 4
2 Theory 5 2.1 An introduction to radiative transfer . . . 5
2.2 The phase function . . . 7
2.3 Normalization of the phase function . . . 9
2.4 The optical depth . . . 10
2.5 Problem formulation . . . 10
2.6 Problem discretization . . . 11
2.7 The scattering function . . . 12
2.8 Evaluation of the H-functions . . . 12
2.9 The characteristic polynomial, Ψ . . . 13
2.10 A matrix formulation . . . 14
2.11 The parametric model . . . 17
3 Numerical implementation 19 3.1 Evaluation of the H-functions in the n:th approximation . . . 19
3.2 An iterative method for evaluating the H-functions . . . 19
3.3 Implementation of the matrix formulation . . . 20
4 A comparison between the two models 23 4.1 An inverse problem formulation . . . 25
4.2 Modifications to the parametric model . . . 26
5 Experiments 29 5.1 The Experimental setup . . . 29
5.2 A comparison between the experimental data and the two models . . . 30
5.2.1 A comparison between RTE solutions and the experimental data . . . 30
5.2.2 A comparison between the experimental data and the modified parametric model . . . 36
6 Discussion and conclusions 38
7 Future work 39
Appendix:
A Convergence of solutions using the Henyey Greenstein phase function I
Chapter 1
Introduction
The theory of light propagation in materials has a huge area of applications. Many applications to this theory use the fact that light propagate differently in different types of materials - That is, the reflected intensity as well as the light distribution inside a material is in some way dependent on material parameters. These material parameters may in turn be estimated using some inverse problem formulation.
For example, these inverse problem formulations may concern finding scattering parameters of human tissue, where the parameters in turn may be used in biomedical imaging applications.
Another possible application is to estimate material parameters such as grain size, humidity and density of snow by just looking at the reflected intensity distribution. These material parame- ters could then be of great help when it comes to avalanche prediction, road safety and friction estimation.
The general theory of the subject that describes light propagation inside materials is the theory of radiative transfer, which gives a description of both absorption and scattering processes in a general medium from a macroscopic perspective. Due to the complexity of the theory, solutions to inverse problem formulations may, however, be very hard to find. The use of parametric models significantly simplifies this matter and one of these parametric models is the model derived by Rosendahl [5].
The Rosendahl parametric model is a two-parameter model, where the intensity is expressed as a linear combination of two cosine functions. The model also has the advantage that the Lambert surface is a special case of the model. That is - by choosing the parameters appropriately, the scattering distribution of a Lambert surface arises.
As of now, no connection between the radiative transfer equation and the parametric model exist. It is the intention of this thesis to investigate whether or not it is possible to find solutions to the radiative transfer equation which matches the parametric distributions. The main idea is to solve the radiative transfer equation for semi-infinite plane parallel cases and find phase functions which matches the parametric distributions.
1.1 Outline
The issue of finding the appropriate phase functions is discussed in four chapters. The first chapter gives a brief theoretical introduction to the radiative transfer equation, as well as the parametric model. The second chapter concerns the issue of implementing the theory and how to solve the equation of radiative transfer numerically. The third chapter gives a first comparison between the parametric model and solutions to the radiative transfer equation, using parameter values from white light measurements on sand grains and glass spheres obtained by Rosendahl [5]. This chapter also presents a possible modification to the parametric model, to better suit radiative transfer approximations.
The last chapter concerns comparisons between frequency dependent intensity distributions obtained from a snow surface, the radiative transfer equation, as well as the parametric model.
1.2 Frequently used variables and abbreviations
RTE - the radiative transfer equation.
HG - the Henyey Greenstein (phase function).
γ - The polar angle in a spherical coordinate system.
φ - The azimuthal angle in a spherical coordinate system.
γ0- The (polar) angle of incidence of the incoming light.
Θ - The angle between the incident and reflected photon.
µ - µ = cos γ.
µ0 - µ0= cos γ0.
µj - The Gauss-Legendre quadrature points.
aj - The Gauss-Legendre quadrature weights.
τ - The optical depth.
p - The phase function.
N - The number of terms in the Legendre expansion of the phase function.
w0- The single scattering albedo.
wi - The i:th coefficient in the Legendre expansion of the phase function.
g - The anisotropy parameter of the Henyey Greenstein phase function.
I - The intensity of interest.
Im- The m:th Fourier component of the intensity.
H(µ) - Chandrasekhar’s H-functions.
Ψ - The characteristic polynomial for a certain H-function.
n, - Model parameters to the parametric model.
λ - Additional model parameter in the modified parametric model.
Chapter 2
Theory
This section is intended to give the reader an introduction to the radiative transfer equation (RTE) and some theory needed for solving the equation using the discrete ordinate method, in the semi-infinite plane parallel case. This section also presents two different implementations of the discrete ordinate method, where the first is the method of expressing the solution in terms of the H- functions, developed by Chandrasekhar [1]. The other method was developed by Stamnes [2, 3] and is a matrix formulation using the same problem formulation and discretization as Chandrasekhar.
The final part of this chapter, presents the parametric model, derived by Rosendahl [5].
2.1 An introduction to radiative transfer
The equation of radiative transfer (or simply, the equation of transfer) is an integro-differential equation and can in its most general form be written as:
−1 a0
dI
ds = I(s) − J(s), (2.1)
where a0 is some, for now undefined constant. The equation gives a relation between the rate of change in intensity (I ) with respect to the spatial variable s, the intensity itself and the source function
J(s) = 1 4π
Z
p(s, s0)I(s, s0)ds0. (2.2) The source function is a summation of intensity contributions from all s0 to s, weighted with some function p. This weighting function is a probability, which gives a macroscopic measure of the probability that an incoming photon from a certain direction s0 is scattered in another direction s.
p(s, s0) is in the following referred to as the phase function.
Before the phase function is properly introduced, a set of ’incoming’ and ’scattered’ angles - that is, a scattering event and a set of spatial variables has to be defined. Chandrasekhar uses a spherical coordinate system to describe the relation between an incoming and scattered photon.
Consider Figure 2.1, and let a surface be defined as the xy-plane of this figure.
Figure 2.1: Definition of angles
A single scattering event is defined as the scattering of an incoming photon, coming from A with the angle γ0relative to the surface normal (the z-axis in Figure 2.1). In a spherical coordinate system, the angle γ is usually referred to as the polar angle. The direction parallel to the surface, φ0, is the angle between the incoming photon relative to some reference direction (in this case, the y-axis). The angle φ is usually referred to as the azimuthal angle. The photon is scattered to B and has another set of angles γ, φ relative to the same normal and reference point as the incoming photon.
The scattering event may also be expressed in terms of the angle Θ between the incoming and scattered photon. The cosine of the angle Θ can in turn be expressed in terms of the angles relative to the surface normal according to
cos(Θ) = cos γ cos γ0+ sin(γ) sin(γ0) cos(φ0− φ). (2.3) It is in this context appropriate to introduce the parameter µ = cos(γ) and µ0 = cos(γ0), which allows the (cosine) angle cos Θ to be expressed as
cos Θ = µµ0+p
(1 − µ2)(1 − µ02) cos(φ0− φ). (2.4) It will later be seen that the solution to the radiative transfer equation itself will be expressed in terms of µ directly and not the angles themselves.
The introduction of µ, may due to the even properties of the cosine function lead to some confusion regarding the directions. Figure 2.2 shows the x-z plane of Figure 2.1, which tells us that when looking in any plane (φ0− φ), all photons heading in the negative z direction have a negative cosine argument (µ < 0). The coordinate system is furthermore chosen such that γ > 0 for x > 0, and γ < 0 for x < 0.
Figure 2.2: the x-z plane of Figure 2.1.
2.2 The phase function
With the definition of the polar and azimuthal angles γ and φ done, a scattering event from a direction (γ0, φ0) to a resulting direction (γ, φ) is occurring with the probability
p = p(γ, φ; γ0, φ0)dφdγdφ0dγ0. (2.5) The probability (density) function, given by 2.5 is a material dependent measure of scattering, known as the phase function. The phase function gives the probability that a photon arriving in the angular window of ([γ0, γ0+ dγ0]; [φ0, φ0 + dφ0]) is scattered into another angular window ([γ, γ + dγ]; [φ, φ + dφ]). Before continuing further, Chandrasekhar expresses the phase function in terms of the angle Θ between the incoming and scattered photon.
p = p(cos Θ)dw0dw,
where dw0 and dw are the appropriate differential element dependent on the angles φ0, γ0 and φ , γ.
As with all probability distributions, the phase function has to be normalized - that is, the phase function has to satisfy
1 4π
Z
Ω
p(cos Θ)d cos Θ = 1. (2.6)
1 4π
Z
Ω
p(cos Θ)d cos Θ = w0≤ 1. (2.7)
The parameter w0 is known as the single scattering albedo, defined as w0= Ps
Ps+ Pa
, (2.8)
where Ps and Pa are the probabilities for a scattering and an absorption event. It can easily be seen from Equation 2.8 that w0 ranges from 0 to 1, where w0 = 1 implies that all photon-atom interactions will lead to scattering. Similarly, w0 = 0 implies that 100% of all interactions will lead to an absorption of the photon. It will later be seen, during the numerical experiments that w0 plays a crucial role when it comes to the behaviour of the solution to the radiative transfer equation.
To solve the radiative transfer equation, the phase function has to be explicitly known, or at least expressed explicitly in some way. Both Chandrasekhar and Stamnes focus on a certain
type of phase functions, namely phase functions which can be expanded as a series of Legendre polynomials in Θ
p =
N
X
l=0
wlPl(cos Θ), (2.9)
where Pl is the Legendre polynomial of order l, and wl is the corresponding weight. It should be noted that w0 = 1 by definition, and that coefficients satisfies wi ≤ 1 i = 1...N . The somewhat strange definition of w0 = 1 follows from that the single scattering albedo is factored out and is already a part of the solution. In theory, the number of terms in the series expansion could reach infinity, but in numerical implementations N will always be finite. Using µ = cos(γ), the argument of the phase function can be expressed in terms of µ and µ0 according to Equation 2.4. The phase function can be developed further by using the addition theorem for spherical harmonics. By doing this, Chandrasekhar writes the phase function in its ’final’ form, where the phase function is expressed in terms of associated Legendre polynomials Plm
p(µ, φ; µ0, φ0) =
N
X
m=0
(2 − δ0,m)[
N
X
l=m
wlmPlm(µ)Plm(µ0)] cos[m(φ0− φ)],
δ0,m=
1 if m = 0 0 otherwise , wlm= wl
(l − m)!
(l + m)!. (2.10)
The introduction of the spherical harmonics expansion of the phase function is necessary and it will later be seen that the expansion in terms of the associated Legendre polynomials, allows for a separation between different Fourier components in the solution. When calculating the terms in the series expansion of the phase function, however, it is easier to only consider the Legendre- polynomial expansion, instead of the associated Legendre-expansion. In the simplest case, the order of the Legendre polynomial is zero, which will result in a phase function consisting of the constant term w0. This case is known as the isotropic case, and corresponds to a uniform spherical probability distribution in µ, φ, µ0, φ0. In the next case, the phase function will be linear with respect to the argument Θ
p1= P0(Θ) + P1(Θ) = w0+ w1cos(Θ) = w0+ w1(µµ0+p
(1 − µ2)(1 − µ02) cos(φ0− φ). (2.11) The linear case, is usually referred to as the linearly anisotropic case. Higher order phase functions will not be shown explicitly here. A common approximation [7, 8] of a more general phase function is the Henyey Greenstein (HG) phase function
p(Θ, g) = 1 − g2
(1 + g2− 2g cos Θ)32. (2.12)
The Henyey Greenstein phase function is a one parameter approximation, where the parameter g determines the degree of anisotropy. A value of g = −1, g = +1 corresponds to complete backward scattering and forward scattering, respectively. A value of g = 0 corresponds to isotropic scattering. It should be noted that the single scattering albedo will still be a part of the solution, it is however factored out from the phase function in this definition. In the corresponding Legendre expansion of the Henyey Greenstein function, the l:th coefficient is given by
wl= gl. (2.13)
Another common phase function is the two-sided Henyey Greenstein phase function. The two- sided HG phase function, is a linear combination of a backward- and forward scattered ’ordinary’
HG phase function
p(Θ, g1, g2, n) = n 1 − g21
(1 + g21− 2g1cos Θ)32 + (1 − n) 1 − g22
(1 + g22− 2g2cos Θ)32. (2.14) The parameter n, determines to what extent the forward and backward scattered part con- tributes to p. Similarily as for the Legendre series expansion of the ’single’ Henyey Greenstein function, the expansion coefficients are given by:
wl= ng1l+ (1 − n)gl2. (2.15)
The idea of having a phase function that is expressible in terms of a backward- and forward scattered part is similar to the concept of the parametric model defined in Section 2.11. The two- sided HG phase function should therefore be a good candidate when trying to unify the Radiative transfer solutions with the parametric model.
2.3 Normalization of the phase function
As stated by Equation 2.6, the phase function has to be normalized to unity. How to do this is not expressed clearly in either the work of Chandrasekhar or Stamnes. I will therefore show how to normalize the phase function, which in turn is equivalent to re-scaling the i:th coefficient wi into a scaled coefficient ˜wi. The unscaled coefficients are assumed to satisfy
0 < w0< 1
−1 ≤ wi≤ 1, i = 1, 2, ..., N. (2.16)
Furthermore, the norm of the phase function can be evaluated by using the fact that the Legendre polynomials satisfies the orthogonality condition
1
Z
−1
Pl(X)Pm(X)dX =
2
2m+1 if l = m,
0 otherwise. (2.17)
By using the orthogonality of the Legendre polynomials, the norm of the phase function can be evaluated according to
< p, p >=
1
Z
−1 N
X
i=0
wiPi N
X
j=0
wjPj =
N
X
j=0
2wj2 2j + 1.
When expressing the phase function as a linear combination of Legendre polynomials, another limitation is that the sum of contributions from each element to p must not exceed unity. Using the orthogonality condition from Equation 2.17 once again, the normalized weights ˜w can be evaluated according to
˜ wi= wi
< Pi, p >
< p, p > = w0
wi3 2i + 1/
N
X
j=0
wj2
2j + 1, i = 0...N. (2.18) That is, if we have chosen a set of weights w0...wN that satisfies Equation 2.16 - the weights has to be normalized according to Equation 2.18. It is then the normalized weights ˜withat should be used in the numerical implementation.
2.4 The optical depth
Another definition necessary, is the definition of the optical depth τ , defined as
τ =
∞
Z
z
κρdz. (2.19)
The parameter ρ is the density of the material and κ is the mass absorption coefficient. The optical depth is defined from the inward normal of Figure 2.3, making it positive in the downward direction. With the definition of the optical depth, we now have three spatial coordinates τ, µ, φ.
Using these we can write the radiative transfer equation on the form µdI(τ, µ, φ)
dτ = I(τ, µ, φ) − J(τ, µ, φ). (2.20)
Note the change in sign when comparing Equation 2.1 with 2.20 due to the definition of the optical thickness. Another difference resulted by the introduction of the optical depth is the removal of the constant a0 from equation 2.1.
2.5 Problem formulation
With the definition of the optical depth, single scattering albedo and phase function done, the radiative transfer problem in question for this thesis can now be formulated more properly.
In the plane parallel case, a beam of intensity I0 excites the system at a (cosine) angle µ0 = cos γ0 relative to the surface normal of the plane, see Figure 2.3. The beam is assumed to be monochromatic, and spatially as well as temporally incoherent. The intensity distribution at an arbitrary optical depth τ will then be the solution to the equation
µdI(τ, µ, φ)
dτ = I(τ, µ, φ) − 1 4π
1
Z
−1 2π
Z
0
p(µ, φ; µ0, φ0)I(τ, µ0, φ0)dµ0dφ0−I0
4e−µ0τ (2.21) where µ = cos γ. The source function has here been separated into two terms. The first term is the ’summation term’, described in the beginning of Section 2.1 and involves light contributions that suffer from one or more scattering processes at an optical depth τ . The other term is a ’direct term’, which corresponds to the directly transmitted light at the same depth τ in the material in question.
The boundary conditions stated by Chandrasekhar [1, p.22] for the plane-parallel case for a slab of optical thickness τ1 are
I(0, µ, φ) = 0, −1 < µ < 0
I(τ1, µ, φ) = 0, 0 < µ < 1.
This means that the above boundary conditions does not take the inward ’excitation’ of the beam into account, a more proper way of writing the boundary conditions are therefore
I(0, µ, φ) =
I0 if µ = µ0,
0 otherwise , −1 < µ < 0 (2.22) I(τ1, µ, φ) = 0, 0 < µ < 1 (2.23) When writing the boundary conditions on this form, it is clear that the boundary conditions are essentially a condition of continuity, saying that the radiation distribution in the downward direction is the same above and under the two planes. A further note is that the second boundary condition (Equation 2.23) assumes that no radiation source is present under the plane at τ = τ1, neither does it take direct reflections into account, since the intensity distribution in the positive
direction at τ = τ1 is forced to be zero. This is however of no importance for the case when τ1→ ∞.
It is in this context also necessary to remind the reader of the difference between the coordinates µ, τ and z. The angles µ are defined such that they are positive in the positive z-direction, negative in the negative z-direction and 0 parallel to the boundary, see Figure 2.3. The optical depth τ is positive in the negative z-direction.
Figure 2.3: A schematic description of a plane-parallel geometry.
2.6 Problem discretization
The work of both Stamnes and Chandrasekhar use the discrete ordinate method for discretizing the problem. The discrete ordinate method is based on the fact that an integral can be approximated (or in some cases, evaluated exactly) according to
b
Z
a
f (µ)dµ =
b
Z
a
g(µ)w(µ)dµ ≈
n
X
j=1
w(µj)f (µj). (2.24)
That is, by choosing appropriate weights, wj and discretization points µj, the integral can be approximated. Chandrasekhar shows [1, ch. 2] that an appropriate choice of points and weights, are the weights given by the Gauss quadrature. In the Gauss quadrature, the discretization points (µj) are given by the roots of a Legendre polynomial. Chandrasekhar argues and finds that dis- cretization points of even order has to be used, and that a discretization in the n:th approximation are given by the roots to a polynomial of order 2n. Furthermore, the weights are given by [1, ch.2 ]:
aj= 1 P2n0 (µj)
1
Z
−1
P2n(µ) µ − µj
dµ. (2.25)
Sykes showed [4] that a better choice of quadrature points is to choose weights and points on the ’half interval’ [0 1], and then mirror the points for the negative part. These were also the points and weights that were adopted by Stamnes and Edstr¨om [3, 7]. In the case of the numerical implementation and theoretical part of this thesis, the points and weights used, are chosen such that they are the same in the both implementations and are therefore evaluated on the full interval in both cases.
With the problem discretized, the solution is assumed to be expressible in terms of a Fourier series expansion
I(τ, µj, φ) =
N
X
m=0
Im(τ, µj) cos[m(φ − φ0)]. (2.26)
Doing this, the integral equation can be reduced to systems of linear differential equations with constant coefficients [3], [1, ch. 3&6] - one for each Fourier component. The advantage with rewriting the system as a system of differential equations is that it can be treated as an eigenvalue problem. The following sections present two methods of solving the corresponding eigenvalue problem, The first is the method by Chandrasekhar, where he derives the characteristic equation of the source function and expresses the solution in terms of H functions. The second method is the matrix formulation by Stamnes.
2.7 The scattering function
The method used by Chandrasekhar (and indirectly Stamnes) when solving equation 2.21, is to rewrite the problem in terms of a scattering- and a transmission function
S = S(τ1; µ, φ; µ0, φ0) (2.27)
T = T (τ1; µ, φ; µ0, φ0) such that the reflected and the transmitted intensities are given by
I(τ = 0, µ, φ) = I0
4µS(τ1; µ, φ; µ0, φ0), 0 < µ < 1 (2.28) I(τ = τ1, µ, φ) = I0
4µT (τ1; µ, φ; µ0, φ0), −1 < µ < 0. (2.29) In the semi-infinite case, however, only the Scattering function is of interest.
In the method used by Chandrasekhar, the solution to Equation 2.21 is obtained by solving a set of independent equations for a series of scattering functions, and then superimposing the result
I(0, φ, µ) = I0 4µ
N
X
l=0
Sl(µ; µ0) cos[l(φ − φ0)]. (2.30)
Chandrasekhar then shows that each Scattering function is in itself an integral equation, which satisfies
(1 µ+ 1
µ0)S(µ, µ0) = H(µ)H(µ0). (2.31)
H are known as H-functions which satisfy the non-linear integral equation
H(µ) = 1 + µH(µ)
1
Z
0
Ψ(µ)
µ + µ0H(µ0)dµ0, (2.32)
where Ψ is a characteristic polynomial in µ. That is - by finding the proper characteristic polyno- mials for a certain phase function, the integral equation 2.32 can be solved, which in turn can be used to express the intensity distribution at the top boundary.
2.8 Evaluation of the H-functions
It has become evident that solving the equation of radiative transfer, (using Chandrasekhar’s method) falls down to Solving the integral Equation 2.32 for the H-functions. The properties of the H-functions and how to solve the integral equation, which determines their behaviour covers an entire chapter in Chandrasekhar’s book but can be condensed into the following:
The H functions can be approximated to an order n, as a rational polynomial in µ [1, p.105]
H(µ) = 1 µ1...µn
n
Q
i=1
(µ + µi) Q
α
(1 + kαµ), (2.33)
where µiare the zeros of P2n(that is, the quadrature points), and kαare the positive, non-vanishing roots to the characteristic equation
1 = 2
n
X
j=1
ajΨ(µj)
1 − k2µ2j. (2.34)
The coefficients aj are the quadrature weights, given by Equation 2.25, µj are the corresponding quadrature points and Ψ is the characteristic polynomial, associated with the corresponding H- function.
Another way of evaluating the H-functions is by using Chandrasekhar’s iterative method [1, p.123-24]. In this method, he makes use of the fact that the H-functions can be written on the form [1, ch.5 eq. 13]
1
H(µ) = [1 − 2
1
Z
0
Ψ(µ)dµ]12 +
1
Z
0
µ0Ψ(µ0)
µ + µ0 H(µ0)dµ0. (2.35) By choosing a starting guess H0 and then evaluate the integral numerically, the values of the H-functions can be evaluated to an arbitrary precision.
2.9 The characteristic polynomial, Ψ
The final step in obtaining the solution to the radiative transfer equation (using Chandrasekhar’s method) is to find the characteristic polynomial Ψ. Using the invariance in the spatial variables together with some clever mathematical tricks, he finds that Ψ is given by Ψ =w20 for the isotropic case. For the linearly anisotropic case, he defines two characteristic functions
Ψ0= 1
2w0[1 + x](1 − w0)µ2, Ψ1=1
4xw0(1 − µ2), (2.36) where w0x = w1, and w0,1 are the Legendre expansion coefficients. The solutions can now, finally be expressed using Equations 2.30 and 2.31 [1, ch. 6]. For the isotropic case, the solution is given by
I(0, µ, φ) = I0µ0
4(µ + µ0)H(µ)H(µ0). (2.37)
For the linearly anisotropic case, the solution is given by I(0, µ, φ) = I0µ0w0
4(µ + µ0)[ψ(µ)ψ(µ0) − xΦ(µ)Φ(µ0)
+(1 − µ2)1/2H1(µ)(1 − µ20)12H1(µ0)] cos(φ0− φ) (2.38) where ψ = H0(µ)(1 − cµ) and Φ(µ) = qµH0(µ), and H0,1is defined by Ψ0,1 in Equation 2.36. The constants q and c are in turn given by
q = 2(1 − w0) 2 − w0α0
, c = xw0(1 − w0)(2 − w0
α1
2 − w0α0
), (2.39)
where α0,1 are moments defined in terms of the H functions, defined as:
αn =
1
Z
0
µnH0(µ)dµ, n = 0, 1. (2.40)
The solution to the radiative transfer equation is visualised in Figure 2.4, for the isotropic (left) and linearly anisotropic case (right). The method used for evaluating the H-functions was the iterative method presented in section 3.2. In both cases, the angle of incidence was µ0= −36◦ and all angles are measured relative to the surface normal. Note that the isotropic case is a special
case of the linearly anisotropic case, by letting x=0 in Equation 2.38, the equation will be equal to 2.37.
Figure 2.4: The solution to the radiative transfer equation for the isotropic (left) and linearly anisotropic case (right), for different single scattering albedo w0, with an illumination angle of γ = −36◦. In the anisotropic case x =0.9 was chosen.
The previous sections has showed how the H-functions can be used to express the solution for the radiative transfer equation, primarily for the isotropic and linearly anisotropic case. The issue of numerically solving the integral equation for H-functions is a different matter, and is treated in section 3.1 and 3.2.
2.10 A matrix formulation
The way of expressing the solution to the Radiative Transfer equation in terms of H-functions has several numerical drawbacks [9]. It is therefore not surprising that the use of Chandrasekhars H- functions is very limited. The theory, based on discretizing the equation and ’rewriting’ the problem as an eigenvalue problem is quite general and has been implemented as a matrix formulation by Stamnes [3]. Note that this section only gives a theoretical introduction to the matrix formulation, a major part of the work in this thesis has been the issue of numerically implementing the matrix formulation for an arbitrary phase function with N (given) coefficients. The implementation itself is discussed in Section 3.3.
Stamnes uses a slightly different notation, compared to Chandrasekhar. The two methods have assumptions regarding the phase function and intensity, namely that the phase function can be expressed in terms of (associated) Legendre polynomials and that the intensity can be expressed as a Fourier cosine series. Furthermore, the two methods are both based on a discretization which uses a numerical quadrature to approximate the integral. With these assumptions, Equation 2.21 can be written as
µdum(τ, µ)
dτ = um−
1
Z
−1
Dm(µ, µ0)um(τ, µ0)dµ0+ Qm(τ, µ), (2.41)
where umis the Intensity for the m:th Fourier component in Equation 2.26. Dm is the associated Legendre expansion of the phase function for the m:th Fourier component
Dm= w0
2
N
X
l=m
(2l + 1) ˜wlmPlm(µ)Plm(µ0), (2.42)
Qm describes the part of light that has reached the depth τ without being subject to any scattering interactions (as discussed briefly in Section 2.5) and is given by:
Qm= X0m(τ, µ)e−µ0τ , (2.43)
X0m=I0w0
4π (2 − δ0,m)
N
X
l=m
(2l + 1) ˜wml Plm(µ)Plm(µ0). (2.44)
˜
wlm= ˜wl(l − m)!
(l + m)!, (2.45)
where wl is the normalized Legendre coefficients (derived in Section 2.3). By discretizing the problem in the same way as described in Section 2.6, the differential equation can be written on the form [3, p.2,3]
du+ dτ du−
dτ
=
−α −β
β α
u+ u−
+
Q˜+ Q˜−
, (2.46)
where
Q˜±= Qm(τ, ±µi)
M = µiδi,j
W = aiδi,j
α = M−1(D+W − I)
β = M−1D−W
D± = Dm(±µi, µj) (2.47)
By formulating the problem in this way, the radiative transfer equation has been reduced to the problem of solving a system of N coupled differential equations. The theory for this type of problem is well known, where the solution can be obtained by solving for the homogeneous (Qm= 0) and particular case separately. For the homogeneous case, the solution will be on the form
u±= G±e−kτ, (2.48)
where k are the eigenvalues, with corresponding eigenvectors G satisfying
α β
−β −α
G+ G−
= k
G˜+ G˜−
. (2.49)
As seen from above, the eigenvectors are ’divided’ into two parts, where G± corresponds to the intensity distribution in the positive and negative z-direction in Figure 2.1, respectively. This formulation has the advantage that existing numerical methods can be used to solve the system of equations that arise.
Furthermore, the particular solution is given by [7, p.456]
Ipm(τ, µi) = Z0m(µi)e−µ0τ , (2.50) where Z0 is evaluated according to
n
X
j=−n
((1 +µj µ0
)δi,j− ajDm(µj, µi))Z0m(µj) = X0(µi) (2.51) and aj are the quadrature weights.
The general solution (for each Fourier component) can then be expressed as a linear combination of the different ui plus the particular solution
Im(τ, µi) =
N
X
j=−N
CjGj(µi) exp(−kjτ ) + Ip(τ, µi), (2.52) where Cj are coefficients to be determined by the boundary conditions. The boundary conditions are formulated in Section 2.5 and imply that
N
X
j=1
CjGj(−µi) + C−jG−j−µi= δµi,µ0− Ip(0, −µi), 0 < µi< 1 (2.53)
N
X
j=1
CjGj(µi)e−τ /µ0+ C−jG−j(µi)eτ /µ0 = 0, 0 < µi< 1. (2.54) By solving the system of equations that arise from these conditions, the coefficients can be determined and the solution can be evaluated which is described in more detail in section 3.3.
The solution to the radiative transfer equation is plotted for the linearly anisotropic case (left) with the same parameters as in Figure 2.4, to illustrate that the two methods only are different implementations of the same theory. The figure to the right, shows the case when the Henyey Greenstein phase function is used with N=40 terms approximating the phase function and w0 = 0.99. Note the significant difference in behaviour when the anisotropy parameter g=-0.9 compared to g=0.9. When g=-0.9, a majority of the light is backscattered in the same direction as the incident beam, while forward scattering parameters gives a more ’smooth’ behaviour. Also note the drop in intensity that is seen when the scattered angle approaches the surface normal.
Even though the matrix formulation is more simple to implement numerically in a programming language such as Matlab, several numerical issues related to the formulation still exist. When the theory above was implemented, four main problems arose. The first problem concerned the evaluation of the eigenvalues and eigenvectors, and the ordering of these. The second problem concerned the evaluation of the associated Legendre polynomials and the third problem was the issue of numerically evaluating the coefficients for the homogeneous solution. The last issue was the problem of finding the quadrature points and weights for the Gauss quadrature. How these problems were handled, and how the solver was implemented, are described in Section 3.3.
Figure 2.5: Solution for the linearly anisotropic case (left) when x=0.9. The figure to the right is the solution for the Henyey Greenstein phase function when w0= 0.99. The illumination angle is in both cases γ0= −36◦.
2.11 The parametric model
As seen from the previous sections, the mathematics behind the radiative transfer equation is far from simple. Another problem with the general theory is that inverse problem formulations become quite complex. With this in mind, Rosendahl formulated a parametric model [5] which more easily could be used for inverse problem formulations. The behaviour of the parametric model has been experimentally verified.
In the parametric model, the intensity I is assumed to be expressible as a linear combination of a forward (C+) and backscattered (C−) intensity distribution
I = (1 − n)C++ nC−, (2.55)
where n is a constant determining the degree of forward, and backward scattering, respectively.
The forward and backscattered intensity distributions are in turn given by
C±= cos α, (2.56)
where the argument α is a skewed version of an ’ordinary’ angle x relative to the surface normal.
This skewness satisfies
α = x + γ cos α, (2.57)
where γ is the angle for the beam of incidence. The angle γ does furthermore coincide with the γ used in the radiative transfer equation. x is the ’ordinary’ angle, see Figure 2.6 and is a model parameter. That is - by solving Equation 2.57 numerically the reflected angle α can be obtained.
An interesting remark is that the model converges to a Lambert surface for n = 0.5 and = 0.
The behaviour of the parametric model is visualized in Figure 2.7 where the parameters are chosen such that they have been experimentally verified [5] to model a reflection distribution from moistened glass spheres.
Figure 2.6: α, x and γ.
Figure 2.7: The intensity I with corresponding forward- and backward scattering functions C±for the case when = 0.84, n = 0.56 and γ0= 40◦
Chapter 3
Numerical implementation
This section covers two possible ways of evaluating the H-functions numerically, one method is by ”directly” solving the characteristic equation 2.34. The other method is an iterative method, based on an iterative scheme by Chandrasekhar, modified by Bosma and Rooij [6]. In addition, this section also covers the implementation of the Matrix formulation. All code has been written in MATLAB.
3.1 Evaluation of the H-functions in the n:th approximation
As stated in Section 2.8, the H-functions in the n:th approximation are dependent on the roots kα
to the characteristic equation
1 = 2
n
X
j=1
ajΨ(µj) 1 − k2µ2j.
The roots can be evaluated by first expanding the characteristic equation 1
2 = a1Ψ(µ1) 1 − k2µ1
+ a2Ψ(µ2) 1 − k2µ2
+ · · · + anΨ(µn) 1 − k2µn
, (3.1)
and then remove the denominator by multiplication. By doing this, the equation can be rewritten as a polynomial equation
n
X
i=1
aiΨ(µi)
n
Y
j=16=i
(1 − k2µj) −1 2
n
Y
j=1
(1 − k2µj) = 0, (3.2)
where µj, aj are the quadrature points and weights. The roots to equation 3.2, can then be found using some numerical method. In the case of the implementation in this thesis, matlabs symbolic engine was used to expand and simplify the expression, where after roots.m could be used to obtain the different kα. One should point out, however, that this method is only applicable when the number of disretization points n is not ’too large’. When the number of discretization points is too large, the inaccuracy in roots.m will lead to complex roots (which in this case gives a sign that something has gone wrong) and other numerical problems.
3.2 An iterative method for evaluating the H-functions
The success in expressing the solution to the RTE in terms of H-functions (in the n:th approxima- tion), lies in the success of evaluating the kαcorrectly. This problem, however, can be avoided, by implementing the iterative scheme suggested by Chandrasekhar and improved by Rooij [6].
The iterative method is based on the relation given by Chandrasekhar
1
H(µ) = [1 − 2
1
Z
0
Ψ(µ)dµ]12 +
1
Z
0
µ0Ψ(µ0)
µ + µ0 H(µ0)dµ0. (3.3) Rooij then suggests an iterative scheme given by
Hn= λHn(µ) + (1 − λ)Hn−1(µ) (3.4)
where
λ = 1 2(1 +p
1 − 2ψ0),
ψ0=
1
Z
0
Ψ(µ)dµ. (3.5)
The iterative scheme can then be summarized into the following
• Make an initial guess Hstart.
• Compute 1/Hnext with Equation 3.3, using numerical integration with Hstart as H.
• Compute new 1/Hstart using the two previous values of H according to Equation 3.4.
• check convergence criterion
• jump to 2nd step if convergence is not met.
The convergence criterion used was the same as the one used by Rooij and is defined as
max[|Hi− Hi−1|] < δ. (3.6)
That is, the absolute value of the maximum difference between two successive iterates should be less than some tolerance δ.
3.3 Implementation of the matrix formulation
The matrix formulation summarized in section 2.10 was implemented in Matlab. Here follows a brief summary on how the numerical issues related to the implementation was handled.
The problem of obtaining the eigenvalues and eigenvectors for the problem formulation in this thesis have been thoroughly examined by several researchers, and the results are summarized by Stamnes in [3]. When writing RTEsolver.m, however the eigenvalue problem given by Equation 2.49 was solved directly using Matlabs eig.m, since the algorithms of this scripts seemed accurate enough.
A remark regarding the eigenvalues and eigenvectors is that the eigenvectors is on the form G = [G+G−]T, where the basis of the eigenvectors are the discretization points µi. The basis and the ordering of the different µi is crucial for obtaining the correct solution, and in this implementation an eigenvector Gi with eigenvalue λi is ordered according to
Gi=G(µn) G(µn−1) · · · G(µ1) G(µ−n) G(µ−n+1) · · · G(µ−1)T
. (3.7) When evaluating the associated Legendre polynomials, the first type of method used was Mat- lab’s built in function legendre.m. This method, however has the (quite severe) disadvantage that it is recursive with respect to m. That is - for a fixed l, all m are evaluated, even though only one term is needed as a part of the solution. Edstr¨om presents a recursive method [7] of evaluating normalized Legendre polynomials
Λml (µ) = s
(l − m)!
(l + m)!Plm(µ) (3.8)
which is recursive in l, where the recursive relation is given by Λ00= 1
Λmm= − r
(1 − µ2)(2m − 1
2m )Λm−1m−1(µ). (3.9)
By expanding Equation 3.9 m times, an explicit expression for Λmmcan be obtained
Λmm= (−1)m(1 − µ2)m2 (2m − 1)!!
2m!!
m2
. (3.10)
The !! operator is known as double factorial and is defined as 2m!! = 2 · 4 · ... · 2m (2m − 1)!! = 1 · 3 · ... · (2m − 1)
for even and odd numbers, respectively. The normalized associated Legendre polynomial can then, for an arbitrary l, m be evaluated by using the following formulas
Λmm+1(µ) = µ√
2m + 1Λmm(µ)
Λml (µ) = µ(2l − 1)Λml−1(µ) −p(l − 1 + m)(l − 1 − m)Λml−2(µ)
p(l − m)(l + m) . (3.11)
The implementation of the normalized Legendre polynomials resulted in a significant speed up when evaluating the sum in Equation 2.42 (roughly 20-200 times faster, depending on the order).
It also made the evaluation of the factorial-expression in Equation 2.10 unnecessary.
The numerical difficulty concerning the evaluation of the boundary conditions arise due to the positive exponential in Equation 2.54. A possible solution to this problem was presented by Stamnes [3] and involves the introduction of a scaling transformation
C−j = ˜C−jexp(−kjτp), (3.12)
where τp is the optical thickness of the medium. Stamnes argues that the scaling transformation only affects the homogeneous solution, which implies that the solution to the homogeneous part can be written as
um(τ, µi) =
N
X
j=1
CjGj(µi)exp(−kjτ ) + ˜C−jG−j(µi)exp[−kj(τp− τ )]. (3.13) By using the scaling transformation given by Equation 3.12, the coefficients can be obtained by rewriting Equation 2.54 as a system of linear equations
Gn(µn) Gn−1(µn) · · · G1 G−n G−n+1 · · · G−1(µn) Gn(µn−1) Gn−1 · · · G1 G−n G−n+1 · · · G−1(µn−1)
..
. . .. · · · ... .. .
.. .
.. .
.. . ..
. . .. · · · ... .. .
.. .
.. .
.. . Gn(µ1) Gn−1 · · · G1 G−n G−n+1 · · · G1(µ1) Gn(−µn) Gn−1 · · · G1 G−n G−n+1 · · · G1(−µn) Gn(−µn−1) Gn−1 · · · G1 G−n G−n+1 · · · G−1(−µn−1)
..
. . .. · · · ... .. .
.. .
.. .
.. . ..
. . .. · · · ... .. .
.. .
.. .
.. . Gn(−µ1) Gn−1 · · · G1 G−n G−n+1 · · · G1(−µ1)
Cn
Cn−1
.. . .. . C1
C˜−n
C˜−n+1
.. . .. . C˜−1
=
0 0 .. . .. . 0
I0δµi, µ0− Z0m(µ−n) I0δµi, µ0− Z0m(µ−n+1)
.. . .. .
I0δµi, µ0− Z0m(µ1)
. (3.14)
In the system of linear equations above, each row consist of eigenvector values for a fixed µi. The I0δµi,µ0 term on the right hand side of Equation 3.14 = I0 when µi = µ0 and 0 otherwise.
With the problem of obtaining the correct coefficients written on this form, the coefficients can be found easily using Matlab’s \-operator.
The quadrature points and weights were first calculated using symbolic expressions for the Legendre polynomials. By doing this, the points and weights could be obtained quite easily, simply by implementing the expression given by Equation 2.25 together with roots.m. As in the case when evaluating the roots to the characteristic polynomial for the H-functions, the accuracy was so low that the solution to the eigenvalue problem of interest (Equation 2.49) will have complex eigenvalues when more than 20 discretization points are used. It is however possible to rewrite the problem of finding the Gauss quadrature points and weights as an eigenvalue problem [10]. The code to this type of problem is available on-line via Matlabs file-exchange, in numerous versions.
In this case, the code used was given by [11].
Chapter 4
A comparison between the two models
The first type of comparisons done was an examination on how the solution for the isotropic and lin- early anisotropic case behaved for different Legendre coefficients w0and w1using Chandrasekhar’s H-functions. This first comparison was done, simply to get a grasp on how the parameters in the equation affected the solution. The idea was then to see if any of the isotropic or linearly anisotropic cases could represent the measured scattering distributions obtained by Rosendahl.
As seen in Section 3.1 and 3.2, the H-functions (and therefore the RTE solution) can be ex- pressed both in terms of a rational polynomial in µ
H(µ) = 1 µ1...µn
n
Q
i=1
(µ + µi) Q
α
(1 + kαµ),
where kαare the roots to the characteristic equation 2.34 and µiare the Gauss quadrature points.
The other option is to evaluate the H-functions iteratively. Both these methods has their own advantages and drawbacks. The first method of expressing the solution in terms of a rational polynomial in µ would make it relatively easy to identify similarities in the coefficients, due to the finite number of terms in the approximation. More exactly, if a low order approximation of the H-functions show the same tendencies as the parametric model, a closed form expression which unifies the two models may be possible. The second method, enables more accurate results since the solution to the RTE can be evaluated with higher accuracy, but makes any equality in ’closed form’ extremely difficult (if not impossible).
The idea of using a low order approximation in the H-functions to approximate the RTE solution and see if the approximate solution show similar behaviour as the parametric model, of course, only works if they actually show similar behaviour. To see whether or not they exhibited similar behaviour, the order in the approximation was varied for both the isotropic and linearly anisotropic case. Figure 4.1 shows the solution for the isotropic case when the maximum value of the distribution has been normalized to unity, together with the parametric model (normalized to unity). In order to increase the visibility, the different curves have been shifted after the normalization was done.
Figure 4.1 shows that the order in the approximation, does not dramatically affect the behaviour of the solution for the isotropic case. There is however a fundamental problem, namely that the isotropic intensity distribution is far too ’flat’ to exhibit similar behaviour.
A similar investigation was done for the linearly anisotropic case, see Figure 4.2. In the linearly anisotropic case, the solution to the radiative transfer equation is more dependent on the order of the approximation. As seen from the figure, there is a significant difference between the first order and higher order approximations. The main issue, however, still remains, namely that the radiative transfer solutions are too flat to describe the same phenomenon as the parametric model.
Figure 4.1: A comparison between different orders in the approximation of the isotropic solution and the parametric model. The parameters used for the isotropic solutions are w0 = 1 and γ = −40◦. The parametric model is plotted for the case when γ = 0.84, n = 0.56 and θ = −40◦
.
Figure 4.2: A comparison between different orders in the approximation of the linearly solution and the parametric model. The parameters used for the isotropic solutions are w0= 1, x = −1 and γ = −40◦. The parametric model is plotted for the case when γ = 0.84, n = 0.56 and θ = −40◦
Even though only two specific cases have been shown above, these plots show that a more general phase function, as well as a more general method of comparing the methods is needed.
Furthermore, if one takes a closer look at Figure 2.4, one sees that the single scattering albedo w0
