A Fast Tomographic Reconstruction Method for Flame Temperature Distribution Measurement Based on Direct Solution Algorithm

The rapid and accurate measurement of the flame temperature distribution is of great significance to the structural design and health diagnosis of the engine. Aiming at the low reconstruction efficiency of traditional flame temperature distribution reconstruction algorithms, a Direct Solution algorithm for flame temperature distribution reconstruction is proposed in this paper based on the structural characteristics of the reconstruction equations. By setting several numerical cases, the performance of the Direct Solution algorithm and some commonly used traditional algorithms, such as Simultaneous Algebraic Reconstruction Technique (SART), Least Squares QR-factorization (LSQR) algorithm, Non-Negative Least Squares QR-factorization (NNLS) algorithm, is compared in the reconstruction of the flame temperature distribution. The results show that the efficiency of the Direct Solution method is 169.4, 7.4, and 3483.3 times higher than that of the SART, LSQR, and NNLS algorithms under the condition of 40 × 40 grids. In addition, with the increase of the number of grids, the growth rate of the reconstruction time of the Direct Solution algorithm is much lower than that of other algorithms. The overall reconstruction accuracy of the Direct Solution algorithm is better than that of SART and LSQR algorithms. This shows that it has an excellent comprehensive performance and has a great application prospect in the rapid reconstruction of the temperature distribution.


INTRODUCTION
Thrust-to-weight ratio is a comprehensive index to measure the performance level and working ability of aero engines. One of the main ways to improve the thrust-to-weight ratio of modern aero engines is to increase the outlet temperature of the combustion chamber . As the temperature rises, turbine blades and vanes need to bear a greater thermal load. The huge thermal gradient causes severe thermal stress and strain on the blades, which greatly reduces the creep life of the blades, and leads to ablation, fracture and failure of the blades (Marahleh et al., 2006). Real-time monitoring and obtaining accurate combustion chamber outlet temperature can respond to local temperature sudden changes in the engine in a timely manner, which is of great significance for the safety of aero-engine operation and prolonged service life (Li et al., 2006;He et al., 2019;Gamil et al., 2020).
Due to the shortcomings of thermocouple and other contact temperature measurement, such as slow dynamic response, interference with combustion field, narrow temperature measurement range, only point measurement, hightemperature oxidation, etc., it is not suitable for cutting-edge combustion chambers under high temperature, high pressure, and high-speed conditions (Yang et al., 2015;Jenkins et al., 2020;Zhang et al., 2016;Lian et al., 2017). In recent years, the flame temperature distribution reconstruction technology based on radiation images has received extensive attention from researchers all over the world (Neuer et al., 2001;Sun et al., 2018;Qi et al., 2019;Zhao et al., 2017;Hossain et al., 2013;Niu et al., 2015). It has the advantages of non-invasiveness and no need for external excitation sources and is more suitable for the measurement of the exit temperature distribution of the combustion chamber of aero engines.
The radiation image method for temperature measurement utilizes flame radiation information from different directions to reconstruct the flame temperature field, and the use of a reasonable reconstruction algorithm is beneficial to improve the efficiency and accuracy of temperature distribution reconstruction. Generally, the number of grids in the reconstruction area will affect the efficiency of the reconstruction algorithm. As the number of grids increases, the difficulty and time to solve the ill-conditioned equations transformed from the inverse radiation problem tend to increase exponentially, which is difficult to meet the current demand for rapid reconstruction of the exit temperature of the combustion chamber of aero engines.
To solve the above-mentioned problems, Zhou et al. (2002) used an improved Tikhonov regularization algorithm to reconstruct a three-dimensional temperature distribution of N x × N y × N z 10 × 10 × 10 from the flame image obtained by the CCD camera, which proved the robustness of the algorithm, but the reconstruction efficiency is not studied. Liu et al. (2010) used the Least Square QR decomposition (LSQR) algorithm to achieve a two-dimensional temperature distribution reconstruction with a discrete area of N x × N y 10 × 10 based on an 8-camera system combined with the Monte Carlo method. However, if the LSQR algorithm is directly used for large-scale illconditioned problems, there will be a "semi-convergence" problem, which will affect the accuracy and reconstruction time. Qiu et al. (2014) combined the Tikhonov regularization algorithm and generalized truncated singular value decomposition hybrid algorithm (TR-GSVD) to reconstruct the temperature distribution in the discrete area of N x × N y × N z 10 × 10 × 12 in the furnace, and the reconstruction time is about 2s. Cai et al. (2013) reconstructed the flame temperature distribution with a grid number of N x × N y × N z 54 × 54 × 10 based on chemiluminescence tomography technology, using Algebraic Reconstruction Technique (ART) and minimization algorithm combined with regularization technology. Numerical simulations and experiments have verified the effectiveness of the method. However, the low efficiency of the "line-by-line correction" of the ART algorithm will affect the reconstruction time. Huang et al. (2018) used the improved Landweber algorithm to reconstruct the three-dimensional temperature distribution of the absorbing and scattering flame based on the light field imaging technology. But it still takes 2.5 s to reconstruct the temperature distribution with the number of grids N φ × N r × N l 1 × 20 × 20. The above-mentioned traditional algorithms have achieved certain results in the reconstruction of the temperature distribution, but the calculation time for the reconstruction of the temperature distribution with a large number of grids is generally too longer. Therefore, it is necessary to develop a flame temperature distribution reconstruction algorithm that greatly improves the reconstruction efficiency under the premise of ensuring accuracy and is less affected by the increase in the number of grids.
In this paper, based on the CCD camera imaging and flame radiative transfer model, the rapid reconstruction method of absorbing and scattering flame temperature distribution is carried out. It takes the generalized source term finite volume method as the model of the forward problem and the proposed direct solution algorithm as the inverse method for reconstructing the temperature distribution. Two numerical cases of flame temperature distribution reconstruction are designed to illustrate the performance of the proposed direct solution algorithm compared with the current main tomographic algorithms, such as Simultaneous Algebraic Reconstruction Technique (SART), LSQR, and Non-Negative Least Squares QR-factorization (NNLS) algorithm.

