基于LSMR算法的火焰三维温度场重建

单良，仰文淇，洪波，周荣幸，孔明

（1 中国计量大学信息工程学院，浙江省电磁波信息技术与计量检测重点实验室，浙江 杭州 310018；2 中国计量大学计量测试工程学院，浙江 杭州 310018）

燃烧和火焰作为自然界和人类生产、生活中普遍存在的现象，在内燃机、电厂锅炉、火箭推进器、燃气轮机和化学工程等工业生产领域扮演着重要角色[1-2]。及时并且准确地获得火焰温度分布，可以更好地探究燃烧的本质和规律，并在此基础上提高燃烧效率，优化燃烧过程[3]。目前图像探测器的分辨率可以达到很高的像素级别，高集成度、高分辨率、高性能、低成本的图像传感器的迅速发展，推动了火焰辐射成像测温技术的发展。火焰辐射成像测温技术具有测量精度高、测量实时性好、测量范围广等优点，逐渐成为目前火焰温度场测量技术领域的研究热点[4-5]。

闫勇等[6]使用单个CCD相机，结合图像处理和双色辐射测温技术，获得了燃烧系统中火焰的三维温度信息。这里使用单个CCD相机拍摄火焰图像，获取的火焰信息有限，导致温度场重建的效果会受到一定影响。如果要提高火焰温度场重建的准确性和效果，则需要使用多台相机从不同的角度拍摄火焰图像。周怀春等[7]使用多个CCD相机获取火焰图像，引入改进的Tikhonov正则化方法重建出三维温度场。岑可法等[8]基于多CCD相机系统，利用代数重建方法（ART）测量火焰断面温度场。这些团队系统研究了求解火焰辐射反问题的算法，并分析了这些方法的重建特性。但是多相机系统的复杂性较大，对不同空间位置的相机进行耦合同步比较困难。光场相机可以通过单台相机在一次曝光中采集四维光场信息[9]。此特点可以解决辐射测温多相机系统光路复杂、同步触发难等问题，在辐射成像三维温度重建时有其独特优势[10-15]。许传龙等[16]将光场成像技术与火焰温度重建结合，然后用LSQR对该数学模型进行求解，重构出火焰的温度分布。LSQR 方法是一种迭代正则化方法，适用于求解大规模不适定问题。但是在求解时无法保证结果的非负性[17]，而且在重建的准确性和效率方面仍有提升空间。孙俊等[17]提出了NNLS-LMBC 混合算法，可同时重建出火焰的温度场和吸收系数，NNLS 算法能够保证求解结果的非负性，且具有较好的计算稳定性。但是计算效率低，计算时间长。

LSMR 算法是一种用于求解稀疏最小二乘问题的迭代算法，其原理是将极小残差法用来求解线性方程组和最小平方问题，已经被应用于图像处理等方面的研究，证明了LSMR 算法在解决不适定问题上具有良好的性能[18-20]。针对LSQR 算法在火焰微元体辐射强度求解中出现负值和NNLS算法的计算效率低等问题，本文提出将LSMR 算法用于火焰光场成像三维温度场重建。构建单光场相机的辐射成像模型，从辐射强度重建结果的非负性、求解精度和求解效率三个方面分析LSMR 算法的求解性能。然后使用LSMR 算法对模拟光场火焰进行温度场重建，验证其进行火焰三维温度场重建的准确性。

1 测量原理与方法

1.1 光场成像

光场相机由主透镜、微透镜阵列和CCD 阵列组成。传统一代光场相机的成像如图1所示。光源点发射出不同方向的辐射光线，经过主透镜后到达微透镜，微透镜根据方向分离这些光线，并在微透镜后面的CCD 阵列上记录下辐射强度信息，同时光线的位置和方向信息可以通过光场相机中CCD阵列和微透镜阵列来确定[21]。

图1 光场相机成像

1.2 火焰辐射光场成像模型

将通过光场相机采集到的火焰辐射信息和重建算法相结合，进行火焰三维温度场的重建。火焰辐射光场图像的辐射强度计算采用的是稳态传输方程，为简化过程，仅考虑火焰的吸收特性[17]，火焰的稳态传输方程可以表示为式(1)。

式中，ILλ为火焰辐射光线L 的光谱辐射强度；κL为光线L 穿过火焰微元体的吸收系数；lL为该微元体的几何长度；Ibλ为火焰微元体的黑体辐射强度，W/(m3·sr)，根据普朗克定律，Ibλ表示为式(2)。

式中，T为微元体的温度，K；c1为第一辐射常数，c1=3.7418×10-16W·m2；c2为第二辐射常数，c2=1.4388×10-2m·K；λ为火焰辐射光线的波长。

