## Abstract

A regularization approach of iterative algorithms was proposed to reconstruct the two-dimensional temperature and concentration distributions based on linear multispectral absorption tomography (MAT). This method introduces a secondary prior into a classical iterative algorithm *via* regularization to improve the reconstruction accuracy. Numerical studies revealed that the regularized iteration outperformed the classical and superiorized versions under various noisy conditions and with different number of spectral lines. The algorithms were also tested with the existing experimental data of a premixed flat flame produced by a McKenna burner. The comparison between the reconstructions and the measured temperature profile using thermocouples confirmed the superiority of our proposed regularized iterative method.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Laser absorption spectroscopy (LAS), due to its species specificity, high sensitivity and ease of implementation, has become one of the standard tools for gas sensing [1–3]. Past decades have seen its development of numerous modalities such as wavelength modulation spectroscopy [4], dispersion spectroscopy [5], cavity enhanced absorption spectroscopy [2], dual comb spectroscopy [6], multi-mode absorption spectroscopy [7], and noise-immune cavity-enhanced optical heterodyne molecular spectrometry (NICE-OHMS) [8]. Due to these advancements, LAS has found extensive applications such as diagnosis of gas turbine engines [9] and rotating detonation engines [10], remote sensing [11], and monitoring of microbial growth [12], to list a few. Despite these successful demonstrations, its line-of-sight (LOS) nature prevents its applications to situations where the physical properties vary along the LOS.

To enable spatial resolution, multiple LOS measurements are typically required and de-convolved either in a linear or a nonlinear tomographic manner [13]. Such combination is called tomographic absorption spectroscopy (TAS) [14]. In the past few decades, TAS based on tunable diode lasers has become mature for combustion and gas dynamics [13–19]. With the development in laser technologies, rapid acquisition of broadband absorption spectrum has provided huge opportunities for multispectral absorption tomography (MAT) as abundant information in the spectral domain can be utilized in the tomographic inversion, making it less sensitive to noise [3,20,21].

Numerous studies have been conducted to discover effective methods to solve the inversion problem of MAT [20,22,23]. Cai *et al.* [20] applied a metaheuristic algorithm to simultaneously reconstruct temperature and concentration fields [20] *via* a nonlinear optimization approach. However, this approach suffers from low computational efficiency as the cost function in the optimization process is repeatedly assessed with a large number of times. In this sense, the linear approach to MAT, which is easier to calculate, has found its advantage [24,25]. The linear approach mainly consists of two steps. In the first step, the distribution of absorption coefficient at each spectral feature is reconstructed using an algorithm such as algebraic reconstruction technique (ART) [26] and maximum likelihood expectation maximization (MLEM) [27]. In the second step, a regression process is applied to fit the corresponding temperature and species for each discretized element. However, since the reconstruction is performed one by one for each absorption feature, the benefit of multispectral information is not fully explored.

In practical applications such as engine measurements, typically only a few LOS measurements are available, making the tomographic inversion problem severely ill-posed. Therefore, a prior reflecting *e.g.* sparseness or smoothness should be added into the reconstruction process to provide more information [28] and improve the reconstruction accuracy. The regularization for nonlinear MAT has been discussed thoroughly by Dai *et al.* [29], and priors for both temperature and absorbing species concentration, including sparseness [30] and model-based priors [31] were assessed. Since the nonlinear solution of MAT is typically modeled as an optimization problem, the implementation of regularization is rather straightforward by adding additional penalty terms to the cost function. However, this is non-trivial for iterative methods as those applied to linear MAT. Superiorization had been found able to handle this issue and the superiorized versions had been proven superior to their respective classical counterpart [32–34].

However, in the superiorized algorithm, the regularization is imposed on the absorption coefficients, which are dependent physical quantities rather than fundamental properties, such as temperature and concentration, of the absorbing gas. In addition, the regularization for each spectral line is independent from each other. Again, the multispectral information is not taken full advantage of since both the reconstruction and the regularization are performed separately for each spectral line. The main contribution of this work is to develop a regularization approach specifically designed for iterative methods. Different from superiorization, the perturbations were directly imposed on the temperature and the concentration fields which would then be fed into the iteration process of each spectral line. By this way, the multispectral information can be utilized simultaneously during the inversion process. Both simulative and experimental studies presented in this work demonstrated its superiority.

## 2. Methods

This section introduces the mathematical formulation of MAT, the classical iterative algorithms, the superiorization and the proposed regularization approach. It should be clarified that Sections 2.1 to 2.3 focus on the existing methods, and Section 2.4 is the major contribution of this work.

#### 2.1 Mathematical formulation of MAT

The principle of TAS has been discussed extensively in numerous papers [25,35,36]. A recent review of both linear and nonlinear TAS has been provided by Cai *et al.* [13]. A brief summary is provided here to facilitate the following discussion.

When a laser beam propagates through a region of interest (ROI) filled with gaseous medium, it will interact with the gas molecules. Phenomena including refraction, scattering and absorption will ensue due to the interaction. In TAS, molecular absorption is assumed to be dominant and other effects are considered to be negligible. The natural logarithm of the ratio between the incident and the transmitted laser intensity at a certain wavelength ($\lambda $) is defined as absorbance according to Beer-Lambert law [37]:

*p*is the absorbance, ${I_t}$ is the transmitted intensity, ${I_0}$ is the incident intensity,

*P*is the pressure, $ T$ is the temperature,

*C*is the mole fraction of the absorbing species,

*e.g.*H

_{2}O;

*S*is the line strength, $\varphi $ is the line shape function, which can be approximated with a Voigt profile for general combustion applications, and $\alpha $ is the absorption coefficient at a specific location denoted by

*l*.

The line strength $S({T,\lambda } )$ has a close connection with temperature and can be expressed as

*Q*is absorber’s partition function; $E^{\prime\prime}$ is the lower-state transition energy of the corresponding spectral line; and $h,c,{k_B}$ are Planck constant, speed of light and Boltzmann’s constant, respectively. The information of $S({{T_0},\lambda } ),\; Q$ and $E^{\prime\prime}$ can be found in the HITRAN database [38]. Besides, when the magnitude of temperature is lower than 2500 K and in the near-infrared range, the last term of Eq. (2) can be neglected [39].

Figure 1 illustrates the general layout of a TAS system. The ROI is discretized into $M = n \times n$ grids; the physical parameters in each grid is assumed to be uniform. Thus, for the *i*-th beam, Eq. (1) can be expressed as

*i*-th beam, $T(j )$ and $C(j )$ are the temperature and concentration in the

*j*-th grid, and $L({i,j} )$ is the length of the

*i*-th beam within the

*j*-th grid.

