红外成像仿真技术通过虚拟场景建模和数字成像计算生成红外相机的仿真探测图像，为红外成像探测系统的设计、研制、评估等工作提供测试与分析工具。随着红外成像探测系统的发展，对红外成像仿真技术的准确度和精细度提出了更高的需求。

红外成像仿真技术的发展主要经历了数值仿真、图形仿真和物理仿真三个阶段：数值仿真方法采用简化模型或经验模型，仅能够计算目标的整体辐射特性，对点或斑目标红外成像仿真时较为有效^[1-2]；图形仿真方法主要关注目标的图像形状和纹理特征，成像模型和仿真结果通常不具有实际辐射物理意义^[3]；物理成像方法对仿真场景的辐射特性进行物理建模，并采用数值方法对成像辐射传输过程进行解算，能够实现辐射能量与图形纹理的一致化高精度仿真，是目前红外成像仿真技术的重要发展方向^[4]。

光线跟踪是一种发展较为成熟的物理成像技术，已应用于国外许多成熟红外成像仿真软件的高精度仿真模块之中。法国OKTAL公司的SE-Workbench软件，在SE-RAY-IR模块中使用光线跟踪技术实现了融合三维目标、大气环境、传感器效应的复杂场景红外图像高精度渲染计算^[5]。美国的DIRSIG^[6]和MuSES^[7]，以及英国的CAMEO-SIM^[8]等红外成像仿真软件也将光线跟踪技术引入其高精度成像仿真计算程序当中。

我国在该领域的研究起步较晚，很多成果是基于国外现有仿真软件开展的应用性分析，缺乏对其中核心成像模型的理论与方法研究成果。西安电子科技大学的张建奇教授团队在OGRE渲染引擎的基础上开发了高真实感红外场景仿真系统，并对光线跟踪红外成像仿真的基础理论和可行性进行了分析^[9-10]；北京理工大学的王霞等人采用Optix并行GPU光线跟踪引擎进行了大场景红外成像仿真，实现了简单30万像素级别红外图像的高速实时渲染^[11-12]。以上研究成果仅能够生成不透光蒙皮目标的红外仿真图像，无法对飞行器目标中常见的尾焰辐射源进行光线跟踪红外物理成像仿真。北京理工大学的郑海晶等人采用蒙特卡洛方法对飞行器尾焰的红外辐射特性进行了仿真计算^[13]，通过尾焰内部的随机光子生成和传输路径采样实现红外辐射强度计算，计算误差优于15%时需要经过10⁷次采样，相较于光线跟踪方法计算效率过低。国内也未见关于光线跟踪红外成像仿真计算优化方法的研究成果，普遍采用的经典辐射传输方程和平均采样方法存在计算效率低、方差收敛速度慢等问题。

文中采用光线跟踪方法对飞行器类目标进行红外物理成像仿真，建立了飞行器红外辐射亮度成像传输模型，并针对蒙皮和尾焰的不同辐射特性提出了相应的光线跟踪计算优化方法。通过仿真实验对模型与方法的准确性和有效性进行了验证，在此基础上开展了隐身飞机红外探测图像仿真和辐射特性分析工作。文中形成的仿真模型和计算方法成果能够为我国自主红外成像仿真平台的发展提供理论支持与技术借鉴。

文中对飞行器类目标的红外入瞳辐射亮度图像进行仿真，仿真图像的像素值以辐射亮度为度量单位。入瞳辐射亮度图像不引入成像探测系统的调制、采样、噪声等影响因素，体现了飞行器红外辐射能量经场景传输后服从几何成像投影关系的空间分布状态。入瞳辐射亮度图像能够反应飞行器的亮度、形态、纹理等图像特征，可作为场景仿真的输入源，也可从中统计辐射强度等辐射特性参数。本节首先建立飞行器的红外辐射亮度成像传输模型，然后针对蒙皮和尾焰的不同辐射特性设计相应的光线跟踪辐射计算优化方法。

入瞳辐射亮度图像中像素值表示像素区域内成像系统的平均探测辐射亮度，可写为：

式中：$L\left( {{\tilde{ p}}} \right)$为图像中的任意像素位置${\tilde{p}}$的平均探测辐射亮度；$A\left( {{\tilde{ p}}} \right)$为像素面积；$L\left( {{{{p}}_0}} \right)$为像素范围内任意点${{{p}}_0}$的探测辐射亮度。