Forward Model
In the radiation image method for temperature measurement, as shown in Figure 1, an ideal thin lens model can be used to describe the imaging process of light from the object space passing through the lens to the camera detector pixels. The light entering the camera is characterized by the line between the center of the CCD pixel and the optical center of the lens. Considering that the flame is a semi-transparent object, the light will pass through the flame. When the deflection effect of the flame refractive index gradient on the light is not considered, the radiation of the flame voxel through which the light passes all contributes to the gray value of the pixel. Therefore, it is necessary to record the number of flame voxels that each detection light passes through, and the distance traveled in the voxel. According to the generalized source term finite volume method (Zhang et al., 2017), the radiation energy received on each pixel can be calculated by the following formula: where I(r p , Ω) is the radiation intensity along the direction Ω at the position r p . n is the total number of the flame voxel that the light passes through. S(r i , Ω) is the generalized source term of the n th voxel along the direction Ω. β is the extinction coefficient of the flame. Δs is the distance that the light travels through the flame voxel, witch can be determined by the reverse tracing of the light.
Since the propagation speed of light is much greater than the equilibrium speed of temperature, the radiative transfer of flame in radiation image method for temperature measurement can be expressed by the steady-state radiative transfer equation : where I(r, Ω) is the radiative intensity along the direction Ω at the flame voxel r. I b (r) is the black body radiative intensity at the flame voxel r. ω is scattering albedo of the flame. Φ(Ω′, Ω) is the scattering phase function, where Ω′ is the incident direction and Ω is the scattering direction. Ω is the solid angle.
The second and third terms on the right side of Eq.
(2) represent emission-enhanced and scattering-enhanced terms, respectively. The sum of the second and third terms can be defined as the generalized source term.
Each detected light of the camera can establish as an equation using Eq. (1), then all detected light of the camera can form reconstruction equations. The generalized source term distribution can be obtained through the inversion calculation of the reconstruction equations. After substituting the generalized source term, the finite volume method can be used to calculate the black body radiative intensity of each flame voxel, and then the temperature of each flame voxel can be obtained by Planck's law.
where A is the coefficient matrix of the reconstruction equations, which is related to the flame extinction coefficient and the optical path of the light passing through the voxels. S is the generalized source term vector. I is the flame radiative intensity distribution received by the CCD camera. M and N are the total number of detected lights and flame voxels.

