Appl. Phys. B 60, 335-340 4(1995)

Applied Physics B and ,a,rs Optics © Springer-Verlag 1995

A stochastic model and a variance-reduction Monte-Carlo method for the calculation of light transport A. V. Starkov z'*, M. Noormohammadian 2, U. G. Oppel z 1 Computing Center, Siberian Division, Russian Academy of Sciences, Pr. Acad. Lavrentieva 6, 630090 Novosibirsk, Russia z Mathematisches Institut der Ludwig-Maximilians-UniversitS.t Mfinchen, Theresienstrasse 39, D-80333 Mfinchen, Germany (Fax: + 49-89/8153-89217, E-mail: [email protected] Received: 16 March 1994/Accepted: 29 March 1994

Abstract. We compare the transport theoretic and a stochastic approach of modeling light propagation in the atmosphere. Computations of LIDAR return signals using algorithms based on the two different approaches show very good agreement of the numerical data. Furthermore, a deeper analysis of the formulas for the LIDAR return signal obtained from two kinds of modeling shows that they are equivalent. Combining these approaches and introducing splitting techniques, variance reduction Monte Carlo methods and codes are designed which are adapted to the configuration of the LIDAR. These codes allow the calculation of multiple scattering effects with high accuracy. Finally, we point out the importance and some perspectives of ultilization of multiple scattering in laser remote sensing of clouds.

PACS: 42.20; 05.60 Different techniques of modeling the process of light propagation have been proposed and used to obtain formulas and algorithms to calculate the multiply scattered parts of the LIDAR return signal. The modeling may be based on transport theory (integro-differential or integral transport equation; stationary or nonstationary), approximations of the transport equations (e.g., small-angle approximation), the physical Monte Carlo simulation, stochastic models (e.g., Rayleigh random walk, transport processes) and iterated applications of Maxwell equations. Using any of these approaches, the multiple-scattering process can be modeled with and without polarization. All these modeling techniques of multiple scattering are essentially based on the idea of iterations of the single scattering process. Again, single scattering may be modeled using different techniques (e.g., linear optics, Maxwell equations, approximations of Maxwell equations). Single scattering is characterized by the scattering distribution (e.g., defined by the phase function or by the * P r e s e n t a d d r e s s : Groupe de Physique Appliqu6e, Universit6 de Gen6ve, CH-1211 Gen6ve, Switzerland

indicatrix) in the models without polarization and by the scattering matrix (e.g., of the Stokes or the Mueller type) in the models with polarization. In April 1991 a cooperation between the Computing Center (CC) of the Russian Academy of Sciences in Novosibirsk and the Institute of Mathematics of the Munich University (LMU) in the field of the calculation of multiply scattered LIDAR return signals was started. This cooperation was supported by the DFG (German Science Foundation) and the DLR (German Aerospace Research Establishment). The modeling of the CC is based on the integral transport equation which is solved by a Neumann series. The iterates of the integral operator of the Neumann series are complicated iterated integrals. These integrals are calculated using variance reduction Monte Carlo methods; e.g., see Kalos et al. [1], Marchuk et al. [2], Kargin et al. [3], and Starkov [4]. Using these methods, it is possible to calculate the iterated integrals of the Neumann series and hence to calculate the return signals of all scattering orders. We shall briefly review the transport theoretical approach, which is heuristic. Light propagation in the atmosphere is considered as a sequence of collisions of photons with the medium. When a collision takes place, absorption or scattering may occur, i.e., the photon may travel no further or it may leave the point of collision in a new direction. The straight path between the two consequent collisions is called a free path, the angle between the previous and the new direction is called the scattering angle. Let x ~ ~a be the position and ~0e S 2 be the direction of the path of the photon (S 2,= {x ~ ~3: i] x I] = 1}), t ~ ~ the time, c the speed of light in the medium, as (x) the scattering coefficient, o-,(x) the absorption coefficient, o-(x) = o-s(x) + ~r,(x) the extinction coefficient, ~t ,= @, ~0') the cosine of the scattering angle, g (/~, x) the indicatrix, Q (x, (0, t) the source function, and I(x, ~o, t) the intensity of light. By definition, I(x, ~o, t) d 5 ° d~o dt is the average number of photons that cross perpendicular a plane element d5 P in the cone of directions [q~,~0+ d~0] during the time interval [t, t + &].

336

The intensity satisfies the integro-differential transport equation:

The integral equation (3) is solved by the Neumann series = ~ of" gto, where Yf"~o may be interpreted as the

1 OI(x, c

0,

t)

3t

+ ¢

