一种火焰三维温度与烟黑体积分数分布同时重建算法

摘要

本发明公开了一种火焰三维温度与烟黑体积分数分布同时重建算法，包括如下步骤：输入火焰辐射光线在不同方向的出射光谱辐射强度信息，离散火焰微元体，建立火焰辐射传输方程；设定目标函数，利用构建的同时重建算法进行迭代计算，重构得到火焰三维温度与烟黑体积分数分布。本发明的算法在迭代过程中，将非负最小二乘算法嵌套在模拟退火算法中，使得SA算法在搜索目标函数过程中搜索效率显著提高。同时由于SA算法具有质量高、初值鲁棒性强、通用易实现等优点，使得NNLS‑SA算法具有良好的全局搜索特性且不易陷入局部最优值，在同时求解火焰三维温度和烟黑体积分数分布过程中具有较高的精度和较好的收敛性。

权利要求

1.一种火焰三维温度与烟黑体积分数分布同时重建算法，其特征在于，包括如下步骤：步骤一、输入火焰辐射光线在不同方向的出射光谱辐射强度信息，离散火焰微元体，建立火焰辐射传输方程；

步骤二、设定目标函数，利用构建的同时重建算法进行迭代计算，重构得到火焰三维温度与烟黑体积分数分布；所述步骤一包括：(11)、输入火焰辐射光线在不同方向的出射光谱辐射强度信息，包括出射光谱辐射强度Iλ(R)(x,y,z,θ,ψ)和出射光谱辐射强度Iλ(G)(x,y,z,θ,ψ)；

其中，Iλ(R)(x,y,z,θ,ψ)为起始点坐标(x,y,z)，方向坐标天顶角为θ，圆周角为ψ的火焰辐射光线在R波长下的出射光谱辐射强度；Iλ(G)(x,y,z,θ,ψ)为起始点坐标(x,y,z)，方向坐标天顶角为θ，圆周角为ψ的火焰辐射光线在G波长下的出射光谱辐射强度；

(12)、离散火焰微元体，对火焰体进行三维网格划分；

(13)、根据步骤(11)输入的辐射强度信息，记录每根光线穿行不同火焰微元体的长度及编号，建立如下方程：

式中，nK代表序号为K的光线穿过的火焰微元体的总个数；K是光线的序号；IλK(s)是序号为K的辐射光线在s方向上λ光谱辐射强度； 是序号为nK的微元体的λ光谱的吸收系数；

ip,jp代表光线沿传输方向穿过的火焰微元体； 是序号为K的火焰辐射光线经过第vn个火焰微元体的λ光谱的黑体辐射强度； 是光线穿过序号为nK的火焰微元体的几何长度；λ是火焰辐射光线的波长；

联立所有方向光线的辐射传输方程，得到如下方程式：Iλ＝AIbλ                          (2)其中，Iλ为火焰各方向的光谱出射辐射强度组成的向量[Iλ1,Iλ2……Iλm……IλK]；A为计算系数组成的系数矩阵，

Ibλ为火焰微元体的黑体辐射强度组成的向量V是火焰控制体的总个数，各微元体编号为1，2，v……V，对于序号为K的火焰辐射光线，穿过的火焰控制体总数为n，沿光线方向将穿过的各微元体编号为v1，v2，vip……vn；计算系数A的计算式如下：所述步骤二包括：

(21)、设置目标函数F，初始化吸收系数向量κ＝κ0，内循环最大迭代次数loopmax，迭代过程是否终止的阈值ε，初始化内循环次数loop＝1，外循环次数Loop＝1，退火计划的初始退火温度Ts,0；

F＝||I′λ(G)‑Iλ(G)||                         (4)其中，Iλ(G)、I′λ(G)分别为测量和计算得到的G波长的光谱辐射强度，吸收系数的初值向量κ0为随机连续分布的向量；

(22)、根据式(3)初始化系数矩阵A＝A0，计算初始目标函数Fg，Fg＝||A0Ibλ(G)‑Iλ(G)||                       (5)式中，Ibλ(G)为当前迭代次数下的G波长的黑体光谱辐射强度；

