CN108225577B - 基于双光谱辐射信息的火焰三维温度分布重建方法 - Google Patents

Abstract

本发明公开了一种基于双光谱辐射信息的火焰三维温度分布重建方法，其特征在于，包括如下步骤：1)、根据采集的火焰双光谱辐射信息，划分火焰微元体，建立火焰辐射传输方程；2)、设定目标函数，根据重建算法迭代更新火焰吸收系数，重建火焰三维温度场。本发明在算法的迭代求解过程中，使用了非负最小二乘算法，能够保证火焰各部分光谱辐射强度的非负性。同时，还将更新的火焰各部分吸收系数映射到设定的区间，保证吸收系数的非负性。算法具有较强的鲁棒性，对初值无依赖性，在输入的数据误差较大，以及无法提供先验的迭代初值而设定随机数值的情况下，算法仍能保持较高的准确度和较好的收敛性。

Description

基于双光谱辐射信息的火焰三维温度分布重建方法

技术领域

本发明涉及一种基于双光谱辐射信息的火焰三维温度分布重建方法，属于火焰温度测量技术领域。

背景技术

温度是表征燃烧火焰一个最主要的参数。火焰三维温度场的测量，对全面研析火焰的燃烧反应速度、三维结构、组分生成、整体特性等有决定性作用，有助于揭示燃烧现象的本质和燃烧过程的规律，也有助于燃烧设备的工程设计改进和运行优化，最终对能源的高效利用和低污染排放产生重大意义。因此，火焰三维温度场的测量是燃烧火焰研究中的一个重要主题。

激光光谱等温度测量技术存在设备使用环境要求高、价格昂贵、操作复杂等缺点，目前还极少用于工业现场的温度测量。图像法火焰温度测量技术，利用图像探测器，拍摄火焰自发辐射图像，探测火焰自身辐射强度进行火焰温度测量，该方法对于工业环境下火焰的三维温度测量具有很大的优势。

辐射图像法主要包括光学分层层析、成像法光学层析法和辐射传输求解法，其中光学层析法，利用多个相机拍摄不同视角的火焰图像，将火焰图像视为空间三维火焰在图像探测器上的投影，利用逻辑滤波反投影(logical filtered back-projection，LFBP)算法和迭代重建和代数重建联合算法(Simultaneous iterative reconstruction andAlgebraic Reconstruction Techniques，SART)重建火焰各横截面的灰度分布，根据双色法计算各截面的温度分布。该研究将传统的双色法火焰二维温度测量技术与光学分层层析技术结合，实现了火焰的三维温度测量，具有重要实践意义。光学分层层析成像法，将火焰视为若干个互相平行的二维断层的组合，对不同断层分别快速聚焦摄像，得到一组辐射图像，采用图像反演算法，重建各断层的原始图像，再根据图像灰度和火焰温度之间的关系，建立火焰的三维温度场。该方法可以利用单相机，快速重建火焰三维温度场。然而，光学层析法和光学分层层析成像法，均未考虑火焰内部光线的实际传播路径和吸收系数分布的影响，难以精确测量存在一定浓度的炭黑等颗粒的火焰温度。

辐射传输求解法，拍摄火焰图像，综合考虑火焰内部的吸收、散射作用，建立关联火焰各部分黑体辐射强度和图像探测器接收到的辐射强度的辐射传输方程，通过一定的反演算法[如吉洪诺夫正则化方法(Tikhonov regularization method)]求解该方程，得到火焰各部分的黑体辐射强度，根据辐射定律计算火焰各部分的温度，实现火焰三维温度分布的重建，主要包括DRESOR(Distribution ofRatios ofEnergy Scattered by the mediumOr Reflected by the boundary surface)法和广义源项多流法[Generalized sourcedmulti-flux method(GSMFM)]。DRESOR法，针对电站锅炉炉膛火焰，综合考虑了炉内燃烧介质的发射散射以及炉墙和水冷壁的对辐射能的发射和反射作用，沿方向

投射到相机图像探测器上的火焰辐射强度通过计算两种DRESOR数来确定的，一种是表示a点为中心的单位体积内发射的能量在以沿方向

方向为中心的单位立体角内直接投射到以b点为中心的单位面积区域的能量份额，此处“直接”是指这部分的发射经过吸收、散射后的剩余部分，另一种是

a点为中心的单位体积内发射的能量在以沿方向