Inversion Principle
The dimension of the coefficient matrix A in the reconstruction equations is determined by the dimensions of the unknowns vector S and the measured value vector I. Because the current camera resolution is relatively large (usually millions), resulting in a huge dimension of the coefficient matrix. This leads to low solution efficiency. However, the solution process can be optimized by analyzing the characteristics of the coefficient matrix. In order to facilitate the demonstration, a twodimensional flame and linear CCD are used to analyze the flame temperature distribution reconstruction in this paper, as shown in Figure 2.
According to Eq. (1), if the light passes through the flame voxel, the elements in the corresponding coefficient matrix are non-zero. If non-zero elements of the coefficient matrix are represented by blue dots, and zero elements of the coefficient matrix are represented by blanks, the coefficient matrix can be represented as Figure 3. It can be seen that the coefficient matrix A is very sparse. There are several rows of coefficient matrix containing only one non-zero element, such as rows 1, 2, / , 6, M-5, M-4, / , M. It means that a specific unknown generalized   source term value can be solved using one equation and one measured value. Due to the proximity of the voxels in space, many rows which contain two non-zero elements and one non-zero element share the same element. For example, row 6 and row 7 share a non-zero element 40. Therefore, the unknowns obtained in the previous calculation can be used to substitute back to a new row for elimination. Using this "solving-substitution" model, only a few equations need to be calculated each time to gradually complete the solution of large-scale linear equations, so that a few unknown equations can be calculated all the time to gradually complete the solution of large-scale coefficient equations. Considering that there are multiple rays of light passing through the same flame voxel at the same time, in order to reduce the influence of measurement error, multiple equations can be solved at the same time to obtain the least square solution of the same non-zero element. For example, rows 1, 2, 3, 4, 5, and 6 can form a sub-matrix to solve the 40th unknown generalized radiation source term.
Since the reconstruction method proposed in this paper is to solve a certain unknown parameter each time directly, this method is called Direct Solution algorithm. The procedure for implementing the Direct Solution algorithm is described as the following steps, and the flowchart is shown in Figure 4.
Step 1. Input the coefficient matrix A (k) and count the number of non-zero elements and corresponding positions in each row.
Step 2. According to the counted number of non-zero elements, arrange the row vectors of the coefficient matrix in ascending order to form a new coefficient matrix A (k) .
Step 3. Extract the row vectors containing the least specific non-zero element a i,j from the matrix to form the solving submatrix A ′(k) m , and the remaining part form the candidate submatrix A ′(k) n .
Step 4. Check whether the solving sub-matrix A ′(k) m is full rank. If it is not full rank, go to Step 5. If it is full rank, go to Step 6.
Step 5. Extract the row vectors containing element a i,j and the least non-zero elements from the candidate sub-matrix A '(k) n . Combine them with the solving sub-matrix A ′(k) m to form a new solving sub-matrix A ′(k) m+t and the remaining candidate sub-matrix constitutes a new candidate sub-matrix A ′(k) m−t . Then go to Step 4.
Step 6. Use the LSQR algorithm to solve the sub-matrix A ′(k) m . Substitute the calculation result back into the candidate submatrix A '(k) n . Eliminate the solved elements in the candidate submatrix to form a new coefficient matrix A (k+1) . Loop to Steps 1 until all unknowns are obtained.

RESULTS AND DISCUSSIONS Evaluation Criterion
In order to evaluate the performance of the Direct Solution algorithm proposed in this paper, several common traditional reconstruction algorithms of flame temperature distribution are selected for comparison, such as SART (Liu et al., 2017), LSQR (Liu et al., 2010), NNLS (Li et al., 2019) algorithms. In terms of reconstruction accuracy, the following criteria are used for evaluation (Resnick et al., 1985): 1) The normalized mean square distance, q; 2) The normalized average absolute distance, r; The worstcase distance, e. The specific expressions are shown in Eqs (5)-(7).
where T i,j is the value of the original value of the temperature at voxel (i, j). R i,j is the reconstructed value of the temperature at voxel (i, j). T is the average value of the temperature at voxel (i, j). The normalized average absolute distance r more sensitively reflects the situation that many points have small errors. The normalized mean square distance q more sensitively reflects the error generated by a few points. The worst-case distance e more sensitively reflects the density difference between the reconstructed value and the original value.
These traditional reconstruction algorithms have been widely used in the reconstruction of flame temperature distribution based on radiation imaging and have been proved to have high reconstruction accuracy. The processes of the temperature reconstruction in these above three algorithms are similar to that of the direct solution algorithm. And the specific principles and steps are also commonly found in other references (Andersen et al., 1984;Teng et al., 1989;Bombara et al., 2011). Thus, this paper will not repeat them here.

Camera and Flame Settings
There is a four-peak 2D flame with a size of L x × L y 60 mm × 60 mm. The flame area is discretized into N x × N y 40 × 40 grids. The temperature distribution is shown  in Figure 5. Considering isotropic scattering, the extinction coefficient of the flame is β 0.02 mm −1 , and the scattering albedo is ω 0.1. The temperature distribution function is defined as follows.
where x and y respectively represent the horizontal and vertical coordinates in the world coordinate system. L x and L y represent the size of the flame in the x and y directions. Eight linear CCD cameras are arranged on the circle 0.55 m from the flame center. As shown in Figure 6, each camera focuses on the plane passing through the flame center. A lens with a magnification ratio of 1:3 is selected. The number of pixels of each camera is 4,096, and the pixel size is 7.04 µmμm, which ensures that the field of view of each camera covers all the flame domains.