随机选取吸收系数向量κ中的第p个分量κp，利用式(6)、(7)进行扰动，得到新的吸收系数向量κnew，并根据式(3)更新系数矩阵A得到Anew：κnew,p＝κp+delta                      (6)其中，delta服从式(7)的标准正态分布：delta～N(0,1)        (7)由更新后的系数矩阵Anew和当前迭代次数下的黑体光谱辐射强度Ibλ(G)，计算当前迭代次数下的火焰G波长下的辐射强度向量的估计值I′λ(G)(κnew)：I′λ(G)(κnew)＝AnewIbλ(G)                  (8)计算当前状态目标函数Fl：

Fl＝||I′λ(G)(κnew)‑Iλ(G)||                    (9)其中迭代过程中的Ibλ(G)通过如下方式获得：首先利用测量得到的R波长的光谱辐射强度向量Iλ(R)、迭代过程中更新的吸收系数向量κnew、更新的系数矩阵Anew，运用非负最小二乘算法求解式(2)，获得火焰各微元体在R波长下的黑体光谱辐射强度Ibλ(R)，通过公式(10)，求解火焰各微元体的温度值Tf；

5

Tf＝c2/λ(R)ln{c1/[λ(R)πIbλ(R)]+1}                   (10)式中，c1为第一辐射常数，c2为第二辐射常数；λ(R)为光线的R波长；

然后由式(10)得到的每个火焰微元体的温度，根据式(11)计算火焰各微元体在G波长下的黑体光谱辐射强度Ibλ(G)；

式中，λ(G)为光线的G波长；

(23)、计算ΔF＝Fl－Fg，更新目标函数值Fg和吸收系数向量κ；

如果ΔF<0，退火过程接受所处状态为重要状态，吸收系数向量κ＝κnew，目标函数Fg＝Fl；

如果ΔF>0，是否接受为重要状态则根据固体处于重要状态的概率函数w来判断；用随机数据发生器产生一个[0,1]区间的随机数ξ，若w>ξ，则新状态为重要状态，更新吸收系数向量κ＝κnew，目标函数Fg＝Fl；否则不更新吸收系数向量和目标函数；

概率函数w为：

w＝exp[‑ΔF/(kbTs)]其中，kb是波尔兹曼常数；Ts是退火计划的退火温度；

(24)、判断内循环是否中止：判断内循环迭代次数loop是否小于内循环最大迭代次数loopmax，如果是，循环到步骤(22)；否则，中止内循环，并进行到下一步骤；

(25)、判断外循环是否中止：根据式(12)判定迭代是否达到收敛条件，如果否，则更新外循环次数Loopnew＝Loop+

1，按照式(13)的退火方式更新退火计划的退火温度Ts，循环到步骤(22)；否则，循环中止，并进行下一步骤；

Loopnew

Ts＝Ts,0·α                       (13)式中，α为温度衰减率；

(26)计算温度、烟黑体积分数，输出结果：根据式(12)计算迭代中止时火焰各微元体的温度Tend及吸收系数κend，输出温度Tend及吸收系数，并根据下式，计算烟黑体积分数；

式中，E(m)是随波长变化的烟黑复折射率m＝a‑bi的函数，其中a，b分别为烟黑复折射率m的实部和虚部。

2.根据权利要求1所述的火焰三维温度与烟黑体积分数分布同时重建算法，其特征在于，温度衰减率α取值为：0.7≤α≤1.0。

3.根据权利要求1所述的火焰三维温度与烟黑体积分数分布同时重建算法，其特征在于，非负最小二乘算法的求解步骤如下：S1：设置序列Γ为空集，序列Z＝[1,2,3……M]，且Ibλ＝0；

T

S2：计算M维向量，η＝A(Iλ‑AIbλ)；

S3：如果序列Z是空集或者对所有的序列i∈Z，均满足ηi≤0，转到步骤S11，否则进行下一步；

S4：找到指数q，使得ηq＝max{ηi:i∈Z}；

S5：将指数q从序列Z移动到序列Γ中；

S6：AΓ记为S×M矩阵；如果i∈Γ，矩阵A的第i列是新矩阵AΓ的第i列；如果i∈Z，那么新矩阵AΓ的第i列是零向量，并计算最小二乘问题 的解Z，解Z只定义满足序列Γ中元素i对应的Zi，序列Z中的元素i对应的Zi为零解；

S7：如果所有满足序列Γ中元素i对应的Zi均大于0，则使Ibλ＝Z，并且转到步骤S2，否则进行下一步；