根据中心透视投影几何成像模型和辐射传输的互易性，将$L\left( {{{{p}}_0}} \right)$写为：

式中：${{\bf{\omega }}_{\rm{o}}}$为${{{p}}_0}$对应的成像视线方向矢量；${{{p}}_1}$为成像视线与飞行器目标的交点；${\tau _{{\rm{atm}}}}\left( {{{{p}}_1},{{{p}}_0}} \right)$为成像视线路径的大气透过率；${L_{\rm{o}}}\left( {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}}} \right)$为飞行器目标的出射辐射亮度。

如图1所示，蒙皮和尾焰是飞行器中的两类主要红外辐射源，需要根据它们的辐射特性差异建立不同的出射辐射亮度方程。飞行器蒙皮为不透光材料，出射辐射主要由自发辐射和反射辐射构成，出射辐射亮度方程可写为：

式中：等号右侧第一项${L_{{\rm{es}}}}\left( {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}}} \right)$表示蒙皮的自发辐射亮度；第二项表示蒙皮的反射辐射亮度，是蒙皮表面法线半球空间$\Omega $内沿${{\bf{\omega }}_{\rm{o}}}$方向的反射辐射亮度积分；${{\bf{\omega }}_{{\rm{is}}}}$为蒙皮反射的入射辐射方向；$L\left( {{{\bf{\omega }}_{{\rm{is}}}},{{{p}}_1}} \right)$为入射辐射亮度；$F\left( {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}},{{\bf{\omega }}_{{\rm{is}}}}} \right)$为反射函数。

图  1  飞行器红外成像辐射亮度传输示意图

Figure 1.  Schematic diagram of aircraft infrared imaging radiance transmission

飞行器尾焰的出射辐射主要由自发辐射和散射辐射构成，其出射辐射亮度方程需要考虑尾焰内部的辐射贡献和辐射传输的衰减，可写为：

式中：等号右侧第一项表示外部辐射穿过尾焰的透射辐射亮度；${{{p}}_2}$为成像视线穿过尾焰的交点；${\tau _{\rm{m}}}\left( {{{{p}}_2},{{{p}}_1}} \right)$为尾焰在${{{p}}_1}$和${{{p}}_2}$路径的透过率；$L\left( {{{\bf{\omega }}_{\rm{o}}},{{{p}}_2}} \right)$为外部辐射亮度；第二项表示尾焰自身的出射辐射亮度，是成像视线穿过尾焰路径上各点的辐射亮度积分；${L_{\rm{m}}}\left( {{{p}},{{\bf{\omega }}_{\rm{o}}}} \right)$为尾焰内部成像视线路径上任意点${{p}}$的出射辐射亮度，可分解为：

式中：等号右侧第一项${L_{{\rm{em}}}}\left( {{{p}},{{\bf{\omega }}_{\rm{o}}}} \right)$表示尾焰的自发辐射亮度；第二项表示尾焰的散射辐射亮度，是${{p}}$点球空间$\Psi $内沿${{\bf{\omega }}_{\rm{o}}}$方向的散射辐射亮度积分；${{\bf{\omega }}_{{\rm{im}}}}$为尾焰散射的入射辐射方向；$L\left( {{{\bf{\omega }}_{{\rm{im}}}},{{p}}} \right)$为入射辐射亮度；$S\left( {{{p}},{{\bf{\omega }}_{\rm{o}}},{{\bf{\omega }}_{{\rm{im}}}}} \right)$为散射函数。

公式（1）~（5）构成了飞行器的红外成像辐射亮度传输方程。其中的蒙皮反射辐射亮度和尾焰散射辐射亮度解算过程存在递归关系，即反射/散射的入射辐射亮度来自于其它位置的出射辐射亮度。通过光线跟踪方法对递归传输方程进行解算，其中的多重积分问题采用蒙特卡洛方法进行采样估计，主要步骤包括：（1）在公式（1）中对像素区域内的像面位置进行随机采样，通过计算采样位置的探测辐射亮度估计像素平均探测辐射亮度；（2）在公式（3）中对蒙皮反射的入射辐射方向进行随机采样，通过计算采样入射辐射产生的反射辐射亮度估计蒙皮的总反射辐射亮度；（3）在公式（4）中对尾焰路径位置进行随机采样，通过计算采样位置的尾焰出射辐射亮度估计尾焰的总出射辐射亮度；（4）在公式（5）中对尾焰散射的入射辐射方向进行随机采样，通过计算采样入射辐射产生的散射辐射亮度估计尾焰的总散射辐射亮度。

