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][11][12][13][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 s denote the start point of the th beam; typically the position of the th source, and r d its end at the corresponding detector such that | | = |r s − r d | is the length of the beam. Further, let p s to be the intensity of the beam leaving r s and p d the intensity of the beam arriving at r d , then by integrating (1.1) over the path of the beam yields where p s ≥ p d ≥ 0, and log p denotes the natural logarithm of p. The imaging problem is then to estimate a bounded function c : Ω → + from a set of m noise contaminated data If p s is known without uncertainty and p m contains additive noise n, then y * = log(p s /p d ) denotes the exact data and y = log(p s /p d + n) 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 d 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 Ac 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][16][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 N 2 m yields an underdetermined system of linear equations y = Ac + n, (1.4) where A ∈ m×N 2 , with m = rank(A) and n is additive zero-mean Gaussian noise of magnitude δ = n , and covariance matrix Γ n . 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 ∈ m×m and V ∈ N 2 ×N 2 are orthogonal matrices and Σ ∈ m×N 2 is a diagonal holding the singular values of A in non-ascending order 1 ≥ σ 2 ≥ . . . ≥ σ m . 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Σ m V m . The columns of U span the range of A, those of V N 2 −m its null space N (A) and V m ∈ N ⊥ (A). Despite being underdetermined, we remark that the m singular values of A reduce at a slow rate and thus the value of σ m 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 ∈ N 2 ×N 2 similar to that considered in [8]. This is equivalent to solving the augmented least-squares problem arg min for a positive parameter λ, whose solution     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 [11,19] 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.

Model-based prior information
This Gaussian assumption is consistent with the experimental observations reported in [10,12], 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 c max 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 u 0 we can approximate the maximum centreline velocity as  indicating that u max decreases inversely proportional with the distance form the nozzle, while the velocity away from the centreline is 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 c 0 is the expected concentration and Γ c 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 * ∈ N 2 is the high-dimensional image and Π ∈ N 2 ×N 2 the orthogonal projection operator for functions in N 2 onto a low-dimensional subspace spanned by a small number of linearly independent basis functions {φ 1 , . . . , φ s } comprising the columns of the tall matrix Φ ∈ N 2 ×s . 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 andr the computed solution of the projected, lowdimensional inverse problem. Our methodology consists of two steps: the projection of the highdimensional image in the low-dimensional subspace inducing r * and the numerical solution of the projected inverse problem yieldingr. This process can be diagrammatically depicted as indicating the two sources of errors affecting the solution. The total image error c * − Φr 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 N 2 can be decomposed as where Π c = Φr for a unique r. Denoting by w c = (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 c . To quantify this error further we project w c onto N (A) using the orthogonal projection operator As Pw c ∈ N (A) then the projected model for the low-dimensional variable r ∈ s simplifies to y = AΦr + Aq c + n, (3.5) where q c = (I − P)w c . In regard to its magnitude, there is little to be said about the projectioninduced data error Aq c , 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 specifyr by formulating an appropriate inverse problem. Let B ∈ m×s a low-dimensional projected model matrix with m ≤ s N 2 , 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σ i , relates to that of the high-dimensional A byσ i ≤ σ i for i = 1, . . . , m and, therefore, we can deduce that 1 >σ 1 ≥σ 2 ≥ · · · ≥σ m > 0. Unfortunately, this condition does not guarantee a small condition number κ(B) despite that 1/σ m is small. In fact, it can be shown that for the projected model in (3.8) and the prior information on the high-dimensional variable c from (2.4). The additive error ε inherits the Gaussian properties of the measurement noise n but it is shifted in mean by an unknown Aq c . Theoretically this term can be estimated by simulation as outlined in [17]; however, in many practical situations this tends to be neglected. The reduced model likelihood then becomes For c 0 ∈ N 2 the mean of the flow-based prior probability density p C (c) in (2.4) then where E[c] denotes the expectation of c, and similarly If Γ c = λ −1 I, for a positive λ then by the orthogonality of Φ we get Γ r = Φ Γ c Φ = λ −1 I, hence in forming the posterior density ofr conditioned on y through Bayes' rule, and tracing its unique maximum a posteriori estimator, yields the low-dimensional problem arg min The computational error r * −r 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σ m ≈ 0 and λ > 0, we can investigate how the estimatorr λ is aided by the prior knowledge of r 0 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 r 0 = 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 λ σ m ≈ 0 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 c and r * − r 0 depend exclusively on the available a priori information and the sufficiency of the reduced basis, the error amplification term max i=1,...,m {σ i /(σ 2 i + λ)} > 1 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 r * −r λ does not imply a small overall error c * − Φr λ , 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 c . 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 c , 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.
Let us now look closer into computing the minimum norm solution for the high-dimensional model (1.4) formulated asĉ = arg min c∈ N 2 c such that y = Ac + n.
Our intent is to show that this is a special case of the subspace projection method with If Σ m ∈ m×m is the non-zero block of Σ that holds the singular values of A, then it is straightforward to show thatr = Σ −1 m U y and thereforeĉ = V mr . Note that in this case the computational error attains its minimum possible value r * −r = τ δ. To see this recall that hence combining with the approximation error we have Choosing the basis Φ = V m is thus optimal for minimizing the computational error, as r − r * → 0 as δ → 0, and no regularization error occurs. Unfortunately, this advantage is diminished by the fact that the basis V m 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 * .

A variant algorithm to impose positivity
Having reduced the dimension of the inverse problem through the projection to a lowdimensional basis, we proceed to suggest a modification that precludes the solution from  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 z where now Π z = Φr, and B(z i ) = K(z i )Φ, 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 i ∈ S and compute c i = e z i , (ii) for i = 1, 2, 3, . . .  Table 1. Three simulated experiments on target concentrations with negligibly small approximation error. Each row corresponds to a different c i from (6.1) and the subscript on c is omitted for clarity in the presentation. In all cases, w c ≈ 0.005 and Aq c ≈ 0.001. The image error c − Φr increases with r * − r 0 , while the value of λ reduces. only the positive part of the function c can be uniquely mapped to a corresponding function z = log c, 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 c 0 . This prior is based on model ( 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 lowdimensional projected model yields a rank-deficient matrix B ∈ 126×144 with κ(B) ≈ 10 16 . For each of the targeted concentration functions c 1 , c 2 and c 3 , we have computed low dimensional estimatorsr 1 ,r 2 andr 3 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 c − Φr ≤ w c + r * −r and r * −r ≤ r * − r 0 in all three cases, where r * = Φ † c is the optimum low-dimensional solution, and r 0 = Φ † c 0 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 * − r 0 , and, therefore, a lower λ value is more appropriate when this discrepancy can be large. Note, however, that the overall image error c − Φr is not proportional to r * − r 0 , despite that w c ≈ 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 ĉ + − c appears to be significantly higher compared with that of the unconstrained solution.
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. 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 * − r 0 . 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 c 4 , c 5 and c 6 show that while λ should be reduced in increasing r * − r 0 in order to relax the influence of the prior guess, the  increasing approximation error Aq c in the data reinstates a high regularization parameter value. This being the primary difference in the two tests, the error bounds c − Φr ≤ w c + r * − r 0 and c − Φr ≤ w c + r * −r were found to hold in this w c 0 case. Moreover, the errors in the constrained and unconstrained solutions are at equivalent levels, as opposed to those for w c ≈ 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ĉ + 5 andĉ + 6 despite the higher overall error.

(a) Choosing a value for the regularization parameter
The choice of the regularization parameter λ has a significant impact on the reconstructed Tikhonov estimatorr λ . 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 c , Lcurve 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 [σ rank(B) ,σ 1 ], whereσ rank(B) 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 {r λ 1 ,r λ 2 , . . . ,r λ M }, from which we compute a residual vector r λ i =r λ i −r λ i−1 , for i = 2, . . . , M and then plot the graph of r λ i on a linear scale, as a function of λ evaluated at the midpoint of the interval [λ i , λ i−1 ]. Our criterion for choosing λ from this curve is that in the neighbourhood of the optimal λ the value of r λ will be minimized indicating a region of stability, in the sense that a small perturbation in λ will only cause a small perturbation on the solutionr λ . More precisely, as λ :σ rank(B) →σ 1 , 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 r λ i curve to converge to zero as λ approachesσ 1 due to over-regularization, which prevents the solutions to deviate from the prior r 0 . The graphs of r λ i for the selection of λ in the reconstruction of c 3 , c 5 and c 6 are shown in figure 9. To demonstrate the validity of this approach we have also computed r * −r λ i 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 r * −r λ i is not known, but the proximity of λ * to the minimum of r λ i can be used as a guide for selecting a near-optimal regularization parameter.

Conclusion
This work proposes a computationally efficient solution of the inverse problem in limiteddata 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.