S8：在序列Γ元素中找出指数u，使其满足以下关系：式中， 是指数为u时的黑体辐射强度；是序列Γ中元素为i时的黑体辐射强度；Zu是指数u对应的序列Z中的值；

S9：设 设Ibλ,new＝Ibλ+β(Z‑Ibλ)；

S10：将所有满足的指数 的指数i，其中i∈Γ，从序列Γ移动到序列Z中，并且转到步骤S6；

S11：计算结束。

说明书

一种火焰三维温度与烟黑体积分数分布同时重建算法

技术领域

[0001] 本发明涉及一种火焰三维温度与烟黑体积分数分布同时重建算法，属于火焰温度测量技术领域。

背景技术

[0002] 燃烧过程是一个极其复杂的物理变化和化学变化相互耦合、相互作用的过程。燃烧火焰中温度、烟黑体积分数及其它相关产物的空间分布情况直接关系到能源转化利用过
程的安全性、有效性和清洁性，这些关键参数的测量工作对设备的安全运行和节能减排具
有重要的意义。

[0003] 近年来，激光，光电子和信息等技术的发展，大大推动了非接触式测温技术的深入研究和应用。非接触式测温技术具有不干扰待测流场的优点，大体上可以分为主动式和被
动式非接触式测温方法两大类。主动式测温方法是对高温燃烧火焰施加光、声等外部激励
信号，通过检测燃烧过程与所施加的外部信号的相互作用结果，进行温度测量。而被动式测
温方法不采用任何外加信号，仅通过检测从燃烧过程中产生的信息进行火焰温度测量，主
要指辐射成像法。燃烧火焰中非侵入式的烟黑颗粒测量方法根据测量过程中光学信号的来
源也可分成两类，一类是基于激光或者其他外部添加光源的测量方法，如消光法和激光诱
导白炽光法；另一类是基于火焰或烟黑颗粒自身的发射光谱的测量方法，如双色法和发射
CT法。基于激光的光学测量方法在火焰温度、烟黑体积分数测量领域己经得到很普遍的应
用，然而由于这些方法都需要外加光源给予初始信号，因此在应用到工业现场进行测量时
会带来很多不便，复杂的现场条件会对布置光源和探测器的位置及校对光路带来困难。通
过测量火焰自身发射光谱信号来直接计算目标火焰中的温度和烟黑颗粒体积分数分布可
以有效地避免这些问题。

[0004] 利用辐射图像法测量火焰三维温度与烟黑体积分数分布，从本质上是根据边界辐射强度分布求解辐射传输体系中的温度和烟黑体积分数，是典型的反问题。目前辐射传输
反问题的求解算法主要包括以下四类：第一类是传统的正则化方法，包括Tikhonov正则化
算法、截断奇异值分解算法以及最小二乘QR分解算法等；第二类是基于梯度的优化算法，包
括共轭梯度方法、Levenberg‑Marquardt(L‑M)算法等；第三类是基于随机搜寻的智能优化
算法，包括遗传算法、模拟退火、微粒群算法、蚁群算法等；第四类是混合优化算法，即将两
种或多种不同类型的算法合成一种混合算法，以提高算法的整体性能。

[0005] 正则化算法是辐射传输反问题求解中应用较为广泛的一类算法，但是该方法仅能处理简单的关于线性方程组求解的反问题。对于火焰温度和烟黑体积分数同时重建问题，
此问题为典型的非线性问题，因此，单纯的采用正则化方法无法实现火焰温度和烟黑体积
分数分布的同时重建。基于梯度的优化算法是一类非线性优化算法。该类算法的优势在于
收敛速度快、结果稳定性好，但是这类算法普遍存在依赖初值、不易处理局部最优值等不
足。相对于传统方法，智能优化算法发展较晚，但具有不依赖初值、不需要对被优化函数求
导等优点。

[0006] 因此，借鉴不同算法的优点，将不同的算法有机的结合起来构成混合优化算法，发挥各自的优势，为温度和烟黑体积分数分布同时求解提供了新的思路。尽管智能优化算法
为不同算法的有机结合提供了更多可能，但是由于遗传算法存在过早收敛、蚁群算法搜索
过程容易出现停滞、以及粒子群算法精度低、易发散等缺点，使得目前混合优化算法在辐射
传输反问题研究中的应用还比较少，特别是火焰三维温度和烟黑体积分数分布同时重建问
题上还鲜有研究，因此新的混合优化算法的提出以及现有混合算法的改进还有待于进一步
的研究。

