- Email: [email protected]
Solution of radiative transfer in anisotropic plane-parallel atmosphere
Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
www.elsevier.com/locate/jqsrt
Solution of radiative transfer in anisotropic plane-parallel atmosphere Tasuku Tanakaa;∗ , Menghua Wangb a
Earth Observation Research Center, National Space Development Agency of Japan, 1-8-10 Harumi, Chuou-ku 1046023, Tokyo, Japan b NASA Goddard Space Center, University of Maryland Baltimore County, Code 970.2, Greenbelt, MD 20771, USA Received 18 September 2002; accepted 27 February 2003
Abstract We solve the radiative transfer problem analytically in the anisotropic, plane-parallel atmosphere. Chandrasekhar formalized the radiative transfer process as a simultaneous, two-variable, non-linear, integral equation and obtained the analytical solution of the second approximation for the isotropic plane-parallel atmosphere. We obtain the second approximation for the anisotropic atmosphere, by integrating the 5rst approximation multiplied by weighting functions which are products of the scattering phase functions. We truncate the second approximation and obtain the radiance at the top of the atmosphere as a quadratic equation with a logarithmic term in the optical thickness. We evaluate the second approximation of the radiance at the top of the atmosphere for Rayleigh scattering and the maritime aerosol atmosphere and compare them with both the exact solution and the single scattering approximation. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Solution of Chandrasekhar integral equation; Radiative transfer process; Anisotropic atmosphere; Iterative integration; Satellite remotesensing; Inversion problem
1. Introduction In the area of satellite remote sensing, we observe the radiance at the top of the atmosphere (TOA) by instruments aboard satellites. And we retrieve the surface re;ectance and the optical thickness from the observed radiance (inversion problem). Since there are two unknown variables for one observation, there must be another observation or a relationship among them. In this paper we ignore the surface re;ectance and seek the solution only for the optical thickness. The retrieval of the surface re;ectance will be discussed in other occasions. ∗
Corresponding author. Tel.: +81-3-6221-9000; fax: +81-3-6221-9192. E-mail addresses: [email protected] (T. Tanaka), [email protected] (M. Wang).
0022-4073/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-4073(03)00105-5
556
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Zenith θ3
i0
i3
Equator
i1
i1
Top of Atmosphere
τ
i2
i3
Azimuth Axis
Earth Surface
i4
i5
ϕ
i’2
3
i’3
i2 i0
Fig. 1. Geometry of radiative transfer.
We show the geometry of the radiative transfer process in Fig. 1. The layer has the vertical optical thickness . The direction of the input solar radiance ˜i0 is the lower bound, the direction of observation by satellites ˜i1 is the upper bound, the direction of the transmitted radiation to the Earth surface ˜i4 is the lower bound, the intermediate direction ˜i2 is the lower bound, and the other intermediate direction ˜i3 is the upper bound. We employ the polar coordinate system in which the zenith angle is measured from the zenith and the azimuth angle ’ is measured from the projection of the solar direction onto the equator. We denote the cosine of the zenith angle and the azimuth angle of the direction ˜in as n and ’n , respectively. We usually employ the generalized re;ectance instead of radiance I , I ; (1) (;˜i1 ;˜i0 ) = F0 |0 | where F0 is the solar irradiance. Since is considered suDciently small to 1, we have the 5rst approximation P(i1 ; i0 ) ; (2) = 4|0 |1 where P(˜i1 ;˜i0 ) is the scattering phase function. The 5rst approximation or the single scattering approximation has been applied to many analyses. Since the advent of Earth observing satellites, in particular the retrieval of the ocean chlorophyll concentration from coastal zone color scanner (CZCS), we need a more accurate formula. Many varieties of calculating methods have been proposed [1,2]. Most of them comprise the forward calculation and the look-up table. They do not give us the explicit form of by . In this paper we seek the quadratic form of in as an explicit form of by , = c1 + c2 2 :
(3)
We employ the iterative integration of Chandrasekhar’s integral equation to solve the inversion problem of the radiative transfer. The radiative transfer process in the atmosphere is governed by
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
557
the integral equation, which was 5rst introduced by Chandrasekhar [3]. He also obtained the explicit solution of radiative transfer in the second approximation for the isotropic plane-parallel atmosphere. We have obtained the explicit solution of the third approximation in the isotropic plane-parallel atmosphere [4]. In this paper we extend the second approximation developed by Chandrasekhar to the anisotropic atmosphere. In Chandrasekhar’s integral equation, there are two unknown functions: the scattering function, S(;˜i1 ;˜i0 ) and the transmitted function, T (;˜i4 ;˜i0 ), for which we seek the solutions (Section 2). We solve this simultaneous integral equation by successive iteration. The 5rst approximation corresponds to the single scattering approximation. Substituting the 5rst approximation into the integration in the original equations, we obtain the second iteration. Adding the 5rst approximation to the second iteration, we obtain the second approximation for S2 (;˜i1 ;˜i0 ) and T2 (;˜i4 ;˜i0 ). We can continue this iterating process and obtain series Sn (;˜i1 ;˜i0 ) and Tn (;˜i4 ;˜i0 ). If the series Sn (;˜i1 ;˜i0 ) and Tn (;˜i4 ;˜i0 ) converge as n approaches to ∞, it satis5es Chandrasekhar’s integral equation (Section 3.1). We decompose the second iteration into a power series expansion in (−) and obtain the quadratic form of the scattering function (Section 3.2). In the process to evaluate S2 (;˜i1 ;˜i0 ), we often encounter a sort of integration, U0 (;˜i0 ) =
0
1
3 |0 | d3 (1 − e−=3 −=|0 | ) : 3 + |0 | 3
(4)
Expanding the term in the integral above into a power series in (−) and integrating inde5nitely those, we obtain 1 (−)2 1 1 d3 −(−) − + − ··· 3 |0 | 2! 3 0 1 log 3 (−)2 1 = −log 3 (−) − − + − ··· : − 3 |0 | 2! 3
U0 (;˜i0 ) =
(5)
Substituting 3 = 1 at the upper integral limit, we can obtain the upper integral value. And the coeDcients of nth power in are separately integrated. However at the lower limit, each integrated term becomes ∞, or has singularities at 3 = 0. For the lower limit, we must sum up all the terms on theinterval (; 1) and then make approach to 0. Brie;y we cannot change the order of “lim” and “ ” at the lower integral limit (Section 3.3). The problem of the radiative transfer is, thus, deduced to evaluate the value of the series expansion at the singularity point. The authors evaluated values of several series expansions at the singularity point, which are necessary to obtain the second and third approximation for the isotropic atmosphere [4]. We show that, for the anisotropic atmosphere, we can exchange the order of “lim” and “ ” for the higher powers of . If we can exchange the order, we can integrate each power (−) separately (Section 3.3). Based on this consideration, we show the coeDcients of the quadratic form of the second iteration expressed as the surface integrations of products of the scattering phase functions on the half unit sphere (Sections 3.4 and 3.5).
558
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Finally, in Section 4, we evaluate the second approximations numerically for the Rayleigh scattering and the maritime aerosol atmosphere. And the results are compared with the single scattering approximation and the exact solutions.
2. Radiation transfer process The scattered radiance at the TOA, I (0;˜i1 ), is expressed with the incident intensity I (0; i0 ) and the scattering function S(; i1 ; i0 ) as 1 S(;˜i1 ;˜i0 )I (0;˜i0 ) d0 ; (6) I (0;˜i1 ) = 4 1 L where 0 is the solid angle subtended around the input direction ˜i0 and integral domain is the lower half of the unit sphere. In the same manner the intensity transmitted to the Earth surface is expressed as 1 I (0;˜i4 ) = I (0;˜i4 ): T (;˜i4 ;˜i0 )I (0;˜i0 ) d0 + exp − (7) 4 |4 | L |4 | Chandrasekhar’s integral equation is expressed [3], 1 1 S(;˜i1 ;˜i0 ) + 1 |0 | (−) (−) P(˜i1 ;˜i0 ) + = 1 − exp 1 |0 | d3 d3 P(˜i1 ;˜i3 )S(;˜i3 ;˜i0 ) − exp − T (;˜i1 ;˜i3 )P(˜i3 ;˜i0 ) + 4 | | 4 3 0 3 U U d2 d2 − exp − + S(;˜i1 ;˜i2 )P(˜i2 ;˜i0 ) P(˜i1 ;˜i2 )T (;˜i2 ;˜i0 ) 4 |2 | 1 4 |2 | L L d2 d3 + S(;˜i1 ;˜i2 )P(˜i2 ;˜i3 )S(;˜i3 ;˜i0 ) 4 |2 | 4 3 U L d2 d3 − T (;˜i1 ;˜i3 )P(˜i3 ;˜i2 )T (;˜i2 ;˜i0 ) ; 4 |2 | 4 3 U L
1 1 − T (;˜i4 ;˜i0 ) |4 | |0 | − exp − P(˜i4 ;˜i0 ) = exp − |0 | |4 |
(8)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2 d3 ˜ ˜ ˜ ˜ − exp − + P(i4 ; i2 )T (; i2 ; i0 ) P(˜i4 ;˜i3 )S(;˜i3 ;˜i0 ) 4 | | | | 4 2 4 3 L U d3 d2 + exp − S(;˜i4 ;˜i3 )P(˜i3 ;˜i0 ) − T (;˜i4 ;˜i2 )P(˜i2 ;˜i0 ) |0 | 4 3 4 |2 | U L d2 d3 + S(;˜i4 ;˜i3 )P(˜i3 ;˜i2 )T (;˜i2 ;˜i0 ) 4 |2 | 4 3 U L d2 d3 − T (;˜i4 ;˜i2 )P(˜i2 ;˜i3 )S(;˜i3 ;˜i0 ) : 4 |2 | 4 3 U U
559
(9)
The cosine of the zenith angle in the lower bound directions, which has a “even index”, is negative. For the convenience of the following evaluation, we replace n by a new parameter, n− which is equal to −n for n = even. Using the new variables, the integration on the lower half sphere is changed to the integration on the upper half unite sphere, shown below −1 2 f(2 ; ··) d2 f(2 ; ·; ·)(−d2 ) d’2 = 4 |2 | 4 (−2 ) 0 0 L 2 1 d2− f(−2− ; ·; ·) d2− − = d’2 = f(− ; ·; ·) : (10) 2 4 2− 4 2− 0 0 U Since we adopt the new variable 2− and the integral domain is changed to the upper half sphere, we omit the U for the integral domain unless the domain has a special meaning. The cosine of the angle between the upper and lower bound directions is thus expressed as below cos(e; o ) = e o + (1 − e2 )1=2 (1 − o2 )1=2 cos(’e − ’o ) = −e− o + (1 − (e− )2 )1=2 (1 − o2 )1=2 cos(’e − ’− o ):
(11)
The phase function is assumed as a function of the angle mn between the directions m and n and is normalized in the whole solid angle 4 . We express the phase function P(im ; in ) as below, P(im ; in ) = P(cos mn ) =
∞
wj Pj (cos mn );
(12)
j=0
where Pj (x) is the Legendre function of the 5rst kind. The input from the upper surface is the solar irradiance, F0 , and is given with the Dirac delta function I (0;˜i) = F0 (˜i − ˜i0 ):
(13)
Substituting the solar irradiance, we obtain the generalized re;ectance from the upper layer S(;˜i1 ;˜i0 ) : (0;˜i1 ) = 41 |0 |
(14)
560
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
3. Second approximation 3.1. Azimuth integration and decomposition of second iteration The simultaneous integral equations of radiative transfer can be solved by successive iteration. The 5rst iterations are the 5rst terms on the right-hand side in Eqs. (8) and (9), expressed below S1 (;˜i1 ;˜i0 ) =
T1 (;˜i4 ;˜i0 ) =
1 1 + − 1 1
−1 (−) (−) 1 − exp P(˜i1 ;˜i0 ); + − 1 0
1 1 − + − 4 0
−1
exp − − 0
− exp − − 4
P(˜i4 ;˜i0 ):
(15)
(16)
A new function U1 (; 0 ; 1 ) is de5ned as −
U1 (; 1 ; 0 ) =
1 1 + − 1 0
−1 (−) (−) 1 − exp : + − 1 0
(17)
Using U1 (; 0− ); 1 ); T1 (;˜i4 ;˜i0 ) is rewritten as
T1 (;˜i4 ;˜i0 ) = exp − − 0
; U1 (; 4− ; (−0− ))P(˜i4 ;˜i0 ):
(18)
Substituting the 5rst iterations into the integrals in the original equations, the second iterations, PS2 (;˜i1 ;˜i0 ) is expressed PS2 (;˜i1 ;˜i0 ) =
1 1 −1 + − 1 1 d3 d3 × P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (0− ) − P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 3 4 3 d2− d2 + P(˜i1 ;˜i2 )P(˜i2 ;˜i0 )U1 (1 ) − P(˜i1 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ) 4 2− 4 2− −1 − exp − − − 1 0 − d3 − d2 ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ + P(i1 ; i2 )P(i2 ; i0 )U1 (−0 ) × P(i1 ; i3 )P(i3 ; i0 )U1 (−1 ) 4 3 4 2− d2− d3 + P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2− P(˜i1 ;˜i3 )P(˜i3 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ; 2− )U1 (−1 ; 3 )) 4 2− −1 − exp − − − 1 0 d2− × P(˜i1 ;˜i3 )P(˜i3 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ; 2− )U1 (−1 ; 3 )) 4 2− −
561
d3 4 3
d3 ; 4 3
(19)
where the arguments in the functions U (;˜i1 ;˜i0 ) are omitted except for the signi5cant ones. Using the addition theorem of the Legendre function, we decompose the phase function into the Bi-Legendre functions, Plm (3 ), as
∞ ∞
(2 − 0; m ) wlm Plm (1 )Plm (3 ) cos m(’3 − ’1 ) ; (20) P(˜i1 ;˜i3 ) = m=0
P(˜i3 ;˜i0 ) =
∞
l=m
(2 − 0; m )
m=0
∞
wkm Pkm (3 )Pkm (−0− ) cos m(’3 − ’0 ) ;
(21)
k=m
where 0; m and wlm are given as below, 0; m = 1 wlm = wl
(m = 0);
0; m = 0
(m = 0);
(l − m)! : (l + m)!
(22) (23)
Due to the orthogonality of the functions Plm (3 ) cos m(’3 − ’1 ) for diQerent m, we obtain the integration 2 1 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) d’3 2 0 ∞ ∞ k m m m m − m m = (2 − 0; m ) (−1) wl wk Pl (1 )Pk (0 )Pl (3 )Pk (3 ) m=0
l; k=m
×cos m(’1 − ’0 ):
(24)
In the same manner, we obtain the integration of the triple product of phase functions 2 2 d’2 d’3 P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) 2 2 0 0
2 2 ∞ ∞ − m m m (2 − 0; m ) wl Pl (1 )Pl (−2 ) cos m(’2 − ’1 ) = 0
0
m=0
l=m
562
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
× ×
∞
(2 − 0; m )
m=0
∞
∞ l=m
(2 − 0; m )
∞
m=0
=
∞
wlm Plm (−2− )Plm (3 )
cos m(’3 − ’2 )
wlm Plm (3 )Plm (−0− ) cos m(’0 − ’3 )
l=m
∞
(2 − 0; m ) cos m(’3 − ’1 )
m=0
d’2 2 d’3 2
(−1)a+b+c
a;b;c=m
{wam wbm wcm Pam (1 )Pam (2 )Pbm (2 )Pbm (3 )Pcm (3 )Pcm (0 )}
:
(25)
We de5ne an integration, Ql;mk (; 0− ), as Ql;mk (; 0− )
=
1
Plm (3 )Pkm (3 )U1 (; 1 ; 0− )
0
d3 : 3
Using Ql;mk (; 0− ), we obtain PS2 ,
1 PS2 (;˜i1 ;˜i0 ) = 2
×
1 1 + − 1 0 ∞
−1 ∞
(2 − 0; m ) cos m(’1 − ’0 )
m=0
wlm wkm {Plm (1 )Pkm (0− )(−1)k }{Ql;mk (; 0− ) − (−1)l+k Ql;mk (; −0− )}
l; k=m
∞
+
l; k=m
wlm wkm {Plm (1 )Pkm (0− )(−1)l }{Ql;mk (; 1 ) − (−1)l+k Ql;mk (; −1 )}
− exp − − − 1 0 ×
∞
−1
wlm wkm {Plm (1 )Pkm (0− )(−1)l }{Ql;mk (−0− ) + (−1)l+k Ql;mk (−1 )}
l; k=m
+
∞ 1 (−1)a+b+c wam wbm wcm Pam (1 )Pcm (0− )Qa;m b (; 1 )Qb;m c (; 0− ) 2 a;b;c=m
−
∞ 1 (−1)b wam wbm wcm Pam (1 )Pcm (0− )Qa;m b (; −1 )Qb;m c (; −0− ) 2 a;b;c=m
(26)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
− exp − − − 0 1
563
−1
∞ 1 − b m m m m m − m m × (−1) wa wb wc Pa (1 )Pc (0 )Qa; b (; −1 )Qb; c (; −0 ) : 2
(27)
a;b;c=m
In the above equation, the sign of each term is determined by the choice of l; k; a; b; c. This determination originates from the even/odd characteristic of the Bi-Legendre functions. 3.2. Zenith integration of second iteration The second iteration PS2 (; i1 ; i0 ) is expressed in Eq. (27) and its component, Ql;mk (; ), has a power series expansion in . In order to obtain the quadratic form of PS2 (; i1 ; i0 ), we decompose Ql;mk (; 0− ) into the power series expansion in (−). First, we integrate U20 (; 0− ) and this is evaluated in the author’s previous work [4] given below, − 1 1 0 d − − d 0 U2 (; i0 ) = 1 − exp − − − U1 (; ; 0 ) = 0 0− + 0 0 ∞ ∞ ∞ (−=0 )r (−=0 )r (−)n − ; (28) = (−) (C + eax()) (r + 1)r! r=0 n=0 (r + 1 + n)n! (r + 1)r! r=0 where the constant C and the function eax(p) are given as C = % + log = 0:5772 · · · + log ;
eax(p) =
∞ (−1)n pn n=1
nn!
;
where % is the Euler constant. The function U20 (; 0 ) is truncated as 1 1 (−)2 − 0 + ···: U2 (; 0 ) = (C − 1)(−) + 1 + C − 2 0− 2! We integrate the moment integration, U2n (; 0− ), as below, 1 d3 − n U2 (; 0 ) = 3n U1 (; 3 ; 0− ) : 3 0
(29)
(30)
(31)
(32)
U2n (; 0− ) is evaluated by partial integration and recurrence relationships (refer to Appendix B). The coeDcients of the series expansion of U2n (; 0− ) in the power (−) are also polynomials of 1=0− .
564
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
The truncated form of U2n (; 0− ) for n ¿ 1 is expressed below. U2n (; 0− ) = (−1)n (−) + u2;n 0
(−)2 (−)2 1 + u2;n 1 + ···; 2! 2! 0−
(33)
where ul;n k is the coeDcient of (−)n =n!(0− )k and is given in Appendix B. It is noted that the coeDcient of the 5rst power (−) does not include 0− . Substituting U2n into Eq. (26), we obtain Ql;mk (; ) as a power series expansion in (−), Ql;mk (; 0 ) =
1 mlk 0
=
− n amlk n 3 U1 (; 3 ; 0 )
n=0 mlk
n amlk n u1; 0
n=0
+
mlk
(−) +
n amlk n u2; 1
n=0
= q1;mlk0 (−) + q2;mlk0
d3 ; 3
mlk
n amlk n u2; 0
n=0
(−)2 2!
(−)2 1 + ···; 2! 0−
(−)2 (−)2 1 + ···; + q2;mlk1 2! 2! 0−
(34)
is the coeDcient of the n degree power of the product of the Bi-Legendre functions where amlk n m mlk Pl (3 ) and Pkm (3 ). The coeDcient qd; e is given as below, mlk qd; e =
ml k
n amlk n ud; e :
(35)
n=0
It is noted that the coeDcient of the 5rst power (−) in Ql;mk (; 0 ) does not include 0− . Substituting Eq. (34) into the 5rst term in the bracket of Eq. (27), we obtain following: Ql;mk (; 0− ) − (−1)l+k Ql;mk (; −0− ) = q2;mlk1
(−)2 2 + ··· 2! 0−
=2q1;mlk0 (−) + 2q2;mlk0
(l + k = even)
(−)2 + ··· 2!
(l + k = odd):
And for the case (l + k = odd), adding the second term to the above, we obtain following, Plm (1 )Pkm (0 )(−1)k {Ql;mk (0 ) − (−1)l+k Ql;mk (−0 )} +Plm (1 )Pkm (0 )(−1)l {Ql;mk (1 ) − (−1)l+k Ql;mk (−1 )}
(36)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
=0
(−)2 + ···: 2!
565
(37)
From the 5rst and second terms in Eq. (27), the 5rst power in (−) vanishes and the second power emerges only for the case (l + k = even). From the third terms in the bracket of Eq. (27), we only need the 5rst power in (−), because it is multiplied by another 5rst power in (−). Substituting Eq. (34) into this term, we obtain Ql;mk (−0 ) + (−1)l+k Ql;mk (−1 ) = 2q1;mlk0 (−) + · · · ; = 0(−) + · · ·
(l + k = even);
(l + k = odd):
(38)
From the fourth and 5fth term, substituting Eq. (34) we obtain following: ∞ 1 (−1)a+b+c wam wbm wcm Pam (1 )Pcm (0 )Qa;m b (1 )Qb;m c (0 ) 2 a;b;c=m
−
∞ 1 (−1)b wam wbm wcm Pam (1 )Pcm (0 )Qa;m b (−1 )Qb;m c (−0 ) 2 a;b;c=m
=0(−)2 + · · · (a + c = even) ∞ b+1 m m m m m mab mbc = (−1) wa wb wc Pa (1 )Pc (0 )q1; 0 q1; 0 (−)2 + · · ·
(a + c = odd):
(39)
a;b;c=m
From the sixth term in the bracket of Eq. (27), no (−)2 term emerges. Thus we obtain the truncated second iterations, PS2 (;˜i1 ;˜i0 ) =
∞
(2 − 0; m )
m=0
(−1)k wlm wkm (q2;mlk1 − 2q1;mlk0 )Plm (1 )Pkm (0 )
(l; k)1
+2
∞
1 1 + − 1 0
−1 ∞
×cos m(’1 − ’0 )
mbc m m (−1)b+1 wam wbm wcm q1;mab 0 q1; 0 Pa (1 )Pc (0 )
(a;b;c)2
(−)2 + ···; 2!
(40)
where (l; k)1 denotes l; k ¿ m and l + k is even and (a; b; c)2 denotes a; b; c ¿ m; a + c is odd. It is noted that the truncated second iteration does not have a 5rst power in (−). It is noted that, the fourth and 5fth terms vanish for wa = 0 and wa = 0 for all possible combinations of a + c = odd such as Rayleigh scattering.
566
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
3.3. Separated integration In this subsection we consider changing order of lim and begin with the integration of U2n (; i7 ), 1 7 n d n 1 − exp − − U2 (; i7 ) = 7 7 + 0
in the integration of the PS2 . We
1 ∞ 1 m−1 n d (−)m 1 + : →0 m! 7 m=1
= −lim
On the integral interval (; 1), we interchange the order of into two parts: one is m 6 n and the other is m ¿ n, m−1 n 1 (−)m n n− m +1 d U2 (; i7 ) = − 7 m! 0 m=1 n
−(−) lim
→0
∞ m=1
1
(−)m (m + n)!
1 1 + 7
(41) 1
m−1
and
+1 7
n
and divide the summation
d :
(42)
The terms of the 5rst part of the summation (m 6 n) are ordinary integrations because the minimum powers in the in the integrand are equal toor greater than 0 and have no singularity at = 0. We can interchange the order of “lim” and “ ” for these terms. If the order is possible to exchange, we can “pick-up” each power in (−) separately. In other words, we can integrate the powers in (−) separately. We call this part as the separated integration part (SIP). The terms of the second part (m ¿ n) involve negative powers in or log , and have singularities at =0. So we cannot interchange the order of “lim” and “ ” and therefore we cannot integrate the terms separately. We call the second part as the non-separated integration part (NIP). It is needless to say that the NIP converges as approaches to 0, because it is equal to the converged U2n (; i6 ) minus the bounded SIP. We apply the above observation to Eq. (27). 3.3.1. The >rst and second terms Due to the discussion in the previous subsection, we can select only the case (l + k = even) and pick up the second power in (−)2 only. Plm (3 )Pkm (3 ) is a polynomial of 3 and its lower powers for l + k = even are given below (refer to Appendix B), Plm (3 )Pkm (3 ) = a0 + a2 32 + · · · ; = a2 32 + a4 34 + · · ·
(m + l = even); (m + l = odd):
(43)
The zeroth power of (−); a0 , is a NIP and the other higher powers of (−) are SIT. Taking account of a0 = Plm (0)Pkm (0), the integration above is then divided into SIP and NIP
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
as below, Ql;mk (; 0− ) = lim Plm (0)Pkm (0) →0 +
1
0
1
U1 (; 3 ; 0− )
567
d3 3
(Plm (3 )Pkm (3 ) − Plm (0)Pkm (0))U1 (; 3 ; 0− )
d3 : 3
(44)
Substituting the second power of (−) in the polynomial U20 (p0 ) into NIP and truncating SIP into the second power of (−), we obtain 1 1 − m Plm (0)Pkm (0) Ql; k (; 0 ) = 1 + C − 2 0− 1 1 1 d3 (−)2 m m m m − (Pl (3 )Pk (3 ) − Pl (0)Pk (0)) + + ···: (45) 3 0− 3 2! 0 For the case (m + l = odd) Eq. (45) holds because Plm (0)Pkm (0) = 0. 3.3.2. The third term We select for the case (l + k = even) and pick up the 5rst power of (−) only. We divide the integration into NIP and SIP, in the same as in Eq. (44). Substituting the 5rst power of (−) in the polynomial U20 (p0 ) into NIP and truncating SIP into the 5rst power of (−), we obtain Ql;mk (; −0− ) = (C − 1)Plm (0)Pkm (0) −
0
1
(Plm (3 )Pkm (3 )
−
Plm (0)Pkm (0))
d3 (−) + · · · + · · · : 3
(46)
For the case (m + l = odd) Eq. (46) holds because Plm (0)Pkm (0) = 0. 3.3.3. The fourth and >fth terms We integrate following for the case (a + c = odd), 1 1 d− d3 Pam (2− )Pbm (2− )Pbm (3 )Pcm (3 )U1 (1 ; 2− ) −2 U1 (3 ; 0− ) : 3 2 0 0
(47)
The lowest powers of 2− and 3 for the product Pam (2− ) · · · are given below (refer to Appendix B), Pam (2− )Pbm (2− )Pbm (3 )Pcm (3 ) =c12 (2− )(3 )2
(m + b = odd; a + b = even);
=c21 (2− )2 (3 )
(m + b = odd; a + b = odd);
=c01 (3 )
(m + b = even; a + b = even);
=c01 (2− )
(m + b = even; a + b = odd):
(48)
568
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
For (m + b = odd) the integration is the product of two SIP: one is for 2− and the other is for 3 . For (m + b = even) the integration is a product of one SIP and NIP. Then the integration is expressed for the four cases in the equation above: 1 1 d2− d3 m − m − m m = (−)2 Pa (2 )Pb (2 )Pb (3 )Pc (3 ) − 3 0 0 2 + · · · (m + b = odd); 1 d3 = −(C − 1)Pam (0)Pbm (0) Pbm (3 )Pcm (3 ) 3 0 1 1 d2− d3 m − m − m m m m (−)2 + {Pa (2 )Pb (2 ) − Pa (0)Pb (0)}Pb (3 )Pc (3 ) − 3 0 0 2 + · · · (m + b = even; a + b = even); 1 d− Pam (2− )Pbm (2− ) −2 = −(C − 1)Pbm (0)Pcm (0) 2 0 1 1 d− d3 (−)2 + Pam (2− )Pbm (2− ){Pbm (3 )Pcm (3 ) − Pbm (0)Pcm (0)} −2 2 3 0 0 + ···
(m + b = even; a + b = odd):
(49)
3.4. Surface integration We can bring back the idea of division of integration to the surface integration on the unit sphere expressed in Eq. (19). This division is possible, because the integrations in Eq. (19) are expressed as the summation of the integrations of the products of Bi-Legendre functions shown in Eq. (27). From the ninth integration in Eq. (19) no second power in (−)2 emerges. 3.4.1. The >rst four integrations The 5rst integration is divided into SIP and NIP expressed as below, d3 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (; 0− ; 3 ) 4 3 d3 = {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )}U1 (; 0− ; 3 ) 4 3 d3 + lim {P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )}U1 (; 0− ; 3 ) : 3 →0 4 3
(50)
where ˜i 3 denotes the projected direction of ˜i3 onto the equator in the unit sphere. Substituting the second power of (−) in the polynomial U20 (p0 ) and the truncated term of the second power into
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Eq. (50), we can obtain d3 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (; 0− ; 3 ) 4 3 d3 (−)2 1 1 =− {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} + 2! 0− 3 4 3 2 2 1 1 d’ 3 (−) ˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) + 1+ C − P( : 2 0− 4 2! 0
569
(51)
Substituting the above equations into the 5rst and second integrations in the bracket in Eq. (19), we obtain d3 − d3 ˜ ˜ ˜ ˜ P(i1 ; i3 )P(i3 ; i0 )U1 (0 ) − P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 3 4 3 (−)2 1 1 d 3 ˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} {P( =− + 4 3 2! 0− 1 2 2 1 1 1 ˜ d’3 (−) ˜ ˜ ˜ : (52) + C− + − P(i1 ; i 3 )P(i 3 ; i0 ) 2 1 0 4 2! 0 We evaluate the third and fourth integrations in Eq. (19) in the same manner. 3.4.2. The >fth and sixth integrations The 5fth integration is divided into SIP and NIP as below, d3 − exp − − − −1 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 0 1 3 2 1 1 d’3 =−2 (C − 1) + − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) 1 0 4 0 (−)2 d 3 + ···: − {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} 4 3 2!
(53)
The equation above holds for the two cases described in (46). We can evaluate the sixth integration in Eq. (19) in the same manner. 3.4.3. The seventh and eighth integrations The seventh integration is divided into SIP and NIP as below, d2− d3 P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3
570
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
=
[P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )
d2− d3 −P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )]U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3 + P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) + =
d2− d3 P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )U1 (1 ; 2− )U1 ((3 ; 0− )] 4 2− 4 3
{P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )
− P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )}
d’2 d3 P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) 4 4 3
− (C − 1) − (C − 1)
d2− d3 4 2− 4 3
d2− d’3 ˜ ˜ ˜ ˜ ˜ ˜ (−)2 : P(i1 ; i2 )P(i2 ; i 3 )P(i 3 ; i0 ) 4 2− 4
(54)
The equation above holds for all the cases described in (49). We can evaluate the eighth integration in Eq. (19) in the same manner. 3.5. The result of the second approximation Expanding the 5rst approximated scattering function into a power series expansion in and truncating up to the second degree, we obtain 2 1 1 ˜ ˜ P(˜i1 ;˜i0 ): + (55) S1 (; i1 ; i0 ) = − 1 |0 | 2! The second iteration is expressed by the integration of the products of the phase functions as below, PS2 (;˜i1 ;˜i0 ) = (3 − 2C)Ie (˜i1 ;˜i0 ) + Iu (˜i1 ;˜i0 ) + Il (˜i1 ;˜i0 ) +2
1 1 + − 1 0
−1
(−)2 : {(C − 1)(Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) + Iul2 (˜i0 ;˜i1 )} 2!
The coeDcients are evaluated on the half unit sphere (we return the variable 2− to 2 ), d3 ˜ ˜ Iu (i1 ; i0 ) = {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} ; 4 3 U
(56)
(57)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2 ; {P(˜i1 ;˜i2 )P(˜i2 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i0 )} 4 |2 | L 2 d’3 ˜ ˜ Ie (i1 ; i0 ) = ; P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) 4 0 2 d’2 d3 Iuu (˜i1 ;˜i0 ) = {P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i1 )} ; 16 2 3 U 0 2 d’3 d2 ˜ ˜ ; Ill (i1 ; i0 ) = {P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i1 )} 2 | | 16 2 L 0 ˜ ˜ Iul (i1 ; i0 ) = [{P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i1 )}; Il (˜i1 ;˜i0 ) =
571
U
(58) (59) (60) (61)
L
− {P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i1 )}; d2 d3 −{P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i1 )}] : 16 2 |2 |3
(62)
Adding the 5rst approximation to the second iteration, we obtain the second approximation, 1 1 (−)2 ˜ ˜ ˜ ˜ P(˜i1 ;˜i0 ) + S2 (; i1 ; i0 ) = P(i1 ; i0 ) − 1 |0 | 2!
−1 1 1 (−)2 + (Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) log + −2Ie (˜i1 ;˜i0 ) + 2 1 |0 | 2! + Iu (˜i1 ;˜i0 ) + Il (˜i1 ;˜i0 ) + (3 − 2%)Ie (˜i1 ;˜i0 ) +2
1 1 + 1 |0 |
−1
{(% − 1)(Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) + Iul (˜i1 ;˜i0 )}
(−)2 : 2!
(63)
From the above equation, we summarize the characteristics of the second approximation. (1) The second approximation is not a quadratic equation of the optical thickness but a quadratic equation with a 2 log(). (2) The second iteration PS2 (;˜i1 ;˜i0 ) does not have the 5rst power of . (3) There exist three kinds of terms which are characterized by the number of multiplication of the scattering phase function: single scattering, double scattering and triple scattering. The single scattering term comes from the 5rst approximation. The surface Ie ; Il ; Iu are the double scattering terms and Iuu ; Ill ; Iul are the triple scattering terms. (4) For the atmosphere with a scattering phase function which does not have the odd order Legendre functions in its series expansion in the cosine of the angle, such as the Rayleigh scattering, there are no triple scattering terms in the second approximation.
572
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 103 Rayleigh Scattering Aerosol: M80 at 443 nm
Phase Function
102
101
100
10-1
10-2 0
20
40
60
80
100
120
140
160
180
Scattering Angle (Deg.) Fig. 2. Scattering phase functions.
4. Calculation of scattering To evaluate the accuracy of the second approximation discussed in the previous sections, we carry out calculations for the Rayleigh and aerosol scattering atmospheres for various cases and compare with the exact solutions. To appreciate the improvement of accuracy from single scattering, we calculate the single scattering re;ectance as (+)!(+)P(i1 ; i0 ) (+;˜i1 ;˜i0 ) = ; 41 |0 |
(64)
where + is the wavelength and !(+) is the single scattering albedo. Fig. 2 illustrates the phase function of air molecules (Rayleigh scattering) as well as for a typical maritime aerosol with relative humidity (RH) of 80.0% (denoted as M80) at wavelength 443 nm. The maritime aerosol model has a single scattering albedo of 0.9929, and the Rayleigh scattering model has a single scattering albedo of 1.0. Note that, in contrast to Rayleigh scattering function which is symmetry in the forward and backward scatterings, the maritime aerosol predominantly scatters in the forward direction. 4.1. The exact solutions The successive-order-of-scattering (SOS) [5] method is used in solving the radiative transfer equation (RTE) for the Rayleigh scattering atmosphere and the aerosol scattering atmosphere. The SOS code for the ocean–atmosphere system was developed for the atmospheric correction algorithm for the ocean color sensors [6,7]. This code computes the upward radiance at the top and the downward radiance at the base of the medium for the ocean–atmosphere system. The code is capable of yielding radiances that are accurate to nearly 0.1% [6]. In the following two subsections, the re;ectance computed by the single scattering and the second approximation are compared with the exact solutions obtained using the SOS code.
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 80 60
20 Rayleigh Atmosphere, θ = 45o 412 nm (S) 443 nm (S)
15
412 nm (D) 443 nm (D)
20 0 -20 -40
490 nm (S) 510 nm (S)
490 nm (D) 510 nm (D)
5 0 -5 -10
-60
-15 ∆φ = 0o -60
(a)
-40
∆φ = 180o -20
0
20
40
∆φ = 180o
∆φ = 0o 60
-20 -80
80
-60
-40
(b)
Solar Zenith Angle (Deg.)
-20
0
20
40
60
80
60
80
Solar Zenith Angle (Deg.) 5
10 Rayleigh Atmosphere, θ = 45 555 nm (S) 670 nm (S)
Rayleigh Atmosphere, θ = 45o
o
765 nm (S) 865 nm (S)
555 nm (D) 670 nm (D)
Error (%)
5
Error (%)
Rayleigh Atmosphere, θ = 45o
10
Error (%)
Error (%)
40
-80 -80
573
0
765 nm (D) 865 nm (D)
0
-5
(c)
-60
-40
-20
0
20
40
Solar Zenith Angle (Deg.)
∆φ = 180o
∆φ = 0o
∆φ = 180o
∆φ = 0o -10 -80
60
-5 -80
80
(d)
-60
-40
-20
0
20
40
Solar Zenith Angle (Deg.)
Fig. 3. Calculation of generalized re;ectances for Rayleigh scattering.
4.2. The Rayleigh scattering atmosphere Fig. 3 provides a quantitative evaluation of accuracy using the second approximation for a Rayleigh scattering atmosphere for the various solar-sensor geometries and for the eight wavelengths at 412, 443, 490, 510, 670, 765, and 865 nm. The Rayleigh optical thicknesses corresponding to the eight wavelengths from 412 to 865 nm are 0.3185, 0.2361, 0.1560, 0.1324, 0.0938, 0.0436, 0.0255 and 0.0155, respectively. In all computations, the Rayleigh atmosphere is treated as a plane-parallel medium without a surface boundary. Figs. 3(a) – (d) show the error (%) in the computed upward re;ectance as a function of the solar zenith angle and for a sensor viewing angle of 45◦ . The left and right parts of each plot in Fig. 3 corresponds to the relative azimuth angles of 0◦ and 180◦ (principal scattering plane), respectively. Results in Fig. 3 show the signi5cant improvement in accuracy using the second approximation compared with the single scattering formula. Using the new formula, the errors are all within 1% for the red and near-infrared wavelengths, while errors are within nearly 10% for the blue bands for
574
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 10
5
5 0
Aerosol, M80 model at 443 nm τa = 0.10, θ = 45o
-5
Error (%)
Error (%)
0
Aerosol, M80 model at 443 nm -10
τa = 0.05, θ = 45o Single-Scattering Approx. Double-Scattering Approx.
-15
-5 -10 -15 Single-Scattering Approx. Double-Scattering Approx.
-20 -25
∆φ = 0o -20 -80
-60
(a)
-40
∆φ = 180o -20
0
20
40
∆φ = 0o 60
-30 -80
80
Aerosol, M80 model at 443 nm
-20
0
20
40
60
80
Solar Zenith Angle (Deg.)
20
Aerosol, M80 model at 443 nm τa = 0.30, θ = 45o
τa = 0.20, θ = 45o
10
Error (%)
0
Error (%)
-40
30
20
10
-60
(b)
Solar Zenith Angle (Deg.)
∆φ = 180o
-10
0 -10 -20
-20
(c)
Single-Scattering Approx. Double-Scattering Approx.
-40 ∆φ = 0o
-40 -80
-30
Single-Scattering Approx. Double-Scattering Approx.
-30
-60
-40
∆φ = 0o
∆φ = 180o -20
0
20
40
Solar Zenith Angle (Deg.)
60
-50 -80
80
(d)
-60
-40
∆φ = 180o -20
0
20
40
60
80
Solar Zenith Angle (Deg.)
Fig. 4. Calculation of generalized re;ectances for maritime aerosol scattering.
zenith angles less than 60◦ . It is important to note that, with the second approximation, the error in a given wavelength (optical thickness) is almost independent of the solar and sensor viewing geometry. 4.3. The maritime aerosol atmosphere A realistic maritime aerosol model as in Fig. 2 is used to study the eDcacy of the double scattering approximation in computing the upward re;ectance. The double scattering approximation is the second approximation without triple scattering. Computations for various coeDcients are much more involved due predominantly to the aerosol forward scattering characteristics. Fig. 4 shows examples of error (%) for sensor viewing angle of 45◦ as a function of the solar zenith angle for the principal scattering plane. Fig. 4(a) – (d) are the results for the aerosol optical thickness of 0.05, 0.1, 0.2, and 0.3, respectively. Again, we include results from the single scattering approximation given by
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
575
Eq. (64) for comparison. Similar to the Rayleigh atmosphere case, the aerosol atmosphere is considered as one plane-parallel layer without the surface boundary. The results in Fig. 4 show signi5cant improvement in accuracy using the double scattering approximation compared with the single scattering formula, in particular, for the part of the results with P- = 0◦ . It is interesting to note that the scattering angles for the P- = 0◦ part (the left part of the plot) vary from 55◦ to 135◦ , (all forward scattering), while the scattering angles for P- = 180◦ (the right part of the plot) change from 135◦ to 180◦ (all backward scattering). We are usually more interested in the results with the large scattering angles because measurements with the large scattering angles are generally accessible for satellite remote sensing. With double scattering approximation and scattering angles ¿ 100◦ , the re;ectance error is usually within about 1% for the aerosol optical thickness of 0.05, while the error is within about 4% for the optical thickness of 0.1. As expected, for the turbid atmosphere, the error increases. The error is proportional to the slant path of the optical thickness, i.e., . =cos . The aerosol re;ectance errors are within about 7% and about 12% for the optical thickness of 0.2 and 0.3, respectively.
5. Conclusion We obtain a new approximated solution in the radiative transfer problem in the anisotropic plane-parallel atmosphere by iterating integration of Chandrasekhar’s integral equation. By the second approximation, the scattering function S(;˜i1 ;˜i0 ) is expressed as a quadratic equation in the optical thickness, , with terms 2 log(). It is noted that the linear term in in the second approximation is identical to the 5rst approximation. This implies that the iterative integration method improves the accuracy from the 5rst approximation and approaches to the solution of radiative transfer. The second approximation is not a pure quadratic equation but a quadratic equation with terms 2 log(). The “log” term comes from the singularity at the lower integral limit. It is possible that, by expanding log around = 0:0, we can obtain the quadratic equation form of the scattering function S(;˜i1 ;˜i0 ). Rather we leave it as “log” to insist the validity of the expression in the range ¿ 0:0. The coeDcients of 2 and 2 log() are evaluated by surface integrations of products of the phase functions on the half unit sphere. The generalized upward re;ectances numerically evaluated by the second approximation for the Rayleigh atmosphere and the maritime aerosol atmosphere are very much closer to the exact solution than the single scattering approximation. In some cases, such as the small optical thickness, they are suDciently close to the exact solution.
Acknowledgements This work is done during Tanaka’s stay at NASA Goddard Space Flight Center as a visiting scientist. The authors expresses their sincere appreciation to the SIMBIOS project people for their assistances.
576
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Appendix A. Derivation of function U2n (; ) We de5ne a new function Wn () as 1 n d3 : 3 exp − Wn () = 3 0
(A.1)
W0 () is evaluated 1 1 1 d3 = 3 exp − exp − − 3 exp − d3 W0 () = 3 3 0 3 0 0 = exp(−) − (−)(C + eax()):
(A.2)
The recurrence relation of Wn () is given as Wn () =
(−) exp(−) + Wn−1 (): n+1 n+1
(A.3)
We can evaluate Wn () by the recurrence relation and the initial function W0 (). Using Wn (), we can derive the recurrence relation for U2n (; ), 1 n− 1 n − exp − Wn−1 − U2 (0 ) : (A.4) U2 (; 0 ) = 0 n 0 Truncating the series expansion up to the second degree in (−), we obtain U21 (; 0 ) = −(−) − U2n (; 0 ) = −
3(−)2 (−)2 1 (−)2 ; +C − 4 2 2 0
(−)2 1 (−)2 (−) − − n 2(n − 1) 2n 0
(n ¿ 2):
Appendix B. Lower powers of the product of the Bi-Legendre function m = e;
l = e;
k = e;
Plm (3 )Pkm (3 ) = a0 + a2 32 ;
m = e;
l = o;
k = o;
Plm (3 )Pkm (3 ) = a2 32 + a4 34 ;
m = e;
l = e;
k = o;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = e;
l = o;
k = e;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = o;
l = e;
k = e;
Plm (3 )Pkm (3 ) = a2 32 + a4 34 ;
m = o;
l = o;
k = o;
Plm (3 )Pkm (3 ) = a0 + a2 32 ;
m = o;
l = e;
k = o;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = o;
l = o;
k = e;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 :
(A.5) (A.6)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
577
References [1] [2] [3] [4] [5] [6]
Hansen JE, Travis LD. Light scattering in planetary atmosphere. Space Sci Rev 1974;16:527–610. Goody RM, Yung YL. Atmospheric radiation. New York, Oxford: Oxford University Press, 1989. Chandrasekhar S. Radiative transfer. Mineola, NY: Dover Publishing, Inc., 1960. Tanaka T. Integration of Chandrasekhar’s integral equation. JQSRT 2003;76:121–44. van de Hulst HC. Multiple light scattering. New York: Academic Press, 1980. Gordon HR, Wang M. Surface roughness considerations for atmospheric correction of ocean color sensors. I: the Rayleigh scattering component. Appl Opt 1992;31:4247–60. [7] Gordon HR, Wang M. Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS preliminary algorithm. Appl Opt 1994;33:443–52.
www.elsevier.com/locate/jqsrt
Solution of radiative transfer in anisotropic plane-parallel atmosphere Tasuku Tanakaa;∗ , Menghua Wangb a
Earth Observation Research Center, National Space Development Agency of Japan, 1-8-10 Harumi, Chuou-ku 1046023, Tokyo, Japan b NASA Goddard Space Center, University of Maryland Baltimore County, Code 970.2, Greenbelt, MD 20771, USA Received 18 September 2002; accepted 27 February 2003
Abstract We solve the radiative transfer problem analytically in the anisotropic, plane-parallel atmosphere. Chandrasekhar formalized the radiative transfer process as a simultaneous, two-variable, non-linear, integral equation and obtained the analytical solution of the second approximation for the isotropic plane-parallel atmosphere. We obtain the second approximation for the anisotropic atmosphere, by integrating the 5rst approximation multiplied by weighting functions which are products of the scattering phase functions. We truncate the second approximation and obtain the radiance at the top of the atmosphere as a quadratic equation with a logarithmic term in the optical thickness. We evaluate the second approximation of the radiance at the top of the atmosphere for Rayleigh scattering and the maritime aerosol atmosphere and compare them with both the exact solution and the single scattering approximation. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Solution of Chandrasekhar integral equation; Radiative transfer process; Anisotropic atmosphere; Iterative integration; Satellite remotesensing; Inversion problem
1. Introduction In the area of satellite remote sensing, we observe the radiance at the top of the atmosphere (TOA) by instruments aboard satellites. And we retrieve the surface re;ectance and the optical thickness from the observed radiance (inversion problem). Since there are two unknown variables for one observation, there must be another observation or a relationship among them. In this paper we ignore the surface re;ectance and seek the solution only for the optical thickness. The retrieval of the surface re;ectance will be discussed in other occasions. ∗
Corresponding author. Tel.: +81-3-6221-9000; fax: +81-3-6221-9192. E-mail addresses: [email protected] (T. Tanaka), [email protected] (M. Wang).
0022-4073/04/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-4073(03)00105-5
556
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Zenith θ3
i0
i3
Equator
i1
i1
Top of Atmosphere
τ
i2
i3
Azimuth Axis
Earth Surface
i4
i5
ϕ
i’2
3
i’3
i2 i0
Fig. 1. Geometry of radiative transfer.
We show the geometry of the radiative transfer process in Fig. 1. The layer has the vertical optical thickness . The direction of the input solar radiance ˜i0 is the lower bound, the direction of observation by satellites ˜i1 is the upper bound, the direction of the transmitted radiation to the Earth surface ˜i4 is the lower bound, the intermediate direction ˜i2 is the lower bound, and the other intermediate direction ˜i3 is the upper bound. We employ the polar coordinate system in which the zenith angle is measured from the zenith and the azimuth angle ’ is measured from the projection of the solar direction onto the equator. We denote the cosine of the zenith angle and the azimuth angle of the direction ˜in as n and ’n , respectively. We usually employ the generalized re;ectance instead of radiance I , I ; (1) (;˜i1 ;˜i0 ) = F0 |0 | where F0 is the solar irradiance. Since is considered suDciently small to 1, we have the 5rst approximation P(i1 ; i0 ) ; (2) = 4|0 |1 where P(˜i1 ;˜i0 ) is the scattering phase function. The 5rst approximation or the single scattering approximation has been applied to many analyses. Since the advent of Earth observing satellites, in particular the retrieval of the ocean chlorophyll concentration from coastal zone color scanner (CZCS), we need a more accurate formula. Many varieties of calculating methods have been proposed [1,2]. Most of them comprise the forward calculation and the look-up table. They do not give us the explicit form of by . In this paper we seek the quadratic form of in as an explicit form of by , = c1 + c2 2 :
(3)
We employ the iterative integration of Chandrasekhar’s integral equation to solve the inversion problem of the radiative transfer. The radiative transfer process in the atmosphere is governed by
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
557
the integral equation, which was 5rst introduced by Chandrasekhar [3]. He also obtained the explicit solution of radiative transfer in the second approximation for the isotropic plane-parallel atmosphere. We have obtained the explicit solution of the third approximation in the isotropic plane-parallel atmosphere [4]. In this paper we extend the second approximation developed by Chandrasekhar to the anisotropic atmosphere. In Chandrasekhar’s integral equation, there are two unknown functions: the scattering function, S(;˜i1 ;˜i0 ) and the transmitted function, T (;˜i4 ;˜i0 ), for which we seek the solutions (Section 2). We solve this simultaneous integral equation by successive iteration. The 5rst approximation corresponds to the single scattering approximation. Substituting the 5rst approximation into the integration in the original equations, we obtain the second iteration. Adding the 5rst approximation to the second iteration, we obtain the second approximation for S2 (;˜i1 ;˜i0 ) and T2 (;˜i4 ;˜i0 ). We can continue this iterating process and obtain series Sn (;˜i1 ;˜i0 ) and Tn (;˜i4 ;˜i0 ). If the series Sn (;˜i1 ;˜i0 ) and Tn (;˜i4 ;˜i0 ) converge as n approaches to ∞, it satis5es Chandrasekhar’s integral equation (Section 3.1). We decompose the second iteration into a power series expansion in (−) and obtain the quadratic form of the scattering function (Section 3.2). In the process to evaluate S2 (;˜i1 ;˜i0 ), we often encounter a sort of integration, U0 (;˜i0 ) =
0
1
3 |0 | d3 (1 − e−=3 −=|0 | ) : 3 + |0 | 3
(4)
Expanding the term in the integral above into a power series in (−) and integrating inde5nitely those, we obtain 1 (−)2 1 1 d3 −(−) − + − ··· 3 |0 | 2! 3 0 1 log 3 (−)2 1 = −log 3 (−) − − + − ··· : − 3 |0 | 2! 3
U0 (;˜i0 ) =
(5)
Substituting 3 = 1 at the upper integral limit, we can obtain the upper integral value. And the coeDcients of nth power in are separately integrated. However at the lower limit, each integrated term becomes ∞, or has singularities at 3 = 0. For the lower limit, we must sum up all the terms on theinterval (; 1) and then make approach to 0. Brie;y we cannot change the order of “lim” and “ ” at the lower integral limit (Section 3.3). The problem of the radiative transfer is, thus, deduced to evaluate the value of the series expansion at the singularity point. The authors evaluated values of several series expansions at the singularity point, which are necessary to obtain the second and third approximation for the isotropic atmosphere [4]. We show that, for the anisotropic atmosphere, we can exchange the order of “lim” and “ ” for the higher powers of . If we can exchange the order, we can integrate each power (−) separately (Section 3.3). Based on this consideration, we show the coeDcients of the quadratic form of the second iteration expressed as the surface integrations of products of the scattering phase functions on the half unit sphere (Sections 3.4 and 3.5).
558
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Finally, in Section 4, we evaluate the second approximations numerically for the Rayleigh scattering and the maritime aerosol atmosphere. And the results are compared with the single scattering approximation and the exact solutions.
2. Radiation transfer process The scattered radiance at the TOA, I (0;˜i1 ), is expressed with the incident intensity I (0; i0 ) and the scattering function S(; i1 ; i0 ) as 1 S(;˜i1 ;˜i0 )I (0;˜i0 ) d0 ; (6) I (0;˜i1 ) = 4 1 L where 0 is the solid angle subtended around the input direction ˜i0 and integral domain is the lower half of the unit sphere. In the same manner the intensity transmitted to the Earth surface is expressed as 1 I (0;˜i4 ) = I (0;˜i4 ): T (;˜i4 ;˜i0 )I (0;˜i0 ) d0 + exp − (7) 4 |4 | L |4 | Chandrasekhar’s integral equation is expressed [3], 1 1 S(;˜i1 ;˜i0 ) + 1 |0 | (−) (−) P(˜i1 ;˜i0 ) + = 1 − exp 1 |0 | d3 d3 P(˜i1 ;˜i3 )S(;˜i3 ;˜i0 ) − exp − T (;˜i1 ;˜i3 )P(˜i3 ;˜i0 ) + 4 | | 4 3 0 3 U U d2 d2 − exp − + S(;˜i1 ;˜i2 )P(˜i2 ;˜i0 ) P(˜i1 ;˜i2 )T (;˜i2 ;˜i0 ) 4 |2 | 1 4 |2 | L L d2 d3 + S(;˜i1 ;˜i2 )P(˜i2 ;˜i3 )S(;˜i3 ;˜i0 ) 4 |2 | 4 3 U L d2 d3 − T (;˜i1 ;˜i3 )P(˜i3 ;˜i2 )T (;˜i2 ;˜i0 ) ; 4 |2 | 4 3 U L
1 1 − T (;˜i4 ;˜i0 ) |4 | |0 | − exp − P(˜i4 ;˜i0 ) = exp − |0 | |4 |
(8)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2 d3 ˜ ˜ ˜ ˜ − exp − + P(i4 ; i2 )T (; i2 ; i0 ) P(˜i4 ;˜i3 )S(;˜i3 ;˜i0 ) 4 | | | | 4 2 4 3 L U d3 d2 + exp − S(;˜i4 ;˜i3 )P(˜i3 ;˜i0 ) − T (;˜i4 ;˜i2 )P(˜i2 ;˜i0 ) |0 | 4 3 4 |2 | U L d2 d3 + S(;˜i4 ;˜i3 )P(˜i3 ;˜i2 )T (;˜i2 ;˜i0 ) 4 |2 | 4 3 U L d2 d3 − T (;˜i4 ;˜i2 )P(˜i2 ;˜i3 )S(;˜i3 ;˜i0 ) : 4 |2 | 4 3 U U
559
(9)
The cosine of the zenith angle in the lower bound directions, which has a “even index”, is negative. For the convenience of the following evaluation, we replace n by a new parameter, n− which is equal to −n for n = even. Using the new variables, the integration on the lower half sphere is changed to the integration on the upper half unite sphere, shown below −1 2 f(2 ; ··) d2 f(2 ; ·; ·)(−d2 ) d’2 = 4 |2 | 4 (−2 ) 0 0 L 2 1 d2− f(−2− ; ·; ·) d2− − = d’2 = f(− ; ·; ·) : (10) 2 4 2− 4 2− 0 0 U Since we adopt the new variable 2− and the integral domain is changed to the upper half sphere, we omit the U for the integral domain unless the domain has a special meaning. The cosine of the angle between the upper and lower bound directions is thus expressed as below cos(e; o ) = e o + (1 − e2 )1=2 (1 − o2 )1=2 cos(’e − ’o ) = −e− o + (1 − (e− )2 )1=2 (1 − o2 )1=2 cos(’e − ’− o ):
(11)
The phase function is assumed as a function of the angle mn between the directions m and n and is normalized in the whole solid angle 4 . We express the phase function P(im ; in ) as below, P(im ; in ) = P(cos mn ) =
∞
wj Pj (cos mn );
(12)
j=0
where Pj (x) is the Legendre function of the 5rst kind. The input from the upper surface is the solar irradiance, F0 , and is given with the Dirac delta function I (0;˜i) = F0 (˜i − ˜i0 ):
(13)
Substituting the solar irradiance, we obtain the generalized re;ectance from the upper layer S(;˜i1 ;˜i0 ) : (0;˜i1 ) = 41 |0 |
(14)
560
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
3. Second approximation 3.1. Azimuth integration and decomposition of second iteration The simultaneous integral equations of radiative transfer can be solved by successive iteration. The 5rst iterations are the 5rst terms on the right-hand side in Eqs. (8) and (9), expressed below S1 (;˜i1 ;˜i0 ) =
T1 (;˜i4 ;˜i0 ) =
1 1 + − 1 1
−1 (−) (−) 1 − exp P(˜i1 ;˜i0 ); + − 1 0
1 1 − + − 4 0
−1
exp − − 0
− exp − − 4
P(˜i4 ;˜i0 ):
(15)
(16)
A new function U1 (; 0 ; 1 ) is de5ned as −
U1 (; 1 ; 0 ) =
1 1 + − 1 0
−1 (−) (−) 1 − exp : + − 1 0
(17)
Using U1 (; 0− ); 1 ); T1 (;˜i4 ;˜i0 ) is rewritten as
T1 (;˜i4 ;˜i0 ) = exp − − 0
; U1 (; 4− ; (−0− ))P(˜i4 ;˜i0 ):
(18)
Substituting the 5rst iterations into the integrals in the original equations, the second iterations, PS2 (;˜i1 ;˜i0 ) is expressed PS2 (;˜i1 ;˜i0 ) =
1 1 −1 + − 1 1 d3 d3 × P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (0− ) − P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 3 4 3 d2− d2 + P(˜i1 ;˜i2 )P(˜i2 ;˜i0 )U1 (1 ) − P(˜i1 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ) 4 2− 4 2− −1 − exp − − − 1 0 − d3 − d2 ˜ ˜ ˜ ˜ ˜ ˜ ˜ ˜ + P(i1 ; i2 )P(i2 ; i0 )U1 (−0 ) × P(i1 ; i3 )P(i3 ; i0 )U1 (−1 ) 4 3 4 2− d2− d3 + P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2− P(˜i1 ;˜i3 )P(˜i3 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ; 2− )U1 (−1 ; 3 )) 4 2− −1 − exp − − − 1 0 d2− × P(˜i1 ;˜i3 )P(˜i3 ;˜i2 )P(˜i2 ;˜i0 )U1 (−0− ; 2− )U1 (−1 ; 3 )) 4 2− −
561
d3 4 3
d3 ; 4 3
(19)
where the arguments in the functions U (;˜i1 ;˜i0 ) are omitted except for the signi5cant ones. Using the addition theorem of the Legendre function, we decompose the phase function into the Bi-Legendre functions, Plm (3 ), as
∞ ∞
(2 − 0; m ) wlm Plm (1 )Plm (3 ) cos m(’3 − ’1 ) ; (20) P(˜i1 ;˜i3 ) = m=0
P(˜i3 ;˜i0 ) =
∞
l=m
(2 − 0; m )
m=0
∞
wkm Pkm (3 )Pkm (−0− ) cos m(’3 − ’0 ) ;
(21)
k=m
where 0; m and wlm are given as below, 0; m = 1 wlm = wl
(m = 0);
0; m = 0
(m = 0);
(l − m)! : (l + m)!
(22) (23)
Due to the orthogonality of the functions Plm (3 ) cos m(’3 − ’1 ) for diQerent m, we obtain the integration 2 1 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) d’3 2 0 ∞ ∞ k m m m m − m m = (2 − 0; m ) (−1) wl wk Pl (1 )Pk (0 )Pl (3 )Pk (3 ) m=0
l; k=m
×cos m(’1 − ’0 ):
(24)
In the same manner, we obtain the integration of the triple product of phase functions 2 2 d’2 d’3 P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) 2 2 0 0
2 2 ∞ ∞ − m m m (2 − 0; m ) wl Pl (1 )Pl (−2 ) cos m(’2 − ’1 ) = 0
0
m=0
l=m
562
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
× ×
∞
(2 − 0; m )
m=0
∞
∞ l=m
(2 − 0; m )
∞
m=0
=
∞
wlm Plm (−2− )Plm (3 )
cos m(’3 − ’2 )
wlm Plm (3 )Plm (−0− ) cos m(’0 − ’3 )
l=m
∞
(2 − 0; m ) cos m(’3 − ’1 )
m=0
d’2 2 d’3 2
(−1)a+b+c
a;b;c=m
{wam wbm wcm Pam (1 )Pam (2 )Pbm (2 )Pbm (3 )Pcm (3 )Pcm (0 )}
:
(25)
We de5ne an integration, Ql;mk (; 0− ), as Ql;mk (; 0− )
=
1
Plm (3 )Pkm (3 )U1 (; 1 ; 0− )
0
d3 : 3
Using Ql;mk (; 0− ), we obtain PS2 ,
1 PS2 (;˜i1 ;˜i0 ) = 2
×
1 1 + − 1 0 ∞
−1 ∞
(2 − 0; m ) cos m(’1 − ’0 )
m=0
wlm wkm {Plm (1 )Pkm (0− )(−1)k }{Ql;mk (; 0− ) − (−1)l+k Ql;mk (; −0− )}
l; k=m
∞
+
l; k=m
wlm wkm {Plm (1 )Pkm (0− )(−1)l }{Ql;mk (; 1 ) − (−1)l+k Ql;mk (; −1 )}
− exp − − − 1 0 ×
∞
−1
wlm wkm {Plm (1 )Pkm (0− )(−1)l }{Ql;mk (−0− ) + (−1)l+k Ql;mk (−1 )}
l; k=m
+
∞ 1 (−1)a+b+c wam wbm wcm Pam (1 )Pcm (0− )Qa;m b (; 1 )Qb;m c (; 0− ) 2 a;b;c=m
−
∞ 1 (−1)b wam wbm wcm Pam (1 )Pcm (0− )Qa;m b (; −1 )Qb;m c (; −0− ) 2 a;b;c=m
(26)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
− exp − − − 0 1
563
−1
∞ 1 − b m m m m m − m m × (−1) wa wb wc Pa (1 )Pc (0 )Qa; b (; −1 )Qb; c (; −0 ) : 2
(27)
a;b;c=m
In the above equation, the sign of each term is determined by the choice of l; k; a; b; c. This determination originates from the even/odd characteristic of the Bi-Legendre functions. 3.2. Zenith integration of second iteration The second iteration PS2 (; i1 ; i0 ) is expressed in Eq. (27) and its component, Ql;mk (; ), has a power series expansion in . In order to obtain the quadratic form of PS2 (; i1 ; i0 ), we decompose Ql;mk (; 0− ) into the power series expansion in (−). First, we integrate U20 (; 0− ) and this is evaluated in the author’s previous work [4] given below, − 1 1 0 d − − d 0 U2 (; i0 ) = 1 − exp − − − U1 (; ; 0 ) = 0 0− + 0 0 ∞ ∞ ∞ (−=0 )r (−=0 )r (−)n − ; (28) = (−) (C + eax()) (r + 1)r! r=0 n=0 (r + 1 + n)n! (r + 1)r! r=0 where the constant C and the function eax(p) are given as C = % + log = 0:5772 · · · + log ;
eax(p) =
∞ (−1)n pn n=1
nn!
;
where % is the Euler constant. The function U20 (; 0 ) is truncated as 1 1 (−)2 − 0 + ···: U2 (; 0 ) = (C − 1)(−) + 1 + C − 2 0− 2! We integrate the moment integration, U2n (; 0− ), as below, 1 d3 − n U2 (; 0 ) = 3n U1 (; 3 ; 0− ) : 3 0
(29)
(30)
(31)
(32)
U2n (; 0− ) is evaluated by partial integration and recurrence relationships (refer to Appendix B). The coeDcients of the series expansion of U2n (; 0− ) in the power (−) are also polynomials of 1=0− .
564
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
The truncated form of U2n (; 0− ) for n ¿ 1 is expressed below. U2n (; 0− ) = (−1)n (−) + u2;n 0
(−)2 (−)2 1 + u2;n 1 + ···; 2! 2! 0−
(33)
where ul;n k is the coeDcient of (−)n =n!(0− )k and is given in Appendix B. It is noted that the coeDcient of the 5rst power (−) does not include 0− . Substituting U2n into Eq. (26), we obtain Ql;mk (; ) as a power series expansion in (−), Ql;mk (; 0 ) =
1 mlk 0
=
− n amlk n 3 U1 (; 3 ; 0 )
n=0 mlk
n amlk n u1; 0
n=0
+
mlk
(−) +
n amlk n u2; 1
n=0
= q1;mlk0 (−) + q2;mlk0
d3 ; 3
mlk
n amlk n u2; 0
n=0
(−)2 2!
(−)2 1 + ···; 2! 0−
(−)2 (−)2 1 + ···; + q2;mlk1 2! 2! 0−
(34)
is the coeDcient of the n degree power of the product of the Bi-Legendre functions where amlk n m mlk Pl (3 ) and Pkm (3 ). The coeDcient qd; e is given as below, mlk qd; e =
ml k
n amlk n ud; e :
(35)
n=0
It is noted that the coeDcient of the 5rst power (−) in Ql;mk (; 0 ) does not include 0− . Substituting Eq. (34) into the 5rst term in the bracket of Eq. (27), we obtain following: Ql;mk (; 0− ) − (−1)l+k Ql;mk (; −0− ) = q2;mlk1
(−)2 2 + ··· 2! 0−
=2q1;mlk0 (−) + 2q2;mlk0
(l + k = even)
(−)2 + ··· 2!
(l + k = odd):
And for the case (l + k = odd), adding the second term to the above, we obtain following, Plm (1 )Pkm (0 )(−1)k {Ql;mk (0 ) − (−1)l+k Ql;mk (−0 )} +Plm (1 )Pkm (0 )(−1)l {Ql;mk (1 ) − (−1)l+k Ql;mk (−1 )}
(36)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
=0
(−)2 + ···: 2!
565
(37)
From the 5rst and second terms in Eq. (27), the 5rst power in (−) vanishes and the second power emerges only for the case (l + k = even). From the third terms in the bracket of Eq. (27), we only need the 5rst power in (−), because it is multiplied by another 5rst power in (−). Substituting Eq. (34) into this term, we obtain Ql;mk (−0 ) + (−1)l+k Ql;mk (−1 ) = 2q1;mlk0 (−) + · · · ; = 0(−) + · · ·
(l + k = even);
(l + k = odd):
(38)
From the fourth and 5fth term, substituting Eq. (34) we obtain following: ∞ 1 (−1)a+b+c wam wbm wcm Pam (1 )Pcm (0 )Qa;m b (1 )Qb;m c (0 ) 2 a;b;c=m
−
∞ 1 (−1)b wam wbm wcm Pam (1 )Pcm (0 )Qa;m b (−1 )Qb;m c (−0 ) 2 a;b;c=m
=0(−)2 + · · · (a + c = even) ∞ b+1 m m m m m mab mbc = (−1) wa wb wc Pa (1 )Pc (0 )q1; 0 q1; 0 (−)2 + · · ·
(a + c = odd):
(39)
a;b;c=m
From the sixth term in the bracket of Eq. (27), no (−)2 term emerges. Thus we obtain the truncated second iterations, PS2 (;˜i1 ;˜i0 ) =
∞
(2 − 0; m )
m=0
(−1)k wlm wkm (q2;mlk1 − 2q1;mlk0 )Plm (1 )Pkm (0 )
(l; k)1
+2
∞
1 1 + − 1 0
−1 ∞
×cos m(’1 − ’0 )
mbc m m (−1)b+1 wam wbm wcm q1;mab 0 q1; 0 Pa (1 )Pc (0 )
(a;b;c)2
(−)2 + ···; 2!
(40)
where (l; k)1 denotes l; k ¿ m and l + k is even and (a; b; c)2 denotes a; b; c ¿ m; a + c is odd. It is noted that the truncated second iteration does not have a 5rst power in (−). It is noted that, the fourth and 5fth terms vanish for wa = 0 and wa = 0 for all possible combinations of a + c = odd such as Rayleigh scattering.
566
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
3.3. Separated integration In this subsection we consider changing order of lim and begin with the integration of U2n (; i7 ), 1 7 n d n 1 − exp − − U2 (; i7 ) = 7 7 + 0
in the integration of the PS2 . We
1 ∞ 1 m−1 n d (−)m 1 + : →0 m! 7 m=1
= −lim
On the integral interval (; 1), we interchange the order of into two parts: one is m 6 n and the other is m ¿ n, m−1 n 1 (−)m n n− m +1 d U2 (; i7 ) = − 7 m! 0 m=1 n
−(−) lim
→0
∞ m=1
1
(−)m (m + n)!
1 1 + 7
(41) 1
m−1
and
+1 7
n
and divide the summation
d :
(42)
The terms of the 5rst part of the summation (m 6 n) are ordinary integrations because the minimum powers in the in the integrand are equal toor greater than 0 and have no singularity at = 0. We can interchange the order of “lim” and “ ” for these terms. If the order is possible to exchange, we can “pick-up” each power in (−) separately. In other words, we can integrate the powers in (−) separately. We call this part as the separated integration part (SIP). The terms of the second part (m ¿ n) involve negative powers in or log , and have singularities at =0. So we cannot interchange the order of “lim” and “ ” and therefore we cannot integrate the terms separately. We call the second part as the non-separated integration part (NIP). It is needless to say that the NIP converges as approaches to 0, because it is equal to the converged U2n (; i6 ) minus the bounded SIP. We apply the above observation to Eq. (27). 3.3.1. The >rst and second terms Due to the discussion in the previous subsection, we can select only the case (l + k = even) and pick up the second power in (−)2 only. Plm (3 )Pkm (3 ) is a polynomial of 3 and its lower powers for l + k = even are given below (refer to Appendix B), Plm (3 )Pkm (3 ) = a0 + a2 32 + · · · ; = a2 32 + a4 34 + · · ·
(m + l = even); (m + l = odd):
(43)
The zeroth power of (−); a0 , is a NIP and the other higher powers of (−) are SIT. Taking account of a0 = Plm (0)Pkm (0), the integration above is then divided into SIP and NIP
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
as below, Ql;mk (; 0− ) = lim Plm (0)Pkm (0) →0 +
1
0
1
U1 (; 3 ; 0− )
567
d3 3
(Plm (3 )Pkm (3 ) − Plm (0)Pkm (0))U1 (; 3 ; 0− )
d3 : 3
(44)
Substituting the second power of (−) in the polynomial U20 (p0 ) into NIP and truncating SIP into the second power of (−), we obtain 1 1 − m Plm (0)Pkm (0) Ql; k (; 0 ) = 1 + C − 2 0− 1 1 1 d3 (−)2 m m m m − (Pl (3 )Pk (3 ) − Pl (0)Pk (0)) + + ···: (45) 3 0− 3 2! 0 For the case (m + l = odd) Eq. (45) holds because Plm (0)Pkm (0) = 0. 3.3.2. The third term We select for the case (l + k = even) and pick up the 5rst power of (−) only. We divide the integration into NIP and SIP, in the same as in Eq. (44). Substituting the 5rst power of (−) in the polynomial U20 (p0 ) into NIP and truncating SIP into the 5rst power of (−), we obtain Ql;mk (; −0− ) = (C − 1)Plm (0)Pkm (0) −
0
1
(Plm (3 )Pkm (3 )
−
Plm (0)Pkm (0))
d3 (−) + · · · + · · · : 3
(46)
For the case (m + l = odd) Eq. (46) holds because Plm (0)Pkm (0) = 0. 3.3.3. The fourth and >fth terms We integrate following for the case (a + c = odd), 1 1 d− d3 Pam (2− )Pbm (2− )Pbm (3 )Pcm (3 )U1 (1 ; 2− ) −2 U1 (3 ; 0− ) : 3 2 0 0
(47)
The lowest powers of 2− and 3 for the product Pam (2− ) · · · are given below (refer to Appendix B), Pam (2− )Pbm (2− )Pbm (3 )Pcm (3 ) =c12 (2− )(3 )2
(m + b = odd; a + b = even);
=c21 (2− )2 (3 )
(m + b = odd; a + b = odd);
=c01 (3 )
(m + b = even; a + b = even);
=c01 (2− )
(m + b = even; a + b = odd):
(48)
568
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
For (m + b = odd) the integration is the product of two SIP: one is for 2− and the other is for 3 . For (m + b = even) the integration is a product of one SIP and NIP. Then the integration is expressed for the four cases in the equation above: 1 1 d2− d3 m − m − m m = (−)2 Pa (2 )Pb (2 )Pb (3 )Pc (3 ) − 3 0 0 2 + · · · (m + b = odd); 1 d3 = −(C − 1)Pam (0)Pbm (0) Pbm (3 )Pcm (3 ) 3 0 1 1 d2− d3 m − m − m m m m (−)2 + {Pa (2 )Pb (2 ) − Pa (0)Pb (0)}Pb (3 )Pc (3 ) − 3 0 0 2 + · · · (m + b = even; a + b = even); 1 d− Pam (2− )Pbm (2− ) −2 = −(C − 1)Pbm (0)Pcm (0) 2 0 1 1 d− d3 (−)2 + Pam (2− )Pbm (2− ){Pbm (3 )Pcm (3 ) − Pbm (0)Pcm (0)} −2 2 3 0 0 + ···
(m + b = even; a + b = odd):
(49)
3.4. Surface integration We can bring back the idea of division of integration to the surface integration on the unit sphere expressed in Eq. (19). This division is possible, because the integrations in Eq. (19) are expressed as the summation of the integrations of the products of Bi-Legendre functions shown in Eq. (27). From the ninth integration in Eq. (19) no second power in (−)2 emerges. 3.4.1. The >rst four integrations The 5rst integration is divided into SIP and NIP expressed as below, d3 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (; 0− ; 3 ) 4 3 d3 = {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )}U1 (; 0− ; 3 ) 4 3 d3 + lim {P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )}U1 (; 0− ; 3 ) : 3 →0 4 3
(50)
where ˜i 3 denotes the projected direction of ˜i3 onto the equator in the unit sphere. Substituting the second power of (−) in the polynomial U20 (p0 ) and the truncated term of the second power into
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Eq. (50), we can obtain d3 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (; 0− ; 3 ) 4 3 d3 (−)2 1 1 =− {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} + 2! 0− 3 4 3 2 2 1 1 d’ 3 (−) ˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) + 1+ C − P( : 2 0− 4 2! 0
569
(51)
Substituting the above equations into the 5rst and second integrations in the bracket in Eq. (19), we obtain d3 − d3 ˜ ˜ ˜ ˜ P(i1 ; i3 )P(i3 ; i0 )U1 (0 ) − P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 3 4 3 (−)2 1 1 d 3 ˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} {P( =− + 4 3 2! 0− 1 2 2 1 1 1 ˜ d’3 (−) ˜ ˜ ˜ : (52) + C− + − P(i1 ; i 3 )P(i 3 ; i0 ) 2 1 0 4 2! 0 We evaluate the third and fourth integrations in Eq. (19) in the same manner. 3.4.2. The >fth and sixth integrations The 5fth integration is divided into SIP and NIP as below, d3 − exp − − − −1 P(˜i1 ;˜i3 )P(˜i3 ;˜i0 )U1 (−1 ) 4 0 1 3 2 1 1 d’3 =−2 (C − 1) + − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) 1 0 4 0 (−)2 d 3 + ···: − {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} 4 3 2!
(53)
The equation above holds for the two cases described in (46). We can evaluate the sixth integration in Eq. (19) in the same manner. 3.4.3. The seventh and eighth integrations The seventh integration is divided into SIP and NIP as below, d2− d3 P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3
570
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
=
[P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )
d2− d3 −P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )]U1 (1 ; 2− )U1 (3 ; 0− ) 4 2− 4 3 + P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )U1 (1 ; 2− )U1 (3 ; 0− ) + =
d2− d3 P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )U1 (1 ; 2− )U1 ((3 ; 0− )] 4 2− 4 3
{P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 )
− P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 )}
d’2 d3 P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) 4 4 3
− (C − 1) − (C − 1)
d2− d3 4 2− 4 3
d2− d’3 ˜ ˜ ˜ ˜ ˜ ˜ (−)2 : P(i1 ; i2 )P(i2 ; i 3 )P(i 3 ; i0 ) 4 2− 4
(54)
The equation above holds for all the cases described in (49). We can evaluate the eighth integration in Eq. (19) in the same manner. 3.5. The result of the second approximation Expanding the 5rst approximated scattering function into a power series expansion in and truncating up to the second degree, we obtain 2 1 1 ˜ ˜ P(˜i1 ;˜i0 ): + (55) S1 (; i1 ; i0 ) = − 1 |0 | 2! The second iteration is expressed by the integration of the products of the phase functions as below, PS2 (;˜i1 ;˜i0 ) = (3 − 2C)Ie (˜i1 ;˜i0 ) + Iu (˜i1 ;˜i0 ) + Il (˜i1 ;˜i0 ) +2
1 1 + − 1 0
−1
(−)2 : {(C − 1)(Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) + Iul2 (˜i0 ;˜i1 )} 2!
The coeDcients are evaluated on the half unit sphere (we return the variable 2− to 2 ), d3 ˜ ˜ Iu (i1 ; i0 ) = {P(˜i1 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 )} ; 4 3 U
(56)
(57)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
d2 ; {P(˜i1 ;˜i2 )P(˜i2 ;˜i0 ) − P(˜i1 ;˜i 2 )P(˜i 2 ;˜i0 )} 4 |2 | L 2 d’3 ˜ ˜ Ie (i1 ; i0 ) = ; P(˜i1 ;˜i 3 )P(˜i 3 ;˜i0 ) 4 0 2 d’2 d3 Iuu (˜i1 ;˜i0 ) = {P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i1 )} ; 16 2 3 U 0 2 d’3 d2 ˜ ˜ ; Ill (i1 ; i0 ) = {P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i1 )} 2 | | 16 2 L 0 ˜ ˜ Iul (i1 ; i0 ) = [{P(˜i1 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i3 )P(˜i3 ;˜i1 )}; Il (˜i1 ;˜i0 ) =
571
U
(58) (59) (60) (61)
L
− {P(˜i1 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i0 ) − P(˜i0 ;˜i 2 )P(˜i 2 ;˜i3 )P(˜i3 ;˜i1 )}; d2 d3 −{P(˜i1 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i0 ) − P(˜i0 ;˜i2 )P(˜i2 ;˜i 3 )P(˜i 3 ;˜i1 )}] : 16 2 |2 |3
(62)
Adding the 5rst approximation to the second iteration, we obtain the second approximation, 1 1 (−)2 ˜ ˜ ˜ ˜ P(˜i1 ;˜i0 ) + S2 (; i1 ; i0 ) = P(i1 ; i0 ) − 1 |0 | 2!
−1 1 1 (−)2 + (Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) log + −2Ie (˜i1 ;˜i0 ) + 2 1 |0 | 2! + Iu (˜i1 ;˜i0 ) + Il (˜i1 ;˜i0 ) + (3 − 2%)Ie (˜i1 ;˜i0 ) +2
1 1 + 1 |0 |
−1
{(% − 1)(Iuu (˜i1 ;˜i0 ) + Ill (˜i1 ;˜i0 )) + Iul (˜i1 ;˜i0 )}
(−)2 : 2!
(63)
From the above equation, we summarize the characteristics of the second approximation. (1) The second approximation is not a quadratic equation of the optical thickness but a quadratic equation with a 2 log(). (2) The second iteration PS2 (;˜i1 ;˜i0 ) does not have the 5rst power of . (3) There exist three kinds of terms which are characterized by the number of multiplication of the scattering phase function: single scattering, double scattering and triple scattering. The single scattering term comes from the 5rst approximation. The surface Ie ; Il ; Iu are the double scattering terms and Iuu ; Ill ; Iul are the triple scattering terms. (4) For the atmosphere with a scattering phase function which does not have the odd order Legendre functions in its series expansion in the cosine of the angle, such as the Rayleigh scattering, there are no triple scattering terms in the second approximation.
572
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 103 Rayleigh Scattering Aerosol: M80 at 443 nm
Phase Function
102
101
100
10-1
10-2 0
20
40
60
80
100
120
140
160
180
Scattering Angle (Deg.) Fig. 2. Scattering phase functions.
4. Calculation of scattering To evaluate the accuracy of the second approximation discussed in the previous sections, we carry out calculations for the Rayleigh and aerosol scattering atmospheres for various cases and compare with the exact solutions. To appreciate the improvement of accuracy from single scattering, we calculate the single scattering re;ectance as (+)!(+)P(i1 ; i0 ) (+;˜i1 ;˜i0 ) = ; 41 |0 |
(64)
where + is the wavelength and !(+) is the single scattering albedo. Fig. 2 illustrates the phase function of air molecules (Rayleigh scattering) as well as for a typical maritime aerosol with relative humidity (RH) of 80.0% (denoted as M80) at wavelength 443 nm. The maritime aerosol model has a single scattering albedo of 0.9929, and the Rayleigh scattering model has a single scattering albedo of 1.0. Note that, in contrast to Rayleigh scattering function which is symmetry in the forward and backward scatterings, the maritime aerosol predominantly scatters in the forward direction. 4.1. The exact solutions The successive-order-of-scattering (SOS) [5] method is used in solving the radiative transfer equation (RTE) for the Rayleigh scattering atmosphere and the aerosol scattering atmosphere. The SOS code for the ocean–atmosphere system was developed for the atmospheric correction algorithm for the ocean color sensors [6,7]. This code computes the upward radiance at the top and the downward radiance at the base of the medium for the ocean–atmosphere system. The code is capable of yielding radiances that are accurate to nearly 0.1% [6]. In the following two subsections, the re;ectance computed by the single scattering and the second approximation are compared with the exact solutions obtained using the SOS code.
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 80 60
20 Rayleigh Atmosphere, θ = 45o 412 nm (S) 443 nm (S)
15
412 nm (D) 443 nm (D)
20 0 -20 -40
490 nm (S) 510 nm (S)
490 nm (D) 510 nm (D)
5 0 -5 -10
-60
-15 ∆φ = 0o -60
(a)
-40
∆φ = 180o -20
0
20
40
∆φ = 180o
∆φ = 0o 60
-20 -80
80
-60
-40
(b)
Solar Zenith Angle (Deg.)
-20
0
20
40
60
80
60
80
Solar Zenith Angle (Deg.) 5
10 Rayleigh Atmosphere, θ = 45 555 nm (S) 670 nm (S)
Rayleigh Atmosphere, θ = 45o
o
765 nm (S) 865 nm (S)
555 nm (D) 670 nm (D)
Error (%)
5
Error (%)
Rayleigh Atmosphere, θ = 45o
10
Error (%)
Error (%)
40
-80 -80
573
0
765 nm (D) 865 nm (D)
0
-5
(c)
-60
-40
-20
0
20
40
Solar Zenith Angle (Deg.)
∆φ = 180o
∆φ = 0o
∆φ = 180o
∆φ = 0o -10 -80
60
-5 -80
80
(d)
-60
-40
-20
0
20
40
Solar Zenith Angle (Deg.)
Fig. 3. Calculation of generalized re;ectances for Rayleigh scattering.
4.2. The Rayleigh scattering atmosphere Fig. 3 provides a quantitative evaluation of accuracy using the second approximation for a Rayleigh scattering atmosphere for the various solar-sensor geometries and for the eight wavelengths at 412, 443, 490, 510, 670, 765, and 865 nm. The Rayleigh optical thicknesses corresponding to the eight wavelengths from 412 to 865 nm are 0.3185, 0.2361, 0.1560, 0.1324, 0.0938, 0.0436, 0.0255 and 0.0155, respectively. In all computations, the Rayleigh atmosphere is treated as a plane-parallel medium without a surface boundary. Figs. 3(a) – (d) show the error (%) in the computed upward re;ectance as a function of the solar zenith angle and for a sensor viewing angle of 45◦ . The left and right parts of each plot in Fig. 3 corresponds to the relative azimuth angles of 0◦ and 180◦ (principal scattering plane), respectively. Results in Fig. 3 show the signi5cant improvement in accuracy using the second approximation compared with the single scattering formula. Using the new formula, the errors are all within 1% for the red and near-infrared wavelengths, while errors are within nearly 10% for the blue bands for
574
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577 10
5
5 0
Aerosol, M80 model at 443 nm τa = 0.10, θ = 45o
-5
Error (%)
Error (%)
0
Aerosol, M80 model at 443 nm -10
τa = 0.05, θ = 45o Single-Scattering Approx. Double-Scattering Approx.
-15
-5 -10 -15 Single-Scattering Approx. Double-Scattering Approx.
-20 -25
∆φ = 0o -20 -80
-60
(a)
-40
∆φ = 180o -20
0
20
40
∆φ = 0o 60
-30 -80
80
Aerosol, M80 model at 443 nm
-20
0
20
40
60
80
Solar Zenith Angle (Deg.)
20
Aerosol, M80 model at 443 nm τa = 0.30, θ = 45o
τa = 0.20, θ = 45o
10
Error (%)
0
Error (%)
-40
30
20
10
-60
(b)
Solar Zenith Angle (Deg.)
∆φ = 180o
-10
0 -10 -20
-20
(c)
Single-Scattering Approx. Double-Scattering Approx.
-40 ∆φ = 0o
-40 -80
-30
Single-Scattering Approx. Double-Scattering Approx.
-30
-60
-40
∆φ = 0o
∆φ = 180o -20
0
20
40
Solar Zenith Angle (Deg.)
60
-50 -80
80
(d)
-60
-40
∆φ = 180o -20
0
20
40
60
80
Solar Zenith Angle (Deg.)
Fig. 4. Calculation of generalized re;ectances for maritime aerosol scattering.
zenith angles less than 60◦ . It is important to note that, with the second approximation, the error in a given wavelength (optical thickness) is almost independent of the solar and sensor viewing geometry. 4.3. The maritime aerosol atmosphere A realistic maritime aerosol model as in Fig. 2 is used to study the eDcacy of the double scattering approximation in computing the upward re;ectance. The double scattering approximation is the second approximation without triple scattering. Computations for various coeDcients are much more involved due predominantly to the aerosol forward scattering characteristics. Fig. 4 shows examples of error (%) for sensor viewing angle of 45◦ as a function of the solar zenith angle for the principal scattering plane. Fig. 4(a) – (d) are the results for the aerosol optical thickness of 0.05, 0.1, 0.2, and 0.3, respectively. Again, we include results from the single scattering approximation given by
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
575
Eq. (64) for comparison. Similar to the Rayleigh atmosphere case, the aerosol atmosphere is considered as one plane-parallel layer without the surface boundary. The results in Fig. 4 show signi5cant improvement in accuracy using the double scattering approximation compared with the single scattering formula, in particular, for the part of the results with P- = 0◦ . It is interesting to note that the scattering angles for the P- = 0◦ part (the left part of the plot) vary from 55◦ to 135◦ , (all forward scattering), while the scattering angles for P- = 180◦ (the right part of the plot) change from 135◦ to 180◦ (all backward scattering). We are usually more interested in the results with the large scattering angles because measurements with the large scattering angles are generally accessible for satellite remote sensing. With double scattering approximation and scattering angles ¿ 100◦ , the re;ectance error is usually within about 1% for the aerosol optical thickness of 0.05, while the error is within about 4% for the optical thickness of 0.1. As expected, for the turbid atmosphere, the error increases. The error is proportional to the slant path of the optical thickness, i.e., . =cos . The aerosol re;ectance errors are within about 7% and about 12% for the optical thickness of 0.2 and 0.3, respectively.
5. Conclusion We obtain a new approximated solution in the radiative transfer problem in the anisotropic plane-parallel atmosphere by iterating integration of Chandrasekhar’s integral equation. By the second approximation, the scattering function S(;˜i1 ;˜i0 ) is expressed as a quadratic equation in the optical thickness, , with terms 2 log(). It is noted that the linear term in in the second approximation is identical to the 5rst approximation. This implies that the iterative integration method improves the accuracy from the 5rst approximation and approaches to the solution of radiative transfer. The second approximation is not a pure quadratic equation but a quadratic equation with terms 2 log(). The “log” term comes from the singularity at the lower integral limit. It is possible that, by expanding log around = 0:0, we can obtain the quadratic equation form of the scattering function S(;˜i1 ;˜i0 ). Rather we leave it as “log” to insist the validity of the expression in the range ¿ 0:0. The coeDcients of 2 and 2 log() are evaluated by surface integrations of products of the phase functions on the half unit sphere. The generalized upward re;ectances numerically evaluated by the second approximation for the Rayleigh atmosphere and the maritime aerosol atmosphere are very much closer to the exact solution than the single scattering approximation. In some cases, such as the small optical thickness, they are suDciently close to the exact solution.
Acknowledgements This work is done during Tanaka’s stay at NASA Goddard Space Flight Center as a visiting scientist. The authors expresses their sincere appreciation to the SIMBIOS project people for their assistances.
576
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
Appendix A. Derivation of function U2n (; ) We de5ne a new function Wn () as 1 n d3 : 3 exp − Wn () = 3 0
(A.1)
W0 () is evaluated 1 1 1 d3 = 3 exp − exp − − 3 exp − d3 W0 () = 3 3 0 3 0 0 = exp(−) − (−)(C + eax()):
(A.2)
The recurrence relation of Wn () is given as Wn () =
(−) exp(−) + Wn−1 (): n+1 n+1
(A.3)
We can evaluate Wn () by the recurrence relation and the initial function W0 (). Using Wn (), we can derive the recurrence relation for U2n (; ), 1 n− 1 n − exp − Wn−1 − U2 (0 ) : (A.4) U2 (; 0 ) = 0 n 0 Truncating the series expansion up to the second degree in (−), we obtain U21 (; 0 ) = −(−) − U2n (; 0 ) = −
3(−)2 (−)2 1 (−)2 ; +C − 4 2 2 0
(−)2 1 (−)2 (−) − − n 2(n − 1) 2n 0
(n ¿ 2):
Appendix B. Lower powers of the product of the Bi-Legendre function m = e;
l = e;
k = e;
Plm (3 )Pkm (3 ) = a0 + a2 32 ;
m = e;
l = o;
k = o;
Plm (3 )Pkm (3 ) = a2 32 + a4 34 ;
m = e;
l = e;
k = o;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = e;
l = o;
k = e;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = o;
l = e;
k = e;
Plm (3 )Pkm (3 ) = a2 32 + a4 34 ;
m = o;
l = o;
k = o;
Plm (3 )Pkm (3 ) = a0 + a2 32 ;
m = o;
l = e;
k = o;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 ;
m = o;
l = o;
k = e;
Plm (3 )Pkm (3 ) = a1 3 + a3 33 :
(A.5) (A.6)
T. Tanaka, M. Wang / Journal of Quantitative Spectroscopy & Radiative Transfer 83 (2004) 555 – 577
577
References [1] [2] [3] [4] [5] [6]
Hansen JE, Travis LD. Light scattering in planetary atmosphere. Space Sci Rev 1974;16:527–610. Goody RM, Yung YL. Atmospheric radiation. New York, Oxford: Oxford University Press, 1989. Chandrasekhar S. Radiative transfer. Mineola, NY: Dover Publishing, Inc., 1960. Tanaka T. Integration of Chandrasekhar’s integral equation. JQSRT 2003;76:121–44. van de Hulst HC. Multiple light scattering. New York: Academic Press, 1980. Gordon HR, Wang M. Surface roughness considerations for atmospheric correction of ocean color sensors. I: the Rayleigh scattering component. Appl Opt 1992;31:4247–60. [7] Gordon HR, Wang M. Retrieval of water-leaving radiance and aerosol optical thickness over the oceans with SeaWiFS preliminary algorithm. Appl Opt 1994;33:443–52.