An efficient approach for limited-data chemical species tomography and its error bounds.

We present a computationally efficient reconstruction method for the limited-data chemical species tomography problem that incorporates projection of the unknown gas concentration function onto a low-dimensional subspace, and regularization using prior information obtained from a simple flow model. In this context, the contribution of this work is on the analysis of the projection-induced data errors and the calculation of bounds for the overall image error incorporating the impact of projection and regularization errors as well as measurement noise. As an extension to this methodology, we present a variant algorithm that preserves the positivity of the concentration image.

Introduction

The concerted effort to create efficient energy technologies with reduced greenhouse gas generation, and to implement them on an industrial scale, has already resulted in the identification of numerous pollutant reduction strategies, such as biomass-derived fuels and their application in aero engines, ammonia-based energy storage, high-efficiency modes of automotive engine operation and fuel cell technology. Many more such options will be identified over the coming years. However, building innovative measurement systems to underpin the development of these technologies depends critically upon the sensitivity and resolving power of the in situ diagnostics. Chemical species tomography (CST) addresses this challenge by imaging quantitatively the concentration of chemical particles in gaseous media and exhaust plumes [1].

In CST, data acquisition typically entails a dense array of collimated laser beams propagating through the medium at various angles spanning the half circle. Measuring the intensity of these beams at their sources and diametrically positioned detectors yields the attenuation information within a noise margin, and that relates linearly, under some assumptions, to the chemical species concentration profile at the measurement plane. When a dense sampling of the domain is feasible, a large number of beams and many projection angles are used, and the image can then be stably reconstructed using Fourier transform-based methods like the popular filtered back-projection algorithm [2]. There are, however, harsh industrial environments where collecting many measurements is deemed impractical, allowing only a sparse beam arrangement and a small number of projection angles [3]. This is the paradigm known as limited data tomography and it is typically addressed in the context of algebraic image reconstruction and compressed sensing algorithms whenever the anticipated solution has limited support, i.e. it has a sparse profile (see e.g. [4,5]). In this work, we focus on the case where the number of projection data is substantially less than the degrees of freedom in the image we seek to reconstruct, rendering the inverse problem severely underdetermined. In addition, the sought chemical species concentration images are expected to be smooth, almost dense as a result of the gas mixing and combustion phenomena [6]. Such imaging problems are known to be ill-posed as they lack uniqueness, although a unique image can still be computed subject to enforcing some form of regularization. As discussed in some detail in [7], the choice of regularization is ultimately linked to whether the resulting inverse problem is discrete ill-posed or rank deficient. As the underlying attenuation model falls within the linear regime of the Beer–Lambert Law, the problem is well suited to the Tikhonov regularization [8], as well as iterative algorithms based on the Landweber iteration, exploiting their convergence properties in solving linear, ill-posed problems [9]. A drawback of these methods is their computational inefficiency in handling large underdetermined high-dimensional systems, as they still rely on inverting dense square matrices in high dimension.