通过增加采样数量能够实现蒙特卡洛方法的积分估计方差收敛，但也存在经典蒙特卡洛方法收敛速度较慢的问题。针对飞行器蒙皮和尾焰的辐射特性，通过拆分出射辐射亮度方程和重要性采样等方法能够实现方差的快速收敛优化。

根据飞行器蒙皮的辐射特性模型可将公式（3）改写为：

式中：等号右侧第一项表示蒙皮的热辐射亮度；$\varepsilon \left( {{{{p}}_1}} \right)$为蒙皮的发射率；$B\left[ {T\left( {{{{p}}_1}} \right)} \right]$为$T$温度的黑体辐射亮度；第二项表示直接光照反射辐射亮度，是外部辐射源直接入射到${{{p}}_1}$点形成的反射辐射，在总反射辐射中占主要部分；${\Omega _{\rm{d}}}$为存在直接光照的入射方向集合；${f_{\rm{r}}}\left( {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}},{{\bf{\omega }}_{{\rm{is}}}}} \right)$为蒙皮材料的双向反射分布函数（BRDF）；${L_{{\rm{de}}}}\left( {{{\bf{\omega }}_{{\rm{is}}}},{{{p}}_{\rm{1}}}} \right)$为外部辐射源的出射辐射亮度；${\theta _{{\rm{is}}}}$为入射辐射方向${{\bf{\omega }}_{{\rm{is}}}}$与蒙皮表面法线方向的夹角；第三项表示间接光照反射辐射亮度，是从辐射源出射后经过两次以上反射形成的反射辐射，在总反射辐射中占比较小，可以采用简化模型进行快速计算^[14]，也可以展开为次级光线的热辐射和直接光照反射辐射进行递归求解。

公式（6）中的直接光照反射辐射项是光线跟踪蒙皮出射辐射亮度计算的主要方差来源，文中采用多重重要性采样方法来提高方差收敛速度。将直接光照反射辐射亮度${L_{{\rm{rd}}}}\left( {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}}} \right)$的表达式转换为辐射源区域积分形式：

式中：${{{p}}_{\rm{L}}}$为外部辐射源上的任意位置；${A_{\rm{L}}}$为外部辐射源面积；${{\bf{\omega }}_{{\rm{is}}}}\left( {{{{p}}_{\rm{L}}}} \right)$为${{{p}}_{\rm{L}}}$到${{{p}}_1}$的入射辐射方向；$\tau \left( {{{{p}}_{\rm{L}}},{{{p}}_1}} \right)$为入射路径环境透过率；$G\left( {{{{p}}_{\rm{L}}},{{{p}}_1}} \right) = {{\left( {\cos {\theta _{{\rm{is}}}}\cos {\theta _{\rm{L}}}} \right)} / {r_{{{{p}}_{\rm{L}}}{{{p}}_1}}^2}}$为几何因子；${\theta _{\rm{L}}}$为入射辐射方向与辐射源${{{p}}_{\rm{L}}}$位置表面法线的夹角；$r_{{{{p}}_{\rm{L}}}{{{p}}_1}}^{}$为${{{p}}_{\rm{L}}}$到${{{p}}_1}$的距离。

根据重要性采样思想，令采样分布概率密度函数与被积函数的形状近似能够提高方差收敛速度。公式（7）的被积函数形状主要受BRDF函数${f_{\rm{r}}}\left[ {{{{p}}_1},{{\bf{\omega }}_{\rm{o}}},{{\bf{\omega }}_{{\rm{is}}}}\left( {{{{p}}_{\rm{L}}}} \right)} \right]$和光照分布函数${L_{{\rm{de}}}}\left[ {{{\bf{\omega }}_{{\rm{is}}}}\left( {{{{p}}_{\rm{L}}}} \right),{{{p}}_1}} \right]$的共同决定，此时直接令入射辐射方向采样分布与被积函数相似较为困难，为解决此问题文中引入多重重要性采样方法^[15]。方法根据BRDF函数和光照分布函数的形状分别进行采样和被积函数计算，并通过自适应函数对积分估计结果进行加权，使采样分布概率密度分段适配被积函数形状来加快方差的收敛。设${{p'}}$为场景中任意的入射辐射起点，将公式（7）的被积函数简写为${l_{{\rm{rd}}}}\left( {{{p'}}} \right)$，多重重要性采样的蒙特卡洛直接光照反射辐射亮度计算方程为：

