基于正则代数重建算法的二维温度场重建

于和放, 胡小梅*

(1.上海大学 机电工程与自动化学院,上海 200444; 2.上海市智能制造及机器人重点实验室,上海 200444)

可调谐半导体激光吸收光谱断层诊断技术(Tunable Diode Laser Absorption Tomography,TDLAT)是可调谐二极管激光吸收光谱技术 (Tunable Diode Laser Absorption Spectroscopy,TDLAS)与计算机断层成像(Computed Tomography,CT)技术相结合的一种新型的流场参数测量技术,可以实现对非均匀温度场的测量。因其具有非侵入式、灵敏度高、抗噪声能力强、实时性好、测量环境适应性强等优势,目前,已广泛应用于燃烧流场的诊断研究[1-3]中。

在基于TDLAT的温度场重建中,重建算法对重建结果的精度起着关键作用。目前常用的重建算法主要包括解析类重建算法和迭代类重建算法。解析类重建算法在完备投影的条件下,可获得较好的重建结果[4],但是由于硬件条件的限制,往往无法实现完备投影。迭代类重建算法可以在迭代过程中加入先验信息,可以应用在不完备投影的温度场重建中[5-6]。夏晖晖等[7]以水汽作为气体吸收分子,使用代数重建算法(Algebraic Reconstruction Technique,ART)对燃烧温度场和浓度场分布进行了仿真重建和实验研究,但是网格划分数量较少,重建结果的分辨率较低。Jeon等[8-11]提出了3种改进的ART,并进行了仿真验证,给出了不同算法的适用范围,并在后续几年中将算法应用在实验研究中。Liu等[12]在扇形投影分布下使用改进的Landweber算法对燃烧温度场进行了仿真重建和实验研究,实验结果表明重建温度场与预期分布吻合较好,表明了该算法检测燃烧火焰的有效性。

上述算法本质上都是求最小二乘解或近似最小二乘解,在重建的过程中会引入噪声和离散误差。为了能够有效抑制噪声,平衡离散误差,本文在ART的基础上提出了基于留一交叉验证的正则ART,并使用典型温度场和典型投影分布对所提出的算法的重建精度和稳定性进行了仿真验证。

1 TDLAS测温和温度场重建原理

1.1 TDLAS测温原理

TDLAS的基本原理是激光穿过流场后的衰减程度不同,根据Beer-Lambert定律[13]可获得激光衰减前后光强的关系[14],即

(1)

式中：v为激光频率(cm-1);a(v)为激光在频率v处的吸光度;l为激光穿过待测区域的路径长度(cm);I0和I1分别为入射光强和出射光强;P为待测区域的压强(atm,1 atm=101 325 Pa);X为气体分子的浓度;T(l)为待测区域的温度(K);S(T) 为气体分子对激光的吸收强度(cm-2/atm);φ(v)为吸收谱线的线性函数。当温度小于2 500 K、激光波长小于 2.5 μm时,线强S(T)可表示为

(2)

式中：T0为参考温度296 K;Q(T)为配分函数;h为普朗克常数(J·S);c为光速(cm/s);E″为低状态能级能量(cm-1);k为玻尔兹曼常数(J/K)。以上这些参数可以通过查询HITRAN 光谱数据库获取。

TDLAS测温通常采用线强比值法,首先根据硬件条件选取2条合适的吸收谱线,然后可采用扫描频率波长直接吸收法获得激光在整个选定吸收谱线的波长范围的积分吸收面积：

(3)

线强比值法就是在同一条激光投影下的一个测量周期中,获取该气体的2个积分吸收面积。一次测量周期的时间很短,所以默认同一测量周期下,温度、浓度和压强不变,因此可以通过求解2个积分吸收面积的比值消除浓度和压强的影响。经过变换后可以获取待测区域的平均温度为

(4)

1.2 非均匀温度场重建原理

上述方法计算的温度结果为投影上气体温度的均值。而实际待测区域多为非均匀温度场,为了能获取非均匀温度场的温度分布,需要将TDLAS和CT技术相结合。TDLAT的原理是获取待测温度场多个方向的投影数据,然后使用相应的重建算法重建温度场。基于ART的温度场重建过程如下。

首先将待测区域离散成N个网格,并且认为每个网格内的温度、压力和浓度都是定值,离散化后积分吸收面积即投影值可以表示为

(5)

式中：i和j为投影和网格的编号;aj为第j个网格内的积分吸光度;Lij为第i条激光束在第j个网格内的路径长度。表示成矩阵的形式为