Various experimental studies [10-14] have affirmed the sensitivity of light attenuation measurements in near and mid infrared radiation to soot and carbon dioxide particles within a jet's exhaust plume and motivate the application of tomography for in situ characterization. Although the gases are in motion, high-speed or simultaneous data acquisition provides static observation conditions, allowing for time-lapse imaging. Assuming negligible optical scattering, the amount of beam energy absorbed at a point is thought to be proportional to the gas density there. Effectively, the attenuation of a beam with density p>0 along an infinitesimal segment dℓ is [1] where c is the two-dimensional chemical species concentration function. Let r denote the start point of the ℓth beam; typically the position of the ℓth source, and r its end at the corresponding detector such that |ℓ|=|r−r| is the length of the beam. Further, let p to be the intensity of the beam leaving r and p the intensity of the beam arriving at r, then by integrating (1.1) over the path of the beam yields where p≥p≥0, and denotes the natural logarithm of p. The imaging problem is then to estimate a bounded function from a set of m noise contaminated data . Before addressing this inverse problem, we make a brief remark on the impact of noise on the synthesized data compared with the actual measurements {p,p}. If p is known without uncertainty and p contains additive noise n, then denotes the exact data and the noisy. Subtracting the one from the other yields indicating that the synthesis of the logarithmic data suppresses the levels of noise in the actual measurements. The presence of p in the denominator suggest using a high-intensity radiation, which is of course consistent with having a high signal to noise ratio in the measurements. Populating the path concentration integrals in (1.2) for many beams over the range of angles in [0,π) yields the linear operator equation where is the Radon transform data of c, and n some additive noise corrupting the data. The mathematical problem of reconstructing c from y based on (1.3) is well studied and analysed in various textbooks, mostly in the context of X-ray computed tomography (e.g. [2,15-17]). Its theory postulates that a reconstruction of c is feasible and stable subject to the sufficiency of the data y. In other words, the Radon operator has a continuous inverse that leads to a unique solution provided that there are enough data to resolve the degrees of freedom in the image. In the limited-data paradigm, however, a stable inversion is no longer feasible without some form of regularization [15]. In this paper, we show that an effective re-parametrization of the unknown in conjunction with a basic regularization scheme can yield a stable, unique solution at a significantly reduced computational cost. Our approach is fundamentally based on projecting onto a low-dimensional image-feature subspace, an operation that inevitably incurs some information loss. This so-called discretization error is a data component that is typically assumed small enough to be neglected. Part of the scope of this work is to quantify its impact on the image error, particularly when the dimension of the physical domain is large and the number of data is very small by comparison. Discretizing the Beer–Lambert law on an N×N grid of square pixels, with N2≫m yields an underdetermined system of linear equations where A∈ℜ, with m=rank(A) and n is additive zero-mean Gaussian noise of magnitude δ=∥n∥, and covariance matrix Γ. Without any loss of generality we assume that the system is normalized so that ∥A∥=1, and that A admits a singular value decomposition A=UΣV ′ where U∈ℜ and V ∈ℜ are orthogonal matrices and Σ∈ℜ is a diagonal holding the singular values of A in non-ascending order 1≥σ2≥…≥σ. In our notation prime denotes transposition. Expressing V and Σ like it is easy to see that A admits an expansion in a truncated basis A=UΣV ′. The columns of U span the range of A, those of V its null space and . Despite being underdetermined, we remark that the m singular values of A reduce at a slow rate and thus the value of σ is maintained well above zero.

An appropriate method for reconstructing a unique solution from the underdetermined model (1.4) is by formulating the Tikhonov problem, using, for example, a smoothness imposing regularization matrix L∈ℜ similar to that considered in [8]. This is equivalent to solving the augmented least-squares problem for a positive parameter λ, whose solution can be computed directly by inverting an N2×N2 matrix. By contrast, our methodology yields image reconstructions without performing any computations in N2 dimensions.

Model-based prior information

At the FLITES experiment [18], 126 light sources and detectors are mounted on a dodecagonal ring structure encompassing the imaging domain of interest Ω=[−0.75,0.75]×[−0.75,0.75] m, which incorporates the gas plume. The measurement plane is normal to the plume propagation axis at a distance of 2 m from the engine's nozzle as shown in figure 1. At the near field, the flow is predominantly axial with a small dispersion, while the detuner positioned immediately after the optical ring vents the engine exhaust, preventing secondary backflows in the measurement plane. In these conditions, we may assume a free turbulent jet model and simulate an ‘expected’ plume trajectory and velocity field at the measurement plane [19,11] or indeed adopt a more sophisticated fluid model if necessary [20]. In a recent study [7], the authors propose using a squared exponential prior that is consistent with turbulent flow mixing models while preserving the smoothness of the concentration profiles. The correlation between the magnitude of the velocity and the concentration of the particles in the plume hints for making a Gaussian assumption on the anticipated concentration profiles. The use of smoothness imposing prior models on this particular problem was originally suggested in [21], and here we extend it to accommodate flow-specific prior information.

Figure 1.

(a) A schematic of the FLITES experimental set-up at INTA testing facilities indicating the position of the engine, the optical ring and the detuner for gas extraction. (b) An oblique view of the engine and the optical plane indicating the beams from one projection angle. (Online version in colour.)

