An inverse radiative transfer problem of simultaneously estimating profiles of temperature and radiative parameters from boundary intensity and temperature measurements

Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

www.elsevier.com/locate/jqsrt

An inverse radiative transfer problem of simultaneously estimating pro)les of temperature and radiative parameters from boundary intensity and temperature measurements Huai-Chun Zhou ∗ , Yu-Bo Hou, Dong-Lin Chen, Chu-Guang Zheng State Key Laboratory of Coal Combustion, School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan, 430074 Hubei, People’s Republic of China Received 3 October 2001; accepted 29 November 2001

Abstract A parallel-plane space )lled with absorbing, emitting, isotropically scattering, gray medium is studied in this paper. The boundary intensity and boundary temperature pro)les are calculated for the inverse analysis. For the simultaneous estimation of temperature, absorption and scattering coe
1. Introduction Inverse radiative transfer problems were classi)ed into two diDerent types, one is sensing (temperature-independent) applications, and the other is energy balance (temperature-dependent) ∗

Corresponding author. Tel.=fax: +86-27-8754-5526. E-mail address: [email protected] (H.-C. Zhou).

0022-4073/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 1 ) 0 0 2 7 4 - 6

606

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Nomenclature ai a1i; j ; a2i; j bi ci C1 ; C2 Dm ; Dc I1 ; I2 J1 ; J2 ka (x) ks (x) L N N0 N1 M Rd S T (x) T1 ; T2 x Y Lx L’ Lyi LV x y  # 1 ; 2  & w

coe
Subscripts 1; 2 L2 m 1 ; 2

refer to boundaries 1 and 2, respectively refers to L2 -type norm measured (containing errors) refer to wavelengths used in two-color method

Superscripts r iteration steps

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

607

applications [1]. For the temperature-independent applications, many research works have been carried out to estimate radiative properties. For example, a method was used to determine the eDective radiative property distributions in cylindrical systems [2], and a simple estimation method for the albedo was studied in [3]. An emerging-Oux )tting (EFF) method for inverse scattering problem was proposed to estimate the albedo and the Henyey–Greenstein phase function of a scattering medium [4], and a method of evaluating high-resolution absorption coe
608

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

which was based on the Monte Carlo method and was combined with a concept of angular factor eDective for image formation. In this paper, a parallel-plane space )lled with absorbing, emitting, isotropically scattering, gray medium is studied in this paper. The bounding surfaces are transparent and there is no externally incident radiation. The boundary intensity and boundary temperature pro)les are calculated for the simultaneous estimation of temperature, absorption and scattering coe
N −1


ai (x=L)i

(1000 K);

i=0

Fig. 1. The parallel-plane medium and its computation condition.

(1a)

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

ka (x) =

N −1


609

bi (x=L)i

(m−1 );

(1b)

ci (x=L)i

(m−1 ):

(1c)

i=0

ks (x) =

N −1
i=0

The span of x=L is [0; 1], and 1000 K is selected as a unit for T (x), which ensures ai , bi , ci approaching each other in the order of magnitude. As a 1-D radiation system, we only calculate the boundary intensity I1 and I2 exiting the boundaries 1 and 2, respectively, that varied symmetrically within the polar angle )eld # ∈ [ −#=2; #=2] within a )xed range of azimuthal angle [0; L’], where L’ is a small angle. Similarly, # ∈ [−#=2; #=2] is also divided into M discrete angles with L#=#=M . M can be large to, for example, 50 or 100, even 1000, which is possible when we use a camera to capture radiation images at the surfaces of the system. The forward problem is just described as below. Given L, T (x), ka (x), ks (x), N0 , N1 , M , L’, 1 and 2 , the boundary intensity vectors I1 and I2 and a kind of boundary temperature vectors T1 and T2 are to be obtained. At )rst, the boundary intensity of I1 and I2 can be calculated from [25,26] I1; i =

N0


Rd1 (j → i)4kaj Tj4 LVj =