方向为中心线的单位立体角内被以b点为中心的单位体积散射的能量份额与4π的乘积，辐射强度与炉膛内各部分的黑体辐射强度的满足线性关系，根据改进的吉洪诺夫正则化方法求解该线性模型，从而得到炉膛内火焰三维温度分布。该方法的提出，一定程度上解决了炉膛内大尺寸火焰多物体(火焰燃烧弥散介质、炉墙、水冷壁等)复杂辐射传输过程(发射、吸收和散射)的求解问题，使得辐射图像法能够运用到恶劣环境下的火焰温度测量，而且DRESOR法对于不同条件(有无高温壁面、散射的各项同性或各向异性等)的辐射传输计算具有较好的适用性。然而，该方法是基于Monte Carlo抽样技术实现的，所需计算次数和时间长，对于具有剧烈变化特性的燃烧火焰的监测是不利的。

广义源项多流法，将计算火焰各方向辐射强度的辐射传输微分方程进行差分，提取各单元的广义源项S，建立投射到图像探测器上各方向辐射强度与火焰各部分的黑体辐射强度的线性方程，利用最小二乘QR分解(1east-square QR decomposition)和共轭梯度(conjugate gradient)混合算法(LSQR-CG)等反演算法求解该线性关系，重建火焰三维温度分布。该方法不依赖蒙特卡洛抽样技术，相比于DRESOR法，具有更高的计算效率。然而，目前所提出的反演算法如LSQR-CG等对求解结果的非负性没有约束力，在求解实际火焰的辐射传输线性方程时，由于测量系统(主要为数码相机)获得的火焰辐射存在较大噪声，求解结果易陷入局部最优甚至产生负值解。火焰的辐射强度和温度值均具有非负性，因而在火焰传输方程的反演算法方面，需要进一步的深入研究。

发明内容

技术问题：根据相机拍摄的火焰辐射图像可以测量火焰温度场，本发明旨在提供一种有效的求解该辐射传输反演问题的重建算法，从而计算获得火焰的三维温度分布。

技术方案：

一种基于双光谱辐射信息的火焰三维温度分布重建方法，其特征在于，包括如下步骤：

1)、根据采集的火焰双光谱辐射信息，划分火焰微元体，建立火焰辐射传输方程；

2)、设定目标函数，根据重建算法迭代更新火焰吸收系数，重建火焰三维温度场。

所述步骤1)包括：

11)、输入采集的火焰辐射信息并进行预处理：

输入火焰不同探测方向的出射光谱辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)；

其中，(x，y，z)为探测线起点坐标，(θ，ψ)为探测线的方向坐标，θ为天顶角，ψ为天顶角，λ(R)和λ(G)表示的两种波长；

根据任一探测线方向的出射光谱辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)建立方程式：

式中，m为探测线的序号，I_λ为沿探测线方向的光谱辐射强度，I_bλ为火焰微元体的自身辐射强度，n为当前探测线穿过的火焰微元体总个数，L为当前探测线穿过的火焰微元体的几何长度，κ为当前探测线穿过的火焰微元体的吸收系数，vn为沿火焰出射方向第n个微元体的编号，i，j分别为当前探测线沿其火焰出射方向穿过的第i，j个微元体；

12)、建立辐射传输方程组：

联立所有探测线方向的光谱出射辐射强度，得到线性方程组式：

I_λ＝A·I_bλ (2)

式中，A为计算系数，I_λ为火焰各方向的光谱出射辐射强度组成的向量[I_λ ¹，I_λ ²，…，I_λ ^m，…，I_λ ^M]，I_bλ为火焰微元体的自身辐射强度组成的向量[I_bλ ¹，I_bλ ²，…，I_bλ ^v，…I_bλ ^V]，A为计算系数A组成的系数矩阵，[A₁ ¹，A₁ ²，...，A₁ ^V；A₂ ¹，A₂ ²，...，A₂ ^V；...；A_M ¹，A_M ²，...，A_M ^V]。

所述步骤2)包括：

21)、设置算法所需的中间参数的初值，吸收系数向量κ＝κ₀，吸收系数的初值向量κ₀可以选择随机向量或单位向量，按照下列各式依次计算对称矩阵Ω，测量值和正向计算值得差值向量r，目标函数的梯度向量g，循环终止指数S_tp和衰减系数μ，

Ω＝J^TJ (3)

r＝I′_λ(G)-I_λ(G) (4)

g＝J^Tr (5)

S_tp＝(||g||≤ε) (6)

