组合优化算法在冲击温度反演中的应用

张宁超, 叶 鑫, 李 多, 解孟其, 王 鹏, 刘福生, 钞红晓

1. 西安工业大学电子信息工程学院, 陕西 西安 710021 2. 西南交通大学高温高压物理研究所, 四川 成都 610031 3. 西北机电工程研究所, 陕西 咸阳 712099

引 言

冲击波物理研究中, 温度是描述冲击压缩下材料状态和性质不可或缺的物理量, 对高压科学、 地球科学、 行星科学和材料科学具有特殊重要性[1-3]。 同时, 特殊材料辐射光谱温度反演研究, 对航空航天、 深空探测和工业制造等领域有重要应用价值[4-5]。 在温度求解过程中, 多光谱的n个通道, 构成了n+1个未知量的欠定方程组, 目标真实温度的获取, 旨在求解该欠定方程组。 当前对光谱温度求解方法主要包括发射率模型假设、 神经网络模型训练等数据处理方法。 前者存在一定局限性, 只有当假定的材料发射率模型与实际发射率相符时才能准确的反演目标温度, 面对未知发射率的材料, 温度反演将束手无策。 早期Flower等[6]利用模型假设法, 开展了较多的实验与应用研究, 虽然得到了一些材料的反演温度, 但模型的不确定性限制了其推广。 杨永军[7]针对高温材料表面温度难以准确测量的问题, 利用光谱仪测量材料的能量信息, 通过多光谱测温获取其表面温度进而计算出光谱发射率, 并利用不锈钢在1 100 K时的发射率进行了验证。 孙红胜等[8]针对弥散介质下的材料温度测量, 分析并总结了几种主要的辐射测温方法的优缺点, 为弥散介质下的高温测量提供了研究方向。 孙晓刚等[9]较早将遗传算法与神经网络相结合, 将其应用于目标温度的测量与反演, 该方法提高了反演的精度, 但需要大量样本数据作为支撑, 求解温度前需对样本进行训练, 所需要的数据集大且训练所需要的时间较长。 然而在冲击温度测量实验中, 材料在高温高压状态下的光谱发射率模型很难获得, 尤其在面对发射率未知的材料或被测材料发生结构相变时, 这种方法将不再适用。 根据约束优化理论, 如果将材料温度求解转化为目标优化问题, 通过建立相应的温度反演模型, 利用优化算法对模型进行求解, 即可实现温度的反演。 根据测量特点及应用范围, 戴景民等[10]对发射率测量方法进行了分类, 分析了发射率测量方面的问题, 将人工神经网络与二次细分方法相结合提高了测温精度, 但该方法所需发射率样本集过多, 存在一定局限性。 刑键等[11]利用约束优化对温度反演进行了系统研究, 利用二次测量法结合迭代递推算法, 提出了无需假设发射率模型的辐射温度测量方法, 利用发射率约束偏差的思想提高了温度反演的速度, 基于广义逆-坐标变换提高了温度反演的精度。 根据应用场景不同, 张福才等[12]使用了不同类型的优化算法, 分别对单目标和多目标模型进行了算法仿真, 并对火箭发动机运行时尾焰温度值进行了测量与反演, 温度反演过程在效率和精度上都取得进展。 冲击温度属于材料内部受到高强度压缩能量在瞬间的释放, 因此冲击温度测量过程更快, 对温度反演的时间效率要求更高, 同时冲击压缩下材料的结构极容易发生相变, 更不能保证发射模型的稳定性与唯一性。 为此本工作基于多光谱测量冲击温度, 将测温方法与优化组合理论相结合, 避免目标发射率模型对温度反演的影响, 提出将平衡优化器算法(equilibrium optimizer, EO)与序列二次规划法(sequential quadratic programming, SQP)算法组合使用, 既提高了温度反演计算的效率, 也确保了温度反演的精度和结果的稳定性。

1 实验部分

1.1 多光谱测温