When there are *Q* beams in total, Eq. (3) can be re-written in a matrix form as

To be clear, the boldfaced symbols refer to matrices (or vectors), while normal symbols refer to scalars or elements in a matrix (or a vector). The vector ${\boldsymbol p}$ can be obtained by LOS measurements, ${\boldsymbol L}$ can be calculated according to the beam arrangement. The absorption coefficient ${\boldsymbol \alpha }(\lambda )$ can be solved with tomographic algorithms. For MAT, each wavelength corresponds to a linear equation system defined by Eq. (4), which has to be solved separately.

#### 2.2 Classical iterative algorithms

The inversion problem of MAT can be solved in a linear manner which mainly consists of two steps. In the first step, the distribution of absorption coefficient at each wavelength, ${\boldsymbol \alpha }(\lambda )$, is reconstructed, for example, with an iterative algorithm. In the second step, a regression process is applied to each grid, fitting temperature and species concentration from the multiple reconstructed absorption coefficients. There are numerous iterative algorithms and their variants for the first step. Testing and comparing them exhaustively is out the scope of this work. Thus, only two most prevalent ones *i.e.* ART and MLEM are investigated here for a demonstrative purpose.

The ART algorithm, originated from the Kacmarz’s method proposed in 1937 [40] was first applied in X-ray tomography by Herman and his coworkers in 1970 [26]. Its iteration process can be mathematically described as [41]:

*k*indicates the current iteration step; ${{\boldsymbol L}_i}$ refers to the

*i-*th row of the weight matrix

**; ${{\boldsymbol p}_i}$ is the**

*L**i-*th element of the absorbance; and ${\lambda _{ART}} \in ({0,2} )$ is the relaxation factor, which controls the convergence rate [42]. In the following studies, ${\lambda _{ART}}$ is set to 0.03.

The MLEM algorithm was first proposed by Dempster in 1977 [27] and has been widely adopted for tomographic reconstructions since 1980s [43], especially in positron emission tomography (PET) [44]. Its iteration step can be expressed as [44]

*i-*th predicted LOS measurement according to the intermediate solution as

In the equations above, *k* again refers to the iteration step; *M* is the total number of pixels; *j* enumerates all the *M* elements in the vector ${\boldsymbol \alpha }{(\lambda )^{({\boldsymbol k} )}}$ within each iteration; *i* indexes the rows in the weight matrix; and *Q* is the total number of LOS measurements.

The termination criterion of an iterative algorithm is usually based on the relative change in the intermediate solution between two consecutive iterations and can be defined as

As for the second step, a regression process [24] is applied to all grids one by one. Suppose a total number of *W* spectral lines were used in the measurement, for the *j*-th grid, ${\alpha _j}({{\lambda_1}} ),{\alpha _j}({{\lambda_2}} )\ldots \; {\alpha _j}({{\lambda_W}} )$ can be reconstructed in the first step. Considering Eqs. (1)–(2) for the *j*-th grid and *r*-th spectral line, the local monochromatic absorbance coefficient ${\alpha _j}({{\lambda_r}} )$ can be expressed as

The logarithmic format of Eq. (9) is

It can be noted that $\ln ( {\frac{{{\alpha_j}({{\lambda_r}} )}}{{S({{T_0},{\lambda_r}} )}}} )$ has a linear relationship with the lower-state energy . Since ${\alpha _j}({{\lambda_r}} )$ can be reconstructed in the previous step and $S({{T_0},{\lambda_r}} )$ and $E_r^{{\prime\prime}}$ are constants given by the HITRAN database [38], a regression process can be implemented for each grid to find the linear relationship between $\ln ( {\frac{{{\alpha_j}(\lambda )}}{{S({{T_0},\lambda } )}}} )$ and $E^{\prime\prime}$, *i.e.* $\ln ( {\frac{{{\alpha_j}(\lambda )}}{{S({{T_0},\lambda } )}}} ) = Y(j )+ Z(j )\cdot {E^{{\prime\prime}}}$. Thus, temperature at the *j*-th grid can be calculated as

*j*-th grid can thence be calculated as

Repeating the regression process for all girds, the distributions of temperature and concentration can be obtained.

#### 2.3 Superiorization of iterative algorithms

The concept of superiorization is originated from the idea of perturbation resilience proposed by Batnariu *et al.* [32] in 2007. This method imposes perturbations to the intermediate solution during the iteration process to obtain a more accurate solution [33,45,46].

To solve Eq. (4), various iterative algorithms including ART and MLEM, have been proposed to find the solution. However, sometimes the solution is under certain constraint, for instance, smoothness. Thus, the reconstruction should take two criteria into consideration, *i.e.* the primary criterion Eq. (4) and the secondary criterion $\varphi ({{\boldsymbol \alpha }(\lambda )} )$, which is usually *a priori* information of ${\boldsymbol \alpha }(\lambda )$. On the condition that the primary criterion has been satisfied, superiorized algorithm applies bounded perturbations to the solution to optimize the secondary criterion at the same time, making the solution *superior*. So far, the most widely adopted secondary criterion is total variation (TV) [30], reflecting piecewise smoothness or sparseness. For a two-dimensional image ** x**, its TV is defined as

Another widely used prior is smoothness [47], which can be quantified by Tikhonov function as:

*n*is the number of grids that are adjacent to the grid $({i,j} )$.

To reduce the secondary criterion, a perturbation in the gradient direction of$\; \varphi ({\boldsymbol x} )$ is added to the intermediate solution in each iteration step. For an original iteration step expressed as ${{\boldsymbol x}^{({k + 1} )}} = {\boldsymbol B}{{\boldsymbol x}^{(k )}}$, the superiorized version can be expressed as

*k*-th iteration; ${\boldsymbol B}$ is the linear operator that describes the iteration; and ${\beta _k}$ is a positive number that guarantees the reduction of $\varphi ({\boldsymbol x} )$.

The superiorization of ART [45] and MLEM [48] has been successfully demonstrated, and has found applications in many fields such as medical computed tomography [33,45,49] and magnetic resonance imaging (MRI) [48,50,51]. However, its application in MAT has not been reported yet. Here, we examine the superiorized version of ART and MLEM as well as other iterative algorithms such as multiplicative algebraic reconstruction technique (MART) [52] and Landweber algorithm [53].

The flow chart of a superiorized algorithm is shown in Fig. 2, and a general iteration step can be summarized as: (1) a perturbation is calculated in the gradient direction of the secondary criterion; (2) judge whether the perturbation will make the secondary criterion decrease; (3) if not, shrink the perturbation and repeat step (2); (4) the solution is adjusted with the perturbation; and (5) an original iteration step is implemented. The stopping rule of a superiorized algorithm can be judged upon residual in the primary criterion or the relative change in the intermediate solution between two consecutive iterations. For example, when the residual is smaller than a threshold *ɛ _{0}*, the iteration terminates.