This Gaussian assumption is consistent with the experimental observations reported in [12,10], although these make reference to variant Gaussian functions modified to have a flatter top, which motivates the use of the so-called super-Gaussian functions (figure 2) where x>0 is the distance along the direction of propagation, ϱ is the cross-jet radial distance from the jet's centreline, θ>0 is the half-angle of the plume cone and is the maximum concentration level at the centreline of the plume that is assumed to scale linearly to the plume velocity there. From the momentum conservation principle [22], if the jet nozzle has a circular shape with diameter d and the gas velocity there is u0 we can approximate the maximum centreline velocity as indicating that decreases inversely proportional with the distance form the nozzle, while the velocity away from the centreline is

Figure 2.

Lateral and vertical cross sections of plume concentration profiles from the super-Gaussian prior model assuming a circular nozzle of diameter d fixed at x=0.1 and the measurement plane at x=2 m (bottom images). For the three plumes, the respective flow parameters are d=0.3 m, θ=0.2 rad (a), d=0.4 m, θ=0.2 rad (b) and d= 0.4 m, θ=0.3 rad (c). In all cases, the fluid velocity at the nozzle is u0=250 m s−1. (Online version in colour.)

Based on this flow model, the concentration image of a certain chemical species at the measurement plane, can be expressed as a random field from a normal probability density function where c0 is the expected concentration and Γ a positive definite covariance matrix.

Subspace projection

The subspace projection relies on the hypothesis that there exists a potentially low-dimensional basis of functions that captures the dominant features of the sought image, and as such this relies explicitly on the available prior information about the solution [16,23]. Suppose c*∈ℜ is the high-dimensional image and Π∈ℜ the orthogonal projection operator for functions in ℜ onto a low-dimensional subspace spanned by a small number of linearly independent basis functions {ϕ1,…,ϕ} comprising the columns of the tall matrix Φ∈ℜ. As we discuss in the following section, there are various options for choosing these functions, and although it is not absolutely necessary, for stability reasons a choice of orthonormal basis is often preferred [16]. Here we consider the case for Φ′Φ=I and ΦΦ′=Π, where Π is idempotent. Further let r*=Φ′c* be the optimal coefficients in the approximation of the targeted image on S and the computed solution of the projected, low-dimensional inverse problem. Our methodology consists of two steps: the projection of the high-dimensional image in the low-dimensional subspace inducing r* and the numerical solution of the projected inverse problem yielding . This process can be diagrammatically depicted as indicating the two sources of errors affecting the solution. The total image error comprises the subspace approximation error, depending on the skill of the basis to express c*, and the computational error which reflects the performance of the image reconstruction algorithm to approximate r* given the properties of the low-dimensional inverse problem and the impact of noise in the data. To estimate this error consider that any bounded image in ℜ can be decomposed as where Πc=Φr for a unique r. Denoting by w=(I−Π)c the subspace approximation error, otherwise put the component of c that does not belong in S, the linear model for the measurements (1.4) becomes yielding the additional data error term Aw. To quantify this error further we project w onto using the orthogonal projection operator P=V V ′, such that As then the projected model for the low-dimensional variable r∈ℜ simplifies to where q=(I−P)w. In regard to its magnitude, there is little to be said about the projection-induced data error Aq, aside the rather obvious upper bound as ∥I−P∥=∥I−Π∥=1. Ultimately, to estimate the overall image error we have and using the triangle inequality and the orthogonality of Φ we obtain

