柴油燃烧烟尘干扰的红外辐射特性研究

2021-12-28马向超蒋佳丽

空天防御 2021年4期

田益，马向超，蒋佳丽

（西安电子科技大学，陕西西安 710071）

0 引言

在现代军事战争中，柴油燃烧所形成的烟尘是非常普遍的。例如：坦克、装甲车等重型武器运行过程中，柴油经过柴油机燃烧后会释放出大量的烟尘；加油站和油井等能源设施也经常被作为军事打击的目标；在某些特殊情况下，为了隐蔽目标，有时也会故意点燃油井、油桶等形成大量燃烧烟尘，以减少目标被命中的几率，如海湾战争中伊拉克故意点燃油井躲避美军红外武器的攻击。

随着红外成像导引头技术［1］的快速发展，虽然采用第三代凝视红外焦平面探测器的导引头系统综合性能大幅提高，但是燃烧烟尘等扩展干扰物依然是制约红外制导导弹作战效果的主要措施之一。当导引头和目标之间存在燃烧烟尘时，由于其消光遮蔽作用，可破坏导引头对目标的发现、识别和跟踪，使得处于搜索状态的导弹不能转入跟踪状态；处于跟踪状态的导弹则会丢失目标，由跟踪转入搜索，或者导致其跟踪误差增大。因此，柴油燃烧所形成的烟尘不仅是己方装备对抗红外成像制导武器的潜在手段之一，也是己方装备的导弹红外导引头所面临的潜在干扰之一。

目前对柴油燃烧烟尘的研究主要在于对烟尘的成分进行物化分析，如柴油完全燃烧时的主要产物是二氧化碳和水，不完全燃烧时还产生一氧化碳甚至碳微小颗粒。

由于柴油的燃烧烟尘干扰普遍存在且容易释放，本文拟通过对这种烟尘的物理化学特性、时空运动特性及其红外辐射的空间、时间和光谱传输特性等进行研究，构建燃烧烟尘的扩散运动和红外辐射传输特性模型，为红外制导武器抗燃烧烟尘干扰能力的评估提供基础。

1 柴油燃烧烟尘的收集和微观形貌分析

柴油燃烧烟尘的组成成分较为复杂，受燃烧条件的影响很大。我们使用煤油灯对柴油进行直接可控燃烧，然后采集其燃烧烟尘。

直接采集烟尘所使用的方法如图1 所示，其目的是为了获得用于微观形貌分析的颗粒物。首先，将导电胶带粘贴到电镜样品台上，靠近燃烧烟尘排放通道，使烟尘颗粒直接流经导电胶带被沾附于样品台，从而获得燃烧烟尘样品。

图1 燃烧烟尘采集示意图Fig.1 Schematic of smoke collection from diesel combustion

图2为所收集的燃烧烟尘样品的扫描电子显微镜（scanning electron microscope，SEM）照片，显示了500 nm 标尺下的表面形貌。从图2 可以看出，燃烧烟尘主要是由小尺寸的颗粒组成，小颗粒物的形状近似为球形，尺寸在20～40 nm之间。

图2 柴油燃烧烟尘的SEM图（标尺500 nm）Fig.2 SEM image of smoke from diesel combustion（scale ruler 500 nm）

2 红外辐射传输模型

2.1 单个粒子消光特性

在建立烟尘的红外辐射传输模型［2］前，需要知道单个烟尘粒子的消光特性参量，如散射截面σsca、消光截面σext、不对称因子g等。为了简化，根据实验分析结果，我们将烟尘粒子视为球形粒子。对于球形粒子，这些参量可由Mie散射理论求解得到，也可采用离散偶极近似（discrete dipole approximation，DDA）算法计算。考虑到通用性，本文采用DDA 算法来计算单个烟尘粒子在特定波长下的消光特性。

当波长不同时，粒子的复折射率会发生改变，因此粒子对辐射的消光特性会随着波长变化。Chang等［3］对烟尘颗粒提出了一组对数拟合函数，复折射率的表达式为n′＝n＋ik，n和k的拟合函数如下所示。

式中：λ为入射光波长，单位μm。

