Single-step matrix operation-based artifact-removing reconstruction method for laser absorption spectrum tomographic image

Abstract

The invention provides a laser absorption spectrum tomographic image de-artifact reconstruction method based on single-step matrix operation. Low-pass Gaussian filtering is fused into an iteration format of image reconstruction, and a single-step matrix operation format for recursive solving of a multi-step iteration result is constructed. According to tomography optical path layout, reconstruction image grid division and Gaussian filter template coefficients, a sensitivity matrix, a Gaussian filter matrix and other matrixes required for solving a multi-step iteration result are calculated and stored in advance, and a solution vector after multi-step iteration is directly obtained through single-step matrix operation. And meanwhile, according to a convergence result of a solution vector of the linear equation set, performing artifact removal on the reconstructed image in combination with non-negativity of the absorptivity and continuity constraint of absorptivity distribution. According to the method, traditional multi-step iteration is replaced with single-step matrix operation, the image reconstruction calculation efficiency is remarkably improved, the image reconstruction precision is effectively improved through low-pass filtering and artifact removal, and a new image reconstruction means is provided for rapid and high-precision imaging monitoring of a combustion field.

Description

Technical Field

[0001] The present invention provides a method for artifact removal and reconstruction of laser absorption spectroscopy tomography images based on single-step matrix operation, which directly obtains the results of multi-step iteration through single-step matrix operation and removes artifacts, belonging to the field of tomography imaging technology. Background Art

[0002] The measurement techniques for the temperature and component concentration in the combustion field can be classified into contact measurement techniques and non-contact measurement techniques according to whether they are in contact with the measured flow field. As a non-contact optical measurement method, the laser absorption spectroscopy technique has the advantages of high accuracy, large measurement field of view, simple structure, strong anti-interference ability, etc., and has been widely used in combustion monitoring. By combining the laser absorption spectroscopy technique with the tomography imaging technique and arranging multiple laser paths to pass through the combustion field for measurement, the inversion of the temperature and component concentration distribution can be realized. The absorption spectrum is non-linearly coupled with flow field parameters such as temperature and concentration, but the non-linear optimization algorithm has problems such as large computational amount and easy to fall into local optimum. Therefore, the calculation process from the absorption spectrum projection measurement value to the temperature and concentration distribution is usually transformed into a linear equation system solving process, mainly including two solving models: the discrete imaging model and the colorimetric imaging model. The discrete imaging model splits the non-linear decoupling process from the absorption spectrum data to the temperature and concentration into two steps: the inversion of the temperature-concentration histogram and the reconstruction of the temperature-concentration pair distribution. The former belongs to the solution of an overdetermined linear equation system, and the latter belongs to the solution of an underdetermined linear equation system. The colorimetric imaging model splits the non-linear decoupling process into two steps: spectral parameter tomography imaging and temperature and concentration decoupling inversion. The former belongs to the solution of an underdetermined linear equation system, and the latter realizes temperature and concentration decoupling through ratio operation.

[0003] The discrete imaging model of laser absorption spectroscopy tomography discretizes the solution domains of temperature and concentration into a finite number of temperature-concentration pairs. By solving the distribution positions of each temperature-concentration pair within the reconstruction region, temperature and concentration distribution imaging is achieved. In 2023, in the paper "A Binary Valued Reconstruction Algorithm for Discrete TDLAS Tomography of Dynamic Flames" published by Qiu et al. in IEEE Transactions on Instrumentation and Measurement, Volume 72, Page 9506614, the discrete imaging model of laser absorption spectroscopy tomography and the corresponding imaging method were first proposed, realizing high-precision imaging of the temperature and water molecule concentration distributions in the combustion field. This method includes two steps: temperature-concentration histogram inversion and temperature-concentration pair distribution reconstruction. First, discrete temperature-concentration pairs are selected based on the prior knowledge of the temperature and concentration values in the combustion field. Using the projection measurement values of the absorption spectrum at multiple wavenumber points, the temperature-concentration histogram on each projection is inversely calculated one by one to obtain the path length of each temperature-concentration pair on each projection. Then, based on the results of the temperature-concentration histogram inversion, the distribution positions of each temperature-concentration pair within the imaging region are reconstructed, thereby obtaining the reconstructed images of the temperature and concentration distributions in the combustion field. The discrete imaging model simultaneously utilizes the multi-point measurement values of the absorption spectrum, avoiding the errors caused by the spectral parameter ratio operation. Compared with the colorimetric imaging model, it has higher noise resistance and can effectively improve the imaging accuracy of temperature and concentration distributions.

[0004] However, the discrete imaging model requires a complete measurement of the absorption spectral profile and an accurate calibration of the wavenumber, which places high demands on the measurement system. The colorimetric imaging model only utilizes single-point information such as the integral area of the absorption rate or the normalized second harmonic peak, and has lower requirements for the measurement system, and is currently more widely used. Implementing temperature and concentration distribution imaging using the colorimetric imaging model involves two steps: one is spectral parameter tomography, and the other is decoupled inversion of temperature and concentration. When performing spectral parameter tomography, spectral parameters such as the absorption rate belong to path integral quantities, and the process of imaging their distribution is transformed into a process of solving a system of linear equations. In the dynamic monitoring of combustion fields, laser absorption spectroscopic tomography sensors usually adopt a fixed optical path structure to avoid the rotation or translation of the laser emission end and the receiving end, thereby obtaining a higher time resolution, but also limiting the number of projections that can be obtained. This results in a strongly underdetermined system of linear equations to be solved, and it is difficult to achieve high-precision image reconstruction only using limited projection measurement data. When performing decoupled inversion of temperature and concentration, the classical colorimetric method uses two absorption spectral lines to achieve temperature inversion by calculating the ratio of the integral absorption rates of the two spectral lines, and calculates the concentration based on the temperature. The integral absorption rate distributions of the two spectral lines each have reconstruction errors, which usually appear as artifacts with relatively small values in non-combustion regions, and the ratio calculation will amplify these errors, further introducing artifacts with relatively large values in the temperature and concentration reconstruction images, reducing the imaging quality of the temperature and concentration distributions.