#### 2.4 Method of regularized iteration

Since the ultimate purpose of MAT is to recover the fields of temperature and species concentration, it is more intuitive to incorporate the priors of temperature and concentration in the iteration process rather than the intermediate physical quantity *i.e.* the distributions of absorption coefficient. To distinguish this approach from superiorization, we call it regularized iteration (RI), in which a constrained optimization problem can be defined as:

The primary criterion in RI contains the measurements of all spectral lines, *i.e.* ${\boldsymbol L\alpha }({{\boldsymbol T},{\boldsymbol C};{\lambda_r}} )= {\boldsymbol p}({{\lambda_r}} ),\; r = 1,2, \ldots ,w.$; and the secondary criterion is the priors of temperature and concentration fields, *i.e.*$\; \varphi ({\boldsymbol T} )$ and $\varphi ({\boldsymbol C} )$. The primary and secondary criteria are optimized alternately in RI: (1) in the *k*-th iteration, the distributions of absorption coefficient $\{{{{\boldsymbol \alpha }^{(k )}}({{\lambda_1}} ),\; {{\boldsymbol \alpha }^{(k )}}({{\lambda_2}} )\ldots {\; }{{\boldsymbol \alpha }^{(k )}}({{\lambda_W}} )} \}$ are reconstructed; (2) temperature and concentration distributions $\{{{{\boldsymbol T}^{(k )}},\; {{\boldsymbol C}^{(k )}}} \}$ are calculated based on the distributions of absorption coefficient in step (1) with regressions; (3) bounded perturbations for ${{\boldsymbol T}^{(k )}}$ and ${{\boldsymbol C}^{(k )}}$ are calculated in the gradient direction of the secondary criterion $\varphi ({\boldsymbol T} )$ and $\varphi ({\boldsymbol C} )$, denoted by ${\boldsymbol v}_T^{(k )}$ and ${\boldsymbol v}_C^{(k )}$ respectively; (4) judge whether the secondary criterion for temperature and concentration fields is reduced; (5) if not, keep reducing the corresponding perturbations; (6) new absorption coefficients ${\boldsymbol \alpha }{^{\prime(k )}}({{\lambda_r}}s )$ are calculated with the perturbed temperature and concentration fields; and (7) the next iteration begins with ${\boldsymbol \alpha }^{({k + 1} )}({\lambda_r} )= {\boldsymbol B \alpha }^{\prime(k )}({{\lambda_r}} )$. Consistent with the other iterative methods, the stopping rule of RI can also be the residual in the primary criterion. The flow chart of a typical RI method is shown in Fig. 3.

As has been mentioned in Section 2.3, the convergence of the superiorized ART and MLEM has been proven in [32] and [48]. As a brief conclusion, as long as the original iterative algorithm satisfies the bounded perturbation resilience (BPR) condition [54], its superiorized version would be convergent with arbitrary bounded perturbations in the iteration process. Superiorization is just one such case of perturbations. In other words, the convergence of an iterative algorithm with perturbation is only related to the original algorithm, as long as the perturbation is bounded. The proof of convergence for RI is similar to that of superiorization. We only need to prove that the perturbation of ${{\boldsymbol \alpha }^{(k )}}$ is bounded. Actually, the only difference is that the perturbation in RI is not applied to ${{\boldsymbol \alpha }^{(k )}}$, but directly to temperature and concentration distributions. According to the Beer-Lambert law, the perturbation in ${{\boldsymbol T}^{(k )}}$ and ${{\boldsymbol C}^{(k )}}$ can be mapped to that of ${{\boldsymbol \alpha }^{(k )}}({{\lambda_r}} )$ *via* the following relationship:

The major difference between superiorization and the RI lies in what priors and how the priors are incorporated in the iteration process. In superiorization, the priors of ${{\boldsymbol \alpha }^{(k )}}$ are introduced in the reconstruction process. Later, a non-regularized regression process is implemented to fit the temperature and concentration for each grid. While in an RI method, regularized iteration and regression are coupled together during the iteration process. In this case, the priors of the temperature and concentration fields rather than those of the distributions of absorption coefficients play a key role in guiding the reconstruction. On the one hand, the secondary criterion of the temperature and concentration fields can be optimized directly. This is more intuitive because temperature and concentration are the fundamental physical properties of substance while absorption coefficient is an artificially defined one.

## 3. Numerical study

#### 3.1 Design of numerical simulations

In this section, simulative studies were conducted with various types of flame phantoms under different “experimental” conditions. The numerical studies can be summarized as: (1) absorption coefficient at each grid and each spectral line was calculated according to Beer-Lambert law with the temperature and concentration phantoms ${{\boldsymbol T}_{act}}$ and ${{\boldsymbol C}_{act}}$; (2) ideal absorbance projections ${{\boldsymbol p}_{cal}}$ were calculated according to Eq. (4); (3) noise was added to projections to mimic real measured projections ${{\boldsymbol p}_{mes}}$; (4) temperature ${{\boldsymbol T}_{rec}}$ and H_{2}O mole fraction ${{\boldsymbol C}_{rec}}$ distributions were reconstructed from the measured projections by a reconstruction algorithm; and (5) analysis and discussion were then performed according to the reconstruction results. It has to be noted that when performing step (2), the pixels should be small enough such that the discretization error could be ignored.

Numerous flame phantoms had been tested in our study. Temperature and H_{2}O mole fraction phantoms mimicking a McKenna flame [55] are shown in Fig. 4. The results based on these two phantoms are presented here as an example. Although the phantoms are axisymmetric, there was no assumption of symmetry and the tomographic reconstruction was still regarded as a two-dimensional inversion problem. For each projection, integrated absorbance of 40 beams was obtained. Uniformly distributed noise was added to the projections to generate noisy data, defined as

*i-*th beam and the

*k-*th spectral line,

*u*is the noise level and $rd$ is a random number which follows a uniform distribution in the range (-1,1). Residual in the primary criterion was set as the stopping criterion, which means the iteration terminates when the criterion becomes less than 0.001. Besides, the maximum number of iterations was set to 5000. All the simulative studies in this section were conducted on a computer with an Intel Core i7-DMI2-X79 PHC 2.60 GHz CPU.

To quantify the reconstruction accuracy, the average errors of temperature and mole fraction were defined as

_{2}O mole fraction, respectively; and

*j*enumerates a total number of

*M*grids.

#### 3.2 Determination of regularization factors

