An explicit formulation for the inverse transport problem using only external detectors: Part II: Application to one-dimensional homogeneous and heterogeneous participating media

An explicit formulation for the inverse transport problem using only external detectors - Part II: Application to one-dimensional homogeneous and heterogeneous participating media

Nancy I. Alvarez AcevedoI,^*; Nilson Costa Roberty^II; Antônio J. Silva Neto^I

^IDepartment of Mechanical Engineering and Energy, Instituto Politécnico, IPRJ Universidade do Estado do Rio de Janeiro, UERJ P.O. Box 97282, 28601-970 Nova Friburgo, RJ, Brazil. E-mail: [e-mail], [e-mail]

^IINuclear Engineering Program, PEN/COPPE, Universidade Federal do Rio de Janeiro, UFRJ, P.O. Box 68509, 21945-970 Rio de Janeiro, RJ, Brazil. E-mail: [e-mail]

ABSTRACT

The steady state inverse radiative transfer problem in one-dimensional participating media is studied as an example of application of the new methodology presented in a accompanying paper by the authors [2]. Spectral methods are used for the appropriate analysis of the direct transport problem. For the inverse problem, we present a matrix that involves only values of the flux intensities at the boundary of the medium. Its columns are built with a set of linearly independent solutions for the system, which is formed when angular half-range prescribed boundary values and the complementary measured half-range boundary values are coupled. The final inverse albedo problem is treated as a full range two point boundary value problem. Test cases results are presented.

Mathematical subject classification: 80A23, 15A29.

Key words: inverse problems, explicit formulation, radiative properties, spectral methods, semi-group theory.

1 Introduction

An important inverse problem with many relevant applications that has been attracting the attention of a growing number of researchers is the one related to the determination of radiative properties in participating media such as the absorption and scattering coefficients. These properties may depend on a number variables: time, position, angle, and energy or wavelength.

This type of inverse problem is applied in non-destructive testing of materials in the field of engineering [7], in biomedicine including the diagnosis of biological tissues and physiologic measurements, in the pharmaceutical industry for on-line measurements of active substance concentration in pharmaceutical formulations, in the food industry to test the quality of the products [1], in remote sensing of marine parameters as chlorophyll concentrations, backscattering coefficient and aerosol parameters such as optical depth [3, 11, 12] or atmospheric pollution research [4], among many other fields of interest.

Inverse radiative transfer problems can be formulated either explicitly or implicitly [10]. We developed an explicit formulation for the inverse radiative transfer problem for estimating the total extinction and scattering coefficients in an absorbing and anisotropic scattering medium only using external detectors, which is presented in an accompanying article [2]. This formulation uses concepts of the discrete ordinates method and matrix operators for modelling the direct problem. The inverse problem is formulated using the relation between the albedo and the collision operator.

In this paper this new formulation is applied to one-dimensional homogeneous and heterogeneous participating media in steady state. The finite dimensional version of the problem is obtained with the discretization of the angular domain using the concepts of the discrete ordinates method, and the expansion of the phase function of anisotropic scattering in a series of Legendre polynomials.

Two matrices are constructed, one for each boundary of the one dimensional domain. Such matrices are built using a proper arrangement of the boundary conditions imposed to the radiative transfer problem, and the experimental data on the radiation that leaves the medium under analysis.

A computational routine was developed and implemented in order to show the feasibility of the formulation, and test cases results are presented.

2 Mathematical formulation of the one dimensional discrete problem

The Linear Boltzmann Equation for the steady state situation, with no spectral dependency [5], usually applied in the modelling of the interaction of radiation with a participating medium is