[0005] To address the problem of the strongly underdetermined system of equations to be solved in spectral parameter tomography, researchers at home and abroad have utilized the prior knowledge of the continuous and smooth distribution of gas parameters in the combustion flow field and introduced continuity constraints during the image reconstruction process to improve the image reconstruction accuracy. The methods of introducing continuity constraints can be mainly divided into two types: one is to construct a basis function reconstruction model and directly add the continuity constraint to the reconstruction model; the other is to use a regularization algorithm and add the continuity constraint to the reconstruction algorithm.

[0006] Introducing continuity constraints through the construction of a basis function reconstruction model means using a linear combination of specific basis functions to represent the continuous absorption rate distribution. The unknowns in the equations to be solved are transformed from the absorption rates in the discrete grid to the coefficients of the basis functions, realizing the mapping of the reconstruction model from a high-dimensional space to a low-dimensional space, thereby reducing the number of unknowns to be solved and improving the image reconstruction accuracy. In 2016, in the paper "An Efficient Approach for Limited-Data Chemical Species Tomography and its Error Bounds" published by Polydorides et al. in Proceedings of the Royal Society A-Mathematical Physical and Engineering Sciences, Volume 472, Issue 2187, Page 20150875, the continuity and smoothness constraints of the gas component concentration distribution were introduced into the reconstruction model, and the reconstruction model was mapped from a high-dimensional space to a low-dimensional space using orthogonal global basis functions, reducing the number of unknowns to be solved. Since the selection of the basis function depends on the prior knowledge of the measured flow field distribution, Polydorides et al. selected a basis function that conforms to the Gaussian distribution according to the simulation results of the carbon dioxide plume distribution at the engine outlet to be measured, analyzed the error introduced by the mapping through theoretical derivation, and verified the algorithm performance through numerical simulation. In 2018, in the paper "Constrained Models for Optical Absorption Tomography" published by Polydorides et al. in Applied Optics, Volume 57, Issue 7, Pages B1-B9, based on the above basis function reconstruction model, by introducing inequality constraints, the boundaries of the carbon dioxide concentration reconstruction results were restricted, and the reconstruction problem was transformed into a regularization problem, which was solved using the Gauss-Newton algorithm. Simulations and experiments have shown that the introduction of inequality constraints not only helps to improve the spatial resolution of the reconstructed image but also can improve the convergence performance of the reconstruction algorithm.In 2022, in the paper "Sparse Zernike Fitting for Dynamic LAS Tomographic Images of Temperature and Water Vapor Concentration" published by Gao et al. on page 5009314 of Volume 71 of the "IEEE Transactions on Instrumentation and Measurement", by fitting the continuous absorption rate distribution with a finite number of Zernike orthogonal polynomials, the quantity to be solved was converted into the coefficients of the Zernike polynomials, effectively reducing the number of unknowns to be solved and improving the accuracy of image reconstruction. However, the fitting error of the Zernike polynomial for the continuous asymmetric distribution is relatively large. To improve the fitting accuracy of the basis function for the local details of the asymmetric distribution, in 2023, in the paper "Radial Basis Function Coupled SART Method for Dynamic LAS Tomography" published by Gao et al. on page 4500110 of Volume 72 of the "IEEE Transactions on Instrumentation and Measurement", a compactly supported radial basis function was used to fit the absorption rate distribution, the quantity to be solved was converted into the fitting coefficients of the compactly supported radial basis function, and the SART algorithm was used to iteratively solve the fitting coefficient vector. Compared with the result of directly reconstructing the absorption rate distribution using the SART algorithm, the image reconstruction error of this method is smaller, but the selection of parameters such as the shape factor and compact support radius of the basis function has a greater impact on the reconstruction result, and it is necessary to predict the shape characteristics of the measured distribution based on prior knowledge and determine the optimal combination of basis function parameters through simulation.In 2024, in the paper "Robust Temperature and Gas Concentration Imaging by LAS Tomography with Adaptive Basis Function Fitting and Artifact Removal" published by Zhang et al. in IEEE Transactions on Instrumentation and Measurement, Volume 73, Page 4505811, a method for adaptively selecting basis function parameters was proposed. The shape factor and basis element center of the basis function were adaptively optimized and selected using the initial image reconstruction result, and the optimized basis function parameter combination was used for secondary reconstruction, reducing the sensitivity of the final reconstruction result to the initial selection of basis function parameters. However, when representing the absorption rate distribution using basis functions, fitting errors are inevitably introduced. It is necessary to combine prior knowledge of the distribution characteristics of the measured flow field and select an appropriate basis function type to control the fitting error within an acceptable range. After selecting the type of basis function, the optimization and selection of parameters such as the shape factor of the basis function also depend on the prior knowledge of the measured flow field, and it is necessary to use computational fluid dynamics simulation to predict the actual flow field or sacrifice a certain computational speed for secondary reconstruction. In addition, although the number of unknowns to be solved is reduced, the basis function reconstruction model is still an underdetermined linear equation system, so it is also necessary to improve the solution accuracy of the reconstruction algorithm.

