Tomography with few data: Use of collocation methods in acoustic pyrometry

Helmut Sielschott^{1} University of Münster Willy Derichs^{2} RWE Energie AG

Published in: ECMI ’94: Student Proceedings Kaiserslautern 1996 pages 251–263

1Institut für Numerische und instrumentelle Mathematik, Einsteinstr. 62, D-48149 Münster. E-mail: helmut.sielschott@uni-muenster.de

2Kraftwerk Niederaußem, Abt. KF-F, Postfach 1461, D-50104 Bergheim

Tomography with few data: Use of collocation methods in acoustic pyrometry

Helmut Sielschott, University of Münster^{∗}

Willy Derichs, RWE Energie AG^{∗∗}

1 Introduction

In mathematics, “tomography” denotes the solution of the so-called “Radon Inversion Problem”, which means the reconstruction of a function from its integrals along lines or planes.

The most famous application is computerized tomography. In this medical use, you always have more then 10^{5} measurements, i.e. lines for which the integral of the function to be reconstructed is measured. In case of acoustic pyrometry the problem is different. The number of data is only 24!

2 Measurement technique

Acoustic pyrometry is a technique for measuring the temperature of a gas. In this case, we are interested in the temperature distribution in one horizontal plane in the combustion chamber of a brown coal fired power station.

Eight transceivers, each consisting of a loudspeaker and a microphon, are distributed at the walls of the combustion chamber in the way shown in figure 1.1. During the measurement one transceiver after the other works as a transmitter, while the others are receiving. The time of flight of the acoustic signal is measured to compute the average speed of sound along the 24 lines between the transceivers.

∗Institut für Numerische und instrumentelle Mathematik, Einsteinstr. 62, D-48149 Münster. e-mail: helmut.sielschott@uni-muenster.de. This work is part of a master thesis supervised by Prof. Dr. F. Natterer.

∗∗Kraftwerk Niederaußem, Abt. KF-F, Postfach 1461, D-50104 Bergheim Fig. 1.1 Scanning geometry with eight transceivers

Knowing the average speed of sound (c), we can calculate a good approximation to the average temperature (T) for each line using Laplace’ equation:

T = \frac{M}{κR}c^{2}.

M denotes the molecular weight, κ = c_{p}/c_{v} is the quotient of specific heats and R is just the universal gas constant. There are some inaccuracies in this way of modeling the problem:

• The gas in the chamber is not an ideal gas, but because of high temperature and low pressure this is a good assumption.

• M and κ are not known exactly.

• The relation of T and c is not linear, but as the temperatures we deal with lie in an interval much higher than zero Kelvin we get a good approximation to the mean temperature on each line. The appropriate way is to reconstruct the values of c(x) first and make use of the Laplace equation afterwards. Tests with realistic examples have shown that the difference in the final reconstructions is less then 3^{o}C and therefore much smaller then the expected accuracy.

• Because of bending phenomena the acustical signal does not travel along the straight line. This is seen as a source of a systematic error at the moment and will be taken into account by the algorithm in future.

• The gas in the combustion chamber is rotating and rising, so the flame is going spiralwise to the top of the furnace. The horizontal part of the movement of gas makes the speed of sound direction-dependend. This is essentially equalized by measuring the speed of sound in both directions for each path and taking the mean value as “true” velocity.

The measurement errors add up to less then one percent as shown in [2]. In addition the error caused by high frequency events in the furnace leads to an overall error of 2.5%.

Two ways for preparing data are implemented: It is possible to take the mean value for each line from 10 successive measurements or for each line and each direction a spline can be computed modeling the value for this specific measurement. Using these splines a reconstruction can be calculated for any given point of time. This is utilised to present the situation in the furnace in a film consisting of consecutive reconstructions.

The first way of exploiting the data leads to an additional error as slow changes in the furnace cannot be modelled appropriate. In this case the overall error is about 4%.

3 Collocation methods

Let f : Ω ⟶ IR denote the temperature at each point of the reconstruction area Ω. As we know the length of the lines, the values

g_{l} := ∫_{L_{l}} f(x) ds, l = 1, …, 24 (1.1)

can be computed for the 24 Lines L_{l} from the measured data.

The simple idea of a collocation method is to build up the reconstruction in the form (n ∈ IN)

f(x) = ∑_{k=1}^{n} ξ_{k}φ_{k}(x) (1.2)

We chose φ_{k}(x) = e^{−λ|x−x_{k}|}. The reconstruction region is devided into 10 × 10 pixels and the x_{k} are the midpoints of these pixels. So we have n = 100 and a clearly underdetermined problem. To get a system of equations for the coefficients ξ_{k}, we combine (1.1) and (1.2):

⇝ g_{l} = ∫_{L_{l}} ∑_{k=1}^{n} ξ_{k}φ_{k}(x) ds

= ∑_{k=1}^{n} ξ_{k} \underbrace{∫_{L_{l}} φ_{k}(x) ds}_{=:a_{lk}} , l = 1, …, 24

⇝ g_{l} = ∑_{k=1}^{n} a_{lk}ξ_{k}  , l = 1, …, 24

With the definitions