图1为冲击温度的实验测量系统。 二级轻气炮发射高速飞片, 产生冲击波通过金属基板, 由于窗口材料一直处于透明状态, 所以光纤束将记录冲击波压缩下金属后界面发光信息。

图1 冲击辐射实验测试系统Fig.1 Shock radiation test system

图2是实验测量得到的金属铜后界面发光强度信息。 飞片撞击材料冲击波传到金属后界面时, 由于透明窗口与金属基板间存在空气间隙,t1时刻出现发光尖峰, 说明此刻光纤开始接收金属基板后界面的辐射信息,t2时刻辐射强度出现明显拐折, 这说明辐射受边侧稀疏波或者追赶稀疏波的影响。

图2 金属铜辐射强度曲线Fig.2 Copper impact radiation intensity curve

1.2 多光谱测温的温度反演模型

由普朗克黑体辐射定律可知,n个通道的辐射高温计中,当材料真实温度为T, 黑体参考温度为T′时, 材料温度与参考温度之间的关系如式(1)

(1)

由式(1)结合辐射测温理论, 理想状态下每个通道求解得到的温度都应相同, 因此各通道反演得到的目标温度理论偏差应趋近于零, 由此可构成目标函数, 如式(2)所示

(2)

式(2)中,Ti为各个通道反演获取的目标温度,E(Ti)为所有通道实际反演温度的平均值。 对式(1)进行变形可得

(3)

(4)

每个通道获取到的温度均值如式(5)

(5)

理想状态下, 各通道反演得到的温度均值与第一个通道温度间的差值也应为0, 由此可以构成等式约束优化模型如式(6)所示

T1-E(Ti)=0

(6)

由发射率0≤ε(λi,T)≤1构成不等式约束优化模型如式(7)所示

0≤ε(λi,T)≤1

(7)

式(7)中,f(x)为各个通道温度值的方差。

2 算法仿真与实例分析

温度反演模型建立了黑体辐射温度、 目标温度、 发射率与波长之间的函数关系, 该模型可利用优化算法求解得出材料的发射率与物理温度值。 传统优化算法求解比较依赖初始点且计算速度较慢, 为了适应冲击温度测量的需求, 首先选用平衡优化器算法, 该算法具有参数少、 耗时少且寻优能力强的优点。

2.1 基于平衡优化器的温度反演

平衡优化器算法(EO算法)是以动态质量平衡为灵感而提出的启发式优化算法。 质量平衡方程可用一阶微分方程表示为

(8)

式(8)中:V为控制容积;c为控制容积内的浓度;Q为控制容积的容积流速;ceq为控制容积平衡状态的浓度;G为控制容积内的质量生成速率。

通过求解式(8)可得

c=ceq+(c0-ceq)F+G(1-F)/λV

(9)

(10)

t=(1-iter/Maxiter)a2iter/Maxiter

(11)

EO中取α用来控制勘探能力,a1越大勘探能力越强, 开采能力将变弱;a2控制开采能力,a2越大开采能力越强, 勘探能力逐渐越弱。 通常a1取值范围为1～3,a2范围为0～2。 由于没有通用的模型, 为了选取最佳的EO算法参数, 选用了四类常见材料发射率模型[13-14]对算法进行仿真验证:

模型A:ε(λ)=ea+bλ,

模型B:ε(λ)=1/[a+bln(λ)/λ],

模型C:ε(λ)=a+bλ,

模型D:ε(λ)=a+bλ+cλ2

式中:ε(λ)为发射率;a、b和c为系数。

上述四种模型涵盖了绝大部分材料的发射率特性, 根据选用材料冲击温度范围, 设材料物理温度为3 000 K, 通过发射率数据对EO算法中的参数进行优化调整, 参考温度选取1 600 K, 具体数据如表1所示。

表1 四种模型的发射率参数Table 1 Emissivity parameters of four models