将火焰划分为V个微元体，依次编号为(1,2, …,v, …,V)，光场相机的图像探测器共有M个像素接收到火焰辐射光线，依次编号为(1, 2, …,m, …,M)。像素m的火焰辐射光线穿过n个微元体，沿着光线的传输方向依次编号为(v1, v2, …,vi, …, vj, …, vn)，将式(1)离散化得到式(3)、式(4)。

1.3 LSMR算法

LSMR迭代算法可以用于求解线性系统Ax-b和最小平方问题min‖Ax-b‖2。LSMR 算法基于Golub-Kahan 双对角化过程，等同于将最小残差法应用于正规方程ATAx=ATb。

对于光场相机图像传感器各像素对应的火焰辐射光线的辐射强度向量Iλ，通过构造系数矩阵A，可建立线性方程Iλ=AIbλ。利用LSMR 算法求解线性方程Iλ=AIbλ的过程等价于求解线性系统Ax-b的过程。

对于线性系统Ax-b，LSMR算法的具体实现是通过Golub-Kahan过程实现双对角化，即式(6)。

选择合适的标量αk≥0和βk≥0，使得uk和vk的范数等于1。式(6)经过k次迭代可得式(7)。

e表示单位矩阵；ek+1表示单位矩阵中第k+1 列向量。结合式(7)可以得到式(8)。

对于正规方程ATAx=ATb，rk=b-Axk为其每一步的迭代残差，每一步迭代过程中通过选择yk寻求方程的近似解，如式(9)所示。

LSMR算法在每一步迭代中选择yk，使‖ATrk‖达到极小值，如式(10)所示。

而LSQR 算法的每一步迭代中选择yk，只使得‖rk‖最小化，如式(11)所示。

在线性系统的求解过程中，LSQR 算法[22]和LSMR 算法都可以得到最小范数解。在LSQR 算法的迭代过程中，只有残差范数‖rk‖单调递减；但是在LSMR算法中，残差范数‖rk‖和‖ATrk‖都是单调递减的。LSMR 算法求解得到的结果比LSQR算法更加精确和稳定。

相比于LSQR算法，LSMR算法求解速度更快，得到的结果更加安全。NNLS 算法[23]可以保证求解的非负性，但是其算法中存在矩阵求逆运算，而LSQR 和LSMR 算法没有这一过程，计算效率会显著提高。

2 数值模拟实验结果与分析

2.1 实验条件

为了对LSMR算法的重建性能进行研究，基于火焰辐射光场成像模型，进行了相关的数值模拟实验。火焰形状设定为底面半径R=20mm 和高度Z=40mm 的圆柱体，吸收系数设置为10m-1。火焰分布模型如式(12)和式(13)所示[17]。式中T为火焰的温度大小，K；z为火焰的轴向坐标；r为火焰的径向坐标。其温度分布满足旋转对称分布，该模拟火焰符合火焰的温度分布特点。火焰温度分布如图2所示，温度范围为900～2100K。

图2 火焰数值模拟模型

以0mm 的位置为火焰的中心层，每隔5mm 取一层火焰温度数据，将被测对象火焰分成7层，每层的大小为72×72。分别用LSQR、NNLS 和LSMR三种算法对火焰的辐射强度进行重建，为了定量分析不同算法的重建效果，在原始火焰辐射强度分布数据上加入1%、5%、10%、15%和20%高斯白噪声，再将重建结果与相对应的原始数据进行对比和分析。其中绝对误差和相对误差的定义如式(14)和式(15)所示。

式中，ΔIabs为重建的绝对误差；ΔIrel为重建结果的相对误差；I为使用算法重建出来的值；Iset为根据式(12)设定的火焰数据。

2.2 火焰辐射强度的非负性分析

在设定的火焰辐射强度分布中叠加1%、5%、10%、15%、20%的噪声，再用LSMR、LSQR 和NNLS 三种算法分别求解火焰的辐射强度。统计7层中所有重建结果中出现负值的数量并计算所占的百分比，结果见表1。

表1 辐射强度重建时出现负值点的数量及百分比

如表1所示，在噪声水平较低的情况下，例如1%噪声时，LSQR、NNLS 和LSMR 三种算法对火焰微元体辐射强度的重建都未出现负值。但是随着噪声的增大，LSQR 算法无法保证求解结果的非负性，对火焰微元体辐射强度求解的结果中会出现越来越多的负值。相比之下，LSMR 和NNLS 求解的结果中未出现负值，保证了求解的非负性。

LSQR 算法对辐射强度的求解会收敛至负值，主要集中在火焰的内焰部分。用LSQR、NNLS 和LSMR三种算法在不同噪声水平下对火焰的辐射强度进行重建后，选取中心切片层大小为21×21的内焰区域，对比和分析此区域的重建效果，如图3所示。