[0007] Introducing continuity constraints by using regularization algorithms means adding a regularization term to the reconstruction algorithm or achieving the regularization effect through low-pass spatial filtering, so that the reconstruction result of the absorption rate distribution conforms to the constraints of spatial smoothness and continuity, and improves the solution accuracy of the underdetermined system of equations. Common reconstruction algorithms include the Tikhonov regularization algorithm, the Landweber algorithm, the Algebraic Reconstruction Technique (ART), the Simultaneous Algebraic Reconstruction Technique (SART), the Modified Adaptive Algebraic Reconstruction Technique (MAART), etc. In 2010, in the paper "Infrared Species Limited Data Tomography through Tikhonov Reconstruction" published by Daun in the Journal of Quantitative Spectroscopy & Radiative Transfer, Volume 111, Issue 1, pages 105-115, the Tikhonov regularization algorithm was used to implement image reconstruction, and the discrete Laplacian operator was introduced into the Tikhonov regularization reconstruction process. This operator simulates the molecular diffusion motion that leads to a smooth distribution of component concentrations and acts as a low-pass filter. In 2015, in the paper "Simultaneous Temperature, Concentration, and Pressure Imaging of Water Vapor in a Turbine Engine" published by Wood et al. in the IEEE Sensors Journal, Volume 15, Issue 1, pages 545-551, the Tikhonov regularization method with the discrete Laplacian operator introduced was applied to the imaging of the temperature, water molecule concentration, and pressure distribution at the annular outlet of a turbine engine. The simulation proved that the average reconstruction error of this algorithm is less than that of the Landweber algorithm.In 2022, in the paper "Configuration Optimization of Optical Tomography Based on Fan-Beam Lasers" published by Li et al. in IEEE Transactions on Instrumentation and Measurement, Vol. 71, page 4505810, according to the prior knowledge of the continuity of the reconstructed image, the result of Tikhonov regularization was used as the initial value of the iteration of the Landweber algorithm to improve the image reconstruction quality and the iteration convergence speed. In 2022, in the paper "Tomographic Imaging of Carbon Dioxide in the Exhaust Plume of Large Commercial Aero-Engines" published by Upadhyay et al. in Applied Optics, Vol. 61, No. 28, pages 8540-8552, the Tikhonov regularization algorithm and the Landweber algorithm were respectively used to reconstruct the carbon dioxide concentration distribution in the exhaust plume of an aero-engine. In the iteration process of the Landweber algorithm, median filtering was introduced. The image reconstruction results obtained by the two algorithms are in good agreement and can reflect the change trend of the carbon dioxide concentration distribution at the exhaust nozzle of the aero-engine under different working conditions. In 2015, in the paper "Tomographic Imaging of the Liquid and Vapour Fuel Distributions in a Single-Cylinder Direct-Injection Gasoline Engine" published by Terzija et al. in International Journal of Engine Research, Vol. 16, No. 4, pages 565-579, the Landweber algorithm was used to realize the imaging of the fuel spray droplet distribution and the fuel vapor concentration distribution of a single-cylinder direct-injection gasoline engine. Median filtering was introduced as a continuity constraint in the iteration, and at the same time, a non-negativity constraint was added to make the concentration imaging results conform to physical meaning.In 2015, in the paper "A Hybrid Tomographic Reconstruction Algorithm for High Speed X-ray Tomography" published by Yang et al. in Computer Physics Communications, Volume 196, pages 27 - 35, median filtering was added to the SART algorithm, and the median filtering result of the previous iteration was used to correct the direct iteration result of the current SART, thereby reducing the salt-and-pepper noise in the reconstructed image. In 2024, in the paper "Laser Beam Optimization for TDLAS Tomography from Asymptotic Point Spread Functions" published by Wen et al. in IEEE Transactions on Instrumentation and Measurement, Volume 73, page 4501311, the above SART algorithm with added median filtering was applied to the field of laser absorption spectroscopy tomography, and the reconstruction of the flame temperature and water molecule concentration distribution was achieved. In 2015, in the paper "Development of a Fan-Beam TDLAS-Based Tomographic Sensor for Rapid Imaging of Temperature and Gas Concentration" published by Liu et al. in Optics Express, Volume 23, Issue 17, pages 22494 - 22511, the modified Landweber algorithm was used to reconstruct the temperature and water molecule concentration distribution of a flat-flame combustion furnace, and low-pass Gaussian filtering was added in each iteration. Compared with the ART algorithm and the joint multiplicative algebraic reconstruction algorithm, this algorithm has higher reconstruction accuracy.In 2020, in the paper "Frequency-Division Multiplexing and Main Peak Scanning WMS Method for TDLAS Tomography in Flame Monitoring" published by Huang et al. in the IEEE Transactions on Instrumentation and Measurement, Vol. 69, No. 11, pp. 9087-9096, dynamic imaging of the temperature and water molecule concentration distribution of an acoustically excited Bunsen flame was carried out. The modified Landweber algorithm was also used to implement image reconstruction, and low-pass spatial filtering was added to the reconstruction algorithm. In 2011, in the paper "Modified Adaptive Algebraic Tomographic Reconstruction of Gas Distribution from Incomplete Projection by a Two-Wavelength Absorption Scheme" published by Li et al. in Chinese Optics Letters, Vol. 9, No. 6, p. 061201, the MAART algorithm was proposed, which adaptively calculates the relaxation factor in each iteration and performs weighted mean filtering on the solution vector obtained in each iteration. In 2021, in the paper "Integral Absorbance Measurement for a Non-Uniform Flow Field Using Wavelength Modulation Absorption Spectroscopy" published by Song et al. in Applied Optics, Vol. 60, No. 17, pp. 5056-5065, the measurement of integral absorbance was realized based on calibration-free wavelength modulation absorption spectroscopy, and the MAART algorithm was used for image reconstruction to achieve imaging of the temperature and concentration distribution of a non-uniform flow field.In 2023, in the paper "Rapid Online Tomograph in Non-Uniform Complex Combustion Fields Based on Laser Absorption Spectroscopy" published by Zhao et al. in Experimental Thermal and Fluid Science, Volume 147, Page 110930, the MAART algorithm was applied to the tomographic imaging of the flame temperature and water vapor concentration distribution at the afterburner exit. In 2017, in the paper "Two-Step Tomographic Reconstructions of Temperature and Species Concentration in a Flame Based on Laser Absorption Measurements with a Rotation Platform" published by Xia et al. in Optics and Lasers in Engineering, Volume 90, Pages 10 - 18, a two-step reconstruction strategy was proposed. First, the MAART algorithm with weighted mean filtering was used to reconstruct the absorption rate distribution. After obtaining the temperature distribution by colorimetry, the reconstructed temperature value was substituted into the calculation model of the absorption rate integral area to obtain a linear equation system with the concentration value as the unknown. Then, the MAART algorithm was used again for iterative solution to obtain the concentration distribution, effectively improving the reconstruction accuracy of the concentration distribution. In 2023, in the paper "Reconstruction Algorithm Optimization Based on Multi-Iteration Adaptive Regularity for Laser Absorption Spectroscopy Tomography" published by Zhao et al. in Applied Sciences-Basel, Volume 13, Issue 21, Page 12083, the MAART algorithm was improved. The reconstruction result of the Landweber algorithm iterated 3 times was used as the initial value of the MAART algorithm iteration, reducing the reconstruction error caused by inappropriate initial value selection of the MAART algorithm.Meanwhile, an algorithm for adaptively adjusting the regularization parameter was proposed, which realized the adaptive iterative calculation of the template coefficients of the weighted mean filter, replaced the fixed weighted mean filter in the original algorithm, and improved the image reconstruction quality. The above reconstruction algorithm effectively improved the image reconstruction accuracy by adding a regularization term or low-pass spatial filters such as median filtering, low-pass Gaussian filtering, and weighted mean filtering. However, considering factors such as the underdetermination of the equation set and the measurement error of spectral data, there are still small numerical artifacts in the reconstructed image of the absorption rate distribution. When using the colorimetric method for temperature and concentration decoupling, the reconstruction errors caused by these artifacts are amplified by the ratio calculation, introducing more serious artifacts in the reconstructed images of temperature and concentration distributions and reducing the reconstruction accuracy of temperature and concentration distributions.