#I(x, ~o, t) 3x

n=0

+ ~(x) I(x, ¢, t)

= a,(x) ~ g(@, ~o' }, x) I(x, p', t) d~o' + Q(x, ~o, t).

contribution of the n-th scattering order. The return signal of the light pulse received by the L I D A R is obtained (1)

Equation (1) is equivalent to the integral transport equation:

~(x, ¢, t) = ~ ~ ~ K(x', ¢', t' ~ x, ~, t) x gt(x', ¢', t') dt' d¢' dx' + ~ ~ Q(x', ~o, t') L ( x ' ~ xl (o) x T(t' ~ t l x, x') dt' dx',

(2)

where ~(x, ¢, t) ,= a(x) I(x, ¢, t) is called the collision density and K(x', (o', t" ~ x, ~o, t) ,= S((p' ~ ~o1x') L ( x ' ~ xl ~o) x T(t' ~ t lx', x) with s(¢

~ ~1 x') ,= ~((~o, ~ ' } , x')/2~,

L ( x ' ~ xi ¢) ,=

a (x) exp [ - v (x' ~ x l ¢)] IX-X'l 2

x5 T(t'--* t l x ' , x ) =

x-x'l

5(t-t'+

{o,

ix-x'

/c),

and r (x' -* x l ~0) the optical distance between x and x'. Again the integral Eq. (2) is equivalent to the integral operator transport equation = ~ 7 " + ~o,

(3)

where a f is defined by the integral kernel K.

Medium of Kind 3

_x(3)

I,x ( l ~ )~1

Mediumof Kind 2

Medium of Kind 1

fr°m the scalar pr°duct < T J ' 0 > = / , = o ~ Y f " g t ° ' 0 >

'

where 0 characterizes the receiver. The iterated integrals (cf"7~o, 0> represent the n times scattered part of the L I D A R return signal. It is difficult, if not impossible, to calculate these iterated integrals analytically or numerically. Variance reduction Monte Carlo methods are used to carry out these calculations. The approach to simulate the transport of light used by the L M U and D L R is based on a stochastic model which describes the process of multiple scattering as a generalized Rayleigh's random walk; it is a general and precise mathematical description of the physical process of multiple scattering; see Oppel et al. [5]. Also this stochastic model allows to derive an exact multiple-scattering L I D A R equation. The total L I D A R return signal may be split up into the parts due to the different orders of scattering. Again, these parts may be expressed by complicated iterated integrals. Using these equations the double-scattering return signal may be calculated fast and easily. We shall briefly review our stochastic model of the process of multiple scattering. The process of multiple scattering of light in the atmosphere is a stochastic process. We assume that light is corpuscular and that photons which are scattered by particles of the atmosphere will travel along straight lines from one point of scattering to the next. Let us consider a photon which is situated in Xo ~ R 3 starting into direction ¢o ~ $2 at the time t= 0 (see Fig. 1). This photon will travel along the halfline H(xo, ~Oo) •.= { X o + t ( O o : t ~ } pointing from Xo in the direction of ~0o. After some time tl it will be scattered at x
~Po

H(xo, CPo;ta, ~ot; ...; t,-1, ¢,,- i) '= {x(,- t)+ ct~o,_l: t ~ R ~ } , Xo

Fig. I. The path of a multiply scattered photon starting at Xo in the direction of ~0o

hence the distribution depends on the "history" (Xo, ~Oo;t 1, (P~; ..-; t,_ 1, ~o,_1)- We call this distribution the

337 "n-th extinction distribution" and denote it by

Halfspace

Slab f I ....

1..0 - 10-t3

i ....

i ....

E,(xo, ~0o;q, (p~; ...; t,-1, f0,-1; "):~(R~) --+ [0, 1],

2.5 - 10-11

where .~(') is the corresponding Borel a-algebra. In many applications we may assume that it is a generalized exponential distribution determined by ~r, i.e., for all t we have

7.5. 10 14

5.0

" 2.0, 10-11 1.5 • 10- n

- I0 -I~

1.0 - i0 -11

En(Xo,

~00; t l , q ) l ; - . . ;

tn-1,

(#n-l;

[0, t])

(4) := ~ ( x ( , _ l ) + c s ~ % _ l ) e x p 0

-

cr(x(,_l)+crg,,_l) dr

ds.

q)o;

t l , (/)1; " - ' ; / n - l ,

~gn-1; l,;-):~(S

2) ~

[0, l].

P,(xo, (Do; tl, Ct; ...; t,-1, ~0,_1; t,; B) is the probability

that a photon with "history" (Xo, fPo; tl, fPl; ...; t,_ 1, fP,- 1; t,) is scattered at x(,) in the direction of ~0, e B, where B is a Borel subset of S z. The stochastic process of LIDAR multiple scattering is determined by the sequence ((E,, P , ) : n e N) of the transition probabilities of the extinction distributions and of the scattering distributions; this sequence determines a probability measure P on the space of the infinite scattering sequences (Xo, ~0o; tl, ~1; ...); for more details see Oppel et al. [5]. Let X~ be the random variable denoting the position of the photon at the time t, let Dt be the random variable denoting the direction of the photon at the time t, and let N, be the random variable denoting the number of scatterings having occured before the time t. The field of view V of the receiver, its inertia, and the length of the laser pulse define a set B ( V ) ~ ~3 x S 2 (which we call the "pulse-pupilla field of view") such that the return signal received at the time t by this receiver is characterized by the probability of the photon to be in the field of view V of the receiver at the time t, i.e., this probability # (t, V) is equal to (see Oppel et al. [5]) Iz(t,

V):=P((X,,Dt)eB(V)) = ~ #,(t, V), where

(5)

n=0

/~.(t, v) ,= p((x,, D,) e B(V); N,=n)

~ s2

5.0 • 10-12

0,0~

0.0 2.0

The distribution of ~, depends on the kind of the scatterer at the n-th scattering point x(,), on the direction ~0,_~, on the polarization, and on the wavelength; hence it depends on the "history" (Xo, (#o; h, ~ol; ...; t,-1, ~o,_1; t,). We denote this "n-th scattering distribution" by

P,(xo,

2.5 - I0 -14

s 2 ~g

i= 1

x E,+l(Xo, 9,o; tl, (01; ...; t,, ~0,; dt,+x) x P,(xo, 9'0; ti, 9,1; ...; t,; dg',) • .. P~(xo, 9'o; tl; dg'~) El(xo, 9'o; dr1).

(6)

2.1

2.2 2.0 Photon Path Length

2.1 [km]

2.2

Fig. 2. Comparision of double-scattering return signals

Equations (5) and (6) together represent an exact multiple-scattering LIDAR equation. The multiple integral /~,,(t, V) describes the n-th scattering order of the return signal, it can be evaluated using numerical, analytical, and Monte Carlo methods or be approximated by properly motivated simplifications of these equations. During the first period of the cooperation the results of the algorithms used by the CC and the L M U / D L R were compared where this was possible. This comparison showed that the two completely different approaches had results which agreed not only qualitatively, but also quantitatively almost perfectly. Figure 2 shows a result of this comparison. Still existing small differences are due to differences in modeling the emitter and the receiver, to different numerical discretization techniques, and to variance produced by the Monte Carlo method. The parameters of the calculations presented in Fig. 2 are the following: monostatic LIDAR geometry, narrow pulsed beam irradiate perpendicular to the C1 cloud; vacuum is outside the cloud; laser pulse length 10 m; distance from LIDAR to the cloud 1 km. Two situations were considered: cloud slab of thickness 100 m with 6 km visibility and 5 mrad receiver FOV; cloud halfspace with 150 m visibility and 10 mrad receiver FOV. During the second period of the cooperation we analyzed the two exact multiple-scattering LIDAR equations which were obtained by the two different approaches. Noormohammadian observed that the two formulas consisting each of a series of iterated integrals could be transformed into each other. Only some few minor agreements concerning the modeling of the emitted beam and the receiver had to be made. Hence, it was possible to show that the transport theoretical approach and the stochastic approach to describe multiple scattering are equivalent. (From the theory of stochastic processes we know that certain stochastic processes are equivalent to differential equations. Also in the case of the multiple scattering process it could be shown that this process is equivalent to a transport equation. But this could be done in full generality only in very simple and specific cases; e.g., see Oppel et al. [5] and Findling [6]. However, in general situations this was not possible.) For a detailed discussion of the equivalence of the stochastic multiple-

338 'l

....

i , , , , i , , , , i

....

i ....

i

. . . .

i , , ' , l ' , , , i , , , , i , ,

1.0 - 10 10 • Sum of all Orders of Scattering O Single Scattering

7.5 , 10 -11

Double Scattering

TrippleScattering 5,0

• 10 -11

2.5 10 11 •

0.0 2.00

2.02

2.04

2.06

Photon Path Length

2.08

[km]

Fig. 3. Multiple-scattering return signals (variance reduction Monte Carlo method)

I ....

I ....

I ....

i ....

[ ....

I ....

i,,''1

....

~....

F''

1.0 10-10 m

/ 7 . 5 . 1 0 -11

/ t

/ 5.0 • i0 - l l

/

~ ~ k

• 30 m VisibiLityinsideCloud [] 120 m Visibility inside Cloud \

× a00 m Visib,ityin,,deCIo0d ~

2.5 " 10 -11

0.0

l l r l [ i r l l l l l l l l l ~ l r l l l l l l l r l l l

2.00

2,02

l l l l l

2.04

2.06

Photon Path Length

....

I ....

~,,

2.08

[km]

Fig. 4. Multiple-scattering return signals (sum of all orders of scattering; variance reduction Monte Carlo method)

um of Kind 3

cp12~~?u

I

¢P13

/

/

.

.

.

.

Medium o{ Kind 2

Mediumof Kind 1

~o Xo

Fig. 5. Cascade process describing the nonphysical process of multiple scattering

scattering process to the transport equation via the comparision of return signals see Kerscher et al. [7]. During the third period of cooperation Starkov developed new variance reduction Monte Carlo methods and a Fortran code for the calculation of multiply scattered LIDAR return signals which allows to calculate the contributions of all orders of scattering. On the basis of this first code Noormohammadian and Starkov developed an improved ANSI-C code. These variance reduction Monte Carlo methods and codes use splitting techniques which are adapted to the specific configuration of the LIDAR (e.g., monostatic or bistatic). Figures 3 and 4 show some results obtained by using these variance reduction Monte Carlo methods. The monostatic LIDAR geometry was used. The narrow pulsed beam irradiates the C1 cloud slab situated in a distance of 1 km from the LIDAR. The receiver FOV is 10 mrad, the laser pulse length is 10 m, and the thickness of the slab is 100 m. Outside of the slab is a Rayleigh atmosphere with the visibility of 10 km. Figure 3 shows the results of the calculation for a value of 30 m visibility inside the slab, and Fig. 4 shows the comparison of the results for the visibilities 30, 120 and 300 m. In our calculations the relative standard deviation was less than 3%. The physical process of light propagation in the atmosphere can be simulated by a Monte Carlo method which is based on our stochastic model. However, using such a physical Monte Carlo simulation for the calculation of LIDAR return signals is very ineffective since only a few photons will arrive at the comparably tiny receiver. At least about 70 million paths of photons would have to be calculated by the physical Monte Carlo simulation to obtain a signal of the quality of the ones shown in the Figs. 3 and 4. Variations of the physical Monte Carlo techniques often lead to a bias. For the construction of unbiased and variance reduced Monte Carlo simulations weight methods are used. They are based on employing "nonphysical" probability distributions. These distributions are constructed properly from the a priori information of the photons importance. Simultaneously weight multipliers are introduced correcting the bias of the estimate. Each trajectory is marked with its weight that represents the number of simulated photons. The variance reduction computationa] schemes are based on the simulation of the transport of a large number of photons with relatively low weights in interesting regions and of only a sufficiently small number of photons with relatively high weights in unimportant regions. The splitting techniques used in the present variance reduction Monte Carlo methods are "point splitting" and "scattering splitting". The point splitting technique is a history dependent combination of "exponential transformation" and splitting in collision points. The scattering splitting technique increases the number of photons scattered towards the receiver. Hence, the applied variance reduction method means a transition from the Monte Carlo simulation of the "physical" multiple scattering process to the Monte Carlo simulation of an "artificial" process which is a pretty complicated cascade process (see Fig. 5). This transition results in consider-

339 able reduction of variance and speeding up of the Monte Carlo simulation. For details see Kerscher et al. [7]. On the basis of the elaborated simulation methods the analysis of the perspectives of the utilization of LIDAR in space for cloud diagnostics was carried out. As compared with ground based and airborn LIDAR systems, which are used for the determination of microphysical parameters of the clouds, the utilization of spaceborn LIDAR for the cloud research is more complicated. This is connected with large distances to the objects under investigation and great moving speed of the space craft. Appearing new questions such as laser pulse energy, optimal optical and geometrical parameters of the experiment, and evaluation of multiple scatterig contributions in the received signal have to be investigated. From this point of view the Monte Carlo computer experiments can give the preliminary useful information about promising sounding schemes.

1.0. t0 -la

7.5 • 1 0 - 1 4

X Total Signal for 0.5 mrad Semiap. Rec. FOV A Total Signal for 0.3 mrad Semiap. Rec. FOV [] Total Signal for 0.l mrad Serniap. Rec. FOV

5.0 • 10-14 t

Q Single Scattering Part of the Return Signal for 0.1, 0.3 and 0.5 mrad Semiap. Rec. FOV, respectively

LI k

2.5 • 10 *4 I

0.0

400.0

400.1

400.2 400.3 400.4 Photon Path Length [km]

400.5

400.6

Fig. 6. Multiple-scattering return signals (variance reduction Monte Carlo method

• 0.5 mrad Semiaperture Rec. FOV ~'-~-

X 0.1 mrad Semiaperture Rec. FOV

~'~

..

100

10

l ~ r

. . . . .

400.0

FI~,,,

400.1

r~

400.2

~1,,,

I

400.3

Photon Path Length

. . . . .

400.4

,'

400.5

r

Since the diameter of the cone of the receiver FOV on the top height of the cloud is compared with the thickness of the cloud, the main part of the multiply scattered light did not go out from the cone, making a lot of contributions to the return signal. Thus, for the LIDAR systems situated at long distances from dense clouds, the level of multiple scattering will be extremely high, so that the physical methods of measurements of single scattering contribution can be complicated or even impossible. The following applications of the space LIDAR for cloud research can be promising: ® The remote sensing of the multilayered structure of the cloud systems. The identification of the cloud layers is carried out on the basis of the analysis and the interpretation of the form of the backscattered signal; ® The polarization characteristics of the return signal allows to determine the phase structure of the cloud. The basis for the utilization of the polarization space LIDAR techniques is the essential differences of the signals and the depolarization degree of the LIDAR returns, for instance for water drop and ice crystal clouds.

1000

[] 0.3 mrad Semiaperture Rec. FOV

Usually, the single scattering return signal is considered as most informative. In principle, from the simple form of LIDAR equation it is possible to obtain the information about the optical properties of the medium. The contribution of multiple scattering can essentially limit this procedure. The problem of the determination of microphysical parameters of the aerosols, optical and geometrical thickness of clouds, and extinction requires high accuracy of measurements and also the correct formula of the received signal. Some of the simulation results are presented in the Figs. 6 and 7. The backscattering LIDAR in space sounds the C1 cloud situated at a distance of 200 km from the LIDAR. The laser wavelength is 1.064 I~m, the semiaperture of the beam is 0.05 mrad, and the pulse length 12 m. The diameter of the receiver telescope is 1 m, the semiaperture of the receiver FOV is supposed varied and equal to 0.1, 0.3, and 0.5 mrad, respectively. The thickness of the cloud is 300 m, the extinction coefficient inside the cloud equals 17.25 km -1. Figure 6 presents the results of the calculations of the total return signal, taking into account all orders of scattering, and the single scattering part of the return signal. Figure 7 shows the ratio of the total return signal to the single scattering part. The maximal value of the ratio here reaches 1000 and even more. The relative standard deviation for the calculations was less than 5 %.

F

400.6

[km]

Fig. 7. R a t i o o f the total r e t u r n signal to the single-scattering p a r t o f the r e t u r n signal ( v a r i a n c e r e d u c t i o n M o n t e C a r l o m e t h o d )

The questions of the utilization of the polarization space LIDAR for the determination of the phase structure of the cloud systems, the discovery of the phase and density inhomogeneities inside of the mixture-type clouds are the promising and require special intensive theoretical and experimental considerations.

340

References 1. M.H. Kalos, P.A. Whitlock: Monte Carlo Methods, Basics, Vol. 1 (Wiley, New York 1986) 2. G.I. Marchuk, G.A. Mikhailov, M.A. Nazaraliev, R.A. Darbinjan, B.A. Kargin, B.S. Elepov (eds.): The Monte Carlo Methods in Atmospheric Optics, Springer Ser. Opt. Sci., Vol. 12 (Springer, Berlin, Heidelberg 1980) 3. B.A. Kargin, K.B. Rakimgulov: Russ. J. Numer. Anal. Math. Modelling 7, 221 (1992)

4. A.V. Starkov: Soviet. Math. Dokl. 43, 92 (1991) 5. U.G. Oppel, A. Findling, W. Krichbaumer, S. Krieglmeier, M. Noormohammadian: A Stochastic Model for the Calculation of Multiply Scattered LIDAR Returns (DLR-FB 89-36, K61n

1989) 6. A. Findling: An Inverse Problem in Radiative Transfer Theory. Dissertation, Universitfit Mfinchen (1993) 7. M. Kerscher, W. Krichbaumer, M. Noormohammadian, U.G. Oppel, A.V. Starkov: Simulation and Inversion of Multiply Scattered LIDAR Signals (DLR, K61n 1993)