where Φ ∈ L²(S ×ω) is the constant velocity radiation intensity, S^N-1 × Ω is the domain of Φ and (x, ω) ∈ Ω × S^N-1, S^N-1≅ [0, 2π) when N = 2; σ_t∈ L^∞(Ω) is the total extinction coefficient (absorption + out scattering); σ_s(x, ω'· ω) is the scattering coefficient; and ω'· ω = cosθ₀, with θ₀ the angle between the direction of incident radiation ω' and the emergent scatteredradiation ω,

σ_s(·, s) ∈ L^∞(Ω), a.e. s ∈ [-1,1]; σ_s(x, ·) ∈ L¹([-1, 1]), a.e. x ∈ Ω.

The one dimensional steady state problem for equation (1) corresponds to the slab problem in the transport theory, and in this situation we expect to reconstruct only spatially homogeneous coefficients considering only experimental data acquired with external detectors. This is, therefore, an important test for the formalism we proposed in the accompanying paper [2].

The one dimensional version of equation (1) is written as

in which we identify

·∇ = ·_x = µ, that is, Φ (ω, x) ≡ Φ(µ, x);

[A₁Φ](µ, x) = σ_t(x) Φ(µ, x) is the extinction operator; and

[A₂Φ](µ, x) = [₌₀σ_sl(x)p_l(µ) p_l(µ') Φ(µ', x)] dµ' is the scattering operator. Here {p_l(µ); l = 0, ..., L} are the normalized Legendre polynomials up to order L.

We introduce here the Hilbert spaces L²((-1, +1) × (0, τ₁)) of square integrable functions in the phase space (-1, +1)×(0, τ₁). Note that τ₁ = τ(1) means the thickness in the direction µ = 1 and that µ = 0 is a singular value for equation (2).

The direct problem for equation (2) has half range boundary value conditions in the two point boundary {0, τ₁}

The parameters to be reconstructed in this case are the vectors of scattering kernel coefficients σ_s = {σ_sl> 0, l = 0, ..., L, ∈ ^L+1} and σ_t∈ ¹. We may suppose σ_sl; l = 0, ..., L, to be real, since these has been built from a truncated Legendre polynomials series of the scattering kernel.

Since we are interested in an absorbing non multiplicative media, we also have σ_t > σ_s0 > σ_s1 > σ_s(l+1); l = 1, ..., L. For a better description of this particular subject see Ref. [5].

The multiplicative operator

is bounded, but have unbounded inverse due to its behavior near µ = 0. Its spectra is continuous in the entire interval [-1, +1] by obvious reasons. It is also self-adjoint.

The operator A₁, here named extinction operator,

is bounded and invertible since σ_t > 0. Also for the two variables µ and x it behaves like the identity operator. It is self-adjoint.

The operator A₂, here named scattering operator

where

[A₂Φ] = -σ_slp_l(µ)p_l(µ')Φ(µ', x) dµ',

is finite rank, and therefore compact. Its set of eigenvalues and eigenfunctions is {(σ_sl, p_l(µ)), for l = 0, ..., L}.

The sum of operators A₁ and A₂ is the collision operator A. This operator

where

[AΦ](µ, x) = σ_tΦ(µ, x) + σ_slp_l(µ)p_l(µ')Φ(ω', x) dµ'

is bounded, strictly positive, self-adjoint and has inverse

where [A^-1f](µ, x) = (1/σ_t)f(µ, x)+[K₁f](µ, x), and

[K₁f](µ, x) =σ_sl/(σ_t - σ_sl)p_l(µ)p_l(µ')Φ(µ', x) dµ'

3 Spectral analysis of the one dimensional problem

In a formal operator notation the equation (2) can be written as

and

or, by defining the current ψ = µΦ

In the two situations we obtain an operator, that for Eq. (11) is defined as

which is the multiplication of two self-adjoint operators A^-1 and T.

They do not commute, so, B is not expected to be itself a self-adjoint operator. But it is not difficult to see that B is similar to a self-adjoint operator, that is

B = A^-1/2 (A^-1/2T A^-1/2)A^1/2

and that the characteristic function in the interval [-1,+1] is a cyclic vector for B, see [5]. So, the spectra of B is simple and there exists a Hilbert space which is equivalent to L²(S^N-1×Ω) and has the appropriate internal product with respect to which the operator B is self-adjoint.

Under this situation there exists a unitary transformation F that diagonalizes A by mapping the original Hilbert space (in which the Legendre polynomials are), in another Hilbert space with a Borel measure (in which the Chandrasekhar polynomials are).

In this transformed space the operator of multiplication by the eigenvalues (the union of a continuous part of spectra that comes from T and a discrete part that comes from operator A) plays the role of the operator B. The semigroup generated by the operator B is holomorphic. For details see [5, 6].

Formally, for

where T_λg = {λg | λ eigenvalue of T^-1A}.

The solution of equation (10) becomes

where

U(x - x₀) = exp [ - (x - x₀)T^-1A]

is the representation of the semigroup generated by T^-1A.

We will now use a set of fundamental solutions of the transport equation (2) and equations (3) and (4) in order to establish the connection between the transport problem with scattering and the equation (15).

4 Chandrasekhar discrete ordinates approximation

This analysis can be done with a distribution approach, but we preferred to adopt the discrete ordinates approximation to work in a finite dimensional phase space and with a matrix operator.

So, let {w_i; i = 1, ..., M} such that w_i = w_M+1-i, i = 1, ..., M/2; M even, are the appropriate weights to be used in the discrete problem. In this case the radiation intensity (or particle flux) Φ is substituted by the vector-valued function,

Φ^T = [Φ₁, ..., Φ_M]

with Φ_i∈ H¹(0, τ₁), i = 1, ..., M.

It is important, due to range related questions, to establish an order in a suchway that the first M/2 components correspond to the finite sequence of quadrature weights {w_i} ∈ (-1, 0). We call this first part of Φ by Φ^-. The remaining sequence {w_i} ∈ (0,+1) will be called Φ⁺ and the vectorial flux will be written

This decomposition is consistent with the half range two point boundary value prescribed in the slab problem and also with the outflux measurement, see Figure (1).

If we denote W the diagonal matrix composed with the quadrature weights, equations (2) and (10) can be written as

where we are using the fact that A is symmetric and commutes with W.

The matrix form of the direct problem is: To find

Φ_i∈ H¹(0, τ₁), i = 1, ..., M;

such that equations (21), (17) and (18) hold.

In order to solve this problem we determine the discrete spectra of the matrix T^-1AW and use the spectral representation of the semi group solution, equation (15), to propagate the two point boundary values throughout the other slab points.

Note that here we have to partition the matrices in order to be consistent with the fact that we only know half range of the angular variable in each boundary position 0 and τ₁.

Note also that the discrete calculated numerical spectra is adherent to the compact union of the continuous and the discrete spectra and that we must be aware of the singularity introduced in the problem by µ = 0. This singularity is avoided with the way the quadrature points are calculated, but this natural regularization is not sufficient when the number of points grows.

Since we are dealing with exponential functions, for very large M will be necessary to filter the smaller values that appear when operator A^-1T is being calculated. The numerical implementation and test cases results are presented in the next section.

In this simplified situation we can derive an explicit formula for the albedo operator, Λ. For this we use the exponential formula (15) and the factorization equation (14) to deduce an explicit expression for the transmission problem.

Let {Φ_i(x), i = 1, ..., M; 0 < x < τ₁} a set of M linearly independent solutions. Then

satisfy

were Δ is a diagonal matrix containing the eigenvalues of [T^-1A].

We now adopt the two range decomposition

decompose the matrix Φ(x) in four parts, Φ_mm(x), Φ_mp(x), Φ_pm(x) and Φ_pp(x) as shown

and prescribe partially the flux in the two surfaces x = 0 and x = τ₁ in the following way

and

where Φ_pm(0) = I and Φ_mp(τ₁) = I are the prescribed intensities for the external sources, Φ_pp(0) = 0 and Φ_mm(τ₁) = 0 are the respective boundary conditions of zero flux in the other side of the domain; and Φ_mm(0), Φ_mp(0), Φ_pm(τ₁), and Φ_pp(τ₁) are the unknown emergent fluxes (or radiation intensities) to be calculated with the albedo operator, or if this is not directly known, measured experimentally.

The transmission operator is

Introducing Eqs. (25) and (26) into Eq. (15) we write

where Tr can also be partitioned in combinations of the two ranges as

Using the definition of the albedo operator one obtains

where

The explicit formula of the albedo that maps the half range influx in the positions {0, τ₁} to the half range outflux in the same positions is then derived.

As we have previously derived, the albedo operator establishes an explicit relation between the influx and the outflux in both sides of the slab. Since in the inverse problem we cannot suppose to know the parameters in the transmission operator used to compose the albedo, we assume the knowledge of some information about its graph. In the ideal case we can suppose that we know completely the albedo's graph. In a real situation, otherwise, we can only suppose we know a linearly independent set of solutions Φ to the problem, derived from some finite set of experiments.

When we have M experiments to generate the columns of the matrix flux in the two sides x = 0 and x = τ₁ of the slab, that is, to form the matrix equation

Φ(τ₁) = exp [ - τ₁T^-1A]Φ(0)

we obtain the explicit formula for A

and inside A we find information about the absorption and the scattering coefficients. Note the resemblance that this equation has with the equation of transmission tomography [8].

In the case of multilayered slab, with transparent contiguous interfaces at positions x = τ₁ < τ₂ < ... < τ_N-1 and boundary surfaces at x = 0 and x = τ_N, we have

Φ(τ_n) = exp [ - τ_{n -1} T^-1A]Φ(τ_{n -1}) ; n = 1, ..., N

and

If inside each layer the operator function x A(x) is continuous in the uniform operator topology and the families {A_n(x), x ∈ [τ_n-1, τ_n]} are stable [9], the formalism for integrals in the measure for x will be preserved,

The quantity inside the integral is the matrix optical thickness of the slab, and the logarithmic expression

is finally obtained.

The question here is about how much information can be extracted from the knowledge of the finite dimensional approximation of the albedo, and consequently in the operator A. Note that in general the number of parameters in A is less than M(M+1)/2, and since in the slab problem the scattering kernel has the symmetry σ_s(x, ω'· ω) = σ_s(x, cos(θ - θ')) which presents a convolutional property, the maximum number in the discrete ordinates method of order M will be L+1, where L is the number of Legendre polynomials in the kernel expansion. Obviously, M > L+1 and we have an excess of M(M+1)/2-(L+1) degrees of freedom in the data matrix -log[Φ(τ_n)Φ^-1(0)].