The parameters including $\gamma $, ${\beta _{T,0}}$ and ${\beta _{C,0}}$ used in the RI method reflect how strong the priors are enforced, which would significantly influence the reconstruction accuracy and efficiency. Thus, a best way should be found to optimize those parameters.

The contraction factor $\gamma $ dictates the speed of convergence and is usually chosen as a positive number lower than but close to 1. To determine an optimal $\gamma $, simulative cases were conducted with $\gamma $ equals 0.9, 0.99 and 0.999, respectively. ${\beta _{T,0}}$ and ${\beta _{C,0}}$ were set to 500000 and 0.003, respectively. The reconstruction domain was discretized into 40×40 grids; two orthogonal projections and ten H_{2}O absorption spectral lines were assumed; and uniform noise with a level of 2% was added to the projections. The regularized MLEM was applied here as an example. The results for different cases were compared in terms of the reconstruction accuracy, the computational time and the final values of the secondary criterion, *i.e.* total variation here. The results are summarized in Table 1.

Although the computational time increased rapidly as $\gamma $ increased, the decrement in reconstruction error still indicated that a $\gamma $ closer to one can achieve better performance. Since our main goal was to minimize the reconstruction error and the total variation differed slightly with different $\gamma $, we set $\gamma = 0.999$ in our following studies.

Besides $\gamma $, the initial values of regularization factors ${\beta _{T,0}}$ and ${\beta _{C,0}}$ need to be determined as well. On the one hand, when they are too large, it is unlikely that the secondary criterion will decrease; on the other hand, if they are too small, then the regularization would have little effect. Thus, ${\beta _{T,0}}$ and ${\beta _{C,0}}$ should be reduced to a point when the algorithm starts to accepts the perturbations, denoted as $\beta _T^{\prime}$ and $\beta _C^{\prime}$. The relationships between the factors ${\beta _{T,0}}$ and ${\beta _{C,0}}$, the criteria $\beta _T^{\prime}$ and $\beta _C^{\prime}$, the errors and the time of reconstruction are illustrated in Fig. 5. It can be seen that for both temperature and concentration, the errors were only influenced by the corresponding ${\beta _0}$ parameter. As either factor increased, the corresponding error decreased at first, and finally kept constant. The trend of the criteria were just opposite to that of the errors. And when the error stopped decreasing, the criterion stopped increasing as well. In addition, the computational time had an overall trend of increasing as the factors increased. Therefore, we set the values of ${\beta _0}$ at the inflecting point of the $\beta ^{\prime}$-${\beta _0}$ curve as the initial factor of regularization, which can reduce computational time cost and control reconstruction errors.

#### 3.3 Results and discussion

Simulative studies on several topics were conducted in this section, including the comparison with classical iterative algorithms, the influence of noise, the effects of number of spectral lines, and the selection of secondary criteria.

### 3.3.1 Comparison between tomographic algorithms

Here, we compare the regularized versions of four representative algorithms denoted as ART-R, MLEM-R, Landweber-R and MART-R. Simulative studies were conducted on the phantoms shown in Fig. 4. The reconstruction domain was discretized into 40×40 grids; two orthogonal projections and ten H_{2}O absorption spectral lines were assumed; uniform noise with a level of 2% was added to the projections; and total variation was selected as the secondary criterion. According to the method detailed in Section 3.2, ${\beta _{T,0}}$ was determined to be 3.6×10^{5} and ${\beta _{C,0}}$ to be 1.8×10^{−5} for both ART-R and Landweber-R. For both MLEM-R and MART-R, ${\beta _{T,0}}$ was set to 4.5×10^{5}, while ${\beta _{C,0}}$ was set to 2.5×10^{−3} for MLEM-R and 5.0 for MART-R, respectively. The corresponding reconstructions are shown in Fig. 6.

According to the termination criteria, ART-R and Landweber-R converged after about 1100 and 1800 iterations respectively, while the total number of iterations for MLEM-R and MART-R were 2400 and 2900, respectively. It can be seen from the figure that the regularized versions had good performance in reconstructing the temperature field. The reconstruction errors were all smaller than 5% and Landweber-R had the smallest error. On the other hand, for the reconstruction of concentration, MLEM-R and MART-R had obvious advantage over ART-R and Landweber-R, and their reconstructions are more visually similar to the phantom than those of ART-R and Landweber-R.

Besides, the relationship between the reconstruction error and the noise level has been illustrated in Fig. 7. As can be seen, the performance of all algorithms were better at low noise levels and deteriorated as the noise level rose. For the reconstruction of temperature, Landweber-R had the smallest error when the noise level was below 4.5%, but its error increased rapidly above that; MART-R had a consistent increasing trend at a relatively steady pace; and ART-R and MLEM-R outperformed Landweber-R and MART-R when the noise level was higher than 5%. For the reconstruction of concentration, MLEM-R and MART-R were consistently better than Landweber-R and ART-R under various noise conditions; MLEM-R had similar performance with MART-R when the noise level was below 4.5% and better performance when the noise level was above 4.5%. Considering the fact that only MLEM-R and MART-R had good performance in reconstructing the concentration distribution under various noise conditions and MLEM-R was superior to MART-T in reconstructing temperature distribution, we will only focus on MLEM-R in the following sections.

### 3.3.2 Comparison between MLEM, its superiorized and RI versions

The performance of MLEM, its superiorized version (MLEM-S) and MLEM-R were compared against each other. Some of the numerical settings such as the phantom, the discretization, the number and arrangement of projections, and the number of spectral lines were the same as those used in Section 3.3.1. 2% uniform noise was added to the projections. MLEM converged within 20 iterations, MLEM-S around 1900, and MLEM-R about 2400. The corresponding reconstructions and the maps of relative errors are shown in Fig. 8. It is shown that MLEM-R had the smallest reconstruction error for both temperature and concentration. The reconstructions of MLEM were over-smoothed, largely deviating from the phantoms; and MLEM-S roughly captured the features of the phantoms but suffered from severe artefacts on the edges. On the contrast, MLEM-R successfully recovered the flatness of both the central and the surrounding areas.

Figure 9 plots the reconstruction error of temperature as a function of noise level for the three methods. The simulation conditions were the same as those used in Fig. 8 except the noise level. As seen, MLEM was insensitive to noise but stabilized at a relative high error. MLEM-S had better performance than MLEM at low noise levels, but its error increased rapidly as noise level rose. The regularized version shows the best performance under all noise levels, suggesting its superior noise immunity.

### 3.3.3 Effect of spectral line numbers

The reconstruction accuracy of MAT depends dramatically on the number and choice of spectral lines. Here, we compare the performance of MLEM, MLEM-S and MLEM-R with different numbers of spectral lines. The spectral lines were chosen randomly from the ten lines listed in Table 2. The results are shown in Fig. 10.