m := 24
A := (a_{lk})_{{l=1,...,m \atop k=1,...,n}} ∈ IR^{m×n}
g := (g_{1}, ..., g_{m})^{t} ∈ IR^{m}
ξ := (ξ_{1}, ..., ξ_{n})^{t} ∈ IR^{n}

we arrive at the problem to find a vector ξ which solves

Aξ = g . (1.3)

The matrix A is ill-conditioned, so we are solving this equation by means of the Singular Value Decomposition (SVD) and digital filtering (Tikhonov– Phillips method).

The idea of the method of Tikhonov–Phillips is to minimize

∥Aξ − g∥^{2} + γ∥ξ∥^{2} (1.4)

The real number γ serves as a parameter. Using the SVD of A we can compute this minimizer directly. The proof for the following theorem can be found in [3]. Theorem 1: Let A ∈ IR^{m×n}. Then there are orthogonal matrices U = (u_{1}, ..., u_{m}) ∈ IR^{m×m} and V = (v_{1}, ..., v_{n}) ∈ IR^{n×n} with:

A = UΣV^{T}

and Σ_{ik} = \{ \begin{matrix}σ_{k} & , i = k \\ 0 & , otherwise \end{matrix} , Σ ∈ IR^{m×n}, σ_{1} ≥ σ_{2} ≥ ... ≥ σ_{min(m,n)} ≥ 0.

With this decomposition, the Moore–Penrose generalized solution to Aξ = g can be written in the following form. Let p′ denote the number of σ_{k} ≠ 0:

ξ_{MP} = ∑_{k=1}^{p′} \frac{1}{σ_{k}} (g, u_{k})v_{k}

The symbol ”(·, ·)” denotes the standard inner product between vectors in IR^{m}. The terms belonging to a high index k (and to a small number σ_{k}) are more sensitive to measurement errors because of the factor \frac{1}{σ_{k}}. To increase stability these terms are weighted using a filter F_{γ}:

ξ = ∑_{k=1}^{p′} F_{γ}(σ_{k})\frac{1}{σ_{k}} (g, u_{k})v_{k}

The method of Tikhonov–Phillips leads to the filter F_{γ}(σ) = \frac{σ^{2}}{σ^{2}+γ} (see [5]), so

ξ_{T P} = ∑_{k=1}^{p′} \frac{σ_{k}^{2}}{σ_{k}^{2} + γ} \frac{1}{σ_{k}} (g, u_{k})v_{k}

is a minimizer of (1.4).

4 Singular Value Decomposition of the discrete Radon Transform

To get a deeper understanding of the described method we now examine the singular value decomposition of the discrete radon transform. Radon’s Transform is a mapping between (weighted) L_{2}-Spaces:

R : L_{2}(Ω, W^{−1}) ⟶ L_{2}(S^{n−1} × IR, w^{−1})
Rf(θ, r) = ∫_{x·θ=r} f(x) ds

We consider a discrete version:

R_{d} : span(φ_{1}, ..., φ_{n}) ⟶ IR^{m}
(R_{d}f)_{l} = ∫_{L_{l}} f(x) ds , l = 1, ..., m.

This is a linear mapping between finite-dimensional hilbert spaces. As we know the SVD for matrices, we can compute the singular value decomposition of R_{d} as well. Here is a more general definition of the SVD, taken from [4]: Definition: Let R be a linear operator between (separable) Hilbert spaces X, Y; R : X ⟶ Y The triple {u_{k}, v_{k}, σ_{k}}_{k≥1} is called a Singular Value Decomposition of the operator R if {u_{k}}_{k≥1} is a complete orthonormal system in X, {v_{k}}_{k≥1} is an orthonormal system in Y, {σ_{k}}_{k≥1} is a set of non-negative real numbers,

Au_{k} = σ_{k}v_{k} and A^{∗}v_{k} = σ_{k}u_{k} .

The singular values σ_{k} are usually ordered such that σ_{1} ≥ σ_{2} ≥ ... ≥ 0.

The SVD of the matrix A in (1.3) now gives us the singular value decomposition of R_{d}: Theorem 2: Let H be a finite dimensional subspace of L_{2} and {φ_{1}, ..., φ_{n}} an orthonormal basis of H. Define A = (a_{lk}) ∈ IR^{m×n} by

a_{lk} := (R_{d}φ_{k})_{l} , l = 1, ..., m, k = 1, ..., n.

Let A = UΣV^{T} denote the SVD of A, notations as in Theorem 1, V = (v_{kj}) and

ψ_{j} = ∑_{k=1}^{n} v_{kj}φ_{k} , j = 1, ..., n.

Then

{(ψ_{j})_{j=1,...,n}, (u_{k})_{k=1,...,m}, (σ_{k})_{k=1,...,min(m,n)}}

is a singular value decomposition of R_{d}. The functions ψ_{j} are called singular functions of R_{d}.

The proof is a straight forward utilization of linearity and orthogonality and can be found in [6].

Using the singular functions of R_{d} as ansatzfunctions of the collocation method, we get

f̃(x) = ∑_{k=1}^{p′} F_{γ}(σ_{k})\frac{1}{σ_{k}} (g, u_{k})ψ_{k}(x) .

The formula for the coefficients is much simpler now and we can give a simple interpretation of regularisation: The Filter F_{γ} directly weights the coefficients of the ansatzfunctions. If we decide to use a truncated singular value decomposition for example, we only drop out some of the ansatzfunctions and solve an overdetermined system. Why is it a form of regularisation to leave out the last 10 singular functions? We can find the answer if we take a look at the singular functions.

Figure 1.2 shows the singular functions for the choice φ_{k}(x) = e^{−λ|x−x_{k}|}. They are ordered in usual reading-order according to decreasing singular values. As the functions φ_{k} are not orthonormal with respect to the usual scalar product in L_{2}, same holds for the singular functions in figure 1.2. The shown scale is valid only for the first singular function. To show the structures, each picture fills the full greyscale range from black to white. With decreasing singular values the singular functions become high frequent. The computation of the coefficient for a high-frequent singular function is less stable then in the low-frequent case, of course. That is why regularisation can be done by reducing the coefficients for singular functions belonging to small singular values.

5 Numerical results

The algorithm will be tested using two phantoms. By “phantoms” we denote artificial temperature distributions produced on a computer for testing the methods. From these phantoms we compute the tomographic data — exact and with simulated measurement errors — and reconstruct the phantoms from these data using the descriped collocation methods.

Figure 1.3 shows the first phantom, the data, and reconstructions. The greyscale lies within 850 and 1250. This is a realistic interval for the temperature in a furnace measured in ^{o}C. In the first row the phantom itself is shown as well as the tomographic data. The data are shown by simply drawing the lines in a greyscale corresponding to the mean temperature on the line. Both, eight and sixteen transceivers are simulated resulting in 24 respectively 96 lines. From these data the phantom is to be reconstructed.

The second row shows reconstructions from the set of 24 lines using γ = 0.1. The pictures are very rough because of the low number of lines. For the second and the third reconstruction the data are corrupted by adding normally distributed noise with a standard deviation of 1% and 2.5% respectively. Although even the reconstruction from badly perturbed data is quite smooth, it has big errors as the right “arm” of the hot part in the phantom is missing completly.

With 96 (16 transceivers) lines we can get much better pictures of course. Here two different regularisations have been used to show their effects. The pictures in the last row fullfill the data in a better way resulting in a slightly better reconstruction from exact data but in serious artefacts with big measurement errors. For the reconstructions in the last two rows the same sets of data have been used. They show that smoothing with this way of regularising is a good way to lower the effects of measurement errors. However, adding normally distributed noise is not the perfect way to model the real measurement error of 2.5% as this error mainly consists of systematic errors. For example errors in κ and M and bending phenomena to not lead to normally distributed noise.

Bad values for κ and M would have serious effects on the absolute values of the final temperature reconstruction but not on structures that can be seen in them. We hope to be able to compensate the effect of bending phenomena in future by using a better model for “travel-paths” of sound. This encourages to use a regularisation for 16-transceivers-measurements that also allows us to reconstruct smaller details. The second phantom (Fig. 1.4) has such details: The blow-in of cold air is simulated by small regions with a lower temperature near the walls of the furnace. These cannot be reconstructed from only 24 lines (8 transceivers) but are clearly found in lowregularised reconstructions using 96 lines (16 transceivers) as can be seen in the last row.

No measurements have been made with 96 lines (16 transceivers), yet. So we can only show a reconstruction from 24 lines (8 transceivers) here (Fig. 1.5). The measured data are taken from [1] and the reconstruction is shown in the way, the film is presented.

References

[1] Derichs, W.; Heß, F.; Menzel, K.; Reinartz, E.: Acoustic Pyrometry: A correlation between temperature distributions and the operating conditions in a brown coal fired boiler. 2^{nd} International Conference On Combustion Technologies For A Clean Environment. Lissabon, July 1993

[2] Derichs, W.; König, J.: Die Schallpyrometrie – Ein Meßverfahren zur Bestimmung der Temperaturverteilung in Kesselfeuerungen. DVV-Kolloquium, September 1990.

[3] Golub, G.H.; van Loan, Ch.: Matrix Computations. Second Edition, North-Holland Verlag 1989.

[4] Maaß, P.: Singular Value Decompositions for Radon Transforms. In: Lecture Notes in Mathematics, Vol. 1497, Herman, Louis, Natterer (Eds.), Mathematical Methods in Tomography. Proceedings. Oberwolfach, 1990.

[5] Natterer, F.: The Mathematics of Computerized Tomography. Teubner, 1986

[6] Sielschott, H.: Tomographie mit wenigen Daten: Bildrekonstruktion aus Daten der Schallpyrometrie mittels Kollokationsverfahren. Master Thesis, Münster, August 1994.

260

Fig. 1.2 Singular functions for φ_{k} = e^{−λ|x−x_{k}|}.

261

Fig. 1.3 Reconstructions of a non-komplex phantom

262

Fig. 1.4 Reconstructions of a komplex phantom Fig. 1.5 Reconstruction from measured data taken from [1]. On the right hand side the

measured data are shown.