Apart from the sensitivity of the data on the coefficients, we expect to have information to treat a number of layers at least equal the integer quotient of (M(M + 1)/(2(L + 1)))-1. Moreover, noise and the extension of data files will probably reduce even more the number of layers that can be estimated.

Also we note that the unitary transformation F is dependent on the magnitude of the scattering and absorption coefficients, and will highly increase the computational cost of non homogeneous problems. Finally, equation (34) is also appropriated for the transient problem but this question will not be discussed in this work.

5 Numerical results

A computer program has been developed and implemented in order to test the ideas that we presented here. In this program the quadrature set and the Legendre polynomials are calculated in real time. When this quadrature set includes very small µ values it is necessary to make a spectral regularization by cutting off the angular directions with small absolute cosine values µ.

In order to establish a very simple reconstruction model we note that by using the Fourier-Legendre series for the Dirichlet kernel approximation in the first part of the operator A we obtain

[AΦ](µ, x) ≅(σ_t(x) - σ_sl(x))p_l(µ)p_l(µ') Φ(ω',x) dx

which is good for a sufficient large L.

In Figure 2 we show a solution for the direct problem with M = 14 discrete ordinate points.

The slab is illuminated in its side x = τ₁ for all ordinates. The values used for properties of this media were: σ_t = .5; σ_sl = [.4, .3, .2, .1, .0]; and τ₁ = 1/2.