根据对燃烧烟尘微观形貌的分析，球形粒子半径取20 nm，依据Chang 等提出的复折射率模型，在红外波段上等间隔选取m个采样点来计算对应波长下单个烟尘颗粒的散射截面、消光截面以及不对称因子，最终得出烟尘粒子在采样点处的光学参量。

2.2 粒子系消光特性

燃烧烟尘是由多个粒子组成的弥散粒子系。根据粒子复折射率及尺度参数计算出单个粒子消光截面及不对称因子后，再结合粒子系的物理参数（如粒子浓度、形状及粒径）可求出粒子系的消光及散射特性。

为了简化模型，本文将粒子系中的粒子视为相同粒径的球体，则单个粒子在某波长处消光截面σext乘以网格区域单位体积内粒子数N的结果即为粒子系在该波长下的消光系数βsys：

由于在单个颗粒的消光中，有一部分是前向散射，依然会进入视场，因此，需要考虑前向散射的贡献。Henyey 和Greenstein 根据多年研究，提出了一个与已知实验数据较为符合的散射相函数经验公式，简称为HG相函数［4］，如式（4）所示。

式中：g是一个无量纲的控制散射相函数分布的常数，为单个烟尘颗粒的不对称因子。当g＜1时，散射是各项异性的，g越大，分布越集中，前向散射占据的比率越大；当g＝1时，散射方向完全集中在前向。

单个网格内粒子系在各角度散射的能量比率为

式中：∑σsca是流场计算中在某波长下单个网格内粒子的累加散射截面；Δs是垂直入射方向的单个网格横截面积；PHG（Θ）为该波长处单个烟尘颗粒HG 的散射相函数。当Θ 取0 时，F（Θ）即为单个网格内粒子系的前向散射比率。

2.3 粒子系光谱消光特性

燃烧烟尘颗粒的光学特性具有强烈的光谱依赖性，粒子系在某个波段内的消光特性通常是对光谱平均获得的，其中常用的平均方法是普朗克平均法，即对之前计算的粒子系在该特定波段内的光谱消光特性结果进行普朗克（Planck）平均，普朗克平均消光系数为

式中：f1为在红外波段λ1～λm上分别等间隔选取m个采样点后计算采样点对应的消光系数βsys，并将βsys数值进行拟合而得的4 阶多项式函数；λ1和λm分别为波段积分下限与上限；Ebλ为特定波长范围内等间距所取的各个波长对应的黑体的光谱辐射照度。

同理，普朗克平均前向散射比率可表达为

式中：f2为利用在红外波段λ1～λm上等间隔选取m个采样点后计算采样点对应的前向散射比率F（Θ），并将F（Θ）数值进行拟合而得的4阶多项式函数。

上述计算过程可用图3所示流程图表达。

图3 粒子系在不同波段的消光特性及前向散射计算流程图Fig.3 Flow chart for calculating the extinction properties and forward scattering of particles in different wavebands

2.4 出射强度

沿厚度方向将烟尘平均分割为Ω份，每份的厚度为Δx，如图4所示。

图4 烟尘平均分割为Ω份的示意图Fig.4 Schematic of dividing the smoke into N parts

出射强度为

其中：I0为入射红外辐射强度。

3 时空扩散模型

3.1 物理模型

燃烧烟尘的扩散运动符合流体的运动规律，因此可以基于流体力学理论对烟尘的时空扩散运动进行建模。基于物理过程的烟尘扩散运动仿真的理论基础是求解其对应的Navier-Stokes（NS）方程，可使用半拉格朗日法与压力泊松方程来高效地求解方程，并使用涡量约束法［5］来降低由半拉格朗日法引起的数值耗散问题，从而准确地模拟出烟尘流场的湍流现象。

一般而言，烟尘的黏性可忽略不计。当烟尘的扩散速度远小于声速时，可压缩性的影响也可以忽略不计，因此视烟尘为无黏性、不可压缩流体［6］。求解烟尘流场的方程可表达为［6］

式中：ρ为密度；p为压强；a为体积力（如重力，惯性力）的加速度；u为速度。

