Inverse scattering problem with two flux methods

0735-1933/93 $6.00 + .00 CopyrightO1993 Pergamon Press Ltd.

INT. COMM. HEAT MASS TRANSFER Vol. 20, pp. 585-596, 1993 Printed in the USA

INVERSE SCATTERING

PROBLEM

W I T H TWO F L U X M E T H O D S

Jen-Hui Tsai Industrial Technology Research Institute Energy & Resources Laboratories System Engineering Lab. Chutung, 310, Taiwan, R.O.C. (Communicated by J.P. Hartnett and W.J. Minkowycz) ABSTRACT The inverse two flux (S---S, S.--P and C-C) approximations are presented for determining single scattering albedo and asymmetry factor from the hemispherical transmittance and reflectance data for an anisotropically scattering slab of finite thickness subjected to isotropic incident radiation at one of its boundary surfaces. Analytical expressions of the reflectance and transmittance are derived with the two flux methods and appfied to determine the parameters with the use of an IMSL subroutine for solving the nonlinear constrained least mean square equation formulated. The results of inversed parameters are compared to those obtained by the P1 method. The accuracy of the predictions of both two-flux and P1 methods are also discussed.

Introduction Inverse analysis of radiative transfer is concerned with the determination of radiative properties from a set of measured radiation quantities. Such a problem has applications in many fields, including, among others, remote sensing for the atmospheric or aerosol properties[i], temperature measurement in glass manufacturing[2], and experiments

585

586

J.-H. Tsai

for obtaining material properties[3,4].

Vol. 20, No. 4

The determination of the complete set of

radiative properties for a scattering medium is generally difficult because it is necessary to solve an inverse scattering problem for radiative transfer. A considerable amount of work has been reported for determining single scattering albedo and the scattering phase function from measurements of surface radiation intensities. Various methods have been utilized.

Saboonchi et al.[5] determined the gray radiative properties of insulating

materials from the

direct

comparison

of the

experimental measurements of

transmittance and reflectance with data established by the isotropic scattering model. Siewert[6,7] has formulated the inverse radiation problem for isotropically and anisotropically scattering, nonemitting media.

He used Chandrasekhar's X and Y

functions to solve for the two--term phase function and employed the FN method to determine the coefficients of a three--term phase function.

Siewert and Dunn[8]

extended the methodology to account for the nonuniform boundary radiance distributions. They used a Fourier transform approach and employed moments of the angular flux in determing the phase function coefficients. Karanjai and Karanjai[O,10] derived expressions for estimating up to three coefficients of the delta-Eddington phase function approximation from angular moments of intensity on the surfaces of the planar media. Ho and Ozisik[ll]presented a solution methodology for the inverse problem to determine the optical thickness and the single scattering albedo of two-layer inhomogeneous plane-parallel media. They adopted a nonlinear least squares approach to solve the inverse problem using detected radiation intensities in the forward and backward directions, which were obtained from the exact intensity expressions given by Cengel et al[12]. Dunn[13] assumed that the optical thickness was know, and tried to determin the single scattering albedo by

the inverse

Monte-Carlo method.

Subramaniam and Menguc[14] applied the same technique to determin the single scattering albedo and asymmetry factor variations within a nonhomogeneous clod of particles.

Wilson and Sen[15] adapted the moment method to determine the single

scattering albedo variation in an nonhomogeneous, semiinfinite parallel layer. Recently,

Vol. 20, No. 4

INVERSE SCATI'ERING P R O B L E M

WITH TWO FLUX METHODS

587

K~miuto and Seki[16] applied the P1 approximation and Tsai and Lin[17] applied the S-P(Sagan-Pollack) two flux method to determine the single scattering albedo and the asymmetry factor of the medium using the hemispherical transmittance data for slabs at various values of thickness. In this work, we apply the simple Sagan-Pollack (S-P), Schuster-Schwartzschild(S-S)

and

Coakley---Chylek(C---C) two

flux methods

to

determine the single scattering albedo and asymmetry factor. Results are compared to PI solutions.

We

consider a homogeneous,

azimuthaUy symmetric, absorbing,