μ＝η·max_{i＝1，...，m}(Ω_ii) (7)

式中，J为函数I′_λ(G)(κ)关于变量κ的一阶偏微分雅各比(Jacobi)矩阵，Ω_ii为矩阵Ω的对角元素，ε是判定迭代过程是否终止的阈值，η为比例系数；

22)计算I′_λ(G)(κ)估计值：

利用非负最小二乘算法求解式(2)，得到火焰各微元体的R波长的光谱黑体辐射强度I_bλ(R)；利用光谱黑体辐射强度I_bλ(R)，根据式(8)计算火焰各微元体的温度值T，

T＝c₂/λ(R)ln{c₁/[λ(R)⁵πI_bλ(R)]+1} (9)

式中，c₁为第一辐射常数，c₂为第二辐射常数，λ(R)为火焰辐射光线的R波长；利用此温度值T，根据式(10)计算火焰各微元体的G波长的光谱黑体辐射强I_bλ(G)：

式中，λ(G)为火焰辐射光线的G波长；再利用当前的系数矩阵A和光谱黑体辐射强度向量I_bλ(G)，根据式(2)计算火焰单色(G)辐射强度向量I′_λ(G)(κ)；

23)、计算迭代步长Δ(κ)：

计算当前吸收系数κ下的一阶偏微分雅各比矩阵J，根据式(3)-(5)计算当前的对称矩阵Ω，测量值和正向计算值的差值向量r及目标函数的梯度向量g：

根据当前的对称矩阵Ω，测量值和正向计算值的差值向量r及目标函数的梯度向量g，求解式(11)得到算法的优化步长Δ(κ)：

(Ω+μI)Δ(κ)＝g (11)

式中，I为单位矩阵；

24)、判断迭代收敛程度：

根据式(12)判断算法的优化步长Δ(κ)是否小于阈值ε：

||Δ(κ)||≤ε||κ|| (12)

如果优化步长Δ(κ)小于阈值ε，循环终止指数S_tp＝true(1)，并且转到步骤28)；如果不小于，继续进行下一步骤；

25)、对更新的参数进行边界约束：

将算法的优化步长Δ(κ)加到吸收系数向量κ上，得到新的吸收系数向量κ_new；吸收系数的初值向量κ₀选择随机向量或单位向量；并且根据式(13)和(14)将新的吸收系数向量κ_new投影到约束区间Q：

P_Q(κ)＝max{min{κ，u}，l} (13)

Q＝{κ∈R^m|l_i≤κ_i≤u_i，i＝1，...，m} (14)

式中，P_Q(κ)为投影函数，R^m为元素个数为m的实数集向量，l([l₁，l₂，...，l_m])为吸收系数约束区间Q下边界，取为([0，0，...，0])，u([u₁，u₂，...，u_m])为吸收系数约束区间Q上边界，取为([1×10⁸，1×10⁸，...，1×10⁸])；

26)、判断迭代指数ρ：

根据式(15)计算迭代指数ρ，

式中，Δ即优化步长Δ(κ)，并判断迭代指数ρ是否大于0，如果大于，按照式(16)和(17)更新衰减系数μ和中间参数υ：

μ＝μυ (16)

υ＝2υ (17)

并且转到步骤28)，否则，继续进行下一步骤；

27)、更新迭代参数：

利用新的吸收系数向量κ_new更新原吸收系数向量κ，利用更新的吸收系数向量κ，根据式(3)-(5)计算当前的Ω，r和g，根据式(18)-(20)计算υ，S_tp和衰减系数μ，

υ＝2 (18)

S_tp＝(||g||_∞≤ε)or(||r||≤ε) (19)

μ＝μ·max[1/3，1-(2ρ-1)³] (20)

28)、判断循环是否终止：

判断循环终止指数S_tp是否等于0并且迭代次数k是否小于最大迭代次数k_max，如果是，循环到步骤22)；否则，终止循环，并且进行到下一步骤；

29)、计算温度求解结果：

根据式(9)计算火焰各微元体的温度值T，输出该温度值及迭代终止的吸收系数值，实现火焰三维温度及吸收系数分布的同时重建。

J为函数I′_λ(G)(κ)关于吸收系数向量κ的一阶偏微分雅各比矩阵，Ω_ii为矩阵Ω的对角元素，ε是判定迭代过程是否终止的阈值。

第一辐射常数c₁为3.7418×10¹⁶W·m³，第二辐射常数c₂为1.4388×10^-2m·K。