j=1

I2; i =

N0


a1i; j Tj4 ;

(2a)

a2i; j Tj4 ;

(2b)

j=1

N0


Rd2 (j →

i)4kaj Tj4 LVj

j=1

=

N0
j=1

where a1i; j = Rd1 (j → i)4kai LVj , a2i; j = Rd2 (j → i)4kai LVj . 4kaj Tj4 LVj refers to the total energy emitted by the jth medium element with absorption coe
C1 1−5 ; eC2 =1 T − 1

S2 (x) =

C1 2−5 eC2 =2 T − 1

(3)

under the same conditions, we can obtain two sets of monochromatic radiative intensity at the two boundaries under the two wavelengths as I1 ;1;i =

N0
j=1

a1i; j (S1 ;j =);

I2 ;1;i =

N0
j=1

a1i; j (S2 ;j =);

(4a)

610

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

I1 ;2;i =

N0


a2i; j (S1 ;j =);

j=1

I2 ;2;i =

N0


a2i; j (S2 ;j =):

(4b)

j=1

Because of the assumption of the gray characteristics, the radiative transfer under monochromatic wavelength obeys the same rules valid for the whole-wavelength thermal radiation transfer, a1i; j and a2i; j , taking the same values used in Eq. (2). Then, two sets of boundary temperature, T1 and T2 , can be derived from Wein’s law of radiation [14,28] as     1 I1 ;1;i 15 1 ln ; (5a) − T1; i = C2 2 1 I2 ;1;i 25  T2; i = C2

1 1 − 2 1

   I1 ;2;i 15 ln : I2 ;2;i 25

(5b)

3. Estimation method Before the estimation method is discussed, the boundary intensity I1 and I2 and the boundary temperature T1 and T2 , all represented by Dc and calculated for the forward problem, are disturbed by measurement error with normal distribution and mean square deviation . Thus, we have Dm = Dc (1 + &);

(6)

where the range of random variables & is chosen as −2:576 ¡ & ¡ 2:576, which represents the 99% con)dence bound for the measured data. For the simultaneous estimation of T (x), ka (x), and ks (x) in the parallel-plane medium, the quantities to be estimated are just ai , bi , ci , which form a vector Y ∈ R(3×N )×1 , Y = [y1 ; : : : ; y3N −1 ]T = [a0 ; : : : ; aN −1 ; b0 ; : : : ; bN −1 ; c0 ; : : : ; cN −1 ]T . Assume that the temperature distribution and the radiative properties are Y(r) = [a0(r) ; : : : ; aN(r)−1 ; b0(r) ; : : : ; bN(r)−1 ; c0(r) ; : : : ; cN(r)−1 ]T after rth iteration. We need (r) (r) (r) (r) T to calculate the updating vector LY(r+1) = [La0(r) ; : : : ; La(r) N −1 ; Lb0 ; : : : ; LbN −1 ; Lc0 ; : : : ; LcN −1 ] for the temperature and radiative properties from the radiative intensity and temperature images Im and Tm from the two boundaries, so that we can obtain an updated vector Y(r+1) . From the )rst approximation,   Ly1   (r)   @I @I(r) .. ;  Im ≈ I(r+1) + (7) ··· .  @y1 @y3N −1  Ly3N −1   Ly1   (r)   @T @T(r) .. ;  (8) wTm ≈ wT(r+1) + w ··· .  @y1 @y3N −1  Ly3N −1

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

611

where w is a factor used to balance the magnitude diDerence between the boundary intensity and temperature. Under the least-square meaning, we have   T  −1  @I(r) @I(r) @I(r) @I(r) Ly1    @y1 · · · @y3N −1   @y1 · · · @y3N −1   .. =       .     @T(r) @T(r)   @T(r) @T(r)  w ···w w ···w Ly3N −1 @y1 @y3N −1 @y1 @y3N −1 