流场中烟尘颗粒受到的浮力与温度有关，温度较高的烟尘颗粒受到浮力的作用上升，通常用如下经验公式描述［7］。

式中：α表示温度对浮力的影响系数；T为烟尘温度；Tair为环境温度；y表示浮力方向上的单位向量。

温度在空间分布中不是均匀的，导致浮力的空间分布也不均匀，所以该浮力对速度场的影响不会被压力项完全消除，会形成对流。

通过对温度的计算可进而求出烟尘颗粒所受的浮力大小与烟尘的自发辐射强度，为此需要引入一个温度标量场。温度场不仅自身具有扩散效应，速度场也会对其产生对流影响。温度场可表示为

式中：cT为温度冷却系数；Tmax为流场最高温度；Tair为环境温度。

风与大气条件和海拔高度有关，流场中用于描述风场的公式为

式中：u0是处于海拔为z0时测量出的风速；z是地面海拔；r用于描述大气中风速随海拔高度的变化情况。在实际仿真计算中取z0为10 m，并设r的默认值为0.3。

3.2 基于网格的烟尘流体仿真方法

不可压缩流体控制方程可化为拉格朗日形式和欧拉形式，两种形式在理论上是等价的。在CFD 中，由于计算精度有限，且两种形式的离散化方式不同，会形成两种截然不同的仿真方法：基于粒子的方法、基于网格的方法。

本文采用的流体计算方法是欧拉法，将流场计算区域划分为离散网格。一个三维空间的规则网格的数据结构就是一个三维数组，只需要记录数组的大小，通过下标就能完成对数组元素的访问。本文的流场计算采用规则的静态结构化网格来描述。依据流场维度将计算区域划分为X×Y×Z的三维立方体网格。

网格将空间均匀地划分，仿真过程中速度、温度、压强、密度这些物理量都记录在网格单元中心。为了简化流场边界处理，在实际仿真过程中需要在流场周围额外加上一层网格［8］。因此在程序中，对每个物理量分配的数组大小为（X＋2）×（Y＋2）×（Z＋2）。

3.3 红外烟尘渲染

红外烟尘，特别是军事上用于遮挡目标而使用的红外烟幕，发烟器都会具有比环境温度高的初始温度，故需要考虑烟尘的红外自发辐射。渲染红外烟尘主要需要考虑的因素是透过率和对辐射的散射。烟尘粒子对光子的相互作用概率决定了其透过率，从视线方向穿透烟尘体的总透过率就是最终被渲染出来的α通道，α通道即在OpenGL 渲染引擎内用于存储透明度值的通道。烟尘的消光系数与消光截面和密度相关，辐射一部分会被外围的烟尘吸收和散射，另一部分会透射出来并进入相机［9］。本文只考虑在视线积分路径上的前向散射。

由于烟尘运动模型采用网格方法构建，这种网格无法通过常规的光栅化方法进行渲染，故采用的是体绘制算法［10］，该算法能够高效地处理烟尘体在视线方向上的透过率和前向散射，以达到实时性的需求。

本文采用了体绘制中最常用的光线投射算法，如图5 所示，展示了从相机视口发出的射线穿过体纹理并进行采样的过程。首先需要计算烟尘占据区域，当从相机引出的射线命中烟尘前表面时开始进行体绘制，光线传输方向与射线方向相反：射线不断地在烟尘体积内步进，不断累计亮度和透过率，当射线离开长方体后表面，体绘制结束，其过程中累计的亮度和透过率就是烟尘的亮度和透过率。

将图5 射线穿越体纹理的过程作为采样合成过程，则合成公式可表达为

图5 射线穿过体纹理［11］Fig.5 Rays passing through the volume texture

式中：Ci和Ai分别是在体纹理上采样所得到的颜色值和不透明度，其实也就是体素中蕴含的数据；和表示累加的颜色值和不透明度。

基于本文所建立的红外辐射传输模型，对体纹理的颜色值Ci与不透明度Ai有

式中：Iself为网格内烟尘的自发辐射强度；Itrans为透射进来的辐射强度；对于烟尘自身红外辐射，将其简化为灰体辐射模型，且使其辐射强度与浓度成正比，该网格单元发出的自身红外辐射强度为