发明内容

[0007] 本发明所要解决的技术问题是针对上述现有求解算法的不足，提供一种火焰三维温度与烟黑体积分数分布同时重建算法。

[0008] 为解决上述技术问题，本发明采用的技术方案是：

[0009] 一种火焰三维温度与烟黑体积分数分布同时重建算法，其特征在于，包括如下步骤：

[0010] 步骤一、输入火焰辐射光线在不同方向的出射光谱辐射强度信息，离散火焰微元体，建立火焰辐射传输方程；

[0011] 所述步骤一包括：

[0012] (11)、输入火焰辐射光线在不同方向的出射光谱辐射强度信息，包括Iλ(R)(x,y,z,θ,ψ)和Iλ(G)(x,y,z,θ,ψ)；

[0013] 其中，Iλ(R)(x,y,z,θ,ψ)为起始点坐标(x,y,z)，方向坐标天顶角为θ，圆周角为ψ的火焰辐射光线在R波长下的出射光谱辐射强度；Iλ(G)(x,y,z,θ,ψ)为起始点坐标(x,y,z)，方
向坐标天顶角为θ，圆周角为ψ的火焰辐射光线在G波长下的出射光谱辐射强度；

[0014] (12)、离散火焰微元体，对火焰体进行三维网格划分；

[0015] (13)、根据步骤(11)输入的辐射强度信息，记录每根光线穿行不同火焰微元体的长度及编号，建立如下方程：

[0016]

[0017] 式中，nK代表序号为K的光线穿过的火焰微元体的总个数；K是光线的序号；IλK(s)是序号为K的辐射光线在s方向上λ光谱辐射强度； 是序号为nK的微元体的λ光谱的吸收
系数；ip,jp代表光线沿传输方向穿过的火焰微元体； 是序号为K的火焰辐射光线经过第
vn个火焰微元体的λ光谱的黑体辐射强度； 是光线穿过序号为nK的火焰微元体的几何长
度；λ是火焰辐射光线的波长；

[0018] 联立所有方向光线的辐射传输方程，得到如下方程式：

[0019] Iλ＝AIbλ                          (2)

[0020] 其中，Iλ为火焰各方向的光谱出射辐射强度组成的向量[Iλ1,Iλ2……Iλm……IλK]；A为计算系数组成的系数矩阵，
Ibλ为火焰微元体的黑体辐射强度组成的向量
V是火焰控制体的总个数，各微元体编号为1，2，v……V，对于序号为K的火
焰辐射光线，穿过的火焰控制体总数为n，沿光线方向将穿过的各微元体编号为v 1，v 2，
vip……vn。

[0021] A的详细计算式如公式(3)，

[0022]

[0023] 步骤二、设定目标函数，利用构建的同时重建算法(公式4‑14)进行迭代计算，重构得到火焰三维温度与烟黑体积分数分布；

[0024] (21)、设置目标函数F，初始化吸收系数向量κ＝κ0，内循环最大迭代次数loopmax，迭代过程是否终止的阈值ε，初始化内循环次数loop＝1，外循环次数Loop＝1，退火计划的
初始退火温度Ts,0；

[0025] F＝||I′λ(G)‑Iλ(G)||                         (4)

[0026] 其中，Iλ(G)、I′λ(G)分别为测量和计算得到的G波长的光谱辐射强度，吸收系数的初值向量κ0为随机连续分布的向量；

[0027] (22)、根据式(3)初始化系数矩阵A＝A0，计算初始目标函数Fg，

[0028] Fg＝||A0Ibλ(G)‑Iλ(G)||                      (5)

[0029] 随机选取吸收系数向量κ中的第p个分量κp，利用式(6)、(7)进行扰动，得到新的吸收系数向量κnew，并根据式(3)更新系数矩阵A得到Anew:

[0030] κnew,p＝κp+delta                      (6)

[0031] 其中，delta服从式(7)的标准正态分布：

[0032] delta～N(0,1)               (7)