Clearly, the first component is the subspace approximation error and it relies entirely on the choice of basis Φ, which in turn reflects the credibility of the prior information one has about c*. To quantify the second term in (3.7), we must first specify by formulating an appropriate inverse problem. Let B∈ℜ a low-dimensional projected model matrix with m≤s≪N2, then the model (3.5) is expressed as with B=AΦ. As ∥A∥=∥Φ∥=1 then we can see immediately that ∥B∥≤1, while from [24], the ith singular value of B, denoted as , relates to that of the high-dimensional A by for i=1,…,m and, therefore, we can deduce that . Unfortunately, this condition does not guarantee a small condition number κ(B) despite that 1/σ is small. In fact, it can be shown that for s17]; however, in many practical situations this tends to be neglected. The reduced model likelihood then becomes For c0∈ℜ the mean of the flow-based prior probability density p(c) in (2.4) then where denotes the expectation of c, and similarly If Γ=λ−1I, for a positive λ then by the orthogonality of Φ we get Γ=Φ′ΓΦ=λ−1I, hence in forming the posterior density of conditioned on y through Bayes’ rule, and tracing its unique maximum a posteriori estimator, yields the low-dimensional problem with solution

The computational error depends on how accurately the solution of the reduced inverse problem (3.11) approximates r* given the necessity of regularization and the various noise components embedded in ε. Assuming the general case where and λ>0, we can investigate how the estimator is aided by the prior knowledge of r0 and how it is corrupted by the measurement noise and approximation errors in ε. Let y=y*+ε and suppose the prior guess on the solution satisfies r0=r*+δr, for an arbitrary δr, then by (3.11) we obtain Rearranging, and taking norms we obtain the computational error upper bound for (3.7) as and therefore Introducing to the general error bound (3.7), we obtain

which indicates that the overall image error in our approach is the sum of the approximation error norm and the norm of the discrepancy between the prior-based expectation and the true image. While the error terms ∥w∥ and ∥r*−r0∥ depend exclusively on the available a priori information and the sufficiency of the reduced basis, the error amplification term can be precomputed for any choice of λ. Given the slow reduction rate in the large singular values of B the amplification may turn out to be small. Note, however, that a small computational error does not imply a small overall error , as the next section demonstrates.

Choice of approximation basis and a special case

We have seen that projecting the unknown high-dimensional image into a low-dimensional subspace reduces the computational complexity by replacing the large matrix A with a smaller matrix B. This computational advantage causes an increase in the ‘noise’ of the data from n to ε, in manifestation of the subspace approximation error w. Moreover, this projection may lead to a rank deficient inverse problem and a new matrix with dispersed singular values, despite the clustering in those of the original large matrix. Choosing the basis Φ appropriately is thus instrumental in the image reconstruction process, in the sense that it controls the approximation error w, but it also has an impact on the computational error. One conclusion that becomes apparent even as early as this stage of the investigation is that the conventional local basis functions with pixel-wise constant support, often adopted by default, might not always be the best basis to model the unknown image when having limited measurements. We propose the use of global basis functions, a smooth basis of orthonormal functions that satisfy the prior (2.4) [16]. Subject to a sufficiently large s, this basis can accommodate a large range of functions, whose features include the expected concentration (see e.g. figure 3). We note that the computational advantage of this approach compared with solving the smoothness imposing high-dimensional Tikhonov regularization problem (1.5) may come at a cost of a higher projection approximation error, as the admissible smoothness of the image is explicitly enforced in terms of this basis. Moreover, other non-smoothness related priors in the context of Bayesian inference, or affine constraints in optimization schemes, may prove challenging to cast in the form of conforming global bases, and in those circumstances it might be easier to revert to the conventional pixel-based basis.

Figure 3.

(a) A high-dimensional image c* and next its projection errors w in s=64 (b) and s=144 (c) discrete cosine transform functions. The relative errors ∥c*−Πc*∥/∥c*∥ are 2×10−3 and 3×10−4, respectively. (Online version in colour.)

Let us now look closer into computing the minimum norm solution for the high-dimensional model (1.4) formulated as Our intent is to show that this is a special case of the subspace projection method with Φ=V (s=m) in which case , Π=P⊥, and thus the basis is optimal in terms of the computational error. Effectively, and hence the measurement error due to the subspace approximation error vanishes leaving ε=n. This unique estimator can be computed rather efficiently, and stably, by inverting only an m×m square matrix B=AV . Given that AA′ is invertible by virtue of κ(A) being small, then can alternatively be computed by solving the low-dimensional least-squares problem If Σ∈ℜ is the non-zero block of Σ that holds the singular values of A, then it is straightforward to show that and therefore . Note that in this case the computational error attains its minimum possible value . To see this recall that r*=V ′c* and hence combining with the approximation error we have Choosing the basis Φ=V is thus optimal for minimizing the computational error, as as , and no regularization error occurs. Unfortunately, this advantage is diminished by the fact that the basis V is unsuitable for reconstructing the anticipated smooth solutions, see for example, the projection of an image into this basis in figure 4, and results in a very large approximation error ∥c*−Πc*∥.