Next we present the results obtained for a inverse problem considering a set of M = 14 experiments and a slab with properties σ_t = .5; σ_sl = [.4, .3, .2, .05, .01, .00]; and τ₁ = 1/2. Figure 3 presents the fluxes matrices whose columns are the incident flux and the emergent flux in each one of M experiments, which are respectively imposed and measured for the complete set of these, as explained in section 4, i.e. Eqs. (25) and (26), respectively. Here, "FLUX MATRIX (A)" corresponding to Φ(0) and "FLUX MATRIX (B)" to Φ(τ₁).

In order to minimize the occurrence of negative fluxes we use a flux of uniform intensity as basis for each incident flux, i.e. the incident flux total has values non negative in all directions, but with a higher value in only one direction.

Figure 4 presents the eigenvalues of these fourteen problems, the exact and the reconstructed values, respectively, for various gaussian noise levels controlled by the variance. We note the good agreement between the exact and the constructed values when the noise level is small.

Figure 5 shows the reconstructed scattering coefficients for several gaussian noise levels and we easily deduce from the results the high sensitivity of the problem with this kind of noise, mainly for the higher order coefficients which has already been reported by Silva Neto and Özisik [11] using an implicit formulation for the inverse problem.

In Figure 6 we show the error in the reconstructed extinction coefficients.