A=La

(6)

式中：A为投影向量;L为系数矩阵;a为网格积分吸光度向量,是未知量。

然后采用ART对上述方程组进行求解。其基本思想是：先给定初始解,在没有先验信息的条件下一般初始解为0;然后计算当前计算投影值和测量投影值的残差,沿着投影的传播路径进行反投影修正;之后不断迭代直到残差值满足重建精度要求为止。ART的特点是每次迭代修正时只使用一个投影,因此只对该投影穿过的网格单元值进行修正。使用第i条投影对第j个网格进行修正的迭代公式为

(7)

式中：k为迭代次数;λ为迭代步长。

求解出每个网格的积分吸光度后,再使用线强比值法求解出每个网格的平均温度,即可获取待测场的温度分布。

2 基于留一交叉验证的正则代数重建算法

2.1 正则代数重建算法

在实际测量投影值时,由于设备和环境的影响,不可避免地会存在一定的噪声。并且根据非均匀温度场的重建原理,在进行重建时,需要将待测区域划分成N个网格,因此基于网格计算出的投影值与测量投影值之间存在一定的离散误差。所以在使用ART对方程组进行求解时,离散误差和投影噪声会随着迭代修正而引入。为了减小重建算法对噪声的敏感度,平衡离散误差对结果带来的影响,本文对ART进行了改进,提出正则代数重建算法。

ART实质上求解的是近似最小二乘解,所以对噪声和离散误差比较敏感,一般的解决办法是在最小二乘准则的基础上增添正则化项,将极小化目标函数改为

(8)

式中：δ为正则化参数。

正则化项的作用就是对解做一些限制,使解更加平滑,δ的作用是控制平滑的程度。采用梯度下降法对目标函数进行求解,将负梯度方向作为迭代方向,可以推导出正则Landweber算法的迭代格式为

ak+1=ak+λ[LT(A-Lak)+δak]

(9)

由式(9)可知,正则Landweber算法对投影的噪声水平和离散误差进行统一衡量,但是在实际情况下,每条投影穿过的重建区域都不相同,因此每条投影的噪声水平和离散误差的大小也不一致。所以有必要对每条投影的噪声水平和离散误差分别进行衡量。ART的特点就是每次迭代修正时只使用一个投影,可以实现对投影噪声水平和离散误差的分别衡量。由式(9)可以推导出正则代数重建算法的格式为

(10)

一条投影穿过待测区域的路径越长,该投影路径上的温度梯度大概率越大,即离散误差越大。所以每条投影的正则化权重由投影穿过的待测区域路径的长度和单位正则化参数共同决定,因此正则代数重建算法的迭代公式可更改为

(11)

式中：μ为单位正则化权参数;li为第i条投影所穿过的待测区域的路径长度。

正则代数重建算法的求解质量依赖于单位正则化参数的选择,当正则化项的权重过小时不能很好地抑制噪声和离散误差,当正则化项权重较大时,会使得重建结果过于平滑。为了选取合适的正则化参数,本文采用留一交叉验证法对正则化参数进行选取。

2.2 基于留一交叉验证的正则化参数选取

交叉验证也称为循环估计,是一种可以评估参数的统计学方法,交叉验证法评估参数的基本思想是将原始样本数据集分割成2个子集,其中一个子集数据用于求解,另一个子集的数据用于验证。对于本文来讲,首先将投影数据分成2个数据集,设置不同的正则化参数,使用其中一个数据集重建温度场,然后使用重建后的结果对另一个数据集的投影值进行计算,将计算投影值和测量投影值进行比较,获取残差值,选取使残差值最小的正则化参数。

由于投影数量较少,在使用普通交叉验证法选取正则化参数时,根据数据集拆分方式的不同,重建结果的差异性较大,使得求解的正则化参数不稳定。

为了减少普通交叉验证方法的不稳定性,可以将数据集划分为m个子集,每次取一个集合用于验证,剩余m-1个集合用于求解,总共须进行m次交叉验证,对这m次交叉验证结果进行累加再取平均值作为最终的验证结果。这种方法称为m折交叉验证法。为了保证在交叉验证中求解的结果与使用整个数据集求解的结果最为接近,本文在划分数据集时,使m等于数据集中投影数据的个数,即每次交叉验证时,验证集只有1个数据用于验证,剩余的数据用于重建温度场,这种方法被称为留一交叉验证法。该方法的极小化目标函数为