[0033] 由更新后的系数矩阵Anew和当前迭代次数下的黑体光谱辐射强度Ibλ(G)，计算当前迭代次数下的火焰G波长下的辐射强度向量的估计值I′λ(G)(κnew)：

[0034] I′λ(G)(κnew)＝AnewIbλ(G)                  (8)

[0035] 计算当前状态目标函数Fl：

[0036] Fl＝||I′λ(G)(κnew)‑Iλ(G)||                   (9)

[0037] 其中迭代过程中的Ibλ(G)通过如下方式获得：

[0038] 首先利用测量得到的R波长的光谱辐射强度向量Iλ(R)、迭代过程中更新的吸收系数向量κnew、系数矩阵Anew，运用非负最小二乘算法求解式(2)，获得火焰各微元体在R波长下
的黑体光谱辐射强度Ibλ(R)，通过公式(10)，求解火焰各微元体的温度值Tf；

[0039] Tf＝c2/λ(R)ln{c1/[λ(R)5πIbλ(R)]+1}              (10)

[0040] 式中，c1为第一辐射常数，c2为第二辐射常数；Tf是每个火焰微元体的温度；λ(R)为光线的R波长；

[0041] 然后由式(10)得到的每个火焰微元体的温度，根据式(11)计算火焰各微元体在G波长下的黑体光谱辐射强度Ibλ(G)；

[0042]

[0043] 式中，λ(G)为光线的G波长；

[0044] 非负最小二乘算法的求解步骤如下：

[0045] S1：设置序列Γ为空集，序列Z＝[1,2,3……M]，且Ibλ＝0；

[0046] S2：计算M维向量，η＝AT(Iλ‑AIbλ)；

[0047] S3：如果序列Z是空集或者对所有的序列i∈Z，均满足ηi≤0，转到步骤S11，否则进行下一步；

[0048] S4：找到指数q，使得ηq＝max{ηi:i∈Z}；

[0049] S5：将指数q从序列Z移动到序列Γ中；

[0050] S6：AΓ记为S×M矩阵；如果i∈Γ，矩阵A的第i列是新矩阵AΓ的第i列；如果i∈Z，那么新矩阵AΓ的第i列是零向量，并计算最小二乘问题 的解Z，解Z只定义满足序列Γ
中元素i对应的Zi，序列Z中的元素i对应的Zi为零解；

[0051] S7：如果所有满足序列Γ中元素i对应的Zi均大于0，则使Ibλ＝Z，并且转到步骤S2，否则进行下一步；

[0052] S8：在序列Γ元素中找出指数u，使其满足以下关系：

[0053]

[0054] 式中， 是指数为u时的黑体辐射强度； 是序列Γ中元素为i时的黑体辐射强度；Zu是指数u对应的序列Z中的值；

[0055] S9：设 设Ibλ,new＝Ibλ+β(Z‑Ibλ)；

[0056] S10：将所有满足的指数 的指数i，其中i∈Γ，从序列Γ移动到序列Z中，并且转到步骤S6；

[0057] S11：计算结束。

[0058] (23)、计算ΔF＝Fl－Fg，更新目标函数值Fg和吸收系数向量κ；

[0059] 如果ΔF<0，退火过程接受所处状态为重要状态，吸收系数向量κ＝κnew，目标函数Fg＝Fl；

[0060] 如果ΔF>0，是否接受为重要状态则根据固体处于该状态的概率函数w来判断；概率函数的计算见式(12)。用随机数据发生器产生一个[0,1]区间的随机数ξ，若w>ξ，则新状
态为重要状态，更新吸收系数向量κ＝κnew，目标函数Fg＝Fl；否则不更新吸收系数向量和目
标函数；

[0061] w＝exp[‑ΔF/(kbTs)]                     (12)

[0062] 其中，kb是波尔兹曼常数，Ts为退火计划的退火温度；

[0063] (24)、判断内循环是否中止：

[0064] 判断内循环迭代次数loop是否小于内循环最大迭代次数loopmax，如果是，循环到步骤(22)；否则，中止内循环，并进行到下一步骤；

[0065] (25)、判断外循环是否中止：

[0066] 根据式(13)判定迭代是否达到收敛条件，如果否，则更新外循环次数Loopnew＝Loop+1，按照式(14)的退火方式更新退火计划的退火温度Ts，循环到步骤(22)；否则，循环
中止，并进行下一步骤；

[0067]