anisotropic scattering, plane-parallel medium with free boundaries. It is exposed to diffuse incident radiattion at the boundary and no entrant radiation at the other. Emission from the medium is negligible. W e applied the S-P, S-S and C---C two flux methods[18] to determine the asymmetry factor and single scattering albedo directly from the hemispherical transmittance and reflectance. The inverse approach involves the use of an IMSL subroutine to solve a nonlinear least square equation constructed from the approximation methods for hemispherical transmittance and reflectance. Results indicate that the inverse C---C two flux solutions of the albedo and asymmetry factor are in better agreement with the exact than the S-P, S-S, PI approximation in most cases. When strongly forward scattering(g=0.85) and high albedo(ur=-0.9) the accuracy of S-P method is better than S-S, C---C and P1 methods.

Fundamental Expressions We

consider a homogeneous,

azimuthally symmetric, absorbing, anisotropically

scattering, plane---parallelmedium of optical thickness r. with free boundaries. It is exposed to diffuse incident radiattion at the boundary surfaces r=0. and no external incident at

r=ro.

Emission from the medium and the boundary surface r'=--rois

negligible. The equation of radiative transfer and the boundary conditions may be written as: OI r

Ir

(

wl

,/~)=~J'_lp(#,#

I(O,#)-- 1 and I ( r . , - # ) = 0

'

)I(r,/~ ' )d# '

O) (2)

588

J.-H. Tsai

Vol. 20, No. 4

where I(r,#) is the radiation intensity, r is the optical variable,# the cosine of the angle between the direction of the radiation intensity and the positive r-axis and w the single scattering albedo. Also, p(#,#') is the azimuthaUy-averaged phase function given by the Henyey--Greenstein phase function. p(#,#')= ~] (2n+l)gnPn(#)Pn(#'), n=0 where g is the asymmetry factor and Pn(#) the Legendre polynomial of order n.

(3)

The two flux models may be formulated by various methods. In general, two flux method arrived at by writing the radiative transfer equation for the forward and backward components seperately, and assuming that the forward intensity is constant over the upper hemisphere and independent of the angle #, and the backward component is a different constant over the lower hemispherical, also independent of #.

The

integro-differential radiative transfer equation would then lead to a set of two linear differential equations.

In this work, we consider three such methods.

respectively the S-S, S-P and C--C approximations.

They are

The differential equations

resulting from these approaches are identical in algebraic form, the only differences appearing in terms making up the constant numerical coefficients.

The equation (1)

appling two-flux method can be written in the form[18]: d dI l r a+~ =v- I

+ (~-)+~(1-B)I + (7-)+oJBI- (r)

(4a)

- d l v a ~ = - I - ( ~ - ) + o~BI+ (,)+w(1-B)I-(r)

(4b)

And the boundary conditions are written as: I+(r)=l/Tr and I-(7". )=0

(5)

where B is the back-scattered fraction and a is a numerical parameter. The definitions of these quantities, which may be different due to different approximations, are listed in Table 1. The B c c ( a ) and BSS are expressed the following formulas: fl0

Bcc(a)=(lp. ) ~: [(2n+l)/(n+l)](-1)ngnPn_l(0)Pn(a )

(Sa)

BSS=(1/2)[1+ ~

(6b)

n=0

n--1

[(2n+l)/(n+l)2](-1)ngnp2n_l(0 )

Vol. 20, No. 4

INVERSE SCATrERING PROBLEM WITH TWO FLUX METHODS

589

Eqns. (4) can be solved by number of standard techniques. We select here the operator approach. Put Eqns. (4) into the form.

{~

[l-a~l-S)]}I+(r)=wSI-(r)

{~ [ 1 - ~ ( 1 - B ) I } I - ( r ) - - - ~ ( 1 - B ) I + Solve equation(7b) for I-(r).

(7a) (r)

(7b)

Substitute into equation(7a) and expand the operator to

get

d2 I+(r) C~I+(r)--0, dr2 where

(8)

C2=~2-}, ~=[I-~(1-B)]/~, and~ B / ~ Equation(8)canbesolvedimmediatelyto give

I +(r)=~:e~p(C~)+I~exp(-C~),

(ga)

where ~ and ]~ are constants of integration. Put Equation (ga) into equation(7a) and solved for I-(r), we then obtain

I-(r)=~vVexp((r)+]~Vexp(-Cr),

(9b)

where

W=R+--~v and V=

~-C I/

•

Finally,apply the boundary conditionsand solvethe resultingexpressionfor ~ and ]~, we get the intensitysolutionsin the finalform

eC(r°-r)-v e-C(ro-r)]

i+(r)=~ w