有益效果：

1.本发明的算法结合了非负约束最小二乘算法NNLS和含边界约束的LMBC算法，从而对火焰各部分黑体光谱辐射强度以及物性参数的重建结果的边界具有约束力，保证解的非负性和大小的合理性。

2.本发明结合了NNLS和LMBC两种算法，能够利用采集的火焰各方向双光谱辐射信息重建火焰三维温度分布，其中LMBC算法具有较强的全局收敛性，对初值无依赖性。因而提出的LMBC-NNLS混合算法具有较强的鲁棒性，在输入的数据误差较大，以及无法提供先验的迭代初值而设定随机数值的情况下，算法仍能保持较高的准确度和较好的收敛性。

附图说明

图1是基于双光谱辐射的火焰三维温度分布重建算法流程图；

图2是增加非负性约束与未增加非负性约束的微元体黑体光谱辐射强度重建结果；

图3是增加非负性约束与未增加非负性约束的微元体吸收系数重建结果；

图4是目标函数的收敛结果；

图5是Case 1火焰的温度及吸收系数分布重建结果的相对误差Δκ和ΔT；

图6是Case 2火焰的温度及吸收系数分布重建结果的相对误差Δκ和ΔT。

具体实施方式

下面，结合附图和具体实施例对本发明作进一步说明：

两种(Case 1和Case 2)不同的温度和吸收系数分布情况的火焰，Case 1的温度T(K)和吸收系数κ(m^-1)为对称分布，分别满足式(1)和(3)，Case 2的温度T(K)和吸收系数κ(m^-1)为对称分布，分别满足式(2)和(4)。火焰的高度Z和底面半径R分别设置为30mm和8mm。

式中，r，x，y和z分别为径向、x轴方向、y轴方向和z轴方向的坐标值，

利用光场相机采集火焰不同方向的辐射强度，光场相机的参数如下，主透镜和微透镜的焦距分别为50mm和600μm，相机主透镜面与微透镜阵列面的距离为53.1mm，微透镜阵列面与图像探测器面的距离为480μm，微透镜阵列面微透镜的排列方式为正四边形排列，微透镜横向(H)和纵向(V)的个数均为60个，各微透镜的直径为95μm，各微透镜覆盖的像素数为12(H)×12(V)，各像素的尺寸为8μm(H)×8μm(V)，光场相机的图像探测器为彩色CCD(Charge Coupled Device)，各像素均能接收火焰的RGB三通道的光谱辐射强度值，火焰RGB通道的波长分别为610nm，530nm，和460nm。

火焰中心处的纵截面与光场相机的主透镜面平行，二者之间的距离为505mm，火焰中心点与主透镜中心点位于同一高度。

利用下式计算火焰各部分的光谱辐射强度，

式中，c₁为第一辐射常数，3.7418×10^-16W·m²，c₂为第二辐射常数，1.4388×10^{- 2}m·K，λ为火焰辐射光线的波长。各像素对应的火焰探测线方向为经过微透镜中心与像素中心的辐射光线，利用辐射传输方程式(6)计算光场相机图像探测面各像素接收的火焰光谱辐射强度I_λ(R)和I_λ(G)。

式中，m为探测线的序号，I_λ为沿探测线方向的光谱辐射强度，[W/(m³·sr)]，I_bλ为火焰微元体的自身辐射强度，[W/(m³·sr)]，n为当前探测线穿过的火焰微元体总个数，L为当前探测线穿过的火焰微元体的几何长度，κ为当前探测线穿过的火焰微元体的吸收系数，A为计算系数，i，j分别为当前探测线沿其火焰出射方向穿过的第i，j个微元体。获得火焰各方向光谱出射辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)，按照式(7)向此强度信号中叠加20dB的高斯白噪声。

式中，SNR表示信噪比(Signal to Noise Ratio)，为光场相机图像探测器上第i个像素的光谱辐射强度，N_i为图像探测器上第i个像素光谱辐射强度的噪声值，k为图像探测器上像素的个数。

下面利用此叠加噪声作为测量误差的火焰各方向光谱出射辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)，根据本发明提出的算法重建火焰的三维温度分布。

基于双光谱辐射的火焰三维温度分布重建算法，图1为该重建算法的详细流程图，包括如下步骤：