[0068] Ts＝Ts,0·αLoopnew                     (14)

[0069] 式中，α为温度衰减率，通常选择0.7≤α≤1.0。

[0070] (27)计算温度、烟黑体积分数，输出结果：

[0071] 根据式(13)计算迭代中止时火焰各微元体的温度Tend及吸收系数κend，输出该温度及吸收系数，并根据下式，计算烟黑体积分数；

[0072]

[0073]

[0074]

[0075] 式中，E(m)是随波长变化的烟黑复折射率m＝a‑bi的函数，其中a，b分别为烟黑复折射率m的实部和虚部。

[0076] 有益效果：

[0077] 1.火焰三维温度与烟黑体积分数分布同时重建问题，是典型的非线性求解问题，在求解过程中既需要对温度和烟黑体积分数进行解耦，同时还需要保证求解结果的准确
性。本发明利用模拟退火算法(Simulated Annealing，SA)，对目标函数值进行全局搜索，同
时在搜索过程中将非负约束最小二乘算法(Non‑Negative Least Squares，NNLS)嵌套在SA
算法中对火焰温度场进行求解，以得到高精度的温度和烟黑体积分数重建结果。

[0078] 2.SA算法在全局搜索过程中，由于每次只对某个或某几个微元体的吸收系数值进行更新，所以全局搜索时间较长，但是嵌套在SA中的NNLS算法能够保证火焰各微元体温度
重建的非负性，使得SA算法在搜索目标函数最小值过程中搜索效率显著提高。同时由于SA
算法具有质量高、初值鲁棒性强、通用易实现等优点，使得NNLS‑SA算法具有良好的全局搜
索特性且不易陷入局部最优值，在同时求解火焰三维温度和烟黑体积分数分布过程中具有
较高的精度和较好的收敛性。

附图说明

[0079] 图1是火焰三维温度与烟黑体积分数分布同时重建算法流程图；

[0080] 图2是分别利用NNLS‑SA与NNLS‑PSO算法重建火焰温度的对比结果；

[0081] 图3是分别利用NNLS‑SA与NNLS‑PSO算法重建火焰烟黑体积分数的对比结果。

具体实施方式

[0082] 下面结合附图和具体施例，进一步阐述本发明。应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种
等价形式的修改均落于本申请所附权利要求所限定的范围。

[0083] 步骤一、设置火焰参数，输入火焰辐射光线在不同方向的出射光谱辐射强度信息，离散火焰微元体，建立火焰辐射传输方程。假设圆柱形火焰模型，高度H＝0.03m，底面半径R
＝0.008m，火焰内部温度和烟黑体积分数分布服从式(1)和(2)。将火焰划分为
(V＝64)个微元体。选取10000根不同方向
光线的出射光谱辐射强度信号Iλ(R)(x,y,z,θ,ψ)和Iλ(G)(x,y,z,θ,ψ)；其中，(x,y,z)为光线
的起始点坐标，x＝‑0.0160，y的变化范围是‑0.0074～0.0074，z的变化范围是0～0.03。(θ,
ψ)为光线的方向坐标，θ为天顶角，ψ为圆周角，θ的变化范围是‑5°～1.2°，ψ的变化范围是
80°～92°，λ(R)和λ(G)分别为R和G通道的波长610nm，530nm；

[0084]

[0085]

[0086] 根据输入的辐射强度信息，记录每根光线穿行不同火焰微元体的长度及编号，建立如下方程：

[0087]

[0088] 式中，nK代表序号为K的光线穿过的火焰微元体的总个数；K是光线的序号；IλK(s)是序号为K的辐射光线在s方向上λ光谱辐射强度； 是序号为nK的微元体的λ光谱的吸收
系数；ip,jp代表光线沿传输方向穿过的火焰微元体； 是序号为K的火焰辐射光线经过第
vn个火焰微元体的λ光谱的黑体辐射强度； 是光线穿过序号为nK的火焰微元体的几何长
度；λ是火焰辐射光线的波长；其中吸收系数κλ和烟黑体积分数fv的关系见式(4)‑(6)；

[0089]

[0090]

[0091]

[0092] 其中，E(m)是随波长变化的烟黑复折射率m＝a‑bi的函数，其中a，b分别为烟黑复折射率m的实部和虚部。

