Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Numerical modeling of the radiative transfer in a turbid medium using the synthetic iteration

Vladimir P. Budak, Gennady A. Kaloshin, Oleg V. Shagalov, and Victor S. Zheltov

Open Access Open Access

Abstract

In this paper we propose the fast, but the accurate algorithm for numerical modeling of light fields in the turbid media slab. For the numerical solution of the radiative transfer equation (RTE) it is required its discretization based on the elimination of the solution anisotropic part and the replacement of the scattering integral by a finite sum. The solution regular part is determined numerically. A good choice of the method of the solution anisotropic part elimination determines the high convergence of the algorithm in the mean square metric. The method of synthetic iterations can be used to improve the convergence in the uniform metric. A significant increase in the solution accuracy with the use of synthetic iterations allows applying the two-stream approximation for the regular part determination. This approach permits to generalize the proposed method in the case of an arbitrary 3D geometry of the medium.

© 2015 Optical Society of America

1. Discretization of radiative transfer equation

The current state of the radiative transfer theory is focused on the development of numerical methods for the solution of the radiative transfer equation (RTE) in a medium with arbitrary geometry. This opens the way to the solution of a significant number of practical problems, among which a special place is occupied by the problem of propagation of laser beams in the atmosphere and the ocean [1–3]. It seems to us that the best way to construct such an algorithm is a generalization of the RTE solutions for planar geometry. There are the developed numerical methods and an exact solution of RTE in this area that allows determining the most effective way to solve the problem in arbitrary geometry.

In the paper [4] comparisons of the main numerical methods of the RTE solution for the slab is carried out. The basic difference of all the algorithms for the solution is connected with the difference methods of the account a strong scattering anisotropy and the presence of the δ-singularity in the radiance angular distribution (RAD) near incidence direction. The obtained coincidence of results of different methods used for the numerical solution of RTE with machine precision shows that the existing algorithms efficiently solve the problem of the strong scattering anisotropy. Moreover, it is obvious that such a coincidence is only possible when using the same RTE solution per se. However, this comparison shows a significant difference (more than two orders of magnitude) in the computation time of these algorithms, indicating a significant distinction in the algorithm implementations.

In [5], a study of numerical methods of solving the RTE and algorithms for their implementation is carried out. To solve RTE by a numerical method, the scattering integral should be replaced by a finite sum that is impossible to the strong anisotropy and the presence of singularities in the radiance angular distribution. The physical basis of RTE is the ray approximation that inevitably generates singularities in the radiance angular distribution that is evident in the sharp boundary between light and shadow. Therefore, the numerical solution of RTE must be done, above all, to eliminate analytically from RAD L(τ,μ,φ) its anisotropic part La(τ,μ,φ), including all the singularities that leads to the formulation of the boundary value problem for the smooth solution regular part [6]:

