An outline of the method of successive order of scattering（MSOS)
Sonoyo Mukai（The Kyoto College of Graduate Studies for Informatics),
Itaru Sano，Makiko Nakata（Kindai University)

Abstract
This short note intends to explain an algorithm for multiple light scattering simulations in the
hazy atmosphere involving polarized radiation eld. It is named MSOS (method of successiveorder of scattering) and available for severe air pollutions from satellite.

1.

Introduction

ground measurements and algorithm improvement are necessary in order to understandIncreased emissions of anthropogenic atmos-

aerosols better. In other words, synthetic analysis

pheric particles (simply called aerosols hereafter)of aerosols from ground [3: Sano et al. 2003], andassociated with economic growth can lead tocomputer [4: Takemura et al. 1999, 5: Mukai et al.

increased concentrations of hazardous air pollu-

2008] as well as space is strongly desired.

tants. Furthermore, yellow dust storms or bio-

Retrieval of optical characteristics of aerosol ismass burnings due to forest re and/or agriculture

achieved by radiation simulation in the Earthburnt eld present serious environmental hazards

atmosphere model. Accordingly various aspects

in East Asia, yet their aerosol properties areshould be taken into account in single and/or

poorly understood. Aerosol retrieval from themultiple light scattering processes concerned withsatellite requires radiation simulations in theaerosols. This short note intends to explain an

Earth atmosphere model. Our research group has

algorithm for multiple light scattering simula-treated with development of efcient algorithms

tions in the hazy atmosphere involving polarized

for aerosol remote sensing including severe hazyradiation eld. It has been already achieved aero-

episodes (dense concentrations of atmospheric

sol retrieval for haze episodes by using satellite

aerosols) [1: Mukai et al. 1992, 2: Mukai et al.

data based on MSOS (method of successive order2015]. It is well known that the space obser-of scattering). Then this efcient algorithm MSOS

vations are particularly available for a systematic

is applied here to vector form radiation simulation

monitoring of the Earth with different spatial,

denoted by Stokes parameters (I, Q, U, V) or (Il, Ir,

spectral, and temporal resolutions. It is of interestU, V) in the optically thick atmosphere. In otherto mention that satellite remote sensing is also

words, this work proposes an algorithm for

useful in a case of heavy air pollution episode.

radiation simulation in the Earth atmosphere

While extreme concentrations of aerosols in theconsisting in the polarized radiation eld. Because

atmosphere can prevent aerosol measurements

the satellite polarimetric sensor POLDER has

with surface-level photometry, satellites can stillshown that the spectro-photopolarimetry is very

be used in such conditions to observe the Earth’s

usefulatmosphere from space. It is highly likely that

Deschamps et al. 1994]. JAXA has been develop-

large scale air pollution will continue to occur. We

ing the new Earth observing satellite GCOM-C,

know with certainty that much more cross-which will board the polarimetric sensor SGLI

validation between the space observations and the(second-generation global imager) in 2017. The

for

the

observationofaerosols[6:

NAIS Journal

45

SGLI has two polarization channels at near-angle, θ, between the outward normal of the plane

infrared wavelengths of 670 and 870 nm. The

parallel medium and the direction of the light, andaerosol models have been developed using theaccumulated measurements during more than tenof the azimuth angle, φ, of the light measured~
from a xed direction. Also, P (Ω, Ω' ) is the matrix

years provided with the world wide aerosolof μ, μ' , φ, and φ' , similarly dened quantities onmonitoring network (AERONET) [7: Dubovik etincident and scattered ray. Hereafter we shall use

al. 2002, 8: Omar et al. 2005].

an abbreviation that Ω stand for (μ, φ). The

The remainder of this article is organized as

standard problem in the functions for a given

follows: Section 2 is the brief interpretation of themedium with constant or vertically-resolved

radiation simulations for remote sensing. Section

single scattering albedo ω. The basic physical3 is the main body of this work on explanation ofconcept is simple. The intensity I (τ, Ω) is reduced

the method of successive order of scatteringby e-dτ after traversing over an optical thickness dτ

(MSOS) for a semi-innite atmosphere.

at a wavelength λ. The equation of radiativetransfer describes this circumstance as [9:2.

Radiation simulations for remote sensing

Algorithms for aerosol retrieval with remote

Chandrasekhar 1960]
μd I(τ, Ω)ωd τ ＝I(τ, Ω)− 4π

EP~(Ω, Ω' )・I(τ, Ω' )d Ω' ,sensing data are divided into three parts in brief(1)

as: 1) satellite data analysis, 2) aerosol modelling,and 3) multiple light scattering calculations in thewhere ω and P~ represent albedo and phase