[0093] 联立所有方向光线的辐射传输方程，得到如下方程式：

[0094] Iλ＝A·Ibλ                             (7)

[0095] 其中，Iλ为火焰各方向的光谱出射辐射强度组成的向量[Iλ1,Iλ2……Iλm……IλK]；A为计算系数组成的系数矩阵，
Ibλ为火焰微元体的黑体辐射强度组成的向量
V是火焰控制体的总个数，各微元体编号为1，2，v……V，对于序号为K的火
焰辐射光线，穿过的火焰控制体总数为n，沿光线方向将穿过的各微元体编号为v 1，v 2，
vip……vn。

[0096] A的详细计算式如公式(8)，

[0097]

[0098] 步骤二：设定目标函数，利用构建的同时重建算法，迭代计算，同时重构火焰三维温度与烟黑体积分数分布。

[0099] (21)设置目标函数F，如式(9)，初始化吸收系数向量κ＝κ0，内循环次数loop＝1，‑5
内循环最大迭代次数loopmax＝50，外循环次数Loop＝1，迭代过程是否终止的阈值ε＝10 ，
退火计划初始温度Ts,0＝2000K；

[0100] F＝||I′λ(G)‑Iλ(G)||                        (9)

[0101] 其中，Iλ(G)、I′λ(G)分别为测量和正向计算得到的光谱辐射强度，吸收系数的初值向量κ0为随机连续分布的向量。

[0102] (22)、根据式(8)初始化系数矩阵A＝A0，计算初始目标函数Fg，

[0103] Fg＝||A0Ibλ(G)‑Iλ(G)||                      (10)

[0104] 随机选取吸收系数向量κ中的第p个分量κp，利用式(11)、(12)进行扰动，得到新的吸收系数向量κnew，并根据式(8)更新系数矩阵A得到Anew:

[0105] κnew,p＝κp+delta                      (11)

[0106] 其中，delta服从式12的标准正态分布：

[0107] delta～N(0,1)               (12)

[0108] 由更新后的系数矩阵Anew和当前迭代次数下的黑体光谱辐射强度Ibλ(G)，计算当前迭代次数下的火焰G波长下的辐射强度向量的估计值I′λ(G)(κnew)：

[0109] I′λ(G)(κnew)＝AnewIbλ(G)                  (13)

[0110] 计算当前状态目标函数Fl：

[0111] Fl＝||I′λ(G)(κnew)‑Iλ(G)||                   (14)

[0112] 利用Matlab的lsqnonneg函数实现NNLS算法来求解式(7)，获得火焰各微元体在R波长下的黑体光谱辐射强度Ibλ(R)，通过公式(15)，求解火焰各微元体的温度值Tf，

[0113] Tf＝c2/λ(R)ln{c1/[λ(R)5πIbλ(R)]+1}               (15)

[0114] 式中，c1为第一辐射常数，c2为第二辐射常数，分别为3.7418×10‑16W·m2和1.4388‑2
×10 m·K；Tf是每个火焰微元体的温度，K；λ(R)为光线的R波长；

[0115] 由式(15)得到的每个火焰微元体的温度，根据Planck定律(式(16))计算火焰各微元体在G波长下的黑体光谱辐射强度Ibλ(G)：

[0116]

[0117] 式中，λ(G)为光线的G波长；根据式(7)，由当前的系数矩阵A和黑体光谱辐射强度Ibλ(G)，计算当前火焰G波长下的辐射强度向量的估计值I′λ(G)(κnew)，根据式(14)计算当前状
态目标函数Fl。

[0118] (24)、计算ΔF＝Fl－Fg，更新目标函数值Fg和吸收系数向量κ；

[0119] 如果ΔF<0，退火过程接受所处状态为重要状态，吸收系数向量κ＝κnew，目标函数Fg＝Fl；

[0120] 如果ΔF>0，是否接受为重要状态则根据固体处于该状态的概率函数w来判断，概率函数的计算见式(17)。用随机数据发生器产生一个[0,1]区间的随机数ξ，若w>ξ，则新状
态为重要状态，更新吸收系数向量κ＝κnew，目标函数Fg＝Fl；否则不更新吸收系数向量和目
标函数；

[0121] w＝exp[‑ΔF/(kbTs)]                     (17)