Theoretically, the more spectral lines are used, the richer information can be incorporated into the inversion, resulting in higher reconstruction accuracy. Nevertheless, adoption of more spectral lines would require more tunable diode lasers or a laser that can span a larger spectral range, leading to both higher experimental and computational costs. Hence, it should be investigated what is the optimal number of spectral lines. The reconstruction error of temperature as a function of number of spectral lines for the three algorithms are shown in Fig. 10(a). The reconstruction error of MLEM only decreased slightly as spectral lines increased. According to the literature, the usage of multispectral information would mainly enhance the noise immunity of the algorithm [3,20,21]. However, as suggested by Fig. 9(a), MLEM was almost insensitive to noise when reconstructing the temperature field. Thus, the benefit of using multiple absorption lines was minimized. On the contrary, MLEM-S was more sensitive to noise. As a result, the reconstruction improved as more spectral lines were used. However, the improvement diminished when more than five spectral lines were used. The behavior of MLEM-R was similar but the error stabilized at a larger number of spectral lines. It is worth noting that the regression and iteration processes were alternately performed in MLEM-R. Thus, the benefit of using multispectral information in the regression process would also be helpful to improve the reconstruction process. This partly contributes to its superiority over the other algorithms.

As can be seen from Fig. 10(b), the computational cost of MLEM was far less compared with MLEM-S and MLEM-R. This was because the repetitive calculation of priors was not required in MLEM. On the other hand, MLEM-R is more computationally efficient than MLEM-S. This can be explained by the fact that MLEM-R only needed to calculate the priors of temperature and concentration within each iteration but MLEM-S had to calculate the priors of the absorption coefficient at all spectral lines. Thus, when two spectral lines were used, the computational costs of MLEM-S and MLEM-R were similar; and more time was required for MELM-S when more than two spectral lines were used. The overall computational cost of MLEM-S had a linear relationship with the number of spectral lines. This can be explained by the fact that the numbers of equations to be solved as well as the priors to be calculated in the reconstruction process and the number of absorption coefficients used in the regression process of each pixel scale linearly with the number of spectral lines. On the contrary, the computational cost of MLEM-R had a positive correlation with number of spectral lines used but not linearly. This is because within each iteration, the number (*i.e*. two) of priors that need to be calculated did not dependent on the number of spectral lines.

### 3.3.4 Different secondary criteria

Except for total variation, other secondary criteria can also be applied in MLEM-R, for example, smoothness. Two phantoms with different characteristics were used to test the suitability of different secondary criteria. One was the McKenna flame, which was piecewise flat; and the other was a smooth phantom with two Gaussian peaks. MLEM-R was applied to reconstruct these two phantoms with different secondary criteria, including total variation (denoted as TV) and smoothness (denoted as Tik). The phantoms and the corresponding reconstructions are shown in Fig. 11. To render more obvious distinctions between the results of different cases and priors, four projections were used here instead of two projections; the orientations of the projections were uniformly arranged within a span of 180 degrees, and the beams in each projection were uniformly distributed. Other conditions were kept the same as in previous simulation cases.

It can be seen from Fig. 11 that the selection of the secondary criterion significantly influences the reconstructions. The piecewise flat phantom was recovered better with total variation, while the smooth phantom with the Tikhonov regularization as the secondary criterion. Therefore, a secondary criterion that is consistent with the target flame field should be selected for the reconstruction. It can be noticed that the reconstruction of Case 1 with total variation as the secondary criterion is not completely piecewise smooth, this was because the optimization of the secondary criterion and the iteration which enforces the constraint were performed alternately, the global optimum was not necessarily guaranteed when the algorithm was stopped. Also, Case 2 with total variation as the secondary is not piecewise smooth at all, this is mainly due to the constraints imposed by the measurement data *i.e.* the projections.

## 4. Experimental studies

#### 4.1 Experimental details

For a demonstrative purpose, we process the experimental data that has been published by Li *et al.* [56] with the algorithms introduced here. Since this article was published in Chinese, for the readers’ convenience, the details of the experiments are briefly summarized here. As shown in Fig. 12, the experimental setup consisted of a McKenna burner, a signal generator (RIGOL DG1062z), laser drivers (Stanford LDC-501), DFB lasers (NEL NLK1E5GAAA), a beam splitter, photoelectric detectors (THORLABS PDA20CS-EC), etalons (PMMZI-1400-300 MHz/ PMMZI-1343-300 MHz), multimode fibers, lenses and a data acquisition module (NI PXI-6110).

The burner was fed with CH_{4}/air and produced a flame which had a flat circular central area with both a high temperature and H_{2}O concentration as well as a flat periphery area with both lower temperature and concentration. The transition between the two regions was rather sharp. The burner was installed on a stepping motor platform with the tomographic system positioned 1 cm above the burner. The target flame was scanned by translating the platform with a step size of 4 mm along two orthogonal orientations. The external diameter of the burner was 76 mm and the ROI was set as an 80 mm×80 mm square and divided into 20×20 grids.

During the experiment, the signal generator would produce a saw-tooth wave to modulate the laser driver. At first, driving current was lower than the threshold and no laser would be produced. The signals detected during this period were recorded as the background noise. Afterwards, as the driving current increased, the four DFB lasers generated laser beams with center wavelengths of 1392 nm, 1343 nm, 1349 nm and 1342 nm, respectively. Each laser beam was then divided into two arms by the beam splitter and used for measurement and reference, respectively. Fiber-coupled detectors were used to record the laser intensity at both arms after passing the etalon and ROI respectively. Finally, the detected signals were acquired by the data acquisition (DAQ) module and then post-processed.

To keep the flame stable during the experiment, the flowrates of CH_{4} and air were kept constant. High speed shroud N_{2} flow was used to eliminate the interference of the ambient air. The equivalence ratio was kept 1 and the flowrate of CH_{4} was kept 1.708 L/min, while the flowrate of air was kept 10.27 L/min. The flowrate of shroud N_{2} flow was set as 15 L/min. The scanning frequency of the signal generator was 100 Hz, and the sampling frequency was 1 MHz. The data was averaged every 50 periods to increase the signal-to-noise ratio.

Ten absorption transitions of H_{2}O were chosen from the HITRAN database [38], which can be covered by the four laser diodes (LD) and are listed in Table 2.

The second column “LD index” indicates a spectral line that is covered by which laser diode. According to Beer-Lambert law and the principle of linear MAT, the properties of spectral lines used in the studies will significantly influence the results. Thus, the spectral lines were selected according to the following principles: (1) to increase the signal-to-noise ratio, intensities of the absorbance at the selected wavelengths should be strong enough in the temperature range of the ROI; (2) to ensure high measurement sensitivity, integrated absorbance of the selected lines should be strongly dependent on temperature in the target range; and (3) there is no interference from the absorption of other species.