图3 中心层内焰区域辐射强度重建对比

图3(a)为该区域的原始辐射强度分布，图3(b)、(c)、(d)分别为LSQR、NNLS 和LSMR 算法的重建结果。可以看出，LSQR 算法在1%噪声时还呈现出较好的重建效果，当噪声增加到5%时，出现了少量的负值点。随着噪声的不断增大，其求解的结果中出现越来越多的负值解。在5%、10%、15%、20%情况下此区域出现负值解的百分比分别为4.3%、14.7%、28.6%、34.7%，中心层内焰区域负值解的百分比远高于所有层全部解的情况。LSMR 和NNLS 算法始终保持着求解结果的非负性。正因为此，在噪声逐渐增大的情况下，LSQR 算法重建的偏差比LSMR 和NNLS 要显著增大。LSMR 和NNLS 算法比LSQR 算法具有更好的求解稳定性。

2.3 火焰辐射强度的误差分析

计算LSMR、LSQR和NNLS三种算法在不同噪声水平下对火焰辐射强度重建的平均相对误差，如图4所示。

图4 LSQR、NNLS 和LSMR算法计算火焰辐射强度的平均相对误差

由图4 可以看出，在噪声很低的情况下，LSQR、NNLS 和LSMR 三种算法的平均相对误差十分接近，这是因为LSQR算法求解的结果中还未出现负值解。随着噪声水平的增加，三种算法的平均相对误差都在增大，但是由于LSQR算法求解结果中出现的负值解变多，平均相对误差也一直高于LSMR 和NNLS 算法。在噪声水平为5%、10%、15%、20%的情况下，LSMR算法重建的平均相对误差比LSQR 算法分别要低11.4%、14.5%、13.3%。LSMR和NNLS算法的平均相对误差都很接近，展现出比LSQR算法更好的抗噪性能。

2.4 火焰辐射强度的计算效率分析

计算效率是评估算法性能的另一个重要指标。在噪声水平为1%、5%、10%、15%、20%的情况下，分别统计LSQR、NNLS和LSMR三种算法多次计算的时间，并求得平均计算时间进行对比，如表2所示。

表2 不同噪声下计算时间的对比

在计算中，系数矩阵会因为所添加的噪声而产生偏差，使得方程组的病态程度加剧。所以在不同的噪声水平下，方程组的计算复杂度也不一样，计算的时间也会受到影响。随着噪声的增加，三种算法的求解时间都随着计算复杂度的变化而发生波动。如表2 所示，得到LSQR 的平均用时为2.90～3.05s，NNLS的平均用时为2201～2278s，LSMR的计算时间仅为0.15s 左右。可以看出，LSMR 算法的计算时间比NNLS 算法降低了4 个数量级。同时LSMR 算法的计算时间只有LSQR 算法用时的二十分之一，LSMR算法展现出在不同噪声水平下进行计算时高效的求解效率。

2.5 温度场重建结果分析

为了验证LSMR 算法对火焰温度场的重建效果，分别叠加1%、5%、10%、15%、20%的噪声，重建后的温度分布和绝对误差分布如图5 和图6所示。

图5 LSMR算法在不同噪声下温度的重建效果

图6 LSMR算法在不同噪声下温度重建的绝对误差分布

如图5 所示，在不同的噪声水平下，LSMR 算法重建出来的火焰温度分布均可以正确反映火焰的温度趋势特征。如图6所示，重建的误差主要分布在火焰的高温区，在1%、5%、10%、15%、20%噪声水平下，LSMR 重建的最大误差分别为4K、29K、68K、141K、215K，对于最高温度在2100K的高温区域，在合理的误差范围之内。

LSMR算法对火焰温度场进行重建时的平均相对误差和最大相对误差见表3，两者都随着叠加噪声的增加而增大，但是在不同的噪声水平下，平均相对误差和最大相对误差都保持在一个较低的水平。LSMR 算法能够准确地重建出火焰的三维温度场，表明了提出的LSMR算法是精确、合理、可靠的。

表3 LSMR算法温度重建的误差

3 结论

基于光场成像技术，采用LSMR方法对火焰的三维温度场重建进行了数值模拟研究，并与LSQR和NNLS方法对火焰辐射强度重建的效果进行了比较，主要结论如下。

（1）使用LSMR 方法对火焰辐射强度求解时，不会出现负值。因此，LSMR 求解的精度和NNLS相当，比LSQR有所提高。

（2）在计算时间方面，LSMR 的平均计算时间仅为0.15s 左右，低于LSQR 方法的2.90～3.05s 和NNLS方法的2201～2278s。

（3）在不同的噪声水平下，LSMR 都能够完成对火焰三维温度场的重建。并且在噪声水平较高的情况下，温度重建的平均相对误差能保持在1.2%以内。