[0008] Aiming at the problem that the ratio calculation in temperature and concentration decoupling inversion amplifies the reconstruction error, researchers at home and abroad have proposed two types of artifact removal strategies: one is to constrain the ratio of the absorption rate distribution and use the characteristics that the ratio of the absorption rate distribution is continuous and smooth in space to remove the artifacts in the reconstructed images of temperature and concentration distributions; the other is to constrain the boundary of the absorption rate distribution, directly remove the artifacts in the reconstructed image of the absorption rate distribution, and then remove the artifacts in the reconstructed images of temperature and concentration distributions.

[0009] Artifact removal is achieved by constraining the ratio of absorption rate distributions, taking advantage of the fact that the absorption rate ratio is only related to temperature. During the image reconstruction process, a constraint of continuous absorption rate ratio is added, or the spatial continuity of the ratio is used to identify and remove artifacts in the reconstructed distributions of temperature and concentration. In 2021, in the paper "Relative Entropy Regularized TDLAS Tomography for Robust Temperature Imaging" published by Bao et al. on page 4501909, Volume 70 of the "IEEE Transactions on Instrumentation and Measurement", the absorption rate ratio of two absorption spectral lines was used as a smoothing constraint. Through the relative entropy regularization term, this constraint was added to the solution model and solved using optimization methods. This method can effectively suppress spike noise in the reconstructed image and reduce some artifacts in the reconstructed images of temperature and concentration distributions, but it cannot remove relatively smooth artifacts. In 2024, in the paper "Robust Temperature and Gas Concentration Imaging by LAS Tomography with Adaptive Basis Function Fitting and Artifact Removal" published by Zhang et al. on page 4505811, Volume 73 of the "IEEE Transactions on Instrumentation and Measurement", the SART algorithm with median filtering was used to reconstruct the distribution of basis function coefficients, obtaining images of the absorption rate distributions of two absorption spectral lines. Temperature and concentration distribution imaging was achieved by colorimetry, and artifacts introduced by ratio calculation were identified and removed based on the similarity between the distribution of the absorption rate ratio of the two spectral lines and the absorption rate distribution of one of the spectral lines. This method achieves artifact removal by removing the basis function components corresponding to the artifacts, but it is only applicable to the basis function model and still has certain limitations.

[0010] Artifact removal by constraining the boundary of the absorption rate distribution refers to adding an edge optimization algorithm or a posteriori topology during the image reconstruction process to obtain a reconstructed image of the absorption rate distribution with clear boundaries, and realizing artifact removal in the reconstructed images of the temperature and concentration distributions by removing the artifacts in the non-burning regions of the reconstructed image of the absorption rate distribution. In 2023, in the paper "Laser Absorption Tomography of Complex Combustion Fields Based on Finite-Element Node Strategy and Adaptive Edge Optimization Algorithm" published by Zhao et al. on page 102251, Volume 46 of Thermal Science and Engineering Progress, the finite-element node strategy was used to replace the traditional discrete grid division strategy, the MAART algorithm with weighted mean filtering was adopted for the reconstruction of the absorption rate distribution, and then the colorimetric method was used to realize the reconstruction of the temperature and concentration distributions. At the same time, aiming at the problem of uneven division of edge nodes, an adaptive edge optimization algorithm was proposed, which effectively reduced the edge reconstruction error caused by the finite-element node strategy and reduced the artifacts in the reconstructed image of the absorption rate distribution. However, the artifact removal strategy of this method can only handle the artifacts caused by the finite-element node division and cannot directly constrain the boundary of the absorption rate distribution in the reconstructed image. In 2022, in the paper "A Novel Parametric Level Set Method Coupled with Tikhonov Regularization for Tomographic Laser Absorption Reconstruction" published by Niu et al. on page 117819, Volume 201 of Applied Thermal Engineering, a Gaussian parametric level set method combined with the regularized Landweber algorithm was proposed. A posteriori topology was added during the solution process to make up for the lack of prior information, and high-precision reconstruction of the temperature and water molecule concentration distributions in the turbulent combustion field was realized. This method constrains the boundary of the absorption rate distribution in the reconstructed image, can obtain a clear boundary contour of the combustion field, effectively removes the artifacts in the reconstructed image of the absorption rate distribution, and thus realizes artifact removal in the reconstructed images of the temperature and concentration distributions.In 2024, in the paper "Nonlinear Deflection Tomography of Inhomogeneous Flame Temperature and Concentration Based on Topology Evolution with Prior Smoothing" published by Niu et al. in Combustion and Flame, Volume 263, Page 113398, this method was used to image the inhomogeneous flame temperature and water molecule concentration, and further considered the influence of the light deflection effect caused by the refractive index change of the inhomogeneous combustion field on image reconstruction. However, the complete solution process of this method requires multiple iterations, and the two sub-steps in each iteration are also iterative solutions, with a complex calculation process and a large amount of computation.