6 Conclusions

In the present work and in the accompanying paper [2] we derive and show the feasibility of a new explicit formulation for the inverse radiative transfer problem. From the knowledge of the imposed incident and measured emergent radiation intensities in a one dimensional homogeneous slab we were able to estimate the anisotropic scattering coefficients, using only external detectors.

We have used both noiseless and noisy synthetic data in our numerical simulations. It was observed that the higher order coefficients are more sensitive to the noise present in the experimental data.

We propose the application of this development for multilayer media and also the possibility of its use in the solution of the inverse problem of radiative transfer in the transient state.

Acknowledgements. The authors acknowledge the financial support provided by CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico, and FAPERJ, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro.

Received: 10/III/09.

Accepted: 31/III/10.

• [1] Ch. Abrahamsson, T. Svensson, S. Svanger and S. Andersson-Engels, Time and wavelength resolved spectroscopy of turbid media using light continuum generated in a crystal fiber Opt. Express, 12(17) (2004), 4102-4112.
• [2] N.I. Alvarez Acevedo, N.C. Roberty and A.J. Silva Neto, An explicit formulation for the inverse transport problem using only external detectors. Part I: Computational Modelling Computational & Applied Mathematics (2010).
• [3] B. Hamre, J.J. Stamnes, W. Li, K. Stamnes and R. Spurr, Fast and accurate retrieval of aerosol and marine parameters for coastal waters Ocean Optics 19^\textth Conference (abstract), Tuscany, Italy (2008).
• [4] M.-J. Jeong and N.C. Hsu, Retrieval of aerosol single-scattering albedo and effective aerosol layer height for biomass-burning smoke: Sinergy derived from "A-Train" sensors Geophys. Res. Lett., 35 (2008), L24801, doi:10.1029/2008GLL036279.
• [5] H.G. Kaper, C.G. Lekkerkjer and J. Hejtmanek, Spectral methods in linear transport theory Operator theory: Advances and applications, Birkhäuser, Verlag, Basel, 5 (1982).
• [6] T. Kato, Perturbation theory of linear operators Classics in Mathematics, Springer-Verlag, Berlin (1995).
• [7] M. Koseki, Sh. Hashimoto, Sh. Sato and H. Kimura, CT Image reconstruction algorithm to reduce metal artifact Journal of Solid Mechanics and Materials Engineering, online ISSN: 1880-9871, 2(3) (2008), 374-383.
• [8] F. Natterer, Algorithms in tomography The State of Art in Numerical Analysis, Clarendon Press, New York, 63 (1997).
• [9] A. Pazy, Semigroups of linear operators and applications to partial differential equations Applied Mathematical Science, Springer-Verlag, New York, 44 (1983).
• [10] A.J. Silva Neto, Explicit and implicit formulations for inverse radiative transfer problems 5^\textth World Congress on Computational Mechanics, Mini-Symposium MS125-Computational Treatment of Inverse Problems in Mechanics, Vienna, Austria (2002).
• [11] A.J. Silva Neto and M.N. Özisik, An inverse problem of simultaneous estimation of radiation phase function, albedo and optical thickness J. Quant. Spectrosc. Radiat. Transfer, 53(4) (1995), 397-409.
• [12] K. Stamnes, H. Eide, W. Li, R. Spurr and J.J. Stamnes, Simultaneous retriefal of aerosols and marine constituents in coastal waters Gayana online ISSN 0717-6538: Pan Ocean Remote Sensing Conference 2004, Concepción, Chile. Proc., Univ. de Concepción, Chile, 68(2) (2004), 543-551.

*

The first author is supported by FAPERJ.

Publication Dates

• Publication in this collection
  22 Nov 2010
• Date of issue
  2010

History

• Received
  10 Mar 2009
• Accepted
  31 Mar 2010