式中：Ibλ为T温度下λ波长的黑体辐射强度；λ1～λ2为所计算辐射的波段；ε为由烟尘粒子的等效发射率；N为网格内的粒子数。

3.4 具体实现方法

针对烟尘的扩散运动模型，基于NS 流体力学控制方程，采用欧拉网格方法求解流场的扩散运动，主要包括：求解压力项、计算对流、降低流场涡量和其他内外力项的解算。将流场计算区域分为40×40×40个网格，对每个网格内的烟尘运动状态进行计算迭代，最终可得出网格内的烟尘的速度、浓度、温度等物理量。

针对烟尘的自发红外辐射，首先建立烟尘的等效发射率模型，然后基于灰体模型建立烟尘粒子系的自身红外辐射模型。针对烟尘红外辐射的时间、空间、波段红外辐射传输特性，采用DDA 算法，且分别在3～5 μm 以及8～12 μm 红外波段等间隔选取40 个采样点，计算单个烟尘粒子在不同波长处的光谱和波段光学参量，如散射与吸收截面、散射相函数等；继而可计算粒子系在各波长采样点处的消光系数和前向散射比例；然后，基于普朗克平均方法获得烟尘粒子系在3～5 μm、8～12 μm 红外波段的消光和散射特性；最后，采用对烟尘仿真区域进行等间隔采样的光线投射算法，结合单个网格内粒子系的消光、散射和自发辐射模型，从成像平面的每个像元引出一条射线，计算途经每一条射线到达视点处的红外辐射，如图6所示。

图6 烟尘辐射传输计算示意图［12］Fig.6 Schematic for calculating the radiation transfer of smoke

依据以上一系列理论和方法，我们使用OpenGL图形API 与CUDA 并行程序搭建整个烟尘的计算和渲染的框架，实现了基于欧拉网格方法的烟尘扩散运动和红外辐射传输特性的联合解算。

4 结果与分析

烟团的运动过程主要由烟尘源（烟尘发生器）类型和边界条件决定，运动导致烟尘的温度和密度随时间变化，所以在某一时刻的温度和密度的分布决定了烟尘的亮度和透过率。下面将结合仿真结果分析扩散时间、风场边界条件以及波段对烟尘运动的影响，并检验其对恒定温度黑体的遮挡效果。

4.1 风场对仿真结果的影响

图7 显示的是在5 m×5 m×5 m 的空间内，无风条件下不同时刻烟的扩散图。发射源持续释放初始速度为1 m/s 且速度方向垂直向上的烟尘，烟尘初始浓度为0.1 kg/m3。发射源温度恒定为350 K，环境温度设为280 K，并且在流场区域放置一个温度与发射源相同的黑体，在每个时刻对烟尘的扩散情况进行计算仿真。

图7 无风情况下不同时刻烟的扩散图Fig.7 Spread of smoke at different times without wind

将烟尘源位置设置为侧面，烟尘初始释放速度设为3 m/s 且方向设为水平方向，其它初始条件保持不变。在有风情况下烟尘的时空运动如图8所示。

图8 有风情况下不同时刻烟的扩散图Fig.8 Spread of smoke at different times in the wind

由仿真结果可以看出，在风力作用下，烟尘的整体形态发生了变化，向着风场方向偏移。一般而言，风速越大，偏移越明显。经与无风情况进行对比还可看出，在有风与无风的条件下，烟尘对350 K 方形黑体遮蔽能力不同：在有风的情况下，烟尘迅速扩散，这导致对黑体的遮蔽时间和效果都明显下降。

4.2 不同波段下的仿真结果

设烟尘初始释放速度为1 m/s，释放浓度为0.1 kg/m3。量化范围上限设置为350 K所对应的黑体辐射亮度，下限设置为280 K。波段分别取3～5 μm 和8～12 μm，并在场景中加入与烟尘初始温度（350 K）相同的黑体方块，用于观察烟尘的遮蔽特性，得到仿真结果如图9所示，左侧是8～12 μm波段，右侧是3～5 μm。