Figure 4.

(a) A target image c* consistent with the prior flow-concentration model (2.4) and (b) its projection V V ′c*. The relative approximation error ∥c*−Πc*∥/∥c*∥ is 0.66 with s=m=126 functions and the condition number of the square matrix B=AV is κ(B)=2.5. (Online version in colour.)

A variant algorithm to impose positivity

Having reduced the dimension of the inverse problem through the projection to a low-dimensional basis, we proceed to suggest a modification that precludes the solution from attaining non-positive values. This is aimed at preserving the feasibility of the images, as chemical concentration is by definition non-negative. In particular, the attenuation model is reformulated in terms of the logarithm of the unknown concentration function [25]. In this context, denoting then the original model (1.4) becomes nonlinear in the new parameters which can then be linearized around a point z∈S to where C=diag(c) and c=e. Denoting K(z)=AC and y(z)=y−Ac+K(z)z, we cast the ith linear model for y at the iterate z as Note that the new unknown variable is still in the high-dimension, so similar to the unconstrained case we can project it onto S to get z=Πz+w where now Πz=Φr, and B(z)=K(z)Φ, yielding the low-dimensional model Assuming once again that the basis functions Φ spanning S are orthonormal then a positive subspace projected estimator of the image can be obtained iteratively using the following algorithm:

(i) Initialize z∈S and compute c=e,

(ii) for i=1,2,3,…

(a) Compute matrix B(z)=K(z)Φ and vector y(z),

(b) Compute using (3.11) with B=B(z) and y=y(z),

(c) Update ,

(iii) end

(iv) .

Although we do not derive error bounds for this positively constrained estimator, we note that as c can no longer admit zero values, some information loss is incurred in the transformation from c to z, and, therefore, the bounds of the previous section no longer apply. In this case, only the positive part of the function c can be uniquely mapped to a corresponding function , while a zero concentration value will typically be assigned a very small yet still positive value.

Numerical results

In order to evaluate the performance of the proposed algorithms and verify the image error bound (3.13), two sets of image reconstruction experiments were performed using simulated data infused with white Gaussian noise of standard deviation equal to 5% of the mean measurement value. In both instances, we consider concentration functions that are similar but not equal to the prior guess c0. This prior is based on model (2.1) with x=2 m, c=5 and θ=0.2 rad, and it appears in figure 3 plotted in a domain Ω=[−0.26,0.26]×[−0.26,0.26] with normalized distance units. The first three target images are of the form where ϕ denotes the jth function in the approximation basis (e.g. the jth column of Φ), β1=0.25, β2=1 and β3=3 is a parameter that controls the deviation of the target function from c0, and ϱ1,ϱ2 are the coordinates on the plane of measurement. The definitions c1, c2 and c3 were chosen to yield a small approximation error, e.g. ∥w∥/∥c∥≈10−4, and consequently the approximation error component in the data Aq is small enough to be ignored, ∥Aq∥/∥y*∥≈10−4. The measurements have been computed on a fine 100×100 square grid model of Ω, while for the imaging problem a coarser 70×70 grid was used. Consistent with the FLITES experiment, 126 attenuation measurements were computed from six projection angles, and a model matrix A∈ℜ126×4900 was assembled. On the coarse grid, we also define Φ∈ℜ4900×144 whose columns represent an orthonormal basis of s=144 discrete cosine basis functions. Forming the low-dimensional projected model yields a rank-deficient matrix B∈ℜ126×144 with κ(B)≈1016. For each of the targeted concentration functions c1, c2 and c3, we have computed low dimensional estimators , and using (3.11). The results obtained are summarized in table 1 and the respective images appear in figure 5. A quick glance at the table confirms the derived error bounds and in all three cases, where r*=Φ†c is the optimum low-dimensional solution, and r0=Φ†c0 the projection of the prior concentration on S. The value of the regularization factor λ used in each case was estimated heuristically as outlined at the end of this section. As anticipated, when the approximation error is very small, the overall reconstruction error scales to the discrepancy between the true image and the prior guess ∥r*−r0∥, and, therefore, a lower λ value is more appropriate when this discrepancy can be large. Note, however, that the overall image error is not proportional to ∥r*−r0∥, despite that ∥w∥≈0.005 remains fixed in all tests. For comparison, a positively constrained solution is also obtained by performing four iterations of the algorithm in section 5. The algorithm exhibits a fast convergence (figure 6), but the value of the error appears to be significantly higher compared with that of the unconstrained solution.