步骤一、输入采集的火焰辐射信息并进行预处理。以火焰底面中心为坐标原点建立笛卡尔右手坐标系统，输入火焰不同探测方向的出射光谱辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)，其中任一探测方向可以视为三维空间中的一根几何光线，本发明中称为探测线。(x，y，z)为探测线起点坐标，(θ，ψ)为探测线的方向坐标，θ为天顶角，ψ为天顶角，λ(R)和λ(G)分别表示RGB三通道滤波片的彩色CCD(Charge Coupled Device)R通道和G通道的波长，分别为610nm，530nm。

将火焰沿周向、径向和轴向划分为4×4×4(V＝64)个微元体，并按一定排列方式将各微元体编号(1，2，v，...，V)，保证每个微元体至少有一根探测线穿过。将各探测线编号(1，2，m，...，M)，任一探测线穿过的火焰微元体总数记为n，则沿火焰出射方向将各微元体的编号依次记为v₁，v₂，vi，(vj)，…，vn。根据任一探测线方向的出射光谱辐射强度信号I_λ(R)(x，y，z，θ，ψ)和I_λ(G)(x，y，z，θ，ψ)建立方程式(1)，联立所有探测线方向的光谱出射辐射强度，得到线性方程组式(2)。

I_λ＝A·I_bλ (2)

式中，m为探测线的序号，I_λ为沿探测线方向的光谱辐射强度，[W/(m³·sr)]，I_bλ为火焰微元体的自身辐射强度，[W/(m³·sr)]，n为当前探测线穿过的火焰微元体总个数，L为当前探测线穿过的火焰微元体的几何长度，κ为当前探测线穿过的火焰微元体的吸收系数，A为计算系数，i，j分别为当前探测线沿其火焰出射方向穿过的第i，j个微元体，I_λ为火焰各方向的光谱出射辐射强度组成的向量[I_λ ¹，I_λ ²，…，I_λ ^m，…，I_λ ^M]，I_bλ为火焰微元体的自身辐射强度组成的向量[I_bλ ¹，I_bλ ²，…，I_bλ ^v，…I_bλ ^V]，A为计算系数A组成的系数矩阵，[A₁ ¹，A₁ ²，...，A₁ ^V；A₂ ¹，A₂ ²，...，A₂ ^V；...；A_M ¹，A_M ²，...，A_M ^V].

步骤二、进行初始化设置：

目标函数设置为：

F(T，κ)＝||I′_λ(G)-I_λ(G)|| (3)

式中，I_λ(G)为火焰各方向的G波长的光谱辐射强度I_λ(G)组成的向量，I′_λ(G)为火焰各方向的G波长的光谱辐射强度向量。

设置算法所需的中间参数的初值，迭代次数k＝0；最大迭代次数k_max＝200；v＝2；η＝2；吸收系数的初值向量κ₀选择单位向量，按照式(4)-(8)依次计算对称矩阵Ω，测量值和正向计算值的差值向量r，目标函数的梯度向量g，循环终止指数S_tp和衰减系数μ，

Ω＝J^TJ (4)

r＝I’_λ(G)-I_λ(G) (5)

g＝J^T _r (6)

S_tp＝(||g||≤ε) (7)

μ＝η·max_{i＝1，...，m}(Ω_ii) (8)

式中，J为函数I′_λ(G)(κ)关于变量κ的一阶偏微分雅各比(Jacobi)矩阵：Ω_ii为矩阵Ω的对角元素；ε是阈值，用于判定迭代过程是否终止的，设置为1×10^-8，比例系数η设置为2。

步骤三、计算I′_λ(G)(κ)估计值。

根据式(1)计算当前吸收系数κ下的系数矩阵A，根据计算得到的吸收系数矩阵A和输入的R波长的光谱辐射强度向量I_λ(R)，利用Matlab的lsqnonneg函数实现NNLS算法来求解式(2)，得到火焰各微元体的R波长的光谱黑体辐射强度I_bλ(R)。利用光谱黑体辐射强度I_bλ(R)，根据式(9)计算火焰各微元体的温度值T，

T＝c₂/λ(R)ln{c₁/[λ(R)⁵πI_bλ(R)]+1} (9)

式中，c₁为第一辐射常数，3.7418×10^-1 ₆W·m²，c₂为第二辐射常数，1.4388×10^{- 2}m·K，λ(R)为火焰辐射光线的R波长。利用温度值T，根据式(10)计算火焰各微元体的G波长的光谱黑体辐射强I_bλ(G)，再利用当前的系数矩阵A和火焰各微元体的G波长的光谱黑体辐射强度向量I_bλ(G)，根据式(2)计算火焰单色(G)辐射强度向量I′_λ(G)(κ)。