(12)

虽然留一交叉验证法所需的计算成本较大,但可以最大限度地利用原始投影数据选出合适的正则化参数。

至此,对正则化参数的选取问题变成了目标函数为单变量函数且函数曲线为单峰的一维最优化问题,但是目标函数无法求导,所以本文采用黄金分割法求解近似最优值,其基本思想是首先确定最优点所在的区间,以固定区间缩短率0.618不断缩小最优点所在的区间范围,设置相应的阈值,区间长度小于阈值时,停止分割,取区间的中点作为近似最优点,将此时正则化权重值作为近似最优正则化参数。

3 仿真验证

基于上述温度场重建原理,应用MATLAB软件编写了仿真重建程序,具体输入参数和仿真重建结果如下。

3.1 仿真数据的设置

3.1.1 仿真温度场的建立

仿真模型采用高斯单峰对称和高斯单峰偏置这2种温度模型,如图1所示,温度范围为500～1 000 K。

图1 仿真温度场模型

3.1.2 投影分布设置

为了验证算法的普适性,分别在扇形投影和平行投影[16]的前提条件下对温度场进行仿真重建。扇形投影分布和平行投影分布的示意图如图2所示,激光器从5个角度发射激光穿过待测区域,每个角度被6个间距相等的探测器探测,共可获取30条有效投影数据。

图2 投影分布示意图

3.1.3 待测区域设置和网格划分

待测区域大小设置为100 mm×100 mm,网格划分为19×19,如图3所示。

图3 待测区域设置和网格划分示意图

从图3可以看出,有部分网格没有被投影穿过,而采用正则代数重建算法进行迭代修正的过程中,只会针对有投影穿过的网格进行修正,为求解没有投影穿过的网格的温度值,本文在迭代的过程中加入了高斯滤波。

3.1.4 气体吸收分子的选取

H2O作为主要的燃烧产物,在1～3 μm波段有大量吸收,并且该波段的激光器较为成熟,所以H2O常作为测量流场温度的媒介,因此在进行仿真温度场重建时采用H2O作为测量媒介,并默认待测区域H2O的浓度为均匀分布。

3.1.5 吸收谱线的选取

在选取吸收谱线时首先要考虑硬件条件,选择的谱线要有利于测量,尽量减少测量误差,并且要保证线强比值R和温度值T是一一对应的关系,即线强比是温度的单调函数。其次要考虑所选择的谱线组对待测区域的温度区间有较高的测温灵敏度[15]。根据仿真温度场的设置,要保证谱线组在500～1 000 K的温度范围内有较高的测温灵敏度。本文仿真测试所使用的一组吸收谱线的参数如表1所示。

表1 谱线参数

3.2 仿真结果及误差分析

为表征温度场的重建精度,对重建结果计算平均相对误差Rave和最大相对误差RE,max。

(13)

(14)

式中：Tm(j)和T(j)分别为模型温度场和重建温度第j个网格的的温度值。

在常压条件下,基于上述参数分别采用ART、正则Landweber算法和正则代数重建算法对温度场进行重建,重建结果的误差如表2和表3所示。

表2 扇形投影分布

表3 平行投影分布

从表2和表3中的数据可以看出,对上述2种温度场而言,相比于ART,正则Landweber算法和正则代数重建算法的重建精度都有一定的提高。其中正则Landweber算法对投影的噪声水平和离散误差进行了统一衡量,优化效果并不明显;正则代数重建算法根据投影穿过待测区域路径的长度和单位正则化参数,动态地调整每条投影的正则化权重,实现了对每条投影的噪声水平和离散误差的分别衡量,使用该算法在扇形投影和平行投影下对上述2种温度场的重建精度都有较明显的提高。

在计算机处理器为AMD Ryzen 7 6800H的硬件条件下,基于正则代数重建算法的温度场重建平均耗时为110.4 s,其中正则化参数的选取耗时占比约为99.68%,原因是采用了留一交叉验证法,计算量较大。

4 结束语

本文依据最优化理论,对ART进行了改进,提出了正则代数重建算法,并使用该算法对单峰对称和单峰偏置温度场进行了仿真重建。在不同投影分布情况下的重建结果表明,相比ART和正则Landweber算法,该算法的重建精度有较为明显的提高,并且稳定性较好。但是该算法的正则化参数选取采用了留一交叉验证法,该方法计算成本较大,不能满足实时重建的要求,因此须进一步研究快速选择正则化参数的方法。