[0011] Based on the above background, the present invention proposes a method for artifact removal and reconstruction of laser absorption spectroscopy tomography images based on single-step matrix operations. According to the prior knowledge of the continuity of the absorption rate distribution, the low-pass Gaussian filtering process is expressed in matrix form and introduced into the iterative calculation format of image reconstruction. Using the recurrence relation, the solution vector after multiple iterations is directly obtained through single-step matrix operations. At the same time, using the convergence result of the solution vector, the non-negativity constraint of the absorption rate, and the continuity constraint of the absorption rate distribution, the artifacts in the reconstructed image of the absorption rate distribution in the non-combustion region are removed. This method adds low-pass Gaussian filtering during the image reconstruction process to achieve regularization, improving the solution accuracy of the underdetermined linear equations in spectral parameter tomography; removing artifacts from the reconstructed image of the absorption rate distribution, constraining the boundary of the absorption rate distribution, and effectively removing the artifacts in the reconstructed images of the temperature and concentration distributions; aiming at the problems of large computational amount and slow calculation speed in the iterative solution process of linear equations, a solution method based on single-step matrix operations is proposed, significantly improving the image reconstruction speed and providing a new image reconstruction means for fast and high-precision imaging monitoring of the combustion field. Summary of the Invention

[0012] The present invention proposes a method for artifact removal and reconstruction of laser absorption spectroscopy tomography images based on single-step matrix operations. The Gaussian filtering process is expressed in matrix form and introduced into the iterative calculation format of image reconstruction. Using the recurrence relation, the solution vector after multiple iterations is directly obtained through single-step matrix operations, and using the convergence result of the solution vector, combined with the non-negativity and continuity constraints of the absorption rate distribution, the artifacts in the reconstructed image are removed. The specific implementation of the method is as follows:

[0013] Step 1: Basic parameter settings, including sensitivity matrix calculation, projection data vector acquisition, and Gaussian filtering matrix calculation.

[0014] The model of laser absorption spectral tomography is as follows:

[0015] Ax = b (1)

[0016] where A is a sensitivity matrix of dimension M×N, which is calculated according to the optical path layout of the laser absorption spectral tomography sensor and the grid division of the reconstructed image. M represents the total number of projection rays, N represents the total number of grids in the reconstructed image, and M < N is satisfied. The element A i，j in A represents the length of the i-th projection ray in the j-th grid; b is a projection data vector of dimension M×1, which is measured by the laser absorption spectral tomography sensor; x is a vector of dimension N×1 rearranged from the reconstructed image; Gx represents Gaussian filtering of the reconstructed image corresponding to x, where G is a Gaussian filtering matrix of dimension N×N, which is calculated according to the template coefficients of the Gaussian filter and the spatial positions of the elements in x in the reconstructed image.

[0017] Step 2: Single-step matrix operation, directly calculating the solution vector after multiple-step iteration using the recurrence relation.

[0018] The iterative format of image reconstruction after introducing Gaussian filtering is as follows:

[0019] x (k) = Gx (k1) + ωV -1 A T W(b - AGx (k1) ) (2)

[0020] where x (k) is the solution vector of the k-th iteration, and x (k-1) is the solution vector of the (k - 1)-th iteration; ω is a relaxation factor, satisfying 0 < ω < 2; V is a diagonal matrix of dimension N×N, and its diagonal elements are A +，j , that is, the column sum of the j-th column of the sensitivity matrix, and there is: W is a diagonal matrix of dimension M×M, and its diagonal elements are , that is, the reciprocal of the row sum of the i-th row of the sensitivity matrix, and there is:

[0021] When k = 1, from formula (2), we get:

[0022] x (1) = Gx (0) + ωV -1 A T W(b - AGx (0) ) (3)

[0023] where x (0) is the initial value of the solution vector. Let B = V -1 A TWA, obtained from Equation (3):

[0024] x (1) =(I - ωB)Gx (0) +ωV -1 A T Wb (4)

[0025] Where I represents the identity matrix. Let D = (I - ωB)G, then from Equation (4):

[0026] x (1) = Dx (0) + IωV -1 A T Wb (5)

[0027] When k = 2, from Equation (2) and Equation (5):

[0028] x (2) = D 2 x (0) +(D + I)ωV -1 A T Wb (6)

[0029] According to the recurrence relation, the expression for the solution vector at the k-th iteration is obtained as:

[0030] x (k) = D k x (0) +(D k-1 + D k-2 +... + D + I)ωV -1 A T Wb (7)

[0031] When k approaches positive infinity, D k approaches the zero matrix. Since in actual calculations, k can only take a finite value, so take all elements of x (0) to be 0 for calculation, excluding the error caused by the value of x (0) .

[0032] Let x (0) = [0 0 … 0] T , then Equation (7) simplifies to:

[0033] x (k) =(D k-1 + D k-2 + … + D + I)ωV -1 A T Wb (8)

[0034] Since D, V, A, and W are only related to the optical path layout of the laser absorption spectroscopy tomography sensor, the grid division of the reconstructed image, and the template coefficients of the Gaussian filter, and are independent of the measured values, for any value of k, D is pre-calculated k-1 +D k-2 +…+D+I and V -1 A T W, and by setting the value of ω, according to formula (8), the calculation result of x (k) can be directly obtained using a single-step matrix operation without performing multi-step iterative operations, thus achieving fast reconstruction.

[0035] Step 3: Remove the artifacts in the reconstructed image. By using the convergence result of the solution vector and combining the non-negativity constraint of the absorption rate and the continuity constraint of the absorption rate distribution, the artifacts in the reconstructed image are removed.

[0036] In laser absorption spectroscopy tomography, the elements of the solution vector represent the absorption rate or the normalized second harmonic peak calculated from the absorption rate, which is non-negative. At the same time, since the temperature and concentration of the target gas molecules are spatially continuous, the reconstructed image related to the absorption rate distribution is also spatially continuous. Using these two properties to constrain x (k) can achieve artifact removal and obtain a high-precision image reconstruction result.

[0037] Let x * be the true value of the solution of the linear equation system shown in formula (1), satisfying Ax * =b. The two-dimensional image obtained after rearranging x * is the discretized true distribution. Transforming formula (2) gives:[[]]

[0038] x (k) =Gx (k-1) +ωB(x * -Gx (k-1) ) (9)

[0039] Define the N×N-dimensional spatial transformation matrix C, and let C=(I - D) -1 ωB. Then, from formula (9), we get:[[]]

[0040] x (k) -Cx * =D(x (k-1) -Cx * ) (10)

[0041] Formula (10) constructs a geometric vector sequence {x (k) -Cx *}, and the modulus of all eigenvalues of its coefficient matrix D is less than 1. Therefore, when k approaches positive infinity, x (k) converges to Cx * .