式中，c₁为第一辐射常数，3.7418×10^-16W·m²，c₂为第二辐射常数，1.4388×10^{- 2}m·K，λ(G)为火焰辐射光线的G波长。

步骤四、计算迭代步长Δ(κ)。

计算当前吸收系数κ下的一阶偏微分雅各比矩阵J，根据式(4)-(6)计算当前的Ω，r和g，求解式(11)得到算法的优化步长Δ(κ)，

(Ω+μI)Δ(κ)＝g (11)

式中，I为单位矩阵。

步骤五、判断迭代收敛程度。

根据下式判断算法的优化步长Δ(κ)是否小于阈值ε(1×10^-8)，如果小于，循环终止指数S_tp＝true(1)，并且转到步骤九，如果不小于，继续进行下一步骤。

||Δ(κ)||≤ε||κ|| (12)

步骤六、对更新的参数进行边界约束。

将算法的优化步长Δ(κ)加到吸收系数向量κ上，得到新的吸收系数向量κ_new，并且根据式(13)将新的吸收系数向量κ_new投影到约束区间Q[如式(14)所示]。

P_Q(κ)＝max{min{κ，u}，l} (13)

Q＝{κ∈R^m|l_i≤κ_i≤u_i，i＝1，...，m} (14)

式中，P_Q(κ)为投影函数，R^m为元素个数为m的实数集向量，l([l₁，l₂，...，l_m])为吸收系数约束区间Q下边界，取为([0，0，...，0])，u([u₁，u₂，...，u_m])为吸收系数约束区间Q上边界，取为([1×10⁸，1×10⁸，...，1×10⁸])。

步骤七、判断迭代指数ρ。

根据下式计算迭代指数ρ，

Δ即优化步长Δ(κ)，并判断迭代指数ρ是否大于0，如果大于，按照下式更新中间参数μ和υ：

μ＝μυ (16)

υ＝2υ (17)

并且转到步骤九，否则，继续进行下一步骤。

步骤八、更新迭代参数。

利用新的吸收系数向量κ_new更新原吸收系数向量κ，利用更新的吸收系数向量κ，根据式(4)-(6)计算当前的Ω，r和g，根据式(18)-(20)计算υ，S_tp和μ，

υ＝2 (18)

S_tp＝(||g||_∞≤ε)or(||r||≤ε) (19)

μ＝μ·max[1/3,1-(2ρ-1)³] (20)

步骤九、判断循环是否终止。

判断循环终止指数S_tp是否等于0并且迭代次数k是否小于最大迭代次数k_max，如果是，循环到步骤三，否则，终止循环，并且进行到下一步骤。

步骤十、计算温度求解结果。

利用利用Matlab的lsqnonneg函数实现NNLS算法来求解式(2)，得到火焰各微元体的R波长的光谱黑体辐射强度I_bλ(R)，根据式(9)计算火焰各微元体的温度值T。

在火焰各方向光谱出射辐射强度信号I_λ(R)的基础上，根据设定的吸收系数初值κ₀，分别利用无非负性约束的现有的LSQR算法和本发明使用的含非负性约束的NNLS算法求解式(2)，结果如图2所示。可以看出，LSQR算法的求解结果中含有较多负值(如微元体6、9和14等)，而LSQR算法的求解结果不含负值，并且其求解结果更接近火焰微元体光谱黑体辐射强度I_bλ(R)的设定值。在火焰各方向光谱出射辐射强度信号I_λ(R)和I_λ(G)的基础上，根据本发明的双光谱辐射的火焰三维温度分布重建算法，分别增加约束(LMBC)和不增加约束(LM)以计算火焰微元体的吸收系数，结果如图3所示。可以看出，LM的求解结果中含有负值(如微元体8)，而LMBC的求解结果不含负值，并且其求解结果更接近火焰微元体吸收系数的设定值。表明了本发明对于火焰各部分黑体光谱辐射强度以及物性参数非负性约束方面的有益效果。

按照式(21)和(22)计算重建结果的相对误差，

式中，κ_rst为重建的吸收系数值，κ_set为设定的吸收系数值，T_rst为重建的火焰温度值，T_set为设定的火焰温度值。