[0122] 其中，kb是波尔兹曼常数，kb＝1.3806488×10‑23，Ts为退火计划的退火温度；

[0123] (25)、判断内循环是否中止：

[0124] 判断内循环迭代次数k是否小于内循环最大迭代次数kmax，如果是，循环到步骤(22)；否则，中止内循环，并且进行到下一步骤；

[0125] (26)、判断外循环是否中止：

[0126] 根据式(18)判定迭代是否达到收敛条件，如果否，则更新外循环次数Loopnew＝Loop+1，按照式(19)的退火方式更新退火计划的退火温度Ts，循环到步骤(22)；否则，循环
中止，并进行下一步骤；

[0127]

[0128] Ts＝Ts,0·αLoopnew                     (19)

[0129] 式中，α为温度衰减率，通常选择0.7≤α≤1.0。

[0130] (27)计算温度、烟黑体积分数，输出结果：

[0131] 根据式(18)计算迭代中止时火焰各微元体的温度Tend及吸收系数κend，输出该温度及吸收系数，并根据下式，计算烟黑体积分数；

[0132]

[0133] 式中，fv为火焰微元体的烟黑体积分数，ppm；κend是火焰微元体的吸收系数，m‑1；E(m)为烟黑复折射率m的函数。

[0134] 为了评价NNLS‑SA算法的性能，本发明引入NNLS‑PSO算法进行比较,利用NNLS‑SA与NNLS‑PSO算法分别重建火焰的温度及烟黑体积分数，其中微粒群优化算法(Particle 
Swarm Optimization,PSO)是根据鸟类的群体觅食行为提出的，也是一种基于随机搜寻的
智能优化算法。

[0135] 按照式(21)计算火焰烟黑体积分数的重建相对误差，式(22)计算火焰温度的重建相对误差，重建结果如图2、3所示。

[0136]

[0137]

[0138] 式中，fv,rst为重建的烟黑体积分数，fv,set为设定的烟黑体积分数，Trst为重建的火焰温度值，Tset为设定的火焰温度值。

[0139] NNLS‑PSO算法主要利用蒙特卡洛方法，每一次迭代，对吸收系数向量κ的全体值进行更新，通过粒子群对解空间进行搜索，并通过比较目标函数的历史最优值得到烟黑体积
分数与温度的最优值。NNLS‑SA算法在每一次迭代过程中，仅对某个或者某几个微元体的吸
收系数值进行更新，进而比较更新前后目标函数的变化，以求得烟黑体积分数与温度的最
优值。

[0140] NNLS‑SA算法将NNLS算法嵌套在SA算法中，保证了温度场求解的温度性，同时利用SA算法对于非线性问题良好的全局搜索特性，对烟黑体积分数进行求解。与NNLS‑PSO算法
相比，算法不盲目的对所有微元体的吸收系数进行随机更新，而是选取某个或者某几个微
元体的吸收系数进行更新，从而算法的精度更高，不易陷入局部极值，稳定性较好，同样的
迭代次数下NNLS‑SA算法的求解精度更高。从图2火焰温度的重建结果中看到，在第46‑48号
网格，温度的真实值是901K，利用NNLS‑SA算法重建的温度分别为894K，894K和901K，利用
NNLS‑PSO算法重建的温度分别为863K，862K，和863K，可以看到利用NNLS‑SA算法重建的火
焰温度更为精确。从图3烟黑体积分数的重建结果也可以看到，在第10‑12号网格，烟黑体积
分数的真实值为0.8031ppm，利用NNLS‑SA算法重建的火焰烟黑体积分数分别为0.8197ppm，
0.8042ppm，和0.8082ppm，利用NNLS‑PSO算法重建的火焰烟黑体积分数分别为0.7283ppm，
0.9224ppm，和1.3466ppm，相比于NNLS‑PSO算法，NNLS‑SA算法对于烟黑体积分数的重建结
果也更为精确。同时，计算了两种算法重建温度和烟黑体积分数的平均相对误差，利用
NNLS‑SA算法重建温度、烟黑体积分数的平均相对误差分别是0.7％和6.05％，利用NNLS‑
PSO算法重建的温度、烟黑体积分数的平均相对误差分别为10.87％和43.85％。综上，NNLS‑
SA算法无论是对于温度还是烟黑体积分数的重建，NNLS‑SA算法均具有更高的重建精度。