Table 1.

Three simulated experiments on target concentrations with negligibly small approximation error. Each row corresponds to a different c from (6.1) and the subscript on c is omitted for clarity in the presentation. In all cases, ∥w∥≈0.005 and ∥Aq∥≈0.001. The image error increases with ∥r*−r0∥, while the value of λ reduces.

i	∥c∥	∥r*−r0∥	λ	∥c−Φr^∥	∥r^∗−r^∥	∥c−c^^+∥
1	15.856	01.060	0.134	0.783	0.783	1.359
2	16.369	04.240	0.016	1.035	1.035	3.566
3	20.272	12.722	0.006	1.659	1.659	9.061

Figure 5.

(a–c) are, respectively, the reconstruction results for c1, c2 and c3 concentration functions from equation (6.1). At the top row, the true concentration targets at fine discretization with N=100, and below their projection Πc, projected image reconstruction and fourth iteration of the positively constrained estimator all in resolution N=70, assuming a subspace of s=144 discrete cosine transform functions. (Online version in colour.)

Figure 6.

Convergence of the positively constrained algorithm on the reconstruction of in (6.1) (a) and in (6.2) (b).

To investigate the impact of a significant approximation error on the image reconstruction, another set of simulated experiments were performed using three different target concentrations. Assuming the same prior for consistency, these new target concentration functions are where (ξ,ϑ) are the polar coordinates of (ϱ1,ϱ2). The coefficients β4=0.01, β5=0.1 and β6=0.5 control the size of the prior guess discrepancy ∥r*−r0∥. Contrary to the previous tests, these incur a significant approximation error in the adopted basis, as it can be seen by comparing the images at first and second rows of figure 7. In this case, we consider a gradually increasing approximation error in conjunction with an increasing prior discrepancy. The results of these experiments are tabulated in table 2, but the reconstructions of c4, c5 and c6 show that while λ should be reduced in increasing ∥r*−r0∥ in order to relax the influence of the prior guess, the increasing approximation error Aq in the data reinstates a high regularization parameter value. This being the primary difference in the two tests, the error bounds and were found to hold in this ∥w∥≫0 case. Moreover, the errors in the constrained and unconstrained solutions are at equivalent levels, as opposed to those for ∥w∥≈0, where the unconstrained estimators were profoundly superior. However, a closer look at the images in figure 7 reveals that quantitatively, the high concentration levels are more accurately reconstructed in the constrained solutions and despite the higher overall error.

Figure 7.

(a–c) are, respectively, the reconstruction results for c4, c5 and c6 concentration functions from equation (6.2). At the top row, the true concentration targets at fine discretization with N=100, and below their projection Πc, projected image reconstruction and fourth iteration of the positively constrained estimator all in resolution N=70, assuming a subspace of s=144 discrete cosine transform functions. (Online version in colour.)

Table 2.

Three simulated experiments on target concentrations with significant approximation error. Each row corresponds to a different c from (6.2). Notice that the overall error is less than the sum of the approximation and computational errors.

i	∥c∥	∥wc∥	∥Aqc∥	∥r*−r0∥	λ	∥c−Φr^∥	∥r^∗−r^∥	∥c−c^^+∥
4	16.633	10.325	0.115	11.050	0.308	10.712	10.633	11.129
5	24.940	13.256	1.149	10.504	0.016	14.086	12.468	16.768
6	67.285	16.283	5.749	52.522	0.064	19.351	10.456	22.574