#### 4.2 Results and discussion

After acquiring the transmitted laser signals, the baseline was fitted with the data points in the non-absorbing region. The time-frequency relationship for each laser was calibrated using an etalon. Later, the absorption spectra could be obtained for each LOS measurements. Voigt fit [57] was then applied and the integrated absorbance for each LOS measurements could be calculated as shown in Fig. 13. Finally, the data was used in the reconstruction of temperature and concentration distributions.

Figure 14(a) shows the distributions of temperature and concentration reconstructed by MLEM-R, in which the total variation was chosen as the target function. As a comparison, the reconstructions from MLEM are presented in Fig. 14(b). It can be seen that the results from MLEM-R have a distinct boundary between low temperature and high temperature regions, which conforms to the features of McKenna flame. On the contrary, the temperature field reconstructed from MLEM exhibits a smoother distribution, indicating a large deviation from the reality; besides, the over-smoothed edge also reflects the inaccuracy of the reconstructions. The distribution of mole fraction should also have a flat central area as well as a flat periphery area but a sharp transition in between. However, the reconstruction by MLEM indicates a much higher concentration on the corners which is considered to be artefacts. In contrast, only minor artefacts can be seen in the reconstruction of MLEM-R, which is mainly due to the incorporation of the smoothness prior of concentration.

To quantify the reconstruction accuracy, B type thermocouples were used to measure the flame temperature profile as a reference. The temperature at 38 uniformly distributed points were measured along the central line y=40 mm on the plane where the height above the burner (HAB) was 1 cm. Considering heat radiation and heat transfer equilibrium, the measurements from thermocouples had been corrected with the method developed by McEnally *et al.* [58]. The measured temperature profile as well as those reconstructed by MLEM, MLEM-S and MLEM-R were all plotted in Fig. 15.

As is shown in Fig. 15, both measured and reconstructed profiles have indicated a high temperature region in the center and shape declining edges. Except MLEM, all other algorithms have generated a flat central region that is consistent with the measurement of thermocouples.

Absolute and relative discrepancies between the temperature reconstructions and the thermocouple data were shown in Fig. 16. It is reasonable that the discrepancy is large at the two flanks, since the B-type thermocouple has a larger uncertainty at lower temperatures. Apart from this, the difference is also partly attributed to inaccurate locating of thermocouples and laser beams, since temperature changes dramatically in this region and a slight misalignment will cause huge difference. Besides, the discretization of ROI is a bit coarse to reflect the temperature change in the flank region. Thus, another indicator, the average relatively discrepancy at the middle ten points, was shown in Fig. 16(b). The definition of relative discrepancy in this region is the same as the relative error defined in Eq. (19). Although MLEM has a smallest overall discrepancy, too large error in the central region implies the failure of this method. Except for MLEM, MLEM-R has the smallest discrepancy with the thermocouple data both for the entire region and the central region. Hence, it can be concluded that MLEM-R has the best performance in solving such an MAT problem.

## 5. Conclusions

This work demonstrated regularization of iterative algorithms for the reconstruction of two-dimensional temperature and H_{2}O concentration fields based on linear MAT. Simulative studies were conducted to evaluate the performance of classical iterative algorithms and the corresponding superiorized and regularized versions. It was demonstrated that the regularized versions outperform the others in terms of the reconstruction accuracy. The regularized version of MLEM was advised for modest noise conditions, and the priors should be selected according to the characteristics of the flame fields as the secondary criterion. The experimental demonstration performed on a McKenna burner with CH_{4}/air premixed flame confirmed the superior reconstruction accuracy of the regularized MLEM. Finally, we would like to point out that the linear MAT problem can also be modeled as an optimization problem and use the priors of the distributions of absorption coefficients as the regularization. The comparison of this approach with the regularized iteration would be one of our future works.

## Funding

National Science and Technology Major Project (2017-III-0007-0033); National Natural Science Foundation of China (51976122, 52061135108).

## Acknowledgement

The authors would like to acknowledge Professor Yair Censor from University of Haifa and Professor Gabor T. Herman from City University of New York, for the useful discussion on superiorization.

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are available in [56].

## References

**1. **Z. H. Wang, P. F. Fu, and X. Chao, “Laser Absorption Sensing Systems: Challenges, Modeling, and Design Optimization,” Appl. Sci. **10**(1), 27 (2019). [CrossRef]

**2. **X. Chao, G. F. Shen, K. Sun, Z. H. Wang, Q. H. Meng, S. K. Wang, and R. K. Hanson, “Cavity-enhanced absorption spectroscopy for shocktubes: Design and optimization,” Proc. Combust. Inst. **37**(2), 1345–1353 (2019). [CrossRef]

**3. **C. Liu and L. J. Xu, “Laser absorption spectroscopy for combustion diagnosis in reactive flows: A review,” Appl. Spectrosc. Rev. **54**(1), 1–44 (2019). [CrossRef]

**4. **G. Mathews and C. Goldenstein, “Near-GHz scanned-wavelength-modulation spectroscopy for MHz thermometry and H2O measurements in aluminized fireballs of energetic materials,” Appl. Phys. B **126**(11), 189 (2020). [CrossRef]

**5. **L. H. Ma, Z. Wang, K. P. Cheong, H. B. Ning, and W. Ren, “Mid-infrared heterodyne phase-sensitive dispersion spectroscopy in flame measurements,” Proc. Combust. Inst. **37**(2), 1329–1336 (2019). [CrossRef]

**6. **A. S. Makowiecki, J. E. Steinbrenner, N. T. Wimer, J. F. Glusman, C. B. Lapointe, J. W. Daily, P. E. Hamlington, and G. B. Rieker, “Dual frequency comb spectroscopy of solid fuel pyrolysis and combustion: Quantifying the influence of moisture content in Douglas fir,” Fire Safety J. **116**, 103185 (2020). [CrossRef]

**7. **S. O’Hagan, J. H. Northern, B. Gras, P. Ewart, C. S. Kim, M. Kim, C. D. Merritt, W. W. Bewley, C. L. Canedy, I. Vurgaftman, and J. R. Meyer, “Multi-species sensing using multi-mode absorption spectroscopy with mid-infrared interband cascade lasers,” Appl. Phys. B **122**(1), 11 (2016). [CrossRef]

**8. **G. Zhao, T. Hausmaninger, W. G. Ma, and O. Axner, “Shot-noise-limited Doppler-broadened noise-immune cavity-enhanced optical heterodyne molecular spectrometry,” Opt. Lett. **43**(4), 715–718 (2018). [CrossRef]