T @I(r) @I(r)

  @y1 · · · @y3N −1  Imeas − I(r+1)   ;  @T(r) @T(r)  wTmeas − wT(r+1) w ···w @y1 @y3N −1

(9)

and J1 = Im − I(r+1) L2 ;

J2 = Tm − T(r+1) L2

(10)

which are the residuals of boundary intensity and temperature, respectively. It is stated that the elements in the updating vector LY = (Ly1 ; : : : ; Ly3N −1 )T calculated from Eq. (9) are very large especially during the early stage of the iteration process. If such a vector LY is used to update the vector Y, then the iteration process will soon diverge. That means the updating magnitude provided by LY is too large to be applied. We adopted a simple method to constrain the updating magnitude. Calculate  as

3N −1  3N −1 

(r+1) = abs(Lyi ) abs(yi ) (11) i=1

then

i=1

 Lyi(r+1)

=

Lyi

if  6 0:5;

(0:5=)Lyi

if  ¿ 0:5:

(12)

This treatment can constrain the updating magnitude provided by LY, as well as maintaining its updating direction. In the middle and )nal stages of iteration,  will become ¡ 0:5, and even lower than 1.0e-6, which means it approaches the end of iteration. The updated vector can be simply calculated from Y(r+1) = Y(r) + LY(r+1) :

(13)

The only issue that remained to be resolved is to calculate  (r)   (r)  @I @T @I(r) @T(r) and ; ··· ··· @y1 @y3N −1 @y1 @y3N −1 the sensitivity functions, which form the Jacobi matrix in Eq. (9). Given a small increment y (it was 0.005 in the present study) for yi , while all the other elements in Y are not changed, repeat the forward calculation, Ii and Ti can be obtained. Then @I(r) ≈ (Ii − I(r) )=y; @yi

@T(r) ≈ (Ti − T(r) )=y: @yi

(14)

612

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

The residuals of boundary intensity and temperature, J1 and J2 expressed in Eq. (10), can be monitored during the updating process. Besides that, the mean square deviations of the temperature T , absorption coe
(1000 K);

ka (x) = 0:15

(m−1 );

ks (x) = 0:17

(m−1 ):

(15)

◦

The thickness was chosen as L = 1 m. N0 = 20; N1 = 200; M = 60; L’ = 0:5 . The two wavelengths for the calculation of boundary temperature using the two-color method were 1 = 700:0 nm and 2 = 546:1 nm. The error of  in Eq. (6) was =0:01. The residuals of boundary intensity and temperature J1 and J2 were shown in Fig. 2(a) and (b), respectively, with diDerent balancing factors, w. When w = 0:0, that means only the boundary intensity measurements are used to estimate uniform T; ka , and ks . J1 decreased during the initial 5 steps of iteration as shown in Fig. 2(a), but it soon increased and )nally stopped at a relatively high level, around 1.0e15, as compared with the other factors. When w = 1:0e2, that means the boundary

Fig. 2. Updating processes of residuals J1 (a) and J2 (b) along with the iterations with diDerent balance factors w for uniform pro)les of T; ka and ks under measurement error of  = 0:01.

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

613

Fig. 3. Updating processes of the uniform pro)les of T (a), ka (b) and ks (c) along with the iterations with diDerent balance factors w under measurement error of  = 0:01.