图4是目标函数的收敛结果，两种温度及吸收系数分布情况的火焰(Case 1和Case2)的目标函数收敛到相同的值，该目标收敛值小于8，表明了在20dB的噪声下，本发明的重建算法依然拥有较好的收敛性。

图5是Case 1火焰温度及吸收系数分布重建结果的相对误差Δκ和ΔT。图6是Case2火焰温度及吸收系数分布重建结果的相对误差Δκ和ΔT。两种情况下的相对误差Δκ的最大值均小于10^-2，而相对误差ΔT的最大值均小于10^-3，表明在20dB的噪声下，火焰温度和吸收系数的重建结果仍具有较高的精确性。

火焰温度重建的精确度高于吸收系数，这是由于火焰的光谱辐射强度对吸收系数的变化不灵敏。火焰情况Case 2的相对误差Δκ(或ΔT)大于火焰情况Case 1的相对误差Δκ(或ΔT)，即吸收系数和温度分布不对称的的火焰重建结果差于对称分布情况下的火焰，表明火焰的温度及吸收系数分布的复杂程度将增加重建结果的误差。

两种火焰情况下的重建结果表明本发明的重建算法对于火焰温度的重建具有较好的精确性和适用性。

Claims (3)

1.一种基于双光谱辐射信息的火焰三维温度分布重建方法，其特征在于，包括如下步骤：

1)、根据采集的火焰双光谱辐射信息，划分火焰微元体，建立火焰辐射传输方程；

2)、设定目标函数，根据重建算法迭代更新火焰吸收系数，重建火焰三维温度场；所述步骤1)包括：

11)、输入采集的火焰辐射信息并进行预处理：

输入火焰不同探测方向的出射光谱辐射强度信号I_λ(R)(x,y,z,θ,ψ)和I_λ(G)(x,y,z,θ,ψ)；

其中，(x,y,z)为探测线起点坐标，(θ,ψ)为探测线的方向坐标，θ为天顶角，ψ为圆周角，λ(R)为火焰辐射光线的R波长，λ(G)为火焰辐射光线的G波长；

根据任一探测线方向的出射光谱辐射强度信号I_λ(R)(x,y,z,θ,ψ)和I_λ(G)(x,y,z,θ,ψ)建立方程式：

式中，m为探测线的序号，I_λ ^m为沿序号m的探测线方向的光谱辐射强度，

为沿火焰出射方向第n个微元体的自身辐射强度，n为当前探测线穿过的火焰微元体总个数，L_vj为当前探测线穿过编号为j的火焰微元体的几何长度，κ_vj为当前探测线穿过编号为j的火焰微元体的吸收系数，vn为沿火焰出射方向第n个微元体的编号，i,j分别为当前探测线沿其火焰出射方向穿过的第i,j个微元体，

为沿火焰出射方向第i个火焰微元体的自身辐射强度；

12)、建立辐射传输方程组：

联立所有探测线方向的光谱出射辐射强度，得到线性方程组式：

I_λ＝A·I_bλ (2)

式中，I_λ为火焰各方向的光谱出射辐射强度组成的向量[I_λ ¹,I_λ ²,…,I_λ ^m,…,I_λ ^M]，I_λ ^m为沿探测线方向序号为m的探测线的光谱辐射强度；I_bλ为火焰微元体的自身辐射强度组成的向量[I_bλ ¹,I_bλ ²,…,I_bλ ^v,…I_bλ ^V]，I_bλ ^v为火焰微元体v的自身辐射强度；A为计算系数A组成的系数矩阵[A₁ ¹,A₁ ²,…,A₁ ^V；A₂ ¹,A₂ ²,…,A₂ ^V；…A_m ^v…；A_M ¹,A_M ²,…,A_M ^V]，A_m ^v为火焰微元体v的自身辐射强度对m方向的火焰光谱出射辐射强度的系数，M为探测线的总数，V为火焰微元体的总数；

所述步骤2)包括：

21)、设置算法所需的中间参数的初值，吸收系数向量κ＝κ₀，吸收系数的初值向量κ₀选择随机向量或单位向量，按照下列各式依次计算对称矩阵Ω，测量值和正向计算值得差值向量r，目标函数的梯度向量g，循环终止指数S_tp和衰减系数μ，

Ω＝J^TJ (3)

r＝I′_λ(G)-I_λ(G) (4)

g＝J^Tr (5)

S_tp＝(||g||≤ε) (6)

μ＝η·max_i＝1,...,m(Ω_ii) (7)