Earth atmosphere model which are called radia-

matrix for single scattering, respectively. The

tive transfer simulations. The radiative transfer

integration is carried out over all solid angle. The

simulations incorporate the Rayleigh scatteringvector I is dened to be the four dimensionalby molecules and Mie scattering by aerosols in the

representing the polarized radiation eld,

atmosphere, and reection by the Earth surface
I＝(Il, Ir, U, V ) or (I, Q, U, V )(2)

into account. Thus the aerosol properties are
estimated by comparing satellite measurementswith the numerical values of radiation simula-

These four quantities represent the Stokes
parameters, and the 4 × 4 matrix P~ represents

tions in the coupled atmosphere-surface model.

the phase relation in the light scattering. The

It is noted that the precise simulation of multiple

boundary conditions for the standard problem arelight-scattering processes is necessary, and needs(i) at the top of the medium, τ ＝ 0, there is noa long computational time especially in the pola-

radiation falling except in direction − μ0. (ii) atrized radiation eld. Therefore, efcient algo-

the bottom of the medium, τ ＝ τ0, for the nite

rithms for radiative transfer problems are indis-

medium, no radiation incident on the medium atpensable to retrieve aerosols from space.

all (namely free surface), or for the semi-innite

The space-borne sensors measure upwellingcase.

radiance at the top of the atmosphere (TOA), andEq. (1) with the boundary conditions above canit is known that incident solar light interacts

be rewritten in terms of source function J (τ, Ω),multiply with atmospheric particles. Now function I (τ, Ω) is dened to be the specic intensityωJ(τ, Ω)＝ 4π

EP~(Ω, Ω' )・I(τ, Ω' )d Ω' ,(3)vector at the optical depth τ in the direction of Ω,
asgiven by Ω＝(μ, φ) and dΩ＝dμ dφ, where μ is the
cosine of the zenith angle θ ( i.e. μ＝cos θ) and φ isthe azimuth angle.

NAIS Journal(4)

That is to say I (τ, Ω) is thefunction of τ, the optical depth, of μ, the cosine of46

J(τ, Ω)＝ωΛτ{J(τ, Ω)}＋e−τ/μ0,

where Λ denotes the usual Λ-operator for semi-

innite or nite medium [10: Busbridge 1960].

terms of ωn, namely Neumann series solution. The

The physical meaning of this equation is clear.

physical meaning of these equations is obvious.

The rst and second terms in Eq. (4) corresponds

Eq. (5) gives the source function which isto the probability of multiple scattering and directcomposed of photons scattered once or more, andtransmission over the distance τ in direction μ0,

Eq. (6) for photons not experienced scattering.

respectively [11: Sobolev 1963]. The reection and

The entire source function given by Eq. (7) is atransmission functions are obtained throughweighted sum of these source functions by ωn,

Laplace transform of J (τ, Ω).

since a photon survives with probability ω at a

A more direct approach to the problem is tosingle act of scattering.

formulate a set of functional equations for the

It is easy to treat Eq. (7) with only rst two orreection and transmission functions. The prin-

three terms [12: van de Hulst 1949]. However, theciple of invariance stated originally by Ambart-

numerical computation must be consulted beyond

sumian and examined extensively by Chandra-

the 3rd order. The direct treatment on thesekhar is elegant to attack the problem. Thereection and transmission functions in thisintegral equation for the reection function in

respect is also possible. We can derive a set ofcase of semi-innite medium which is related torecurrence relation between the n-th orderH-function is given by Chandrasekhar. And, infunction add the lower ones from some mathe-case of nite medium, the reection and trans-

matical manipulation starting from Eq. (6) or

mission functions are related to X and Y function

from some physical reasoning. The problem in thisand the equations for them are given also by

successive approximation method is the slowness

Chandrasekhar.

of the convergence in Eq. (5). When ω is small
compared to one, or when the thickness of the

3. The method of successive order of scattering

medium is small, the computation is easy.
However, when ω→1 and medium becomes semi-

3-1

The successive scattering approximationinnite, the labour is formidable. This problem

The idea of successive order of scattering iswill be considered later.

fundamentally related to the probabilistic approach to the theory of radiative transfer. Many

3-2

Brief concept of MSOS

papers have been written in this eld and much

It seems not trivial to examine if the successive

progress has been made. However, since it is too

scattering technique mentioned in the previous

heavy to review all of these here, only a narrowsection is useful to the radiation eld reectedeld which concerns directly to the method of

from the optically semi-innite atmosphere model

successive order of scattering (MSOS) is simply[13: van de Hulst 1963, 14: Sovolev 1968, 15, 16:reviewed in this section.

Uesugi&Irvine 1969, 1970]. This technique is