设定a2不变, 在a1取不同数值时, 利用EO算法对目标温度进行反演, 结合常见材料的发射率特征, 将发射率的范围限定为0.1≤ε(λ,T)≤0.9, 结果如表2所示。

表2 不同a1值的四种模型的温度反演结果(K)Table 3 Temperature inversion results of four models with different a1 values (K)

对表2中的温度数据进行求解, 获取不同a1取值下的发射率数值。 结果如图3所示。 模型A中当a1取值为1.8和2.0时发射率曲线与理论值最为接近, 当a1取值为2.8时, 由于算法的开采能力大幅减弱, 导致反演的发射率与理论值误差较大; 模型B中当a1取值为1.4、 1.8、 2.0和2.4时得到的发射率曲线与理论值最为接近, 其中当a1为1.4时的发射率低于理论值; 模型C中当a1取值为2.0和2.4时得到的发射率曲线与理论值最为接近,a1为2.0时的发射率低于理论值,a1为2.4时的发射率高于理论值; 模型D中a1取值为1.8、 2.0、 2.2和2.4时的发射率与理论值更为接近, 其中a1为1.8和2.4时的发射率高于理论值。

图3 不同a1值时四种模型的发射率趋势Fig.3 Emissivity trend of four models at different a1 values

设定a1不变, 选取a2范围为0.2～1.8, 步长为0.2, 对四种模型的温度进行求解, 结果如表3所示。

表3 不同a2值的四种模型的温度反演结果(K)Table 3 Temperature inversion results of four models with different a2 values (K)

对表3中的温度数据进行求解, 获取不同a2取值下的发射率数值, 如图4所示。

图4 不同a2值时四种模型的发射率趋势Fig.4 Emissivity trends of four models at different a2 values

结合图4和表3的数据, A模型中a2取值为0.8、 1.0与1.8时, 发射率变化曲线与理论曲线接近; B模型中a2取值为0.4或1.0时, 发射率值曲线与理论值接近; C类模型a2取值为1.0或1.4时, 温度反演效果较佳; D类模型a2取值为1.0或1.8时发射率与理论值更为贴近。 结合图1可得, 当a1取1.8,a2取1.0时, EO算法的温度反演效果最佳。 然而, 在多次求解寻找最优值时发现, EO算法得到的温度值波动较大, 稳定性较低。 这是由于EO算法在优化过程中, 每一次迭代会等概率的随机从平衡池中选择一个当前的最优解进行迭代。 由于从平衡池中每次选取的当前最优解存在随机性, 使得迭代的寻优方向变得随机, 导致每次计算出的温度值之间的浮动较大, 稳定性偏低。

2.2 基于平衡优化器结合多起点序列二次规划的温度反演

为解决平衡优化器求解温度值稳定性偏低, 且误差较大的问题, 将平衡优化器算法(equilibrium optimizer, EO)与序列二次规划法(sequential quadratic programming, SQP)相结合, 对温度反演方法进行完善, 解决EO算法反演求解稳定性不足的问题。 EO算法运算速度快, 收敛性强, 但最终求解结果的稳定性较差, 解的准确性较低, 易陷入局部最优。 SQP法边界搜索能力强, 解的稳定性高, 但初始点的不同, 会导致运算量增大而影响收敛速度, 增加算法的计算效率。 结合两种算法各自的优点, 用EO算法求得的解作为SQP算法的初始点之一, 利用SQP方法展开多起点迭代寻优, 实现多起点寻优的方式求解, 保证精度和稳定性高于单一算法, 实现精度高、 速度快且稳定性强的温度反演。

利用SQP方法来完成温度求解, 随机在可行域范围内初始化多个起始点, 将求得的解进行对比并输出满足式(7)的最优解, 即为目标精度最高的温度值。 仿真结果表明起始点个数n直接影响温度值精度与计算时间。 分别取3, 6, 9, 12和15为起始点, 对四种发射率模型进行求解, 当起始点的个数为3或15时, 求得的四种模型对应温度值误差都很高, 起始点个数取9时, 四种模型误差较低。 不同起始点个数四种模型温度反演结果如表4所示。