[0042] When k is large enough, x (k) ≈ Cx * , that is, the solution vector of the k-th iteration is the spatial transformation of the true distribution. Based on the calculation result of x (k) , artifact removal is performed, and x est is used to represent the solution vector after artifact removal. By calculating the difference between cx est and x (k) , the difference between cx est and Cx * is estimated, so as to estimate the proximity between x est and x * . According to the non-negativity constraint, the elements of x est are set to 0 one by one. When the difference between Cx est and x (k) decreases, the corresponding element is an artifact, and its value is kept as 0 to make x est approach x * . For the salt-and-pepper noise introduced by the zeroing operation, according to the continuity constraint, median filtering is performed on the two-dimensional image rearranged from x est . The zeroing operation and median filtering are repeated until the difference between Cx est and x (k) no longer changes. At this time, the two-dimensional image rearranged from x est is the reconstructed image after artifact removal.

[0043] The artifact removal process is specifically implemented according to the following steps:

[0044] First step, according to the non-negativity of the absorption rate, the reconstructed image is globally constrained. The elements in x (k) that are less than 0 are set to 0, denoted as and the difference between and x (k) is calculated: Set the threshold ε0.

[0045] Second step, according to the non-negativity of the absorption rate, the reconstructed image is pixel-by-pixel constrained. According to the non-negativity constraint of the absorption rate, each element of is operated on one by one. The 1st, 2nd,..., Nth elements in are set to 0 one by one, and the difference between and x (k) is calculated: 6 1，1 , δ 1，2 , …, δ 1，N . If δ 1，j < δ 1，j-1 , then the value of the th element is kept as 0; if δ 1，j ≥ δ 1，j-1, then set the value of the ,-th element to the corresponding element value in After operating on each element of

[0046] Step 3: Determine and remove the artifacts in the reconstructed image according to the continuity of the absorption rate distribution. Update according to the continuity constraint of the absorption rate distribution Set rearrange it into a two-dimensional image, and perform median filtering. Determine the elements that are still 0 after median filtering as artifacts and keep their values as 0. For the elements that are not 0 after median filtering, set their values to the corresponding element value in Calculate the difference between (k) and x And calculate 6 2，0 relative to 6 1，0 the change of: ε1 = |δ 2，0 - δ 1，0 | / |δ 1，0 |.

[0047] Step 4: Judge whether the artifact removal is completed according to the relative change of the reconstructed image before and after artifact removal, and output the final result of the reconstructed image. Compare ε1 with the given threshold ε0. If ε1 < ε0, then output the artifact removal result and rearrange it to obtain the final result of the reconstructed image; if ε1 ≥ ε0, then repeat Step 2 and Step 3 until ε l < ε0, output the artifact removal result and rearrange it to obtain the final result of the reconstructed image.

[0048] Advantages of the present invention: The image reconstruction method proposed by the present invention constructs an iterative calculation format in matrix form, directly obtains the solution vector after multiple-step iteration through single-step matrix operations, and realizes the fast reconstruction of tomographic images; uses the convergence result of the solution vector, the non-negativity constraint of the absorption rate, and the continuity constraint of the absorption rate distribution to remove artifacts, realizes the high-precision reconstruction of laser absorption spectroscopy tomographic images, and provides a new image reconstruction means for the fast and high-precision imaging monitoring of the combustion field. Description of the Drawings

[0049] Figure 1 is the flowchart of the artifact removal and reconstruction method for laser absorption spectroscopy tomographic images based on single-step matrix operations.

[0050] Figure 2 is the schematic diagram of the beam layout of the laser absorption spectroscopy tomographic imaging sensor.

[0051] Figure 3 is the original distribution of temperature.

[0052] Figure 4 is the original distribution of the mole fraction of water molecules.

[0053] Figure 5 is the reconstructed image of temperature.

[0054] Figure 6 is the reconstructed image of the mole fraction of water molecules. Detailed implementation manners

[0055] The present invention will be further described below with reference to examples.

[0056] In this example, according to the principle of laser absorption spectroscopic tomography, an artifact removal reconstruction method for laser absorption spectroscopic tomography images based on single-step matrix operation is adopted for image reconstruction. The schematic diagram of the beam layout of the laser absorption spectroscopic tomography sensor used in this example is as shown in the appendix Figure 2 shown. Among them, the frame of the sensor is circular, and 9 sector laser sources are evenly distributed on the circumference at equal angular intervals. There are 7 detectors between two adjacent laser sources, which are also arranged on the circumference at equal angular intervals. The laser emitted by each sector laser source can be received by 3 groups of detectors on the opposite side, that is, 21 detectors, and a total of 189 projection rays are generated, that is, M = 189. Image reconstruction is performed on the area inside the circular sensor, that is, the diameter of the reconstructed image is equal to the diameter of the sensor. The reconstruction grid along the diameter is set to 60, and the total number of reconstruction grids is 2828, that is, N = 2828. Set a group of circular original distributions of temperature and mole fraction of water molecules as shown in the appendix Figure 3 and the appendix Figure 4 shown. Their diameters are also equal to the diameter of the circular sensor. The discrete grid along the diameter is set to 100, so as to simulate the continuous measured distribution in actual measurement.

[0057] In this example, two water molecule absorption spectral lines with central wave numbers at 7185.6 cm -1 and 7444.4 cm -1 are selected, and the temperature and the mole fraction of water molecules are reconstructed by using the principle of colorimetry. The size of the Gaussian filter template is 3×3, and the template coefficients are: G1 = 9.847×10 -1 , G2 = 3.807×10 -3 , G3 = 1.472×10 -5 , and it satisfies G1 + 4G2 + 4G3 = 1.

[0058] The specific implementation steps of the image reconstruction process are as follows:

[0059] Step 1: Set basic parameters.

[0060] According to the optical path layout of the laser absorption spectral tomography sensor and the grid division of the reconstructed image, a sensitivity matrix A of 189×2828 dimensions is calculated. For the water molecule absorption spectral line with a central wavenumber of 7185.6 cm -1 Based on the principle of laser absorption spectrum measurement, the wavelength modulation method is adopted. From the original distributions of temperature and water molecule concentration, the projection values of the normalized second harmonic peaks of two spectral lines on 567 laser paths are calculated respectively, obtaining a projection data vector b of 189×1 dimension. According to the Gaussian filter template coefficients and the grid division of the reconstructed image, a Gaussian filter matrix G of 2828×2828 dimensions is calculated.

