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: I ( r p , Ω ) = S ( r n , Ω ) [ 1 − exp ( − β n Δ s n ) ] + ∑ i = 1 n − 1 S ( r i , Ω ) [ exp ( − ∑ j = i + 1 n β j Δ s j ) − exp ( − ∑ j = i n β j Δ s j ) ] (1) 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 (Qi et al., 2015): dI ( r , Ω ) β ds = − I ( r , Ω ) + ( 1 − ω ) I b ( r ) + ω 4 π ∫ 4 π I ( r , Ω ′ ) Φ ( Ω ′ , Ω ) d Ω (2) 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. S ( r , Ω ) = ( 1 − ω ) I b ( r ) + ω 4 π ∫ 4 π I ( r , Ω ′ ) Φ ( Ω ′ , Ω ) d Ω (3)

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. A ⋅ S = I , A ∈ ℝ M × N , S ∈ ℝ N , I ∈ ℝ M (4) 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.

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 two-dimensional 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 sub-matrix A m ′ ( k ) , and the remaining part form the candidate sub-matrix A n ′ ( k ) .

Step 4. Check whether the solving sub-matrix A m ′ ( k ) 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 n ' ( k ) . Combine them with the solving sub-matrix A m ′ ( k ) to form a new solving sub-matrix A m+t ′ ( k ) and the remaining candidate sub-matrix constitutes a new candidate sub-matrix A m-t ′ ( k ) . Then go to Step 4.

Step 6. Use the LSQR algorithm to solve the sub-matrix A m ′ ( k ) . Substitute the calculation result back into the candidate sub-matrix A n ' ( k ) . Eliminate the solved elements in the candidate sub-matrix to form a new coefficient matrix A ( k + 1 ) . Loop to Steps 1 until all unknowns are obtained.

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 worst-case distance, e . The specific expressions are shown in Eqs (5)–(7). r = ∑ i = 1 N ∑ j = 1 N | T i , j − R i , j | / ∑ i = 1 N ∑ j = 1 N | T i , j | (5) q = ∑ i = 1 N ∑ j = 1 N ( T i , j − R i , j ) 2 / ∑ i = 1 N ∑ j = 1 N ( T i , j − T ¯ ) 2 (6) e = max 1 ≤ i ≤ [ N 2 ] 1 ≤ j ≤ [ N 2 ] | M i , j − E i , j | , { M i , j = ( T 2 i , 2 j + T 2 i + 1,2 j + T 2 i , 2 j + 1 + T 2 i + 1,2 j + 1 ) / 4 E i , j = ( R 2 i , 2 j + R 2 i + 1,2 j + R 2 i , 2 j + 1 + R 2 i + 1,2 j + 1 ) / 4 (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.

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. T ( x , y ) = 1200 × exp [ − 0.05 × ( x + 0.2 L x ) 2 − 0.075 × ( y + 0.2 L y ) 2 − 1.375 ] + 1200 × exp [ − 0.075 × ( x + 0.2 L x ) 2 − 0.05 × ( y − 0.2 L y ) 2 − 0.925 ] + 1200 × exp [ − 0.075 × ( x − 0.2 L x ) 2 − 0.05 × ( y + L y ) 2 − 0.925 ] + 1200 × exp [ − 0.05 × ( x − 0.2 L x ) 2 − 0.075 × ( y − L y ) 2 − 1.375 ] + 300 ( K ) (8) 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.

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 ill-conditioned 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.

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.