**9. **P. J. Schroeder, R. J. Wright, S. Coburn, B. Sodergren, K. C. Cossel, S. Droste, G. W. Truong, E. Baumann, F. R. Giorgetta, I. Coddington, N. R. Newbury, and G. B. Rieker, “Dual frequency comb laser absorption spectroscopy in a 16 MW gas turbine exhaust,” Proc. Combust. Inst. **36**(3), 4565–4573 (2017). [CrossRef]

**10. **A. P. Nair, D. D. Lee, D. I. Pineda, J. Kriesel, W. A. Hargus, J. W. Bennewitz, S. A. Danczyk, and R. M. Spearrin, “MHz laser absorption spectroscopy via diplexed RF modulation for pressure, temperature, and species in rotating detonation rocket flows,” Appl. Phys. B **126**(1), 20 (2020). [CrossRef]

**11. **G. B. Rieker, F. R. Giorgetta, W. C. Swann, J. Kofler, A. M. Zolot, L. C. Sinclair, E. Baumann, C. Cromer, G. Petron, C. Sweeney, P. P. Tans, I. Coddington, and N. R. Newbury, “Frequency-comb-based remote sensing of greenhouse gases over kilometer air paths,” Optica **1**(5), 290–298 (2014). [CrossRef]

**12. **J. Shao, J. Xiang, O. Axner, and C. Ying, “Wavelength-modulated tunable diode-laser absorption spectrometry for real-time monitoring of microbial growth,” Appl. Opt. **55**(9), 2339–2345 (2016). [CrossRef]

**13. **W. W. Cai and C. F. Kaminski, “Tomographic absorption spectroscopy for the study of gas dynamics and reactive flows,” Progress in Energy and Combustion Sci. **59**, 1–31 (2017). [CrossRef]

**14. **C. Liu, L. J. Xu, and Z. Cao, “Measurement of nonuniform temperature and concentration distributions by combining line-of-sight tunable diode laser absorption spectroscopy with regularization methods,” Appl. Opt. **52**(20), 4827–4842 (2013). [CrossRef]

**15. **L. Ma, X. S. Li, S. T. Sanders, A. W. Caswell, S. Roy, D. H. Plemmons, and J. R. Gord, “50-kHz-rate 2D imaging of temperature and H2O concentration at the exhaust plane of a J85 engine using hyperspectral tomography,” Opt. Express **21**(1), 1152–1162 (2013). [CrossRef]

**16. **Z. C. Qu, R. Ghorbani, D. Valiev, and F. M. Schmidt, “Calibration-free scanned wavelength modulation spectroscopy-application to H2O and temperature sensing in flames,” Opt. Express **23**(12), 16492–16499 (2015). [CrossRef]

**17. **T. Yu, B. Tian, and W. W. Cai, “Development of a beam optimization method for absorption-based tomography,” Opt. Express **25**(6), 5982–5999 (2017). [CrossRef]

**18. **F. Stritzke, S. van der Kley, A. Feiling, A. Dreizler, and S. Wagner, “Ammonia concentration distribution measurements in the exhaust of a heavy duty diesel engine based on limited data absorption tomography,” Opt. Express **25**(7), 8180–8191 (2017). [CrossRef]

**19. **J. Foo and P. A. Martin, “Tomographic imaging of reacting flows in 3D by laser absorption spectroscopy,” Appl. Phys. B **123**(5), 160 (2017). [CrossRef]

**20. **W. W. Cai, D. J. Ewing, and L. Ma, “Application of simulated annealing for multispectral tomography,” Comput. Phys. Commun. **179**(4), 250–255 (2008). [CrossRef]

**21. **J. Q. Huang, H. C. Liu, J. H. Dai, and W. W. Cai, “Reconstruction for limited-data nonlinear tomographic absorption spectroscopy via deep learning,” J. Quant. Spectrosc. Radiat. Transfer **218**, 187–193 (2018). [CrossRef]

**22. **D. Verhoeven, “Limited-data computed tomography algorithms for the physical sciences,” Appl. Opt. **32**(20), 3736–3754 (1993). [CrossRef]

**23. **H. H. Xia, R. F. Kan, Z. Y. Xu, Y. B. He, J. G. Liu, B. Chen, C. G. Yang, L. Yao, M. Wei, and G. L. Zhang, “Two-step tomographic reconstructions of temperature and species concentration in a flame based on laser absorption measurements with a rotation platform,” Opt. Lasers Eng. **90**, 10–18 (2017). [CrossRef]

**24. **X. L. An, A. W. Caswell, J. J. Lipor, and S. T. Sanders, “Determining the optimum wavelength pairs to use for molecular absorption thermometry based on the continuous-spectral lower-state energy,” J. Quant. Spectrosc. Radiat. Transfer **112**(14), 2355–2362 (2011). [CrossRef]

**25. **S. J. Grauer, J. Emmert, S. T. Sanders, S. Wagner, and K. J. Daun, “Multiparameter gas sensing with linear hyperspectral absorption tomography,” Meas. Sci. Technol. **30**(10), 105401 (2019). [CrossRef]

**26. **R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography,” J. Theoretical Biol. **29**(3), 471–481 (1970). [CrossRef]

**27. **A. P. Dempster, “Maximum likelihood from incomplete data via the EM algorithm,” Journal of the Royal Statal Society **39**(1), 1–38 (1977).

**28. **L. Ma and W. W. Cai, “Determination of the optimal regularization parameters in hyperspectral tomography,” Appl. Opt. **47**(23), 4186–4192 (2008). [CrossRef]

**29. **J. H. Dai, T. Yu, L. J. Xu, and W. W. Cai, “On the regularization for nonlinear tomographic absorption spectroscopy,” J. Quant. Spectrosc. Radiat. Transfer **206**, 233–241 (2018). [CrossRef]

**30. **L. I. Rudin and S. Osher, “Total variation based image restoration with free local constraints,” Proceedings ICIP-94 (Cat. No.94CH35708), 31–35 vol. 31 (1994).

**31. **K. J. Daun, S. J. Grauer, and P. J. Hadwin, “Chemical species tomography of turbulent flows: Discrete ill-posed and rank deficient problems and the use of prior information,” J. Quant. Spectrosc. Radiat. Transfer **172**, 58–74 (2016). [CrossRef]

**32. **D. Butnariu, R. Davidi, G. T. Herman, and I. G. Kazantsev, “Stable Convergence Behavior Under Summable Perturbations of a Class of Projection Methods for Convex Feasibility and Optimization Problems,” IEEE J. Sel. Top. Signal Process. **1**(4), 540–547 (2007). [CrossRef]