Choosing a value for the regularization parameter

The choice of the regularization parameter λ has a significant impact on the reconstructed Tikhonov estimator . In abstract form, this problem has been studied extensively in the context of Bayesian estimation for linear inverse problems and various approaches have been proposed on optimizing this parameter in conjunction to the noise levels in the data and the confidence one has on the prior information [16]. The generalized cross validation (GCV) and L-curve methods are among the most popular tools for choosing λ, for example Ma & Cai [26] on their implementation in chemical species tomography. These schemes refer predominantly to linear Gaussian models assuming some knowledge on the noise's first and second statistical moments. Alternatively, one may use criteria based on the singular value decomposition in the cases where the model matrices are rank deficient [27]. As in our setting, regularization is applied to the projected problem, which includes an unknown data error component Aq, L-curve and GCV are not straightforward to apply. Instead, we resort in a heuristic criterion for choosing a sub-optimal λ as follows. We generate a dense grid of M logarithmically equally spaced values in the interval , where is the smallest, non-zero singular value of B, which can be easily identified from the ‘jump’ in the singular values as illustrated in figure 8. Assigning in turn, each value in that interval to λ and solving the problem (3.11) yields a set of solutions , from which we compute a residual vector , for i=2,…,M and then plot the graph of on a linear scale, as a function of λ evaluated at the midpoint of the interval [λ,λ].

Figure 8.

The singular values of matrices A∈ℜ and B=AΦ with Φ∈ℜ, where m=126, s=144 and N=70. Note the slow rate in the reduction of the singular values of A and the effective rank deficiency of B. Φ consists of a basis of discrete cosine functions.

Our criterion for choosing λ from this curve is that in the neighbourhood of the optimal λ the value of will be minimized indicating a region of stability, in the sense that a small perturbation in λ will only cause a small perturbation on the solution . More precisely, as , we anticipate the residual norm to gradually reduce as regularization sets in to alleviate the instabilities caused by the rank deficiency of matrix B. This reduction will lead to a minimum near the optimal λ*, before the residual norm increases again in response to the regularization beginning to affect some of the larger singular values that are clustered closely together. Thereafter, we expect the curve to converge to zero as λ approaches due to over-regularization, which prevents the solutions to deviate from the prior r0. The graphs of for the selection of λ in the reconstruction of c3, c5 and c6 are shown in figure 9. To demonstrate the validity of this approach we have also computed using r*=Φ†c for each value of λ in the same interval, regarding as optimal parameter choice λ* the argument that minimizes this error discrepancy. In reality of course is not known, but the proximity of λ* to the minimum of can be used as a guide for selecting a near-optimal regularization parameter.

Figure 9.

Plots of the and for 1000 values of λ equally spaced within [10−4,1]. Note the proximity of the optimal regularization parameter λ* to the argument that minimizes . The graphs refer to the experiments in reconstructing (a), (b) and (c). (Online version in colour.)

Conclusion

This work proposes a computationally efficient solution of the inverse problem in limited-data chemical species attenuation tomography with flow model-based prior information. We have showed that projection of the unknown function in a smooth low-dimensional basis reduces drastically the dimensionality of the problem and probably yields rank-deficient linear systems that are suitable to Tikhonov regularization. This image reconstruction approach is computationally efficient as it avoids inverting high-dimensional matrices; however, it incurs errors due to the subspace projection and the need for regularization in obtaining a stable solution of the projected inverse problem. To quantify these errors and their impact on the reconstructed image, we provide an upper bound on the overall image error which incorporates the subspace approximation and regularization errors and can be used in the interpretation of the reconstructed images. In this bound the approximation error appears both as an offset as well as a component of the data error, however, the amplification factor of these errors in the image reconstruction is maintained at small levels. This framework is complemented by an variant algorithm for solving the positively constrained imaging problem and a heuristic method for selecting the required regularization parameter. The numerical simulations affirm the performance of the algorithms proposed and the validate the error bounds.