式中：${{{p'}}_m}$和${{{p'}}_n}$分别为服从BRDF函数和光照分布函数形状的随机采样入射辐射起点；${P_{\rm{f}}}\left( {{{{{p'}}}_m}} \right)$和${P_{\rm{l}}}\left( {{{{{p'}}}_n}} \right)$为两种采样分布的概率密度函数；${M_{\rm{f}}}$和${N_{\rm{l}}}$为两种采样分布的采样数量；${w_{\rm{f}}}\left( {{{{{p'}}}_m}} \right)$和${w_{\rm{l}}}\left( {{{{{p'}}}_n}} \right)$为多重重要性采样权重函数。

文中使用的二次自适应权重函数可表示为：

直接光照反射辐射亮度的多重重要性采样计算流程如图2所示，主要步骤包括：（1）对外部辐射源的出射位置进行随机采样，在入射光线无遮挡的情况下进行光源采样反射辐射亮度积分估计；（2）若外部辐射源服从Delta分布（点光源或平行光源），直接将光源采样估计结果作为直接光照反射辐射亮度；（3）若外部辐射源不是Delta分布，在反射点进行BRDF入射光线方向采样，通过递归计算入射光线的辐射亮度并进行BRDF采样反射辐射亮度积分估计；（4）基于光源采样估计结果和BRDF采样估计结果进行多重重要性采样积分估计解算，作为直接光照反射辐射亮度。

图  2  直接光照反射辐射亮度多重重要性采样计算流程

Figure 2.  Flow chart of direct lighting reflection radiance calculation using multi-importance sampling

忽略外部透射辐射贡献，根据飞行器尾焰的辐射特性模型将公式（4）和公式（5）改写为：

式中：$d$为成像视线穿过尾焰的路径长度；${{{p}}_t}$为尾焰内成像视线路径上与${{{p}}_1}$距离为$t$的点；${\tau _{\rm{m}}}\left( t \right) = {\tau _{\rm{m}}}\left( {{{{p}}_t},{{{p}}_1}} \right)$为路径的透过率；${\mu _{\rm{e}}}$为尾焰的发射系数；${\mu _{\rm{s}}}$为散射系数；${\mu _{\rm{a}}}$为吸收系数；${\mu _{\rm{t}}} = {\mu _{\rm{a}}} + {\mu _{\rm{s}}}$为消光系数；${s_{\rm{p}}}\left( {{{{p}}_t},{{\bf{\omega }}_{\rm{o}}},{{\bf{\omega }}_{{\rm{im}}}}} \right)$为尾焰的散射相函数。

尾焰的出射辐射亮度计算需要对散射辐射亮度积分和路径辐射亮度积分两处进行蒙特卡洛计算优化。散射辐射亮度的计算优化方法与蒙皮反射辐射亮度相似，将散射入射辐射分解为直接光照和间接光照两部分分别计算，直接光照散射辐射亮度计算采用多重重要性采样方法进行优化，间接光照散射辐射亮度使用简化模型或递归方程求解。下面重点介绍路径辐射亮度积分计算的优化方法。

尾焰的出射辐射主要由自发热辐射组成，为方便描述将公式（10）中的散射辐射项省略，并定义${L_{{\rm{me}}}}\left( t \right) = $$ {\mu _{\rm{e}}}\left( {{{{p}}_t}} \right)B\left[ {T\left( {{{{p}}_t}} \right)} \right]$，则有：

采用透过率指数函数分布对成像视线路径距离进行采样，采样概率密度函数为${P_{\rm{m}}}\left( t \right) = {\mu _{\rm{t}}}\left( t \right){\tau _{\rm{m}}}\left( t \right)$，公式（12）的蒙特卡洛计算方程可写为：