Performance Comparison
According to the arrangement of the camera system and the flame, the flame image of each camera, as shown in Figure 7, is obtained through the calculation of the radiative transfer equation and the reverse tracing of the light. The light that does not pass through any flame voxel is removed to obtain the number of effective rays is 28,304. Thus the scale of the reconstruction coefficient matrix is 28,304 × 1,600, in which the proportion of non-zero elements to the total elements is 2.46%, the rank of the coefficient matrix is 1,600, and the condition number is 526.13. It can be seen that the reconstruction equations are a typical large-scale sparse illconditioned equation system. Since the traditional general solution algorithms do not take into account the characteristics of the reconstruction coefficient matrix, they cannot give full play to their advantages.
The case was implemented using MATLAB code, and the developed program was executed on an Intel(R) Core(TM) i9-10940X CPU @ 3.30 GHz PC. The temperature distribution of the above-mentioned four-peak flame is reconstructed according to the measured flame image using SART, LSQR, NNLS, and Direct Solution algorithm, respectively. The calculation efficiency and accuracy results are shown in Table 1. It can be seen from the table that the reconstruction accuracy of the Direct Solution algorithm is better than the SART and LSQR algorithms and slightly worse than that of the NNLS algorithm in the three reconstruction quality evaluation indicators. This is because the Direct Solution algorithm uses the LSQR algorithm to calculate the small-scale equations, and the NNLS algorithm performs non-negative constraints on the calculation results based on the physical characteristics of the flame so that the method can avoid some local optimal solution that do not conform to the laws of physics.
In terms of reconstruction efficiency, the Direct Solution algorithm is 169.4 times higher than the SART algorithm, 7.4 times higher than the LSQR algorithm, and 3483.3 times higher than the NNLS algorithm, which has obvious advantages. This is because the Direct Solution algorithm gradually decomposes the equation system into multiple sub-equation systems and solves them one by one. Each sub-equation system has a few unknowns to solve, avoiding the direct calculation of large-scale equations, which can greatly reduce the calculating time.

Influence of the Grids Number
The size and temperature of the flame area are kept unchanged, and the discrete grids of the flame area is set as 10 × 10, 20 × 20, 40 × 40, and 80 × 80, as shown in Figure 8. In this case, all the camera and lens parameters are kept unchanged. The total number of the detected light is 32,768. After ray tracing, the rays that do not pass through any flame control body are eliminated, and the number of effective rays under the four grids is 28,304.
The flame temperature distribution is reconstructed using the SART algorithm, LSQR algorithm, NNLS algorithm, and Direct Solution algorithm proposed in this paper, and the overall error distribution and calculation time change of the four algorithms under different grids are obtained. The accuracy and efficiency of these algorithms are shown in Figure 9 and Figure 10. It can be seen from the figure that as the number of grids increases, the scale of the equation system expands, and the error presents an upward trend. The reconstruction accuracy of the Direct Solution algorithm is less than that of the NNLS algorithm but better than the SART and LSQR algorithms.
It can be seen from Figure 10 that as the number of grids gradually increases, the calculation time of the four algorithms increases. The calculation time of the Direct Solution algorithm increases slowest, indicating that it is least affected by the increase in the scale of the equations. The direct solution algorithm takes less time than the other three algorithms, indicating that when the number of grids increases, the direct solution algorithm takes less time to extract and solve the partial extraction of the overall equation system. This reduces the overall time-consuming to solve all the unknowns, greatly shortens the reconstruction time compared to other algorithms. The above numerical results show that the Direct Solution algorithm can quickly reconstruct the flame temperature distribution under the premise of ensuring accuracy. In occasions where real-time reconstruction is required, this algorithm has greater application prospects.

CONCLUSION
In this paper, a Direct Solution algorithm is proposed. It can quickly reconstruct the flame temperature distribution. According to the characteristics of the reconstruction coefficient matrix, it first performs statistical classification according to the minimum non-zero element criterion and the full rank criterion of the sub-equation system. And then, it solves the small-scale equations group by group according to the method of "solving-substitution", which greatly improves the reconstruction efficiency of large-scale sparse equations. The performance of the proposed method is compared with that of the traditional tomographic reconstruction algorithms. Under the condition of 40 × 40 grids, the normalized average absolute distance criterion r, normalized mean square distance criterion q, and worst-case distance criterion e of the proposed algorithm are all smaller than that of SART and LSQR algorithms. Compared with the SART, LSQR, and  Frontiers in Energy Research | www.frontiersin.org November 2021 | Volume 9 | Article 790581 NNLS algorithms, the reconstruction efficiency is increased by 169.4, 7.4, and 3483.3 times, respectively. On this basis, the influence of the grids number on the reconstruction algorithm is analyzed. The results show that as the number of grids increases, the reconstruction time of the Direct Solution algorithm increases the slowest. The overall reconstruction accuracy is higher than SART and LSQR methods, which can meet the accuracy requirements on many occasions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
BZ: conceptualization, writing -original draft, and methodology. W-JP: coding and data processing. JL: Investigation and Visualization. Z-HL: writing editing. C-LX: paper reviewing and supervision. All authors contributed to the article and approved the submitted version.