表4 不同起始点数n的四种模型的温度反演结果Table 4 Temperature inversion results of four models with different starting points n

为了验证本温度反演方法的有效性, 分别利用EO-SQP组合算法、 基于内点罚函数(interior penalty function, IPF)的对四种模型温度进行验证, 将两种方法计算的精度和时间进行对比, 表5中a1,a2起始点的个数取值分别为1.8, 1.0, 9。

表5 两种方法温度反演结果对比Table 5 Comparison of temperature inversion results between two methods

结果表明IPF方法在求解B类模型的温度时, 误差较小, 仅为0.39%。 但在求解A、 C和D三类模型时, 温度误差普遍较大。 尤其在求解C类模型时, 误差达到了2.37%, 且IPF在反演温度时, 程序迭代运算所需要的时间较长。 本文算法求解得到的四类模型温度值, 误差均小于1%, 且算法的运行时间小于3 s, 相较于IPF, 在提高精度的前提下, 反演效率大幅提升, 求解精度和速度都要优于IPF。

由亮度温度模型, 各通道反演出的光谱发射率可表示为

(12)

根据式(12), 对表5中两种方法反演出的温度结果进行计算, 获取其对应的发射率曲线, 并与目标发射率理论值进行对比, 如图5所示。

图5 3 000 K时四种模型发射率变化曲线Fig.5 Emissivity curves of four models at 3 000 K

图5中IPF方法和本组合反演算法对模型A、 B和D都较适用且发射率反演误差较小, 但IPF应用于C类材料模型的温度反演时, 发射率曲线与目标发射率理论值偏差较大, 结合表5中数据可知, 本文提出的组合优化算法适用范围广。

2.3 实例分析

利用典型的辐射高温计, 获取8波长目标冲击发光强度值。 测试系统各通道的中心波长分别为: 0.809、 0.779、 0.702、 0.650、 0.589、 0.533、 0.509和0.488 μm。 分别利用最小二乘法(least square method, LSM)、 内点罚函数法(IPF)以及本算法(EO-SQP)对冲击过程各时刻温度进行求解, 温度反演结果如图6所示。

图6 三种方法温度反演结果Fig.6 Temperature inversion results of three methods

依据金属测温理论, 实验直接获取的辐射对应界面温度, 并不代表在该压力下金属的卸载温度。 由“理想界面模型”, 金属卸载温度、 界面温度以及窗口冲击温度三者间的关系如式(13)。

(13)

利用式(13)计算得到本实验冲击温度、 卸载温度以及界面温度数值, 将其与本算法、 最小二乘法以及内点罚函数法反演得到的温度值进行对比。

表6是本次实验相关量的理论值与实验计算结果。 最小二乘法计算出的界面温度绝对误差约为16.44%, 内点罚函数法计算出的界面温度绝对误差约为11.37%, 本算法求解出的界面温度绝对误差约为4.41%, 得到的结果更接近理论计算值。

表6 三种方法的计算结果比较Table 6 Comparison of calculation results of three methods

3 结 论

要获得材料真实的冲击温度值, 除了精密的测试手段与方法, 温度反演模型与温度反演算法的选取至关重要。 针对以往冲击温度测量方法适用性差, 温度反演速度慢的问题, 将普朗克理论与目标优化相结合, 建立基于约束优化的温度反演模型, 并利用平衡优化器算法对模型进行求解以获取材料温度值, 提高了温度反演的运算效率以及温度反演方法的适用性; 针对平衡优化器算法求解精度较差且稳定性偏低的问题, 利用多起点序列二次规划法进一步完善求解, 避免平衡优化器算法陷入局部最优, 以提高温度反演的精度和效率。 结果表明, 本文提出的算法是一种有效的冲击波压缩下材料温度反演算法, 实现了高效、 高精度以及适用性广的温度反演。