公式（13）是透光介质路径辐射亮度积分计算中的经典点积分方法，通过计算采样路径点处的被积函数值来实现积分结果的估计^[16]。点积分方法的计算方程形式简单，在自发射能力不强介质的路径辐射亮度解算时效果较好，但对飞行器尾焰类介质的辐射计算存在明显的不足。飞行器尾焰类介质具有自发热辐射能力强、消光能力弱的特点，此时公式（12）的被积函数形状主要受${L_{{\rm{me}}}}\left( t \right)$的影响，通常与透过率采样概率密度函数形状存在较大的差异。如图3所示，当介质消光系数较小时，成像视线的路径透过率随距离增加下降较慢，透过率采样获得的采样路径长度通常较远，积分估计过程容易忽略路径距离较近区域的辐射亮度贡献，在介质厚度较薄的情况下可能存在大量采样路径距离完全穿透发射介质的极限情况，造成较大的计算方差。

图  3  成像视线透过率采样的采样距离分布示意图

Figure 3.  Sampling distance distribution of line-of-sight transmittance sampling

针对点积分方法在飞行器尾焰路径辐射亮度积分计算时方差较大的问题，文中提出一种线积分计算思想。在每个采样位置处对整个路径距离的尾焰辐射亮度进行积分计算，并引入权重函数避免多次路径积分的重复计算问题，如图4所示。将路径辐射亮度积分的被积函数替换为$\displaystyle\int_0^t {{\tau _{\rm{m}}}\left( s \right){L_{{\rm{me}}}}\left( s \right){\rm{d}}s} $的形式，定义线积分权重函数为${w_t}\left( s \right)$，公式（12）可改写为：

采用透过率路径采样时公式（14）的蒙特卡洛计算方程为：

图  4  线积分方法路径辐射亮度积分计算原理

Figure 4.  Schematic diagram of path radiance integral calculation based on line integration method

线积分方法在薄介质辐射亮度计算时方差较小，但对消光系数分布均匀的厚介质辐射亮度计算时方差会随采样路径长度的增加而变大，而点积分方法对此种情况仅存在较小的恒定方差。文中参考多重重要性采样思想结合两种积分方法的优势提出了一种适用性更强的飞行器尾焰路径辐射亮度的混合积分计算方法，通过设置透过率相关权重分配函数使点积分计算结果在消光系数较大区域具有更大的贡献，使线积分计算结果在消光系数较小区域具有更大的贡献。混合积分方法的计算方程为：

公式（17）中等号右侧第一项采用点积分方法进行计算，第二项采用线积分方法进行计算，采用透过率路径采样的混合积分蒙特卡洛计算方程为：