temperature begins to participate in the iteration process, J1 could sharply drop to a lower level around 1.0e11, and then oscillate obviously. As w increases to 1.0e3, 1.0e4, and 1.0e5, J1 could also sharply drop to the level around 1.0e11, and fortunately the oscillation reduced largely. Similar results were obtained for J2 as shown in Fig. 2(b). Fig. 3(a)–(c) displayed the updating process of T; ka , and ks , respectively. The results obtained as w = 0:0 are also unacceptable, but the results obtained with other non-zero w are more reasonable. This simple example indicates that the boundary intensity measurement alone is not enough to estimate simultaneously the temperature (source) and the radiative properties (both absorption and scattering coe
614

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Furthermore, with suitable balance factor w around 1.0e4, it is obvious that the oscillation of scattering coe
(1000 K);

(16a)

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Fig. 4. Estimation results of 16 times of calculations for measurement errors of  = 0:01; 0:02; 0:03, respectively. (a) Temperature vs. absorption coe
615

Fig. 5. Estimation results of 16 times of calculations for the geometric thickness of 1.0, 4.0, and 8:0 m with measurement error of  =0:01. (a) Temperature vs. absorption coe
ka (x) = 0:24 + 0:18(x=L) − 0:20(x=L)2

(m−1 );

(16b)

ks (x) = 0:21 + 0:22(x=L) − 0:25(x=L)2

(m−1 ):

(16c)

In order to observe the inOuence of optical thickness on the estimation, three diDerent values of geometrical thickness, L = 1; 4 and 8 m, are adopted in the simulation study, which correspond to optical thicknesses of 0.5, 2.0 and 4.0, respectively. The boundary intensity and temperature are shown in Fig. 6(a) and (b) with diDerent optical thicknesses. When the optical thickness becomes large, the emittance of the medium will also increased due to the larger volume of the medium with the same radiative property pro)les. But under this circumstance, the high-temperature zone will also shift from the boundaries, so that the boundary intensity and temperature shown in Fig. 6(a)

616

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Fig. 6. The boundary intensities (a) and temperatures (b) at the two boundaries of the medium with diDerent optical thicknesses of 0.5, 1.0 and 2.0 when the three parameters T; ka and ks have parabolic pro)les inside the medium.

and (b) does not vary very much. Contrary to the boundary intensity, the boundary temperatures display even more pro)les. This is because the temperature was calculated from the two-color method by Eq. (5) in which the ratio of the monochromatic intensity under two wavelengths is dominated by the emittance from the high-temperature zone of the medium. The estimation results are shown in Figs. 7, 8, and 9 for T; ka , and ks , respectively, with a measurement error of  = 0:001. The calculation for the diDerent optical thicknesses is repeated 10 times with diDerent random value series, and the estimated results for the three parameters are derived from the average of the last 15 times of calculation over the total of 30 times of calculation. In these )gures the shadow areas are just covered by all the results from the 10 times of calculation. For the temperature pro)les shown in Fig. 7, when the optical thickness increases, the estimation does not change more obviously than those observed in Figs. 8 and 9 for absorption and scattering coe
H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Fig. 7. The estimation results for T with measurement error of  = 0:001 when the three parameters have parabolic pro)les inside the medium. The optical thickness is 0.5 (a), 1.0 (b) and 2.0 (c), respectively.

617

Fig. 8. The estimation results for ka with measurement error of  = 0:001 when the three parameters have parabolic pro)les inside the medium. The optical thickness is 0.5 (a), 1.0 (b) and 2.0 (c), respectively.

618

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

Fig. 9. The estimation results for ks with measurement error of  =0:001 when the three parameters have parabolic pro)les inside the medium. The optical thickness is 0.5 (a), 1.0 (b) and 2.0 (c), respectively.

the least-squares method. To avoid over-updating for the quantities, the relative updating magnitude during the iteration process is constrained not to be ¿ 0:5. In order to compare the eDects with and without boundary temperature measurement on the simultaneous estimation, uniform pro)les of temperature, absorption and scattering coe
H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

619

than measurement errors. More interestingly, once the absorption coe
620

H.-C. Zhou et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 74 (2002) 605 – 620

[14] Zhou HC, Yuan P, Sheng F, Zheng CG. Simultaneous estimation of the pro)les of the temperature and the scattering albedo in an absorbing, emitting, and isotropically scattering medium by inverse analysis. Int J Heat Mass Transfer 2000;43(23):4361–4. [15] Zhou HC, Han SD, Zheng CG. Simulation on simultaneous reconstruction of temperature distribution, absorptivity of wall surface and absorption coe