{μ∂Lr(τ,l^)∂τ+Lr(τ,l^)=Λ4π∮x(l^,l^′)Lr(τ,l^′ ) dl^′+Δ(τ,l^);Lr(τ,  l^ )|τ=0, μ>0=0;Lr(τ, l^ )|τ=τ0, μ<0=−La(τ,l^),
where L(τ,l^)=La(τ,l^)+Lr(τ,l^) is the radiance of the light field in the medium at the optical depth τ in the direction l^. We use Cartesian coordinate system OXYZ where axis OZ is directed down perpendicularly to the slab border:, l^0={1−μ02,  0,  μ0}, where μ=(l^,z^), μ0=(l^0,z^), z^ is the unit vector along OZ. x(l^,l^′) is the phase function; Λ is the single scattering albedo. The unit vectors have a symbol “^”, the columns have the right arrow, the rows have the left arrow, the matrices have the double arrow above the symbol. All notations were taken from [5].

The source function on the right side of the Eq. (1) is determined by the residual of the approximate solutions La(τ,μ,φ) that describes the anisotropic part of the RTE solution [5]. The second boundary condition is introduced because the approximation La(τ,μ,φ) may not satisfy the exact boundary conditions.

The source function is given by

Δ(τ,μ,φ)=−μdLa(τ,μ,φ)dτ−La(τ,μ,φ)+Λ4π∮x(l^,l^′)La(τ,μ′,φ′)dl^′.

The direct estimation of expression (2) is impossible since the solution anisotropic part contains the δ-singularity of the angular distribution. To eliminate this problem, we isolate obviously the angular singularity – the direct non-scattered radiation:

La(τ,μ,φ)=e−τ/μ0δ(l^−l^0)+L˜a(τ,μ,φ),
that reduces Eq. (2) to the expression

Δ(τ,μ,φ)=−μdL˜a(τ,μ,φ)dτ−L˜a(τ,μ,φ)+Λ4π∮x(l^,l^′)L˜a(τ,μ′,φ′)dl^′+Λ4πe−τ/μ0x(l^0,l^).

L˜a(τ,μ,φ) is the regular anisotropic part of the solution. It is anisotropic by angle, but the smooth function, therefore the source function in Eq. (4) is also a smooth function and does not present any problem for the numerical calculation.

Note that since Lr(τ,μ,φ) is the quasi-isotropic part of the solution, and could be represented as a finite set of elements (basis functions). For example, if the discrete ordinates method (DOM) is used for the RTE discretization this representation takes a form:

L→r(τ,μi±,φ)=∑m=0M(2−δ0,m)cos(mφ)C→m(τ,μi±),  C→±(τ)≡C→m(τ,μi±),
where C→(τ)≡[C→−(τ),  C→+(τ)]T are the columns of the expansion coefficients of radiance into Fourier series in azimuth, μj±=0.5(ζj±1), ζj are zeroes of the Gaussian quadrature of the N/2 order.

This allows replacing the scattering integral by the finite sum. Thus the boundary value problem leads to the system of inhomogeneous linear differential equations of the first order with constant coefficients [5]:

dC→(τ)dτ=−B↔C→(τ)+M↔−1Δ→(τ),  B↔≡M↔−1(1↔−A↔W↔),
where M↔≡diag(μi+,μi−), W↔≡0.25Λdiag(wi,wi) are the Gaussian quadrature weights, x(γ)=∑k=0K(2k+1)xkPk(cosγ), Pkm(μ), Pk≡Pk0(μ) are the associated Legendre polynomials, A↔≡[∑k=0K(2k+1)Pkm(μi±)xkPkm(μj±)].

One can see that any expression of RAD over a finite basis, as DOM, and the method of spherical harmonics is reduced to the representation Eq. (5). Note that any numerical solution method of RTE leads to this resulting representation. In this sense, the method of successive scattering order is an iterative method of solving the eigenvalue problem, and the Monte Carlo is a statistical one. For definiteness, we will use DOM everywhere further.

2. Solution regular part

The solution of Eq. (5) is represented [5, 6] by the sum of the general solution of homogeneous equation and the particular solution of the inhomogeneous one:

−C→(0)+P↔(0,τ0)C→(τ0)=∫0τ0P↔(0,τ)M↔−1Δ→(τ,μ0)dτ
where P↔(t,τ)≡eB↔(τ−t) is a solution of the homogeneous equation – a propagator.

The problem of the solution Eq. (6) is connected with negative and positive exponents in the expression that leads to fast worsening of the system matrix conditions. To eliminate this effect, we used scale transformation [5, 6]:

−S↔U↔−1C→(0)+H↔U↔−1C→(τ0)=S↔∫0τ0eΓ↔tU↔−1M↔−1Δ→(t,μ0)dt,
where the matrix exponent is represented in the form
eB↔τ0=U→e−Γ↔τ0U→−1,
U→ is the matrix of eigenvectors of the matrix B↔; Γ↔=diag(Γ↔−,  Γ↔+) is the diagonal matrix of eigenvalues sorted by ascending ordering so that Γ↔+=−Γ↔−≡[γi], S↔=[1↔00e−Γ↔+τ0],  H↔=[eΓ↔−τ0001↔].

The columns C→+(0),  C→−(τ0) in Eq. (7) are the streams incident on the slab and are determined by the boundary conditions; C→−(0),  C→+(τ0) are the desired solutions that determine the reflected and transmitted streams outgoing from the slab. We solve Eq. (7) relatively to the streams outgoing from the slab and get the solution in the form of scatterers [5, 6]:

[C→−(0)C→+(τ0)]=[F→−F→+]+[R↔T↔T↔R↔][C→+(0)C→−(τ0)]
where [F→−F→+]=h↔J↔,  J↔≡S↔∫0τ0eΓ↔τU↔−1M↔−1Δ→(τ,μ0)dτ,  [R↔T↔T↔R↔]=h↔[u↔11−eΓ↔−τ0u↔12e−Γ↔+τ0u↔21−u↔22]h↔≡[−u↔12eΓ↔−τ0u↔11−e−Γ↔+τ0u↔22u↔21]−1,  U↔−1≡[u↔11u↔12u↔21u↔22].

The solution in the form of scatterers (9) gives the relation between radiation emerged by the slab with incident one on its borders. The advantage of this solution representation is its physical clarity: the radiation outgoing from one slab boundary is connected with the reflection of the radiation incident on this boundary, the transmittance of the radiation incident on another boundary and the slab self-irradiation. Since Eq. (9) relates only to the regular part, the sources are connected with the residual of the anisotropic part. In this sense, R↔ is a matrix of discrete values of the radiance factor for the reflection, and T↔ is for the transmission. However, the light field at points inside the slab cannot be calculated by using this solution.

From Eq. (9) it is obvious that the complexity of its solution depends on the number of azimuthal harmonics M, ordinates N and harmonics in the representation of the anisotropic part K that, in general, are approximately equal to each other. However, with the successful elimination of the anisotropic part it can be achieved that K >> N >> M, and the angular distribution of the smooth part is quasi-isotropic in this case. In [5] it is shown that the best method of the anisotropic part elimination is the small-angle modification of the spherical harmonics method (MSH) [5, 6] – program MDOM. The MSH solution has a form

L˜a(τ,μ,φ)=∑k=0∞∑m=−kk2k+14πZ˜k(τ)Qkm(μ0)Qkm(μ)eimφ,
where Z˜k(τ)=e−dkτ/μ0−e−τ/μ0,  dk=1−Λxk,  Qkm(μ)≡(k−m)!(k+m)1Pkm(μ) are the half-normalized Schmidt polynomials.

We substitute the expansion of the solution anisotropic part into the expression (4) of the source function. Let’s analyze the item with derivation separately, using the ratio as defined above

μdL˜a(τ,μ,φ)dτ=μμ0∑k=0K∑m=−kk2k+14π(Λxke−τ/μ0−dkZ˜k(τ))Qkm(μ0)Qkm(μ)eimφ.

After transferring in the last expression to the azimuth harmonics, one can get

Δm(τ,μ)=(1−μμ0)∑k=mK2k+14π(Λxke−τ/μ0−dkZ˜k(τ))Qkm(μ0)Qkm(μ).

Then we perform discretization of the obtained Eq. by the sighting angle using Sykes scheme:

Δ→m(τ)=(1↔−M↔/μ0)∑k=mK2k+14π(Λxke−τ/μ0−dkZ˜k(τ))Qkm(μ0)Q→km.

In accordance with Eqs. (7) and (10) the source function in the Eq. (7) in this case is given by

J↔=∑k=mK2k+14πQkm(μ0)((1−Λxk)J↔k−J↔0)U↔−1(1↔−μ0M↔−1)Q→km,
where

J↔k≡diag(1γiμ0−dk)(H↔e−dkτ0/μ0−S↔),  J↔0≡diag(1γiμ0−1)(H↔e−τ0/μ0−S↔).

3. Number of discrete ordinates for representation of solution regular part

However, when comparing different algorithms an interesting feature is revealed [7]. At the accuracy of less than 1% MDOM excels in speed by orders of magnitude, when the mean convergence is used. Nevertheless, MDOM requires approximately the same time as all the other algorithms in case of the uniform metric. In the uniform metric computations require approximately the same number of ordinates N and the azimuthal harmonics M for all the methods of the solution anisotropic part elimination.

A reasonable question is raised here: how many discrete ordinates N are necessary to represent the smooth part. Since a regular part of the solution is a smooth function, its series of Legendre polynomials has a finite number of terms N:

C(μ)=∑k=0N2k+12ckPk(μ).

Legendre polynomials are correctly represented by the Lagrange interpolation formula

Pk(μ)=∑i=1N+1Pk(μi)PN+1(μ)(μ−μi)P′N+1(μi),
where μi is the root of polynomial PN+1(μ).

Accordingly, this leads to the expression

C(μ)=∑i=1N+1PN+1(μ)(μ−μi)P′N+1(μi)∑k=0N2k+12ckPk(μi)=∑i=1N+1C(μi)PN+1(μ)(μ−μi)P′N+1(μi),
that corresponds to Lagrange interpolation formula for the function С(μ).

Equations (13)-(15) are the analog of the Nyquist – Shannon sampling theorem for the angular spectrum of the radiance distribution in Legendre polynomials. Consequently if we would like to reproduce details of the structure of the angular distribution of the solution regular part (for example, in the rainbow, the glory), the sampling interval should be comparable with the size of desired details of the radiance angular distribution that determines uniquely the number of ordinates N. А smaller number of ordinates in the solution will lead to smoothing the fine-structure details of the angular distribution. Therefore, the given level of accuracy determines the number of ordinates N for the definition of the solution regular part. Accordingly, the presence of singularities in the angular distribution in the backward hemisphere eliminates the differences in the methods of the anisotropic part elimination.

In nuclear physics [7] a method of synthetic iterations (SI) is suggested. The accuracy of calculations for smaller N can be increased significantly using ideas of SI. By SI each step of iterations is divided into two steps [7]: at the first step the approximate angular distribution of the RTE solution is sought, and at the second step it is specified by the algorithm of standard iteration. Since MDOM provides the fast means convergence by power, we can expect a substantial acceleration of the RAD calculation when using synthetic iterations. For calculating the iteration of the solution smooth part it is necessary to know its angular distributions at an arbitrary point of the slab, but the ratio of Eq. (9) expresses the emergent radiation through the incident one.

4. Radiance angular distribution inside slab

After the solution of the discrete RTE, taking into account the boundary conditions, the two-points boundary value problem is transformed into two one-point boundary value problems (with initial conditions) that allows using each of them to determine the radiance field inside the slab:

  • through the field on the upper boundary
    U↔−1C→(τ)=e−Γ↔τU↔−1C→(0)+e−Γ↔τ∫0τeΓ↔tU↔−1M↔−1Δ→(t,μ0)dt,
  • through the field on the lower boundary
    U↔−1C→(τ)=eΓ↔(τ0−τ)U↔−1C→(τ0)−e−Γ↔τ∫ττ0eΓ↔tU↔−1M↔−1Δ→(t,μ0)dt.

Note that the same function is on the left side of both Eqs. (16) and (17), and the expressions on the right sides are different ones. In this case, the expression (16) contains only positive exponents for the upper hemisphere of sighting directions, and expression (17) contains the same for the lower one. Since left sides are equal, then right sides are also equal. That allows us to take the expression for the lower hemisphere from Eq. (16), and for the upper hemisphere from Eq. (17):

U↔−1C→(τ)=[eΓ↔−(τ0−τ)u↔11C→−(τ0)+eΓ↔−(τ0−τ)u↔12C→+(τ0)e−Γ↔+τu↔21C→−(0)]+e−Γ↔τ[−(∫ττ0eΓ↔tU↔−1M↔−1Δ→(t,μ0)dt)−(∫0τeΓ↔tU↔−1M↔−1Δ→(t,μ0)dt)+].

The obtained expression contains only negative exponents. Therefore, the radiance field calculation does not present any problems for the arbitrary slab optical depth. Substituting Eq. (10) for the source function and calculating the corresponding integrals in Eq. (18) similarly to Eq. (11), we get the resultant expression for the radiance field inside the slab:

U↔−1C→(τ)=[eΓ↔−(τ0−τ)00e−Γ↔+τ]G→(τ0)+F→(τ),
where

F→(τ)≡∑k=mK2k+14πQkm(μ0)diag(dke−dkτ/μ0γiμ0−dk−e−τ/μ0γiμ0−1)U↔−1(1↔−μ0M↔−1)Q→km≡∑k=mKdiag(dke−dkτ/μ0γμ0−dk−e−τ/μ0γiμ0−1)f→k,
G→(τ0)≡[u↔11C→−(τ0)+u↔12C→+(τ0)u↔21C→−(0)]−[(F→(τ0))−(F→(0))+].

The radiance angular distributions at some points of the turbid medium are presented in Fig. 1 as polar diagrams.

 figure: Fig. 1

Fig. 1 Normalized radiance angular distribution inside the slab.

Download Full Size | PDF

5. Method of synthetic iterations

For the iteration calculation, we take RTE in the Peierls integral form, which for the plane unidirectional source is given by [6]:

L(τ,l^)={  e−τ/μ0δ(l^−l^0)+Λ4πμ∫0τe−(τ−t)/μ∮x(l^,l^′)L(t,l^′)dl^′dt,  μ≥0;−Λ4πμ∫ττ0e−(τ−t)/μ∮x(l^,l^′)L(t,l^′)dl^′dt,  μ<0.

To express the radiance of the first iteration, we substitute the approximate solution, obtained by the Eq. (19), into the right side of Eq. (22) that gives us a formula for the scattered component:

L(1)(τ,l^)={  Λ4πμ∫0τe−(τ−t)/μ∮x(l^,l^′)(La(t,l^′)+Lr(t,l^′))dl^′dt,  μ≥0;−Λ4πμ∫ττ0e−(τ−t)/μ∮x(l^,l^′)(La(t,l^′)+Lr(t,l^′))dl^′dt,  μ<0.

The full field inside the slab takes a form:

L(τ,μ,φ)=La(τ,μ,φ)+Lr(τ,μ,φ)=∑k=0K2k+14πZk(τ)Pk(ν)+∑m=0M(2−δ0m)Cm(τ,μ)cos(mφ),
where Zk(τ)=Z˜k(τ)+e−τ/μ0≡e−dkτ/μ0,  ν=(l^,l^0).

On the basis of the addition theorem for the Legendre polynomial, the scattering integral in Eq. (24) can be written in the form:

∮x(l^,l^′)L(t,l^′)dl^′≡Sa(τ,l^)+Sr(τ,l^)=∑k=0K(2k+1)xkZk(t)Pk(l^⋅l^0)+2π∑m=−MM(2−δ0m)cos(mφ)∑k=0K(2k+1)xkQkm(μ)∫−11Cm(t,μ′)Qkm(μ′)dμ′.

The first iteration of the solution anisotropic part is the integral of MSH along the path:

L+a(τ0)=Λ4πμ∫0τ0e−(τ0−t)/μ∮x(l^,l^′)La(t,l^′)dl^′dt=Λμ04π∑k=0K(2k+1)xkμ0−μdk(e−dkτ0/μ0−e−τ0/μ)Pk(l^⋅l^0),
L−a(0)=−Λ4πμ∫0τ0et/μ∮x(l^,l^′)La(t,l^′)dl^′dt=Λμ04π∑k=0K(2k+1)xkμ0−μdk(1−e(1/μ−dk/μ0)τ0)Pk(l^⋅l^0).

The integral of the solution regular part in Eq. (25) is calculated using Gaussian quadrature that gives the expression in the matrix form

S→m(τ)≡[S→−m(τ)S→+m(τ)]=M↔−1A↔W↔[C→−(τ)C→+(τ)]=M↔−1A↔W↔C→(τ),
where S→r(τ)=∑m=−MM(2−δ0m)cos(mφ)S→m(τ).

The last expression (28) in accordance with Eq. (19) results in

S→m(τ)=X↔{[eΓ↔−(τ0−τ)00e−Γ↔+τ]G→(τ0)+F→(τ)},  X↔≡M↔−1A↔W↔U↔,
from which we have the expression for the integrals of the solution regular part along the path that is the radiance of the first iteration of the smooth part

[L→˜−(0)L→˜+(τ0)]=∑m=0M(2−δ0m)cos(mφ)diag(−1↔, e−τ0/μi)∫0τ0diag(et/μi)S→m(t)dt.

It appears from Eq. (30) that the calculation of the first iteration of the solution regular part is reduced to the following integrals:

∫0τ0diag(et/μi)S→m(t)dt=∫0τ0diag(et/μi)X↔[eΓ↔−(τ0−t)00e−Γ↔+t]dtG→(τ0)+∫0τ0diag(et/μi)X↔F→(t)dt.

Note that after some simple algebraic transformations one can obtain relations:

diag(et/μi)X↔[eΓ↔−(τ0−t)00e−Γ↔+t]=(X↔⊗[e(1/μi−γj)t])[eΓ↔−τ0001↔],
diag(et/μi)X↔F→(t)=∑k=mKX↔⊗[et/μi(dke−dkt/μ0γjμ0−dk−e−t/μ0γjμ0−1)]f→k,
where we use the MATLAB notation: ⊗ denotes the element-by-element multiplication of two matrices, [aij] denotes the matrix with a common element aij.

Evaluating integrals in (31) over dt leads to expressions:

[L→˜−(0)L→˜+(τ0)]=∑m=0M(2−δ0m)cos(mφ)(P→m(1)(τ0)+∑k=mKP→k(2)(τ0)),
where

P→m(1)(τ0)=(X↔⊗[e−τ0/μi+−e−γj+τ01/μi+−γj+e−(1/μi++γj+)τ0−11/μi++γj+−e−(1/μi++γj+)τ0−11/μi++γj+−e−τ0/μi+−e−γj+τ01/μi+−γj+])G→(τ0),
P→k(2)(τ0)=μ0(X↔⊗[−dk(e(μ0/μi−dk)τ0/μ0−1)(γjμ0−dk)(μ0/μi−dk)+e(μ0/μi−1)τ0/μ0−1(γjμ0−1)(μ0/μi−1)dk(e−dkτ0/μ0−e−τ0/μi)(γjμ0−dk)(μ0/μi−dk)−e−τ0/μ0−e−τ0/μi(γjμ0−1)(μ0/μi−1)])f→k,Γ↔±=[γj±].

Since the anisotropic part of the solution includes the direct radiation, the iteration will exactly contain the first orders of scattering, and the other orders of scattering will be greatly improved in accuracy. The advantage of this approach is the recalculation of the angular distributions of the solutions for any arbitrary grid angles. A numerical comparison of the first synthetic iteration with MDOM is shown in Fig. 2, where Δt is the time of calculation.

 figure: Fig. 2

Fig. 2 Comparison of synthetic iteration (SI) with MDOM for the radiation reflected from the slab.

Download Full Size | PDF

From Fig. 2, one can see that the greatest difficulty for computing represents the neighborhood of the glory (Fig. 2(b)). For the calculation by MDOM, it is required N = 401 and M = 256 in this area that corresponds to the description of the glory peak of about 1°. In order to achieve the same accuracy, SI requires only N = 41, M = 4, which reduces the calculation time almost by 100 times. Accordingly SI from MDOM allows calculating the RAD with the convergence in the uniform metric with 1% accuracy and computational speed of less than 1 sec.

The simple form of the RTE solution with the anisotropic part elimination on the basis of MSH [5, 6] allows obtaining an analytical expression for the RAD inside the slab and calculating its first iteration. The first iteration of the solution [5, 6] substantially refines RAD that allows using MDOM with low values of N and M as an initial approximation that speeds up greatly the computation.

6. Quasi two-stream approximation

It allows making a conjecture that the most effective way of the RTE solution is reduced to the synthetic iterations using the two-stream approximation together with the solution anisotropic part elimination by MSH as a method of determining the smooth part. To move to the two-stream approximation, it is necessary to integrate RTE over the upper and the lower hemispheres

{∫Ω−μ∂L˜(τ,l^,l^0)∂τdl^+∫Ω−L˜(τ,l^,l^0)dl^=Λ4π∫Ω−∮x(l^′,l^)L˜(τ,l^′,l^0)dl^′dl^+∫Ω−Δ(τ,l^,l^0)dl^,∫Ω+μ∂L˜(τ,l^,l^0)∂τdl^+∫Ω+L˜(τ,l^,l^0)dl^=Λ4π∫Ω+∮x(l^′,l^)L˜(τ,l^′,l^0)dl^′dl^+∫Ω+Δ(τ,l^,l^0)dl^,
where Ω+ and Ω− represent the upper and lower hemispheres, respectively.

Let’s introduce new variables:

E˜↓↑(τ)≡∫Ω±L˜(τ,μ)dl^,  Δ↓↑(τ)≡∫Ω±Δ(τ,l^,l^0)dl^,  μ↓↑=12π∫Ω±μdl^ ,
βc≡14π∫Ω±x(l^′,l^)dl^|l^′∈Ω±,  βo≡14π∫Ω±x(l^′,l^)dl^|l^′∈Ω∓
allowing rewriting the system (37) in matrix form
M↔dE→(τ)dτ=−A↔E→(τ)+Δ→(τ),
where
E→(τ)≡[E˜↓(τ)E˜↑(τ)],  M↔≡diag(μ↑,μ↓),  A↔=[1−Λβc−Λβo−Λβo1−Λβc]≡1↔−Λ[βcβoβoβc],  Δ→(τ)=[Δ↑(τ)Δ↓(τ)].
The boundary conditions for the system (39) follow from the boundary conditions (1) taking into account their integrations over the upper and lower hemispheres:

E˜↓(τ)|τ=0=0, E˜↑(τ)|τ=τ0=−∫Ω−La(τ0,μ,φ)dl^≡−Ea↑(τ0)

Similarly Eq. (5), the solution of Eq. (39) has the form

[−u11eγ1τ0u12−u21eγ2τ0u22][E˜↑(0)E˜↓(τ0)]=S→(0,τ0)+[eγ1τ0u11eγ2τ0u21]Ea,
where U↔−1=[u11u12u21u22],  B↔≡M↔−1A↔,  eB↔τ=U↔eΓ↔τU↔−1,  S→(τ1,τ2)≡∫τ1τ2eΓ↔tU↔−1M↔−1Δ→(t)dt.

Accordingly, for the field inside the medium we obtain expressions

E→(τ)=U↔(e−Γ↔τU↔−1[E˜↑(0)0]+e−Γ↔τS→(0,τ)),  E→(τ)=U↔(e−Γ↔(τ−τ0)U↔−1[−EaE˜↓(τ0)]−e−Γ↔τS→(τ,τ0)),
combining which similarly Eq. (18) one can obtain a stable expression for the field inside the medium.

The expression for the radiance at any point of the field in the medium can be written as

L(τ,l^)=La(τ,l^0,l^)+12πE˜↓(τ)θ(μ)+12πE˜↑(τ)θ(−μ),  θ(μ)={1, μ≥0,0, μ<0.

Further calculations are carried out in full accordance with synthetic iterations. Because of the linearity RTE the iteration of the total solution is the sum of the iterations of the MSH and the regular part. The expression for the first iteration of MSH was obtained above. Accordingly, the expression for the iteration of the regular part have the form for the radiance of the transmitted radiation μ≥0:

L˜(τ0,μ)=Λe−τ0/μ4πμ[∫0τ0(E˜↓(τ)+E˜↑(τ))et/μdt+l(μ)∫0τ0(E˜↓(τ)−E˜↑(τ))et/μdt],
for the radiance of the reflected radiation μ≤0:
L˜(0,μ)=−Λ4πμ[∫0τ0(E˜↓(τ)+E˜↑(τ))et/μdt+l(μ)∫0τ0(E˜↓(τ)−E˜↑(τ))et/μdt],
where

l(μ)=∑j=1(K+1)/2(4j−1)x2j−1P2j−1(μ)αj,  αj≡∫01P2j−1(μ)dμ,  αj+1=−2j−12(j+1)αj,  α1=0.5.

Calculations based on the obtained expressions (44)-(45) are different from corresponding calculations by the synthetic iterations shown in Fig. 2 no more than 2%.

The two-stream approach has long been widely used in research on the theory of radiative transfer. However, in the case of turbid media with strongly anisotropic scattering this solution allows to describe only the reflection coefficient, for example [8], Meador and Weaver. The elimination of the solution anisotropic part allows the two-stream approach to find a quasi-isotropic angular distribution of the solution regular part. Using iteration to the obtained complete solution permits us to refine the details of the radiance angular distribution, almost without increasing the calculation time. Thus, we are able to describe not the reflection coefficient of a turbid medium slab but its radiance factor.

7. Conclusion

The solution of RTE in MSH is possible in an analytical form for any medium geometry [6]. There are no analytical problems for the derivation of the expression for the solution iteration in an arbitrary 3D geometry. However, the two-stream approach is based on the plane symmetry of a slab. The two-stream approximation is the simplest case of the discrete ordinates method. The solution in a P1 approximation of the spherical harmonics method or another the diffusion approximation is equivalent in accuracy to the two-stream approach [9], Davis and Marshak.

The diffusion approximation and MSH admit simple solutions of RTE for an arbitrary medium geometry. Therefore, the use of the synthetic iterations with RTE solution at the first step in the diffusion approximation eliminating the anisotropy based on MSH allows solving the problems of transport theory in the arbitrary medium geometry.

MSH describes the radiance angular distribution of concentrated sources more accurately. Therefore, we can expect that these estimates of the algorithm accuracy will be valid in the case of three-dimensional geometry too.

Acknowledgment

The Ministry of Education and Science of the Russian Federation (Grant no. 14.604.21.0042) supported this work.

References and links

1. F. Hanson and I. Bendall, “Off-axis laser beam imaging and characterization with two cameras,” Appl. Opt. 52(22), 5342–5347 (2013). [CrossRef]   [PubMed]  

2. J.-P. Cariou, “Off-axis detection of pulsed laser beams: simulation and measurements in the lower atmosphere,” Proc. SPIE 5086, 129–138 (2003).

3. N. Roy and F. Reid, “Off-axis laser detection model in coastal areas,” Opt. Eng. 47(8), 086002 (2008). [CrossRef]  

4. A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C-Labonnote, B. Mayer, Q. Min, T. Nakajima, Y. Ota, A. S. Prikhach, V. V. Rozanov, T. Yokota, and E. P. Zege, “Benchmark results in vector atmospheric radiative transfer,” J. Quantum Spectrosc. Radiat. Transf. 111(12-13), 1931–1946 (2010). [CrossRef]  

5. V. P. Budak, D. S. Efremenko, and O. V. Shagalov, “Efficiency of algorithm for solution of vector radiative transfer equation in turbid medium slab,” J. Phys. Conf. Ser. 369, 012021 (2012). [CrossRef]  

6. V. P. Budak, D. A. Klyuykov, and S. V. Korkin, “Convergence acceleration of radiative transfer equation solution at strongly anisotropic scattering,” in Light Scattering Reviews 5: Single Light Scattering and Radiative Transfer, A.A. Kokhanovsky eds. (Springer Praxis Books, 2010), pp.147–204.

7. M. L. Adams and E. W. Larsen, “Fast iterative methods for discrete-ordinates particle transport calculations,” Prog. Nucl. Energy 40(1), 3–159 (2002). [CrossRef]  

8. W. E. Meador and W. R. Weaver, “Two-stream approximations to radiative transfer in planetary atmospheres: A unified description of existing methods and a new improvement,” J. Atmos. Sci. 37(3), 630–643 (1980). [CrossRef]  

9. A. B. Davis and A. Marshak, “Multiple scattering in clouds, insights from three-dimensional diffusion/P_1 theory,” Nucl. Sci. Eng. 137, 251–288 (2001).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (2)
Fig. 1
Fig. 1 Normalized radiance angular distribution inside the slab.

View in Article | Download Full Size | PDF

Fig. 2
Fig. 2 Comparison of synthetic iteration (SI) with MDOM for the radiation reflected from the slab.

View in Article | Download Full Size | PDF

Equations (52)

Equations on this page are rendered with MathJax. Learn more.

{ μ ∂ L r (τ, l ^ ) ∂τ + L r (τ, l ^ )= Λ 4π ∮ x( l ^ , l ^ ′ ) L r (τ, l ^ ′  ) d l ^ ′ +Δ(τ, l ^ ); L r (τ,   l ^  ) | τ=0, μ>0 =0; L r (τ,  l ^  ) | τ= τ 0 , μ<0 =− L a (τ, l ^ ),
Δ(τ,μ,φ)=−μ d L a (τ,μ,φ) dτ − L a (τ,μ,φ)+ Λ 4π ∮ x( l ^ , l ^ ′ ) L a (τ, μ ′ , φ ′ )d l ^ ′ .
L a (τ,μ,φ)= e −τ/ μ 0 δ( l ^ − l ^ 0 )+ L ˜ a (τ,μ,φ),
Δ(τ,μ,φ)=−μ d L ˜ a (τ,μ,φ) dτ − L ˜ a (τ,μ,φ)+ Λ 4π ∮ x( l ^ , l ^ ′ ) L ˜ a (τ, μ ′ , φ ′ )d l ^ ′ + Λ 4π e −τ/ μ 0 x( l ^ 0 , l ^ ).
L → r (τ, μ i ± ,φ)= ∑ m=0 M (2− δ 0,m )cos(mφ) C → m (τ, μ i ± ) ,   C → ± (τ)≡ C → m (τ, μ i ± ),
d C → (τ) dτ =− B ↔ C → (τ)+ M ↔ −1 Δ → (τ),   B ↔ ≡ M ↔ −1 ( 1 ↔ − A ↔ W ↔ ),
− C → (0)+ P ↔ (0, τ 0 ) C → ( τ 0 )= ∫ 0 τ 0 P ↔ (0,τ) M ↔ −1 Δ → (τ, μ 0 )dτ
− S ↔ U ↔ −1 C → (0)+ H ↔ U ↔ −1 C → ( τ 0 )= S ↔ ∫ 0 τ 0 e Γ ↔ t U ↔ −1 M ↔ −1 Δ → (t, μ 0 )dt ,
e B ↔ τ 0 = U → e − Γ ↔ τ 0 U → −1 ,
[ C → − (0) C → + ( τ 0 ) ]=[ F → − F → + ]+[ R ↔ T ↔ T ↔ R ↔ ][ C → + (0) C → − ( τ 0 ) ]
L ˜ a (τ,μ,φ)= ∑ k=0 ∞ ∑ m=−k k 2k+1 4π Z ˜ k (τ) Q k m ( μ 0 ) Q k m (μ) e imφ ,
μ d L ˜ a (τ,μ,φ) dτ = μ μ 0 ∑ k=0 K ∑ m=−k k 2k+1 4π ( Λ x k e −τ/ μ 0 − d k Z ˜ k (τ) ) Q k m ( μ 0 ) Q k m (μ) e imφ .
Δ m (τ,μ)=( 1− μ μ 0 ) ∑ k=m K 2k+1 4π ( Λ x k e −τ/ μ 0 − d k Z ˜ k (τ) ) Q k m ( μ 0 ) Q k m (μ) .
Δ → m (τ)=( 1 ↔ − M ↔ / μ 0 ) ∑ k=m K 2k+1 4π ( Λ x k e −τ/ μ 0 − d k Z ˜ k (τ) ) Q k m ( μ 0 ) Q → k m .
J ↔ = ∑ k=m K 2k+1 4π Q k m ( μ 0 )( (1−Λ x k ) J ↔ k − J ↔ 0 ) U ↔ −1 ( 1 ↔ − μ 0 M ↔ −1 ) Q → k m ,
J ↔ k ≡diag( 1 γ i μ 0 − d k )( H ↔ e − d k τ 0 / μ 0 − S ↔ ),   J ↔ 0 ≡diag( 1 γ i μ 0 −1 )( H ↔ e − τ 0 / μ 0 − S ↔ ).
C(μ)= ∑ k=0 N 2k+1 2 c k P k (μ) .
P k (μ)= ∑ i=1 N+1 P k ( μ i ) P N+1 (μ) (μ− μ i ) P ′ N+1 ( μ i ) ,
C(μ)= ∑ i=1 N+1 P N+1 (μ) (μ− μ i ) P ′ N+1 ( μ i ) ∑ k=0 N 2k+1 2 c k P k ( μ i ) = ∑ i=1 N+1 C( μ i ) P N+1 (μ) (μ− μ i ) P ′ N+1 ( μ i ) ,
U ↔ − 1 C → ( τ ) = e − Γ ↔ τ U ↔ − 1 C → ( 0 ) + e − Γ ↔ τ ∫ 0 τ e Γ ↔ t U ↔ − 1 M ↔ − 1 Δ → ( t , μ 0 ) d t ,
U ↔ − 1 C → ( τ ) = e Γ ↔ ( τ 0 − τ ) U ↔ − 1 C → ( τ 0 ) − e − Γ ↔ τ ∫ τ τ 0 e Γ ↔ t U ↔ − 1 M ↔ − 1 Δ → ( t , μ 0 ) d t .
U ↔ −1 C → (τ)=[ e Γ ↔ − ( τ 0 −τ) u ↔ 11 C → − ( τ 0 )+ e Γ ↔ − ( τ 0 −τ) u ↔ 12 C → + ( τ 0 ) e − Γ ↔ + τ u ↔ 21 C → − (0) ]+ e − Γ ↔ τ [ − ( ∫ τ τ 0 e Γ ↔ t U ↔ −1 M ↔ −1 Δ → (t, μ 0 )dt ) − ( ∫ 0 τ e Γ ↔ t U ↔ −1 M ↔ −1 Δ → (t, μ 0 )dt ) + ].
U ↔ −1 C → (τ)=[ e Γ ↔ − ( τ 0 −τ) 0 0 e − Γ ↔ + τ ] G → ( τ 0 )+ F → (τ),
F → (τ)≡ ∑ k=m K 2k+1 4π Q k m ( μ 0 )diag( d k e − d k τ/ μ 0 γ i μ 0 − d k − e −τ/ μ 0 γ i μ 0 −1 ) U ↔ −1 ( 1 ↔ − μ 0 M ↔ −1 ) Q → k m ≡ ∑ k=m K diag( d k e − d k τ/ μ 0 γ μ 0 − d k − e −τ/ μ 0 γ i μ 0 −1 ) f → k ,
G → ( τ 0 )≡[ u ↔ 11 C → − ( τ 0 )+ u ↔ 12 C → + ( τ 0 ) u ↔ 21 C → − (0) ]−[ ( F → ( τ 0 ) ) − ( F → (0) ) + ].
L(τ, l ^ )={    e −τ/ μ 0 δ( l ^ − l ^ 0 )+ Λ 4πμ ∫ 0 τ e − (τ−t) /μ ∮ x( l ^ , l ^ ′ )L(t, l ^ ′ )d l ^ ′ dt ,  μ≥0; − Λ 4πμ ∫ τ τ 0 e − (τ−t) /μ ∮ x( l ^ , l ^ ′ )L(t, l ^ ′ )d l ^ ′ dt ,  μ<0.
L (1) (τ, l ^ )={    Λ 4πμ ∫ 0 τ e − (τ−t) /μ ∮ x( l ^ , l ^ ′ )( L a (t, l ^ ′ )+ L r (t, l ^ ′ ) )d l ^ ′ dt ,  μ≥0; − Λ 4πμ ∫ τ τ 0 e − (τ−t) /μ ∮ x( l ^ , l ^ ′ )( L a (t, l ^ ′ )+ L r (t, l ^ ′ ) )d l ^ ′ dt ,  μ<0.
L(τ,μ,φ)= L a (τ,μ,φ)+ L r (τ,μ,φ)= ∑ k=0 K 2k+1 4π Z k (τ) P k (ν) + ∑ m=0 M (2− δ 0m ) C m (τ,μ)cos(mφ) ,
∮ x( l ^ , l ^ ′ )L(t, l ^ ′ )d l ^ ′ ≡ S a (τ, l ^ )+ S r (τ, l ^ )= ∑ k=0 K (2k+1) x k Z k (t) P k ( l ^ ⋅ l ^ 0 ) +2π ∑ m=−M M (2− δ 0m )cos(mφ) ∑ k=0 K (2k+1) x k Q k m (μ) ∫ −1 1 C m (t, μ ′ ) Q k m ( μ ′ )d μ ′ .
L + a ( τ 0 )= Λ 4πμ ∫ 0 τ 0 e − ( τ 0 −t) /μ ∮ x( l ^ , l ^ ′ ) L a (t, l ^ ′ )d l ^ ′ dt = Λ μ 0 4π ∑ k=0 K (2k+1) x k μ 0 −μ d k ( e − d k τ 0 / μ 0 − e − τ 0 /μ ) P k ( l ^ ⋅ l ^ 0 ) ,
L − a (0)=− Λ 4πμ ∫ 0 τ 0 e t/μ ∮ x( l ^ , l ^ ′ ) L a (t, l ^ ′ )d l ^ ′ dt = Λ μ 0 4π ∑ k=0 K (2k+1) x k μ 0 −μ d k ( 1− e (1/μ − d k / μ 0 ) τ 0 ) P k ( l ^ ⋅ l ^ 0 ) .
S → m (τ)≡[ S → − m (τ) S → + m (τ) ]= M ↔ −1 A ↔ W ↔ [ C → − (τ) C → + (τ) ]= M ↔ −1 A ↔ W ↔ C → (τ),
S → m (τ)= X ↔ { [ e Γ ↔ − ( τ 0 −τ) 0 0 e − Γ ↔ + τ ] G → ( τ 0 )+ F → (τ) },   X ↔ ≡ M ↔ −1 A ↔ W ↔ U ↔ ,
[ L → ˜ − (0) L → ˜ + ( τ 0 ) ]= ∑ m=0 M (2− δ 0m )cos(mφ)diag( − 1 ↔ ,  e − τ 0 / μ i ) ∫ 0 τ 0 diag( e t/ μ i ) S → m (t)dt .
∫ 0 τ 0 diag( e t/ μ i ) S → m (t)dt = ∫ 0 τ 0 diag( e t/ μ i ) X ↔ [ e Γ ↔ − ( τ 0 −t) 0 0 e − Γ ↔ + t ]dt G → ( τ 0 )+ ∫ 0 τ 0 diag( e t/ μ i ) X ↔ F → (t)dt .
diag( e t/ μ i ) X ↔ [ e Γ ↔ − ( τ 0 −t) 0 0 e − Γ ↔ + t ]=( X ↔ ⊗[ e (1/ μ i − γ j )t ] )[ e Γ ↔ − τ 0 0 0 1 ↔ ],
diag( e t/ μ i ) X ↔ F → (t)= ∑ k=m K X ↔ ⊗[ e t/ μ i ( d k e − d k t/ μ 0 γ j μ 0 − d k − e −t/ μ 0 γ j μ 0 −1 ) ] f → k ,
[ L → ˜ − (0) L → ˜ + ( τ 0 ) ]= ∑ m=0 M (2− δ 0m )cos(mφ)( P → m (1) ( τ 0 )+ ∑ k=m K P → k (2) ( τ 0 ) ) ,
P → m (1) ( τ 0 )=( X ↔ ⊗[ e − τ 0 / μ i + − e − γ j + τ 0 1/ μ i + − γ j + e −(1/ μ i + + γ j + ) τ 0 −1 1/ μ i + + γ j + − e −(1/ μ i + + γ j + ) τ 0 −1 1/ μ i + + γ j + − e − τ 0 / μ i + − e − γ j + τ 0 1/ μ i + − γ j + ] ) G → ( τ 0 ),
P → k (2) ( τ 0 )= μ 0 ( X ↔ ⊗[ − d k ( e ( μ 0 / μ i − d k ) τ 0 / μ 0 −1 ) ( γ j μ 0 − d k ) ( μ 0 / μ i − d k ) + e ( μ 0 / μ i −1) τ 0 / μ 0 −1 ( γ j μ 0 −1)( μ 0 / μ i −1) d k ( e − d k τ 0 / μ 0 − e − τ 0 / μ i ) ( γ j μ 0 − d k ) ( μ 0 / μ i − d k ) − e − τ 0 / μ 0 − e − τ 0 / μ i ( γ j μ 0 −1)( μ 0 / μ i −1) ] ) f → k , Γ ↔ ± =[ γ j ± ].
{ ∫ Ω − μ ∂ L ˜ (τ, l ^ , l ^ 0 ) ∂τ d l ^ + ∫ Ω − L ˜ (τ, l ^ , l ^ 0 )d l ^ = Λ 4π ∫ Ω − ∮ x( l ^ ′ , l ^ ) L ˜ (τ, l ^ ′ , l ^ 0 )d l ^ ′ d l ^ + ∫ Ω − Δ(τ, l ^ , l ^ 0 )d l ^ , ∫ Ω + μ ∂ L ˜ (τ, l ^ , l ^ 0 ) ∂τ d l ^ + ∫ Ω + L ˜ (τ, l ^ , l ^ 0 )d l ^ = Λ 4π ∫ Ω + ∮ x( l ^ ′ , l ^ ) L ˜ (τ, l ^ ′ , l ^ 0 )d l ^ ′ d l ^ + ∫ Ω + Δ(τ, l ^ , l ^ 0 )d l ^ ,
E ˜ ↓↑ (τ)≡ ∫ Ω ± L ˜ (τ,μ)d l ^ ,   Δ ↓↑ (τ)≡ ∫ Ω ± Δ(τ, l ^ , l ^ 0 )d l ^ ,   μ ↓↑ = 1 2π ∫ Ω ± μd l ^   ,
β c ≡ 1 4π ∫ Ω ± x( l ^ ′ , l ^ )d l ^ | l ^ ′ ∈ Ω ± ,   β o ≡ 1 4π ∫ Ω ± x( l ^ ′ , l ^ )d l ^ | l ^ ′ ∈ Ω ∓
M ↔ d E → (τ) dτ =− A ↔ E → (τ)+ Δ → (τ),
E → (τ)≡[ E ˜ ↓ (τ) E ˜ ↑ (τ) ],   M ↔ ≡diag( μ ↑ , μ ↓ ),   A ↔ =[ 1−Λ β c −Λ β o −Λ β o 1−Λ β c ]≡ 1 ↔ −Λ[ β c β o β o β c ],   Δ → (τ)=[ Δ ↑ (τ) Δ ↓ (τ) ].
E ˜ ↓ (τ) | τ=0 =0,  E ˜ ↑ (τ) | τ= τ 0 =− ∫ Ω − L a ( τ 0 ,μ,φ)d l ^ ≡− E a↑ ( τ 0 )
[ − u 11 e γ 1 τ 0 u 12 − u 21 e γ 2 τ 0 u 22 ][ E ˜ ↑ (0) E ˜ ↓ ( τ 0 ) ]= S → (0, τ 0 )+[ e γ 1 τ 0 u 11 e γ 2 τ 0 u 21 ] E a ,
E → (τ)= U ↔ ( e − Γ ↔ τ U ↔ −1 [ E ˜ ↑ (0) 0 ]+ e − Γ ↔ τ S → (0,τ) ),   E → (τ)= U ↔ ( e − Γ ↔ (τ− τ 0 ) U ↔ −1 [ − E a E ˜ ↓ ( τ 0 ) ]− e − Γ ↔ τ S → (τ, τ 0 ) ),
L(τ, l ^ )= L a (τ, l ^ 0 , l ^ )+ 1 2π E ˜ ↓ (τ)θ(μ)+ 1 2π E ˜ ↑ (τ)θ(−μ),  θ(μ)={ 1, μ≥0, 0, μ<0.
L ˜ ( τ 0 ,μ)= Λ e − τ 0 /μ 4πμ [ ∫ 0 τ 0 ( E ˜ ↓ (τ)+ E ˜ ↑ (τ) ) e t/μ dt +l(μ) ∫ 0 τ 0 ( E ˜ ↓ (τ)− E ˜ ↑ (τ) ) e t/μ dt ],
L ˜ (0,μ)=− Λ 4πμ [ ∫ 0 τ 0 ( E ˜ ↓ (τ)+ E ˜ ↑ (τ) ) e t/μ dt +l(μ) ∫ 0 τ 0 ( E ˜ ↓ (τ)− E ˜ ↑ (τ) ) e t/μ dt ],
l(μ)= ∑ j=1 (K+1)/2 (4j−1) x 2j−1 P 2j−1 (μ) α j ,   α j ≡ ∫ 0 1 P 2j−1 (μ)dμ ,   α j+1 =− 2j−1 2(j+1) α j ,   α 1 =0.5.
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use