式中，J为函数I′_λ(G)(κ)关于变量κ的一阶偏微分雅各比矩阵；I′_λ(G)为火焰单色(G)辐射强度向量的计算值；I_λ(G)为焰单色(G)辐射强度向量的测量值；Ω_ii为矩阵Ω的对角元素，ε是判定迭代过程是否终止的阈值，η为比例系数；max_i＝1,…,m(Ω_ii)为矩阵Ω对角元素的最大值；

22)计算I′_λ(G)(κ)估计值：

利用非负最小二乘算法求解式(2)，得到火焰各微元体的R波长的光谱黑体辐射强度向量I_bλ(R)；利用光谱黑体辐射强度向量I_bλ(R)，根据式(8)计算火焰各微元体的温度值T，

T＝c₂/λ(R)ln{c₁/[λ(R)⁵πI_bλ(R)]+1} (8)

式中，c₁为第一辐射常数，c₂为第二辐射常数；利用此温度值T，根据式(9)计算火焰各微元体的G波长的光谱黑体辐射强度向量I_bλ(G)：

式中，再利用当前的系数矩阵A和光谱黑体辐射强度向量I_bλ(G)，根据式(2)计算火焰单色(G)辐射强度向量I′_λ(G)(κ)；

23)、计算优化步长Δ(κ)：

计算当前吸收系数向量κ下的一阶偏微分雅各比矩阵J，根据式(3)-(5)计算当前的对称矩阵Ω，测量值和正向计算值的差值向量r及目标函数的梯度向量g：

根据当前的对称矩阵Ω，测量值和正向计算值的差值向量r及目标函数的梯度向量g，求解式(10)得到算法的优化步长Δ(κ)：

(Ω+μI)Δ(κ)＝g (10)

式中，I为单位矩阵；

24)、判断迭代收敛程度：

根据式(11)判断算法的优化步长Δ(κ)是否小于阈值ε：

||Δ(κ)||≤εκ|| (11)

如果优化步长Δ(κ)小于阈值ε，循环终止指数S_tp＝true(1)并且转到步骤28)；如果不小于，继续进行下一步骤；

25)、对更新的参数进行边界约束：

将算法的优化步长Δ(κ)加到吸收系数向量κ上，得到新的吸收系数向量κ_new；吸收系数的初值向量κ₀选择随机向量或单位向量；并且根据式(12)将新的吸收系数向量κ_new投影到约束区间Q：

P_Q(κ)＝max{min{κ,u},l} (12)

Q＝{κ∈R^m|l_i≤κ_i≤u_i,i＝1,...,m} (13)

式中，P_Q(κ)为投影函数；R^m为元素个数为m的实数集向量；l([l₁,l₂,…,l_m])为吸收系数约束区间Q下边界，取为([0,0,…,0])；u([u₁,u₂,…,u_m])为吸收系数约束区间Q上边界，取为([1×10⁸,1×10⁸,…,1×10⁸])；

26)、判断迭代指数ρ：

根据式(14)计算迭代指数ρ，

并判断迭代指数ρ是否大于0，如果大于，按照式(15)和(16)更新衰减系数μ和中间参数υ：

μ(new)＝μυ (15)

υ(new)＝2υ (16)

并且转到步骤28)，否则，继续进行下一步骤；

27)、更新迭代参数：

利用新的吸收系数向量κ_new更新原吸收系数向量κ，利用更新的吸收系数向量κ，根据式(3)-(5)计算当前的Ω,r和g，根据式(17)-(19)计算υ，S_tp和衰减系数μ，

υ＝2 (17)

S_tp＝(||g||_∞≤ε)or(||r||≤ε) (18)

μ(new)＝μ·max[1/3,1-(2ρ-1)³] (19)

28)、判断循环是否终止：

判断循环终止指数S_tp是否等于0并且迭代次数k是否小于最大迭代次数k_max，如果是，循环到步骤22)；否则，终止循环，并且进行到下一步骤；

29)、计算温度求解结果：

根据式(8)计算火焰各微元体的温度值T，输出该温度值及迭代终止的吸收系数值，实现火焰三维温度及吸收系数分布的同时重建。

2.根据权利要求1所述的火焰三维温度分布重建方法，其特征在于：第一辐射常数c₁为3.7418×10^-16W·m²。

3.根据权利要求2所述的火焰三维温度分布重建方法，其特征在于：第二辐射常数c₂为1.4388×10^-2m·K。

Images