Eq. (4) suggests an iterative solution of a formin n times scattering case [11: Sobolev 1963],

available for the aerosol episodes. Then we treat
with the radiation simulation algorithms in asemi-innite polarized radiation eld.

−τ/μ0

Jn(τ, Ω)＝ωΛτ{Jn−1(τ, Ω)}＋e
J0(τ, Ω)＝e−τ/μ0,

nB1,

(5)

Supposethere is an incident radiation of ux F in direction(μ0, φ0) falling on the top of a semi-innite(6)
atmosphere, and the ux is supposed to bepolarized as

and
∞

Σ

ωnJn(τ, Ω),J(τ, Ω)＝F＝(Fl, Fr, FU, FV)

(8)

(7)

n＝0

which represents the expansion of J (τ, Ω) in

Let the diffusely reected radiation intensity
vector be I (Ω).

Decomposing I (Ω) into the

NAIS Journal

47

successively scattered intensity I (n: Ω) from thesame reasoning of Eq.(7), i.e.,
∞

I(Ω)＝

E P~ (Ω, Ω')・F(Ω')d Ω'
1 Δ ~
Δ
−
R
d Ω'
4π E ( μ − μ' ) (1:Ω, Ω')・F(Ω')μ'
1
＋
R~
d Ω'.
4π E (1 : Ω, Ω')・F(Ω')μ'

Δ
＝ 4πμ

−

−

Σ

ωn I(n : Ω).

(9)

n＝1

−

(15)

In other words, the above diffusely reected

Since the thickness (Δ) of added layer is very

intensity denotes the emergent intensity fromthin, the intensity I (1: Ω) in Eq. (13) must be

TOA (τ＝0) and it is derived by using reection~
matrix R

equal to I (1: Ω, Δ) given by Eq.(15). As a result,

I(Ω)＝

1
4π

E R~(Ω, Ω')・F(Ω')μ' d Ω',

(10)

−

we get the following relation

E ( μ − μ' )R~ (1 : Ω,Ω')・F(Ω')μ' d Ω'1r μ E P~(Ω, Ω')・F(Ω')d Ω'.
1Δ−

where the integration covers inward hemisphere

−

(16)

and the reection matrix is expressed by theNow, letting
Stokes parametersF(Ω')＝δ ( μ−μ0 )・δ ( φ'−φ0 )E,
Rll

Rlr

RlU

(17)

RlV

Rrl Rrr RrU RrVR~＝
RUl RUr RUU RUVRVl RVr RVU RVV

where E is the unity vector (1,1,1,1), we obtain

(

)

(11)
from eq. (16)~~(μ＋μ0 )R (1 : Ω, Ω0 )＝P (Ω, Ω0 ).

(18)

Also we can dene the n-th order reection in
In a similar manner as above, we have theeq.(11) as:
second order one as∞Σω R (n : Ω, Ω ),R~(Ω, Ω0)＝

n ~

(12)

0

n＝1

I(2 : Ω)＝ 1
4π

E R~(2 : Ω, Ω')・F(Ω')μ' d Ω'−

(19)

where n is the number of scattering. The matrix R~
Under the same reasoning as above, we get(n: Ω, Ω0) is the n-th order reection describingtheradiation

emerging

at

the

TOA

after

scattering n times within the atmosphere. That is~
to say reection R (Ω, Ω0) can be calculated as the~
Neumann series of the nth-order of reectance R＝(n: Ω, Ω0).

＋E d4πΩ'' E d4πΩ' (1− μ )Δ−−

R~ (1 : Ω, Ω')・ΔP~ (Ω', Ω'')・F(Ω'')

Now, when a thin layer of thickness Δ is added at

E d4πΩ'' E d4πΩ' μ P~ (Ω, Ω')・
Δ
R~ (1 : Ω', Ω'')・(1− μ'' ) F(Ω'')μ''Δ
Δ＋E d Ω'' (1− μ )R~ (2 : Ω,Ω'')・(1− μ'' ) F(Ω'')μ''4π＝ΔE d Ω'' E d Ω' R~ (1 : Ω, Ω')P~(Ω', Ω'')・F(Ω'')
4π
4π
d
Ω''
d
Ω' P~ (Ω, Ω')・
＋ΔE
4π E 4π

the TOA of this semi-innite medium, the radia-

R~ (1:Ω', Ω'')・F(Ω'')μ''tion eld becomes−Δ

The equation for the n-th order of reection will
be derived here as follows. We have a dening~
equation for R (1: Ω, Ω0) as
1I(1 : Ω)＝ 4πμE P~(Ω, Ω')・F(Ω')d Ω'−

E R~ (1 : Ω, Ω')・F(Ω')μ' d Ω'
Δ
Δ
1
＋ 4π
E (1− μ )R~ (1:Ω, Ω')・(1− μ' ) F(Ω')μ' d Ω'
1I(1 : Ω, Δ)＝ 4π

−

−

(14)
or

NAIS Journal

Δ

−

−

−

(13)

where the integration is over inward half space.

48

I(2 : Ω, Δ)

−

＋

−

＋

E d4πΩ'' ( μ ＋ μ'' )R~ (2 : Ω, Ω'')・F(Ω'')μ''
1
＋ 4π
(20)
E R~(2 : Ω, Ω')・F(Ω')μ' d Ω'.
Δ

Δ

−

−

From Eqs. (19) and (20), and by applying Eq.
(17),

we have~(μ＋μ0 )R (2 : Ω, Ω0 )

μ
4π
μ0
＋ 4π
＝

E R~ (1 : Ω, Ω')・P~ (Ω', Ω )d Ω'
E P~(Ω', Ω )・R~(1 : Ω, Ω')d Ω'.

~
The single scattering matrix P (Ω, Ω0) can be0＋

expanded in a Fourier series with respect to the
(21)

0

＋

azimuth angle as

Generally, the equation for the n-th order R~ (n)P~ (Ω, Ω0 )≡P~ ( μ, φ, μ0, φ0 )(nB3) can be derived in a similar manner asshown above. The n-th order emergent intensity I∞Σ~＝P 0( μ, μ0 )＋2

l＝1(n, Ω) is dened as
1
I(n : Ω)＝ 4π[P~c(l )( μ, μ0 )cosl ( φ−φ0 )E R~(n : Ω, Ω')・F(Ω')μ' d Ω'
−

(22)

Then

(25)

Similarly the reection matrix (R~ ) can be also
expanded in a Fourier seriesI(n : Ω, Δ)＝~＋P s(l )( μ, μ0)sinl ( φ−φ0 )]E d4πΩ' E d4πΩ'' (1− μ )Δ−

R~ (n : Ω, Ω0 )≡R~ (n : μ, φ, μ0, φ0 )

−

∞

Σ

R~ (n−1 : Ω, Ω')ΔP~ (Ω', Ω'')・F(Ω'')~＝R 0(n : μ, μ0 )＋2

E[R~c(l )(n : μ, μ0)cosl ( φ−φ0 )

l＝1

E

d Ω'
d Ω'' Δ P~ (Ω, Ω')・
μ
＋ 4π
− 4π
Δ
R~ (n−1 : Ω', Ω'')・ 1− μ'' F(Ω'')μ''
n−2
d Ω'
d Ω''
d Ω''' 1− Δ ・
＋
μ
4π
4π
4π
−
＋
−
n '＝1
R~ (n' : Ω, Ω')

(

ΣE

E

~
＋R s(l )(n : μ, μ0 )sinl ( φ−φ0 )]

)

E

(

)

(26)

The equations for Fourier coefcients corresponding to equations (18), (21) and (24) are
derived. It is clear that computation time is saved~
×ΔP (Ω', Ω'')・

by using the Fourier expansion. However Eq.(12)Δ~
ΔR (n−n'−1 : Ω'', Ω''') 1− μ''' F(Ω''')μ'''Δ＋ d Ω'' 1− μ R~ (n : Ω, Ω''')・F(Ω''')μ'''.
− 4π

(

E

(

shows the problem of slowness convergence of))MSOS in ω→1. After several improvements have(23)

been made on the present form, the MSOS isavailable for data analysis of the satellite data in

Thus we havea global scale.
~( μ＋μ0 )R (n : Ω, Ω0 )E R~ (n−1 : Ω, Ω')・P~ (Ω', Ω )d Ω'μ＋ 4π E P~ (Ω, Ω')・R~ (n−1 : Ω', Ω )d Ω'μ μ＝ 4π 4π ΣE E R~ (n' : Ω, Ω')・P~ (Ω', Ω'')

μ
＝ 4π

−

0

Acknowledgements

0

S. Mukai really thanks Dr. Akira Uesugi for his

0

＋

n−2

previous works and his valuable suggestions and

0

−

＋

discussions. Several applications of this workn '＝1

×R~ (n−n'−1 : Ω'', Ω0 )d Ω' d Ω''.

(24)

were supported in part by the Global Change
Observation Mission-Climate project by JAXA

Thus the polarized radiation eld reected from(no. JX-PSPC-434796), the Global Environment
the optically semi-innite atmosphere is calcuResearch Fund of the Ministry of Environment,
lated based on Eq. (12) by utilizing Eqs. (18), (21)Japan (S-12), and JSPS KAKENHI Grant Numand (24). The space-borne sensors measure upber 15K00528.
welling radiance at the top of the atmosphere~(TOA), namely the reectance (R )) from the
Earth’ s atmosphere. We name this technique
MSOS (method of successive order of scattering)as a successive scattering method for semi-innite atmosphere.

【References】[1] S. Mukai, I. Sano, K. Masuda, and T. Takashima,
Atmospheric correction for ocean color remote sensing:
Optical properties of aerosols derived from CZCS
imagery, IEEE Transactions on Geoscience and

NAIS Journal49

Remote Sensing, 30, 818-824, 1992.
[2] S. Mukai, M. Nakata, M. Yasumoto. I. Sano I and A.
Kokhanovsky, A study of aerosol pollution episode due
to agriculture biomass burning in the east-central
China using satellite data, Frontiers in Environmental
Science, section Environmental Informatics, 3:57. doi:
10.3389/fenvs.2015.00057, 2015.
[3] I. Sano, S. Mukai, Y. Okada, B.N. Holben, S. Ohta, andT. Takamura, Optical properties of aerosols during
APEX and ACE-Asia experiments, Journal ofGeophysical Research, 108, 8649, doi: 10. 1029/
2002JD003263, 2003.
[4] T. Takemura, T., T. Nozawa, S. Emori, T. Y. Nakajima,
and T. Nakajima, Simulation of climate response toaerosol direct and indirect effects with aerosol
transport-radiation model, Journal of GeophysicalResearch, 110, D02202, doi: 10. 1029/2004JD005029,
2005.
[5] M. Mukai, T. Nakajima, and T. Takemura, A study ofanthropogenic impacts of the radiation budget and thecloud eld in East Asia based on model simulations
with GCM, Journal of Geophysical Research, 113,
D12211, doi:10.1029/ 2007JD009325, 2008.
[6] P. Y. Deschamps, F. M. Bréon, M. Leroy, A. Podaire, A.
Bricaud, J. C. Buriez, and G. Seze, The POLDER
mission: Instrument characteristics and scienticobjectives,” IEEE Transactions on Geoscience andRemote Sensing, 32, 598-615, 1994.
[7] O. Dubovik, B. N. Holben, T. F. Eck, A. Smirnov, Y. J.
Kaufman, M. D. King, D. Tanré, and I. Slutsker,
Variability of absorption and optical properties of key
aerosol types observed in worldwide locations, J.
Atmos. Sci. 59, 590-608, 2002.
[8] A. H. Omar, J,-G. Won, D. M. Winker, S.-C. Yoon, O.
Dubovik, and M. P. McCormick, Development of global
aerosol models using cluster analysis of Aerosol
Robotic Network (AERONET) measurements, Journalof Geophysical Research. 110, doi: 10. 1029/
2004JD004874, 2005.
[9] S. Chandrasekhar, Radiative transfer, Dover publica-

50

NAIS Journal

tions, inc., New York, 1960.
[10] I. W. Busbridge, The mathematics of radiative transfer, Cambridge, 1960.
[11] V. V. Sobolev, A treatise on radiative transfer, (translated by S. I. Gaposchkin), van Nostlsnd, Princeton,
NJ, 1963.
[12] H.C. Van de Hulst, The atmosphere of the Earth andplanets, (ed. G. P. Kuiper), Chicago Univ. Press,
Chicago, 1949.
[13] H. C. Van de Hulst, A new look at multiple scattering,NASA Institute for Space Studies, New York, 1963.
[14] V. V. Sobolev, Anisotropic Scattering of Light in a
Semiinnite Atmosphere, SOVIET Astronmy-A. J.,
12,(3), 1968.
[15] A. Uesugi and William M. Irvine, Computation of
Synthetic Spectra for a Semi-Innite Atmosphere. J.
Atmos. Sci., 26, 973-978, 1969.
[16] A. Uesugi, A. and Irvine, W. M. Multiple scattering ina plane-parallel atmosphere I. Successive scatteringin a semi-innite medium The Astrophysical Journal,
159, 127-135, 1970.

……………………………………………………………………
◆著者紹介

向井

苑生

Sonoyo Mukai

京都情報大学院大学
応用情報技術研究科
環境リモートセンシングセンター(KCGI-REESIT）教授

佐野到

Itaru Sano

近畿大学理工学部情報学科
大学院総合理工学研究科教授

中田真木子

Makiko Nakata

近畿大学総合社会学部准教授