飞行器尾焰的辐射特性空间分布状态采用三维轴向均匀网格进行建模，每个网格基元内部的辐射特性参数相同，此时混合积分方法的计算流程如图5所示，主要步骤包括：（1）随机采样成像视线的路径截止透过率${\tau _{{\rm{m}}\xi }}$，然后沿视线方向逐网格基元进行辐射亮度与透过率计算；（2）计算视线穿过当前网格基元后的透过率${\tau '_{\rm{m}}}$，若${\tau '_{\rm{m}}} > {\tau _{{\rm{m}}\xi }}$则累加当前网格基元的线积分辐射亮度贡献，然后重复进行下一网格基元计算；（3）若${\tau '_{\rm{m}}} \leqslant {\tau _{{\rm{m}}\xi }}$则在该网格基元内进行散射位置采样，计算尾焰的采样路径混合积分辐射亮度和透过率，并采用多重重要性采样方法对散射辐射亮度进行递归计算；（4）成像视线穿透尾焰或在尾焰内发生散射的情况下，路径采样截止，返回尾焰的出射辐射亮度计算结果。

图  5  尾焰出射辐射亮度混合积分计算流程

Figure 5.  Flow chart of exhaust plume radiance calculation using combining integration method

设计了数字仿真实验对文中提出的飞行器红外物理成像模型与光线跟踪计算方法进行正确性和有效性验证。

仿真实验的辐射强度计算相对误差（Relative Error, RE）表示为：

式中：$\hat I$为文中方法辐射强度计算结果；$I$为辐射强度理论计算结果。

仿真实验的辐射亮度仿真图像计算相对均方根误差（Relative Root Mean Square Error, RRMSE）表示为：

式中：${\hat L_n}$为图像中第$n$个像素的文中方法辐射亮度计算结果；${L_n}$为图像中第$n$个像素的辐射亮度理论计算结果；$N$为辐射亮度仿真图像的总像素数。

在数字仿真场景中构建面积为1 m²的正方形平板进行蒙皮出射辐射强度仿真计算实验，从仿真辐射亮度图像中计算成像视线方向的平板等效辐射强度，通过与理论计算结果的对比验证文中方法的计算准确性。辐射亮度仿真图像尺寸为1000×1000，如无特殊说明像素的采样光线数量为100条。

在仅考虑平板自发射特性的情况下，将平板的表面发射率$\varepsilon $依次设为0.5和0.9，表面温度设为400 K，对短波（1~3 μm）、中波（3~5 μm）和长波（8~12 μm）三个红外谱段的平板辐射强度随成像视线与平板法线夹角${\theta _{\rm{o}}}$的变化情况进行计算，文中方法与理论方法的计算误差如图6所示。

图  6  蒙皮自发辐射强度实验计算相对误差

Figure 6.  Calculation relative errors of emissive radiation intensity experiments of skin

在仅考虑平板反射特性的情况下，将其平板的表面材料BRDF设为朗伯模型，其总反射率$\rho $依次设为0.5和1.0，入射平行辐射源的辐射亮度设为100 W/(sr·m²)，入射辐射方向与平板法线的夹角${\theta _{\rm{i}}}$依次设为0º、30º、60º，文中方法与理论方法的平板辐射强度计算误差随${\theta _{\rm{o}}}$的变化情况如图7所示。

图  7  蒙皮反射辐射强度实验计算相对误差

Figure 7.  Calculation relative errors of reflective radiation intensity experiments of skin

文中方法的自发辐射强度计算绝对误差最大值为0.043 W/sr，相对误差最大值为0.048%；反射辐射强度计算绝对误差最大值为0.005 W/sr，相对误差最大值为0.023%。正方形平板辐射强度仿真计算结果说明了文中蒙皮出射辐射亮度计算方法的正确性。

将数字仿真场景中的平行辐射源替换为半径1 m的漫发射圆盘辐射源进行蒙皮辐射亮度图像仿真计算实验。圆盘辐射源放置于平板的右上方位置，圆盘法线与平板法线平行，辐射源出射辐射亮度设为100 W/(sr·m²)，平板的表面材料BRDF设为朗伯模型，总反射率$\rho $设为0.5。采用文中直接光照多重重要性采样方法（简称直接光照采样）和经典反射方程BRDF采样方法（简称BRDF采样）的仿真图像随像素采样光线数的变化情况如图8所示，仿真图像的相对均方根误差如表1所示。

图  8  蒙皮反射辐射亮度仿真图像

Figure 8.  Reflective radiance simulation images of skin

圆盘辐射源对平板不同位置的入射辐射方向和辐射亮度存在差异，仿真图像呈现为自右上至左下逐渐变暗的渐变纹理，如图8中参考图像所示。反射辐射亮度方程的被积函数形状主要由光源分布函数决定，与BRDF函数形状的差异较大。表1说明此时直接光照采样方法比BRDF采样方法具有更快的方差收敛速度，图8中直接光照采样方法的仿真图像更为精细，验证了文中飞行器蒙皮辐射亮度图像仿真方法的有效性。

在数字仿真场景中分别构建XYZ三轴方向长度均为2 m的厚介质和XY轴方向长度为2 m、Z轴方向长度为0.5 m的薄介质进行尾焰出射辐射强度仿真计算实验。介质的正方体网格基元边长为5 mm，介质内部辐射特性参数均匀分布且不考虑散射效应。厚介质的吸收系数和发射系数均为0.1 cm⁻¹、黑体辐射亮度为10 W/(sr·m²)，薄介质的吸收系数和发射系数均为0.0001 cm⁻¹、黑体辐射亮度为10000 W/(sr·m²)。分别采用点积分方法及文中提出的线积分方法和混合积分方法对两种介质的等效辐射强度进行计算，仿真图像尺寸为1000×1000，如无特殊说明像素的采样光线数量为100条。设成像视线方向与介质+Z轴方向的夹角为${\theta _{\rm{o}}}$，三种积分方法与逐网格积分理论方法的计算结果对比如图9和图10所示。

图  9  厚介质自发辐射强度实验计算相对误差

Figure 9.  Calculation relative errors of emissive radiation intensity experiments of thick medium

图  10  薄介质自发辐射强度实验计算相对误差

Figure 10.  Calculation relative errors of emissive radiation intensity experiments of thin medium

厚介质的吸收系数较大，辐射能量的传输衰减较快，透过率采样的路径距离不会超出介质厚度。此时点积分方法的计算误差能始终保持在较低水平，辐射强度计算绝对误差的最大值为0.004 W/sr，相对误差最大值为0.007 %。线积分方法的计算方差与采样距离有关，因此在厚介质中计算误差较大，辐射强度计算绝对误差的最大值为0.022 W/sr，相对误差最大值为0.050 %。混合积分方法能够在一定程度上弥补线积分方法的不足，计算误差介于点积分和线积分之间，辐射强度计算绝对误差的最大值为0.013 W/sr，相对误差最大值为0.027 %。

薄介质的吸收系数较小，辐射能量的传输损失较小，同时介质厚度较薄也使透过率采样存在大量路径距离超过介质厚度的情况。采用点积分方法会出现大量误差较大的零值估计结果，辐射强度计算绝对误差的最大值为1.264 W/sr，相对误差最大值为0.634 %。线积分方法能够考虑采样路径范围内的所有辐射贡献，有效降低采样路径距离超出介质厚度造成的计算误差，辐射强度计算绝对误差的最大值为0.220 W/sr，相对误差最大值为0.110 %。混合积分方法的计算误差与线积分方法相当，但误差随探测角度变化的起伏更小，辐射强度计算绝对误差的最大值为0.135 W/sr，相对误差最大值为0.068 %。

在数字仿真场景中构建X轴方向长度为4 m、Y轴方向长度均为1 m、Z轴方向长度为0.2 m的渐变介质进行尾焰辐射亮度图像仿真计算实验。介质的正方体网格基元边长为5 mm，介质的黑体辐射亮度均匀分布为100 W/(sr·m²)，吸收系数和发射系数在0~0.004 cm⁻¹范围内沿+X轴方向逐渐增大，在YZ轴方向均匀分布。视线沿+Z轴方向对渐变介质进行成像仿真，三种积分方法仿真图像随像素采样光线数的变化情况如图11所示，仿真图像的相对均方根误差如表2所示。

图  11  尾焰自发辐射亮度仿真图像

Figure 11.  Emissive radiance simulation images of exhaust plume

随着介质吸收系数和发射系数的逐渐增大，仿真图像的辐射亮度自左向右逐渐增加，图9中线积分方法和混合积分方法的仿真图像噪声更小，仅需要少量的采样光线数即可实现较为精细的图像仿真效果，与参考图像基本无目视差异。表2说明文中提出的混合积分方法具有最快的方差收敛速度，在像素采样光线数量为1000时仿真图像的相对均方根误差仅为0.020 %，验证了文中飞行器尾焰辐射亮度图像仿真方法的有效性。

采用文中飞行器红外物理成像仿真方法对F35隐身飞机进行红外辐射成像仿真计算，并对其红外辐射特性进行分析。采用Fluent软件对高速巡航飞行状态的飞机流场进行仿真，仿真计算的环境参数如表3所示。

飞机蒙皮表面温度场和尾焰流场仿真计算结果如图12所示。飞机蒙皮的长度为15.67 m、高度为2.19 m、翼展为10.72 m，蒙皮表面的最大温度约为360 K。飞机尾焰区域的最高温度约为2000 K，垂直喷射方向最大直径约为1.2 m，沿喷射方向最大长度约为30 m，垂直喷射方向的最大投影面积约为32 m²。

图  12  隐身飞机流场仿真计算结果

Figure 12.  Simulation calculation results of stealth aircraft flow field

分别对采用红外隐身技术和未采用红外隐身技术两种情况的飞机辐射特性参数进行配置：未采用红外隐身技术条件下，蒙皮的红外发射率为0.8、反射率为0.2，尾焰仅由燃料完全燃烧后的高温气体构成；采用红外隐身技术条件下，蒙皮的红外发射率为0.2、反射率为0.8，在发动机喷口四周添加固体粒子喷射装置，采用气溶胶红外隐身技术对尾焰红外辐射进行抑制，粒子直径为10 μm、粒子注入流量为2.0 kg/s、粒子材料复折射率为2+5×10⁻⁵ i。尾焰中气体分子的辐射特性参数基于HITEMP数据库的分子谱线数据进行计算^[17]，固体粒子的辐射特性参数计算采用Mie散射模型^[18]，并采用逐线方法进行谱段积分^[19]。

基于上述飞机红外辐射特性建模与仿真计算结果，选择短波（1~3 μm）、中波（3~5 μm）、长波（8~12 μm）三个常用红外探测宽谱段进行隐身飞机入瞳辐射亮度图像仿真计算。分别设计了地基探测和空间探测两种仿真环境场景，场景中的大气环境参数均采用MODTRAN标准大气模型；地基探测场景中红外相机高度为0 km，飞机高度为10 km，探测距离为11.547 km，如图13（a）所示；空基探测场景中红外相机和飞机高度均为10 km，探测距离为5 km，如图13（b）所示。红外相机的探测器像素尺寸为500×500，每像素采样光线数为1000条，反射/散射递归深度为10次，地基探测相机的角分辨率为0.005 mrad，空间探测相机的角分辨率为0.01 mrad，使两种探测场景下仿真图像中的飞机最大投影面积均约为35000个像素。红外相机沿目标本体坐标系三个基准面绕坐标轴对目标进行成像仿真计算，探测角度间隔为5º，并从辐射亮度图像中计算其等效红外辐射强度，如图13（c）所示。

图  13  隐身飞机仿真成像几何关系。（a）地基探测场景；（b）空基探测场景；（c）探测视线方向

Figure 13.  Geometric relationship of stealth aircraft simulation imaging. (a) Ground-based detection scene; (b) Air-based detection scene; (c) Line of detection sight

不同谱段、不同隐身条件、不同探测场景的隐身飞机红外辐射强度仿真计算结果如图14所示。低层大气环境的气体压强和吸收气体分子浓度较高，地基探测场景的大气透过率小于空基探测场景，因此三个谱段的地基探测飞机辐射强度比空基探测更小。红外隐身技术能够抑制飞机的自发热辐射能量，但也会在一定程度上增加飞机的反射辐射能力，短波谱段在能够接收太阳光照的探测方向出现了隐身条件飞机辐射强度高于未隐身条件的现象；中波和长波谱段的太阳光照辐射能量可以忽略，采用红外隐身技术后飞机的辐射强度均有所下降。总的来说，未采用红外隐身技术的情况下，飞机的辐射强度在长波最大、中波次之、短波最小；采用红外隐身技术后，飞机的辐射强度在长波最大、短波次之、中波最小。

图  14  隐身飞机红外辐射强度仿真计算结果

Figure 14.  Simulation calculation results of infrared radiation intensity of stealth aircraft

不同条件的隐身飞机红外入瞳辐射亮度仿真图像如图15所示，相同谱段图像的量化系数相同。红外隐身材料抑制了飞机蒙皮的自发辐射能力并增强了其反射能力，短波谱段中隐身情况下的飞机蒙皮图像更为清晰，中波和长波谱段中未隐身情况下的飞机蒙皮图像更为清晰。固体粒子气溶胶增强了尾焰的散射能力，短波和中波谱段仿真图像中隐身情况下的飞机尾焰区域的像素辐射亮度有所降低，但增加了尾焰在图像中的投影面积。长波谱段不覆盖尾焰高温气体的发射谱线，尾焰的辐射能量贡献很小，隐身和未隐身情况下仿真图像中均无法观察到清晰的尾焰轮廓。

图  15  隐身飞机红外辐射亮度仿真图像

Figure 15.  Infrared radiance simulation images of stealth aircraft

文中建立了飞行器红外成像的辐射亮度传输模型，并采用光线跟踪方法实现了飞行器红外入瞳辐射亮度图像的仿真计算。提出了直接光照多重重要性采样蒙皮辐射亮度计算方法和混合积分尾焰辐射亮度计算方法，有效提高了光线跟踪计算的方差收敛速度。通过仿真实验对文中的模型与算法成果进行了正确性和有效性验证，蒙皮辐射强度计算相对误差小于0.048 %，尾焰辐射强度计算相对误差小于0.068 %，每像素采样光线数为1000条时蒙皮和尾焰仿真图像的相对均方根误差分别为1.785 %和0.020 %。最后采用文中方法对隐身飞机红外辐射特性进行仿真分析，仿真结果表明红外隐身技术能够降低中波和长波谱段的飞机辐射强度，但在短波谱段增加飞机的反射辐射强度。

下一步将在文中理论方法基础上开展红外物理成像仿真引擎的研发工作，并结合相关场景辐射特性建模成果与红外探测系统仿真模型逐步开展我国自主红外成像探测场景仿真软件系统的研制工作。