weCr ° _ Ve_(ro

,

(10a)

and

1 [ eC( ~o-~ ) -e-C(~o-~)] I-(r)=Ir

.

.

Wet, to _ Ve-t~ro

W V.

(lOb)

Once the equation of transfer is solved, the hemispherical transmittance T and reflectance R based on the two-flux model are readily obtained as

T(ro)=l/~, and

(lla)

R( ro)=~S[sinh(~Iro/a)]/(#~7),

(11b)

where

590

J.-H. Tsai

Vol. 20, No. 4

~=¢(1-w)[1-~1-2B)J and ~--'-cosh(~/ro/a)+[1-~(1-B)][sinh(~/ro/a)]/~/. Similar results of the hemispherical transmittance and reflectance with the use of P1 method are given by Kamiuto and Seki[16] as

(12a) R(ro)=(-l+4w-3wg)sinh(~ro)/¢

(12b)

where ¢=(7-4w-3ug)sinh(~ro)+4~cosh(~-o) and

~=43(l-w)(l-~g)

Comparing P1 and two-flux methods, we note that both approximations yield the same expression of the hemispherical transmittance except for numerical factors.

Inverse Method It is seen from Eqs.(ll) and (12) that the right hand side of this equation involves the two parameters w and g. The main purpose of the present study is to evaluate w and g in Eqs.(11) and (12) from a set of hemispherical transmittance and reflectance. Whose values are obtained experimentally in actual situations. In this paper we use the Differential-Discrete--Ordinate(quadrature points=32) numerical value instead of the experiment data and use the inverse method to estimate the parameters w and g. T w o approaches are utihzed in the present study to estimate the w and g. (a) w and g are determined from 0nly transmittance to minize the squared error 7 over various sample thickness N 7= Z [Tn-Tn ]2 (13) 1 (b) w and g are determined from both transmittance and reflectance so as to minize the

following sum of squared error 7 over thickness N Z {[Tn-Tn]2+ [Rn-Knl2 }. 1 In above expressions, T n and 1~n

(14) are'measured'

transmittance

and

reflectance,respectively. While T and R n are 'simulated' transmittance and reflectance n

Vol. 20, No. 4

INVERSE SCATI'ERING PROBLEM WITH TWO FLUX METHODS

591

which are functions of unknown parameter of c# and g and are expressed by Eqs. (11) and (12). In Eqs. (13) and (14), N is the number of samples at various values of thickness. In the following calculations, we consider the optical thickness ranged from 1 to 10 and N equal to 10. The central problem is then in locating the minimum, say 7rain, of the g-c# surface. The conditions of this minimum provide the required inverse solution in each case. Expressions for T n and R n are replaced by P1 approximation(16), which are given in Eqs.(12) when P1 solutions are searched. It is apparent from the definitions of c# and g that c# takes a value between 0 and 1, while g lies between -1 and 1. The optimal values of c# and g must be determined within these regions when minimizing 7. To estimate the values of these parameters, a FORTRAN

SUBROUTINE

program B C O N F

in the IMSL

determine the optimal values of c# and g. The subroutine B C O N F minizes a function of N

package was used to gives solutions which

variables subject to bounds on the variables using a

quasi-Newton method.

Results and Discussions Numerical results are presented to illustrate the applications of the foregoing inverse analysis for determining

~o and g.

DDO(quadrature

points-32)

solutions of

hemispherical transmittance and reflectance are used to simulate the experimental measurements and parameters w and g obtained with the subroutine program BCONF. In the inverse computations, parameters involved in Eqs (1) and (2) were varied in the following regions -0.85
Table 2 shows the comparison of inverse

results obtained from Eq.(13) for P1, S - P and C-C approximations.

For these

methods, agreement between calculated and assumed experimental values of c# and g is mediocre only at high albedo or at strongly back scattering and is considerably poor at low values of albedo. Numerical computations with the use of transmittance only also indicate that the inverse solutions are highly dependent on the initial guesses. This is due to the nonlinear nature of the inverse analysis and different initial guess would

592

J.-H. Tsai

Vol. 20, No. 4

result in different local minimum. The results given in this table are ones whose squared error is minimum among those obtained with various initial guess.

The accuracy of

solutions can be improved with the use of both hemispherical transmittance and reflectance, that is, from Eq. (14), as shown in Table 3. Compared with Table 2, Table 3 shows the results which are much more accurate with the exact. In table 3, it also shows that the use of C--C two flux method(for a=2/3) to estimate the albedo and asymmetry factor gives more reasonable results over whole g-w domain than other approximation methods. The maximum error for all cases are below 8%. While for S-P two flux method is better only for strongly forward scattering(i.e., large asymmetry factor) and high albedo. scattering.

For P1 method is poor predicition for strongly forward

This problem is caused by the fact that calculation of hemispherical

transmittance and reflectance by the S-S two flux method yields reasonable results for most cases.

When strongly forward scattering(g=0.85) and high albedo(w=0.9) the

accuracy of S-P for transmittance and reflectance is more agreement exact value than C-C and S-S two flux.

While strongly forward scattering and low albedo the

reflectance for P1 method has even become negative and hence meaningless.

The

above-mentioned dependence of the solutions with the initial guess can be changed as well with the use of both hemispherical transmittance and reflectance for P1, and two flux methods.

With the use of both transmittance and reflectance for determing the

radiative parameters, results become spacewise independent of the initial guesses.

Conclusion The inverse two flux approximation is presented for determining single scattering albedo and asymmetry factor from the hemispherical transmittance and/or reflectance. The formulated nonlinear constrained least mean square equation solved by the use of the subroutine program BCONF in IMSL package. The S-S, S-P, C-C two flux and P 1 methods

yield

similar

analytical

expressions

form

for

the

transmittance and reflectance but with different numerical factors.

hemispherical Numerical

Vol. 20, No. 4

INVERSE SCAqTERING PROBLEM WITH TWO FLUX METHODS

Table 1. Summary of v a r i o u s t v o - f l u x approximations

I+(T) ]gI(r,#)#d#

I"(T) a C-C ]gI(T,-#)#d# 0
1

B BCC(a)" Bsp BSS

• Bcc(.):-½ilop(•* • **

BsP=½(l-g) 110 BSS=-~O~_Ip(#,# )d#'d#

Table 2.

Comparison of i n v e r s e p a r a m e t e r s obtained by the P1, C-C

and S-P two f l u x methods from Equation (13).

g

-0.85 P1

0.1

w=0.568

0

S-P

C-C

0,609

0,169

P1 0.615

0.85

S-P

C-C

0.655

0,328

P1 0.681

S-P

C-C

0.719

0.476

g=-1.00 -1.000 -1,000 -1.000 -1.000 -0,977 -1.000 -1.000 -0.975 0.3

w=0.614

0.653

0,322

0.725

0.738

0.560

0.821

0.793

0.761

g = - l . 0 0 -1.000 -0,984 -1.000 -0.915 -0.970 -0.820 -0.566 -0.858 0.5

~=0,690

0,728

0,494

0.779

0.746

0.706

0,855

0.832

0.806

g=-1.00 -1,000 -0.974 -0,767 -0,520 -0.817 -0.245 -0.057 -0.256 0.7

0.9

w=0,795

0.826

0.758

0.896

0,880

0,862

g = - I . 0 0 - 0 . 9 9 9 -0,983 -0,388 -0.179 -0.479

0,249

0.372

0,256

~=0.924

0,894

0.951

0.943

0.934

0,029 -0.139

0.669

0,726

0.647

0.917

0.682

0.887

0,819

0.920

~ = - 1 . 0 0 - 0 . 7 7 1 -0.983 -0.135

0.791

0.908

593

594

J.-H. Tsai

Table 3.

Vol. 20, No. 4

Comparison of inverse parameters obtained by the P1,

S-S, S-P and C-C two fluxs from Equation(14).

w

0.1 0.3 0.5 0.7 0.9

0.1 0.3 0.5 0.7 0.9

0.1 0.3 0.5 0.7 0.9

P1

S-S

S-P

C-C

w=0.321 g=-0.243 w=0.460 g=-0.544 ~=0.602 g=-0.710 w=0.780 g=-0.816 w=0.909 g=-0.888

g=-0.85 w=0.321 #=0.216 g=0.723 g=0.307 w=0.460 ~=0.377 g=O.180 g=-0.220 ~=0.602 w=0.540 g=-0.156 g=-0.437 w=0.750 w=0.711 g=-0.364 g=-0.557 w=0.909 w=0.895 g=-0.501 g=-0.632

w=0.095 g=-0.784 w=0.280 g=-0.870 #=0.469 g=-0.894 w=0.666 g=-0.905 w=0.878 g=-0.906

w=0.312 g=O.O01 w=0.437 g=-0.024 w=0.576 g=-0.037 w=0.730 g=-0.042 w=0.903 g=-0.042

g=O.O ~=0.312 w=0.205 g=0.900 g=0.654 w=0.437 w=0.350 g=0.681 g=0.356 w=0.576 w=0.510 g=0.5196 g=0.227 w=0.730 w=0.688 g=0.404 g=0.156 w=0.903 w=0.888 g=0.322 g=O.114

w=0.082 g=0.0064 ~=0.250 g=-0.043 w=0.434 g=-0.049 w=0.640 g=-0.048 w=0.871 g=-0.043

w=0.320 g=0.264 w=0.458 g=0.532 w=0.601 g=0.679 w=0.751 g=0.768 w=0.910 g=0.822

g=0.85 w=0.322 w=0.215 g=l.O00 g=0.965 w=0.458 w=0.374 g=0.991 g=0.923 w=0.601 w=0.540 g=0.980 g=0.903 w=0.751 w=0.712 g=0.970 g=0.889 w=0.910 w=0.896 g=0.960 ~=0.873

~=0.093 g=0.886 w=0.277 g=0.852 w=0.469 g=0.842 w=0.668 g=0.833 w=0.880 g=0.818

calculations show that with the use of transmittance only, the results are strongly dependent on the initial guess for these methods.

While with the use of both

transmittance and reflectance results become spacewise independent of the initial guess and are more accurate than that obtained with use of transmittance only. The results of

Vol. 20, No. 4

INVERSE SCATFERING PROBLEM WITH TWO FLUX METHODS

595

the albedo values and asymmetry factor with C - C two flux method for a=-2/3 are more reasonable results over whole g-• domain. But for the strongly forward scattering and high albedo case the S-P two flux has a excellent predictions for aJ,g.

Nomenclature B : Back-scattered fraction defined by Eqs.(6) g : asymmetry factor I : radiation intensity Pn : Legendre polynomial of degreen n p : phase function R : hemispherical reflectance T

: hemispherical transmittance

7 : quantity defined by equations(13) and (14) w : single scattering albedo 7- : optical thickness #

:cos0

a

: Numerical parameter

i.

References E.A. Ustinov, Cosmic Res. 15(5), 667(1977).

2.

R.Van Laethem, M. Leger and E.Plumat, J. Am. Ceram. Soc. 44,321(1961).

3.

J.A. Roux and A.M. Smith, AIAA 16th Thermophysics Conf., Paio Alto, California, 23(1981).

4.

J.F. Sacadura, G.Uny and A.Venet, Proc. 8th Heat Transfer Conf. 2, 565(1986).

5.

A. Saboonchi, W.H. Sutton and T.J. Love, J. Thermophysics and Heat Transfer 2, No. 2, 97(1988).

6.

C.E. Siewert, J. Appl. Math. Phys. 30, 517(1979).

7.

C.E. Siewert, J.Quant.Spectrosc.Radiat.Transfer 22, 441(1979).

8.

C.E. Siewert and W.L. Dunn, J. Math. Phys. 23, 1376(1982).

596

J.-H. Tsai

Vol. 20, No. 4

9.

S. Karanjai and M. Karanjai, Astrophys Space Sci. 117, 151(1985).

10.

S. Karanjai and M. Karanjai, Astrophys Space Sci. 117, 387(1985).

11.

C.H. Ho and M.N. Ozisik, J.Quant.Spectrosc.Radiat.Transfer 40, 553(1988).

12. Y.A. Cengel M.N. Ozisik and Y. Yener, Int. J. Heat Mass Transfer 27, 1919(1984). 13. W.L.Dunn, J.Comput. Phys. 41, 154(1981). 14.

S. Subramanian and M.P. Menguc, Int. Heat Mass Transfer 34, No.l, 253(1991).

15.

S.J.Wilson and K.K. Sen, Astrophys Space Sci. 119, 93(1986).

16.

K.Kamiuto and J.Seki, J.Quant.Spectrosc.Radiat.ransfer 37, 455(1987).

17. J.D. Lin and J.H. Tsal, A S M E Paper 91-WA-HT-13, presented at the 1991 A S M E Winter Annual Meeting December, 1991, Atlanta, Georgia. 18. J.J.,Buglia, Introduction to the theory of atmospheric Radiative Transfer, N A S A Reference Publication No 1156(1986). 19.

M.N. Ozisik,Radiative Transfer, John Wiley and Sons, Inc.,(1973). Received March 30, 1993