**33. **R. Davidi, G. T. Herman, and Y. Censor, “Perturbation-resilient block-iterative projection methods with application to image reconstruction from projections,” Int. Trans. Oper. Res. **16**, 505–524 (2009). [CrossRef]

**34. **C. Y. Shui, J. Q. Huang, and W. W. Cai, “On the superiorization of inversion algorithms for tomographic absorption spectroscopy,” Phys. Gases **5**, 28–37 (2020).

**35. **W. W. Cai and C. F. Kaminski, “A tomographic technique for the simultaneous imaging of temperature, chemical species, and pressure in reactive flows using absorption spectroscopy with frequency-agile lasers,” Appl. Phys. Lett. **104**(3), 034101 (2014).10.1063/1.4862754

**36. **W. W. Cai and C. F. Kaminski, “Multiplexed absorption tomography with calibration-free wavelength modulation spectroscopy,” Appl. Phys. Lett. **104**(15), 154106 (2014). [CrossRef]

**37. **L. Ma, X. S. Li, W. W. Cai, S. Roy, J. R. Gord, and S. T. Sanders, “Selection of Multiple Optimal Absorption Transitions for Nonuniform Temperature Sensing,” Appl. Spectrosc. **64**(11), 1274–1282 (2010). [CrossRef]

**38. **L. S. Rothman, I. E. Gordon, Y. Babikov, A. Barbe, D. C. Benner, P. F. Bernath, M. Birk, L. Bizzocchi, V. Boudon, L. R. Brown, A. Campargue, K. Chance, E. A. Cohen, L. H. Coudert, V. M. Devi, B. J. Drouin, A. Fayt, J. M. Flaud, R. R. Gamache, J. J. Harrison, J. M. Hartmann, C. Hill, J. T. Hodges, D. Jacquemart, A. Jolly, J. Lamouroux, R. J. Le Roy, G. Li, D. A. Long, O. M. Lyulin, C. J. Mackie, S. T. Massie, S. Mikhailenko, H. S. P. Muller, O. V. Naumenko, A. V. Nikitin, J. Orphal, V. Perevalov, A. Perrin, E. R. Polovtseva, C. Richard, M. A. H. Smith, E. Starikova, K. Sung, S. Tashkun, J. Tennyson, G. C. Toon, V. G. Tyuterev, and G. Wagner, “The HITRAN2012 molecular spectroscopic database,” J. Quant. Spectrosc. Radiat. Transfer **130**, 4–50 (2013). [CrossRef]

**39. **M. G. Allen, “Diode laser absorption sensors for gas-dynamic and combustion flows,” Meas. Sci. Technol. **9**(4), 545–562 (1998). [CrossRef]

**40. **S. Kaczmarz, “Angenäherte Auflösung von Systemen linearer Gleichungen,” Bull. Intern. Acad. Polonaise Sci. Lett **35**, 355–357 (1937).

**41. **Y. Censor and A. Lent, “An iterative row-action method for interval convex programming,” J. Optim. Theory Appl. **34**(3), 321–353 (1981). [CrossRef]

**42. **G. T. Herman, A. Lent, and P. H. Lutz, “Relaxation methods for image reconstruction,” Commun. ACM **21**(2), 152–158 (1978). [CrossRef]

**43. **K. Lange and R. Carson, “EM reconstruction algorithms for emission and transmission tomography,” J. Comput. Assist. Tomogr. **8**, 306–316 (1984).

**44. **L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Trans. Med. Imaging **1**(2), 113–122 (1982). [CrossRef]

**45. **Y. Censor, R. Davidi, and G. T. Herman, “Perturbation resilience and superiorization of iterative algorithms,” Inverse Probl. **26**(6), 065008 (2010). [CrossRef]

**46. **E. Garduno and G. T. Herman, “Superiorization of the ML-EM Algorithm,” IEEE Trans. Nucl. Sci. **61**(1), 162–172 (2014). [CrossRef]

**47. **K. J. Daun, “Infrared Species Limited Data Tomography through Tikhonov Reconstruction,” Ht2009: Proceedings of the Asme Summer Heat Transfer Conference 2009, Vol 1, 187–196 (2009).

**48. **S. S. Luo and T. Zhou, “Superiorization of EM algorithm and its application in single-photon emission computed tomography(SPECT),” Inverse Probl. Imaging **8**(1), 223–246 (2014). [CrossRef]

**49. **T. Nikazad, R. Davidi, and G. T. Herman, “Accelerated perturbation-resilient block-iterative projection methods with application to image reconstruction,” Inverse Probl. **28**(3), 035005 (2012). [CrossRef]

**50. **Y. Censor and A. J. Zaslavski, “Convergence and perturbation resilience of dynamic string-averaging projection methods,” Comput. Optim. Appl. **54**(1), 65–76 (2013). [CrossRef]

**51. **W. M. Jin, Y. Censor, and M. Jiang, “Bounded perturbation resilience of projected scaled gradient methods,” Comput. Optim. Appl. **63**(2), 365–392 (2016). [CrossRef]

**52. **D. Verhoeven, “Multiplicative algebraic computed tomographic algorithms for the reconstruction of multidirectional interferometric data,” Opt. Eng. **32**(2), 410–419 (1993). [CrossRef]

**53. **L. Landweber, “An iteration formula for Fredholm integral equations of the first kind,” Am. J. Math. **73**(3), 615–624 (1951). [CrossRef]

**54. **Y. Censor, “Weak and Strong Superiorization: Between Feasibility-Seeking and Minimization,” Analele Stiint. Univ. Ovidius C. **23**, 41–54 (2015).

**55. **W. W. Cai, X. S. Li, F. Li, and L. Ma, “Numerical and experimental validation of a three-dimensional combustion diagnostic based on tomographic chemiluminescence,” Opt. Express **21**(6), 7050–7064 (2013). [CrossRef]

**56. **K. Li, B. Zhou, J. Z. Jia, Y. Du, Y. Lu, H. Y. Cheng, S. M. Wang, and H. Wang, “Alternating Iterative Optimization Algorithm Based on Multispectral Absorption,” J. Eng. Thermophysics **39**, 1598–1608 (2018).

**57. **J. J. Olivero and R. L. Longbothum, “Empirical fits to the Voigt line width: a brief review,” J. Quant. Spectrosc. Radiat. Transfer **17**(2), 233–236 (1977). [CrossRef]

**58. **C. S. McEnally, U. O. Koylu, L. D. Pfefferle, and D. E. Rosner, “Soot volume fraction and temperature measurements in laminar nonpremixed flames using thermocouples,” Combustion and Flame **109**(4), 701–720 (1997). [CrossRef]