[0061] Step 2: Single-step matrix operation.

[0062] According to the sensitivity matrix A, calculate V -1 A T W; take k = 2 15 , from B = V -1 A T WA and D = (I - ωB)G, calculate D k-1 +D k-2 +…+D+I. According to formula (8), take ω = 1, and obtain x (k) through single-step matrix operation.

[0063] Step 3: Reconstructed image artifact removal, including the following five sub-steps:

[0064] First step: According to the non-negativity of the absorption rate, perform overall constraint on the reconstructed image. Set the elements in x (k) that are less than 0 to 0, obtaining Calculate the 2828×2828-dimensional spatial transformation matrix C from C = (I - D) -1 ωB. Calculate the difference between and x (k) : Set the threshold ε0 = 10 -6 .

[0065] Second step: According to the non-negativity of the absorption rate, perform pixel-by-pixel constraint on the reconstructed image. According to the non-negativity constraint of the absorption rate, perform one-by-one operation on the elements of , and set the 1st, 2nd, …, 2828th elements in to 0 one by one, and calculate the difference between and x (k) : δ 1，1 , δ 1，2 , …, δ 1，2828 . If δ 1，j < δ 1，j-1 , then keep the value of the, -th element as 0; if δ1，j ≥δ 1，j-1 , then set the value of the \(j\)-th element to the corresponding element value in, where \(j = 1, 2, \ldots, 2828\). After operating on each element of , obtain

[0066] Step 3: Determine and remove the artifacts in the reconstructed image according to the continuity of the absorption rate distribution. Update according to the continuity constraint of the absorption rate distribution Arrange rearranged into a two-dimensional image and perform median filtering. The template size of the median filter is \(3\times3\). Determine the elements that are still 0 after median filtering as artifacts and keep their values as 0. For the elements that are not 0 after median filtering, set their values to the corresponding element value in. Calculate the difference between (k) and \(x\): and calculate \(\delta\) 2，0 relative to \(\delta\) 1，0 : \(\varepsilon_1 = |\delta\) 2，0 - \(\delta\) 1，0 | / |\(\delta\) 1，o |.

[0067] Step 4: Determine whether the artifact removal is completed according to the relative change of the reconstructed image before and after artifact removal, and output the final result of the reconstructed image. Compare \(\varepsilon_1\) with the given threshold \(\varepsilon_0\). If \(\varepsilon_1 \lt \varepsilon_0\), then output the artifact removal result and rearrange to obtain the final result of the reconstructed image; if \(\varepsilon_1 \geq \varepsilon_0\), then repeat Step 2 and Step 3 until \(\varepsilon\) l \lt \varepsilon_0, output the artifact removal result and rearrange to obtain the final result of the reconstructed image.

[0068] For the water molecule absorption spectral line with a central wave number of 7185.6 cm -1 , when \(l = 20\), obtain the artifact removal result and rearrange to obtain the final result of the reconstructed image. For the water molecule absorption spectral line with a central wave number of 7444.4 cm -1 , repeat the above process. When \(l = 18\), obtain the artifact removal result and rearrange to obtain the final result of the reconstructed image.

[0069] According to the principle of colorimetry, the reconstructed distribution of temperature is calculated as shown in Appendix Figure 5 , and the reconstructed image of the water molecule mole fraction is as shown in Appendix Figure 6 . Perform nearest-neighbor interpolation on the original distributions of temperature and water molecule mole fraction so that the vector rearranged from them has the same number of elements as the vector \(x\)(k) Same dimension. Let x ori represent the vector rearranged after interpolation of the original distribution of temperature or mole fraction of water molecules. Let x rec represent the vector rearranged from the reconstructed image of temperature or mole fraction of water molecules. Then the image reconstruction error is defined as: e = ||x rec - x ori ||₂ / ||x ori ||₂. The calculated Figure 5 image reconstruction error of the corresponding temperature reconstructed image is 0.1792, Figure 6 the image reconstruction error of the corresponding mole fraction of water molecules reconstructed image is 0.3907, and the average time to obtain the solution vector through single-step matrix operation is 0.029 s. While under the same other conditions, when using the joint algebraic reconstruction algorithm for reconstruction, the image reconstruction error of the obtained temperature is 1.0334, and the image reconstruction error of the mole fraction of water molecules is 1.3832, and the average time to obtain the solution vector through iterative operation is 0.254. It can be seen that compared with the classical joint algebraic reconstruction algorithm, the image reconstruction method proposed by the present invention effectively reduces the image reconstruction error and improves the operation speed.

[0070] The above description of the present invention and its implementation manners is not limited thereto. What is shown in the drawings is only one of the implementation manners of the present invention. Without departing from the gist of the present invention, structures or embodiments designed without creativity and similar to the technical solution all fall within the protection scope of the present invention.

Claims

1. A method for artifact removal reconstruction of laser absorption spectroscopic tomography images based on single-step matrix operation, characterized in that, The Gaussian filtering process is expressed in matrix form and introduced into the iterative calculation format of image reconstruction. Using the recurrence relation, the solution vector after multiple-step iterations is directly obtained through single-step matrix operations. Then, by using the convergence result of the solution vector and combining the non-negativity of the absorption rate and the continuity constraint of the absorption rate distribution, the artifacts in the reconstructed image are removed. The specific implementation of the method is as follows: Step 1: Set basic parameters, including calculating the sensitivity matrix, obtaining the projection data vector, and calculating the Gaussian filtering matrix; The model of laser absorption spectroscopy tomography is: Ax = b (1) Among them, A is a sensitivity matrix of M×N dimensions, which is calculated according to the optical path layout of the laser absorption spectroscopy tomography sensor and the grid division of the reconstructed image. M represents the total number of projection rays, and N represents the total number of grids in the reconstructed image, satisfying M < N. The element A in A i,j represents the length of the i-th projection ray in the j-th grid. b is a projection data vector of M×1 dimensions, which is measured by the laser absorption spectroscopy tomography sensor. x is an N×1-dimensional vector rearranged from the reconstructed image. Let Gx denote the Gaussian filtering of the reconstructed image corresponding to x, where G is an N×N-dimensional Gaussian filtering matrix, which is calculated according to the template coefficients of the Gaussian filter and the spatial positions of the elements in x in the reconstructed image; Step 2: Perform single-step matrix operations to directly calculate the solution vector after multiple-step iterations using the recurrence relation; The iterative format of image reconstruction after introducing Gaussian filtering is: x (k) = Gx (k-1) + ωV -1 A T W(b - AGx (k-1) ) (2) where x (k) is the solution vector of the k-th iteration, x (k-1) is the solution vector of the (k - 1)-th iteration, ω is the relaxation factor, satisfying 0 < ω < 2, V is an N×N diagonal matrix, whose diagonal elements are A +,j , that is, the column sum of the j-th column of the sensitivity matrix, and there is: W is an M×M diagonal matrix, whose diagonal elements are , that is, the reciprocal of the row sum of the i-th row of the sensitivity matrix, and there is: When k = 1, it is obtained from formula (2): x (1) = Gx (0) + ωV -1 A T W(b - AGx (0) ) (3) where x (0) is the initial value of the solution vector; let B = V -1 A T WA, obtained from formula (3): x (1) =(I - ωB)Gx (0) + ωV -1 A T Wb (4) where I represents the identity matrix; let D = (I - ωB)G, and it is obtained from formula (4): x (1) = Dx (0) + IωV -1 A T Wb (5) When k = 2, it is obtained from formula (2) and formula (5): x (2) = D 2 x (0) + (D + I)ωV -1 A T Wb (6) According to the recurrence relation, the expression of the solution vector for the k-th iteration is: x (k) = D k x (0) +(D k-1 + D k-2 +…+ D + I)ωV -1 A T Wb (7) As k approaches positive infinity, D k approaches the zero matrix. Since in actual calculations, k can only take on finite values, we take all elements of x (0) to be 0 for the calculation, excluding the error caused by the value of x (0) . Let x (0) = [0 0 … 0] T , then formula (7) is simplified to: x (k) =(D k-1 +D k-2 +…+D+I)ωV -1 A T Wb (8) For any value of k, according to formula (8), the calculation result of x is directly obtained by using a single-step matrix operation. (k) ​ Step 3: Remove the artifacts in the reconstructed image. By using the convergence result of the solution vector and combining the non-negativity constraint of the absorption rate and the continuity constraint of the absorption rate distribution, the artifacts in the reconstructed image are removed; In laser absorption spectroscopy tomography, the elements of the solution vector represent the absorption rate or the normalized second harmonic peak calculated from the absorption rate, which is non-negative. At the same time, since the temperature and concentration of the target gas molecules are spatially continuous, the reconstructed image related to the absorption rate distribution is also spatially continuous; For the x calculated using formula (8) (k) , when k approaches positive infinity, x (k) converges to Cx * , where C is an N×N-dimensional space transformation matrix, and there is: C = (I - D) -1 ωB, x * is the true value of the solution of the linear equation system shown in formula (1), satisfying Ax * = b, and the two-dimensional image obtained after rearranging x * is the discretized true distribution; When k is large enough, x (k) ≈ Cx * , that is, the solution vector of the k-th iteration is the spatial transformation of the true distribution; Based on the calculation result of x (k) , artifact removal is performed, and the solution vector after artifact removal is represented by x est . By calculating the difference between Cx est and x (k) , the difference between Cx est and Cx * is estimated, so as to estimate the closeness between x est and x * ; According to the non-negativity constraint, the elements of x est are set to 0 one by one. When the difference between Cx est and x (k) decreases, the corresponding element is an artifact, and its value is kept as 0 to make x est approach x * ; For the salt-and-pepper noise introduced by the zeroing operation, according to the continuity constraint, median filtering is performed on the two-dimensional image rearranged from x est ; The zeroing operation and median filtering are repeated until the difference between Cx est and x (k) no longer changes. At this time, the two-dimensional image rearranged from x est is the reconstructed image after artifact removal.

2. The artifact removal reconstruction method for laser absorption spectroscopic tomography images based on single-step matrix operation according to claim 1, wherein The artifacts are removed by using the convergence result of the solution vector, the non-negativity constraint of the absorption rate, and the continuity constraint of the absorption rate distribution. The specific implementation is as follows: Step 1: According to the non-negativity of the absorption rate, perform an overall constraint on the reconstructed image; set the elements in x (k) that are less than 0 to 0, denoted as Calculate the difference from x (k) : Set the threshold ε0; Step 2: According to the non-negativity of the absorption rate, perform pixel-by-pixel constraint on the reconstructed image; perform operations on each element of one by one, set the 1st, 2nd, …, Nth elements in to 0 one by one, and calculate the difference between and x (k) : δ 1,1 , δ 1,2 , …, δ 1,N . If δ 1,j < δ 1,j-1 , then keep the value of the jth element as 0. If δ 1,j ≥ δ 1,j-1 , then set the value of the jth element to the value of the corresponding element in , where j = 1, 2, …, N. After operating on each element of , the obtained vector is denoted as Step 3: Determine and remove the artifacts in the reconstructed image according to the continuity of the absorption rate distribution; update according to the continuity constraint of the absorption rate distribution Rearrange as a two-dimensional image and perform median filtering. Determine the elements that are still 0 after median filtering as artifacts and keep their values as 0. For the elements that are not 0 after median filtering, set their values to the corresponding element values in; calculate the difference between (k) and x and calculate δ 2,0 with respect to δ 1,0 : ε1 = |δ 2,0 - δ 1,0 | / |δ 1,0 |; Step 4: Determine whether the artifact removal is completed based on the relative change in the reconstructed images before and after artifact removal, and output the final result of the reconstructed images; compare ε1 with the given threshold ε0. If ε1 < ε0, output the artifact removal result and perform rearrangement to obtain the final result of the reconstructed image; if ε1 ≥ ε0, repeat Step 2 and Step 3 until ε l < ε0, output the artifact removal result and perform rearrangement to obtain the final result of the reconstructed image.

Images