HTML
-
随着红外探测器件性能的提高以及大规格、多波段红外焦平面阵列技术的日益成熟，红外探测技术在目标探测预警、跟踪识别等方面得到广泛应用，而通过数值计算的方法，快速准确的仿真飞机、导弹等空中目标蒙皮和尾焰辐射经大气传输后到达探测器形成的红外图像，是开展目标红外特性研究、评估红外探测设备性能的有效手段[1]。
蒙皮和尾焰是空中目标红外辐射的主要来源，蒙皮辐射计算采用普朗克公式计算，尾焰辐射计算相对复杂，又分为辐射组分光谱辐射特性计算和辐射传输计算两部分互相关联的内容。在HITEMP、CDSD等光谱数据库的基础上，国内外学者采用逐线计算、谱带模型等辐射特性计算方法和离散坐标、蒙特卡洛等辐射传输方法开展了计算研究和验证工作[2-4]，基于OpenGL等软件包，张建奇团队[5-6]、黄建峰[7]等实现了红外图像的生成。
CUDA是NVIDIA公司提供的GPU并行编程环境，国外许多学者已成功用于大气、云层等的辐射传输计算[8-9]，国内西安电子科技大学的Guo Xing[10]实现了空中目标对太阳和环境反射辐射的并行计算，中国科学院大学的蔡方伟[11]实现了目标和背景在探测器入瞳处的辐射照度的并行实时计算，以上研究均表明：CUDA的计算速度可提高百倍以上，是一种有效提高计算效率的并行方法。
文中以空中目标蒙皮和大尺度三维尾焰的红外辐射特性为研究对象，以目标蒙皮温度分布和三维尾焰流场作为输入，根据蒙皮和尾焰不同的几何成像关系并结合红外辐射算法建立了目标二维成像算法，解决了蒙皮对尾焰辐射的遮挡、以及蒙皮辐射穿过尾焰流场时的辐射传输问题，开发了红外辐射计算软件，最终实现目标辐射经大气传输后到达探测器入瞳前的红外图像仿真，并对计算结果进行了分析。此外，在程序中的蒙皮投影计算和尾焰光谱辐射计算部分采用CUDA并行计算方法，提高了计算速度，并分析了不同运算量下的并行计算加速比。
-
已知目标蒙皮温度后，可由普朗克公式计算蒙皮辐射亮度：
式中：
${L_\lambda }$ 为蒙皮光盘辐射亮度；${\varepsilon _\lambda }$ 为蒙皮光谱发射率；${c_1}$ 、${c_2}$ 为第一、第二辐射常数；$\lambda $ 为波长；$T$ 为蒙皮温度。采用单线组谱带模型（SLG）[12]计算尾焰气体辐射，微元中气体的光学厚度计算公式为：
式中：
${X_{\lambda,{\rm{tot}}}}$ 为气体组分的总光学厚度；下标i表示尾焰中的第i种辐射气体，文中主要考虑CO2、H2O和CO气体；${Y_i}$ 为谱线展宽作用下的光学厚度参量；${X_i}^*$ 为弱线近似的光学厚度参量；${X_{C,i}}$ 为碰撞展宽影响的光学厚度参量值；${X_{D,i}}$ 为多普勒展宽影响的光学厚度参量值。采用视在光线法（LOS）计算尾焰红外辐射传输[13]问题，LOS方法的基本思想是将三维非均匀介质中的辐射传输简化为一维多层介质的辐射传输问题：射线沿探测方向与流场相交，将相交部分沿射线方向分层，并假定同层介质均匀、等温，则沿着射线方向逐层积分即可计算该探测方向上的红外辐射与透过率。
如图1所示，LOS射线上的红外辐射传输方程的离散表达式为：
式中：
${I_\lambda }$ 为LOS出射光谱辐射；${I_\lambda }^0$ 为LOS入射光谱辐射，如壁面、环境等背景的光谱辐射，若不考虑背景辐射，则该项取零；$I_{b,\lambda }^{n - 1}$ 为第n-1层气体的黑体光谱辐射；$\sum\limits_{j = n}^N {X_{\lambda,{\rm{tot}}}^j} $ 为第n层至LOS射线入射点的气体总光学厚度。
Figure 1. Schematic of the LOS method
LOS射线上的光谱透过率计算公式可离散为：
-
由于蒙皮和尾焰遮挡关系复杂，如图2所示，尾焰流场可能会对蒙皮辐射存在遮蔽影响，因此，根据探测器成像几何关系以及目标蒙皮、尾焰成像特性，采用不同的成像方法，实现探测器入瞳处的目标红外图像计算。其中，大气辐射传输通过调用Modtran计算。
Figure 2. Schematic of target infrared imaging calculation
-
如图2中A、B线所示，目标蒙皮成像采用正向光线追迹的计算方法，即以蒙皮各面元为计算起始点，计算蒙皮红外辐射成像。能量上，计算面元光谱辐射经尾焰和大气辐射传输后，到达探测器入瞳处的表观辐射，空间上，由面元方向和中心坐标计算在探测器上的成像位置，并将落入同一像元的蒙皮辐射累加，从而实现目标蒙皮的红外辐射成像计算，此外，由于蒙皮边缘不能充满探测器的单个像元视场，同时计算蒙皮在探测器各像元上的成像面积用于后续成像计算。
面元红外辐射经尾焰和大气辐射传输后到达探测器入瞳前的表观辐射亮度
${L'_{t,\lambda }}$ 计算公式为：式中，
${L_{t,\lambda }}$ 为面元的真实光谱辐射亮度；${\tau _{p,\lambda }}$ 和${L_{p,\lambda }}$ 分别为面元至探测器路程上的尾焰透过率和红外辐射，当路程不穿过尾焰时，${\tau _{p,\lambda }}$ =1，${L_{p,\lambda }}$ =0，路程穿过尾焰时，由公式（4）、（5）计算路径上的红外辐射和透过率；${\tau _{\rm atmo,\lambda }}$ 和${L_{\rm atmo,\lambda }}$ 分别为面元至探测器路程上的大气透过率和红外辐射，将公式（4）代入公式（6）中，此时面元的表观辐射亮度计算公式可以表示为：若简单地以面元中心点或顶点坐标计算面元落入哪一探测器像元，当探测器象元分辨率不能远大于面元尺度时，如图3中（a）例所示，采用简单投影算法会导致应分布在10个像元内的辐射能量集中于1个像元，导致计算误差很大，因此，采用将面元细分为更小面元的方法，如图3中（b）所示，使细分面元远小于探测器的像元分辨率，再采用简单投影的算法计算面元成像，提高计算准确性。
Figure 3. Schematic of subdivision surfaces projection imaging
-
如图2中C线所示，发动机尾焰成像采用反向光线追迹的计算方法，即以探测器各像元为计算起始点，以像元探测方向为指向生成LOS射线，能量上，首先计算该射线上的尾焰真实光谱辐射，然后计算尾焰辐射经大气传输后到达探测器入瞳处的表观辐射，空间上，对蒙皮位置的成像像元，当蒙皮在该像元的成像面积占像元应成像面积的比例高于一定值时，由于计算蒙皮辐射时已考虑尾流辐射传输影响，不再重复计算，当蒙皮成像面积占比小于一定值时，计算尾焰的表观光谱辐射，并根据蒙皮辐射及蒙皮与尾焰面积的各自占比，采用加权平均的方式计算该像元的红外辐射亮度：
式中：
$\eta $ 为蒙皮的成像面积占比；${L'_{t,\lambda }}$ 为蒙皮的表观光谱辐射亮度；${L'_{p,\lambda }}$ 为尾焰的表观光谱辐射亮度。其中，蒙皮成像面积占比的定值在程序中取0.99，即当目标机体的成像面积占整个成像像元面积的99%以上时，不再考虑剩余1%部分的辐射影响，由此引起的计算误差一般小于1%。
此外，为捕捉尾焰核心区的精细辐射分布，根据核心区贴近于本体的位置特征，对相应位置的探测像元进行加密处理以提高计算精度。
-
CUDA是CPU和GPU混合计算的平台，在CUDA架构下，一个程序分为主机端和设备端两部分，其中主机端代码在CPU上执行，设备端代码则在GPU上执行，一个CUDA程序中可以依次执行多个设备端代码。由于显卡通过PCI-E接口与主机互联，并不能直接访问主机内存，因此设备端代码运行前需先将计算参量拷贝由主机内存拷贝至显卡内存，计算结束后再把计算结果由显卡内存拷贝至主机内存。
目标红外成像计算需要在面元网格、光谱等层面上进行大量计算，适合采用硬件调度模式为SPMD（Single Program Multiple Data，即每个线程执行相同的代码，但每个线程调用的数据并不相同）的GPU并行方法提高计算速度。文中基于VS 2013（C++）+CUDA 8.0软件平台开发了目标红外成像计算软件，程序流程如图4所示。
Figure 4. Flow chart of target infrared imaging program
程序中，在蒙皮投影计算和尾焰光谱辐射计算部分采用CUDA并行，其并行层面分别为面元细分数和光谱数，即分别建立与面元细分数和光谱数量相同的GPU线程进行计算。
-
基于建立的目标红外辐射成像计算软件，以某型号飞机为例计算了其红外图像，对计算结果进行了验证，并对CUDA并行加速性能进行分析。
-
飞机蒙皮温度分布和尾焰流场作为计算输入，如图5所示，图5中左图为飞机蒙皮温度分布，右图上为尾焰流场的温度分布，右图下为尾焰流场CO2浓度分布。
Figure 5. Temperature of aircraft skin and flowfield of engine plume
设定观测方位角、俯仰角均为40°，软件计算结果如图6所示。整体上，蒙皮和尾焰的辐射分布与其温度分布一致，另取机身长度和翼展长度为特征值，目标成像的尺度与特征值大小误差在一个像素的尺度范围内，表明投影成像算法准确。
Figure 6. Aircraft infrared radiance image of mid-wave
依据Avital等人对BEM发动机尾焰的地面红外辐射实验数据[14]对尾焰辐射计算进行了验证，选取与实验相同的边界条件计算了发动机流场，使用建立的红外成像软件计算1.5~5.5 μm波段内的尾焰光谱辐射强度，计算结果与实验结果的对比如图7所示，整体上计算结果与实验结果吻合较好，其中，3.5 μm为HCl气体辐射带，由于辐射计算所用数据库光谱分辨率较测量仪器更高，因此，计算结果较实验结果起伏变化较大，但辐射总量基本符合。
Figure 7. Comparison of the plume radiation data between calculated and measured
-
如图4所示，由于程序中既有串行代码部分，又有并行代码部分，对软件的整体计算速度进行统计分析并不能准确体现CUDA并行的加速能力，因此，以本体网格投影子程序为对象，统计该子程序的运行时间，编写相同功能、相同算法的串行程序并统计运行时间，以此分析CUDA的加速性能。测试电脑CPU为i7-7700，主频为3.6 GHz，显卡为NVIDIA GTX 1060，频率1.6 GHz，显存6 G，CUDA核心数为1280。
将面元细分为不同数量，记录相应串行程序和并行程序的运行时间，并定义加速比为二者运行时间之比，结果如表1所示。
从表1中可以看出，并行程序计算耗时始终小于串行程序，且随着面元细分数的增加，加速比增大，当面元细分数达50 000时，加速比达563.91，CUDA并行加速效果明显。此外，当提分面元数从500增大至50 000时，串行程序计算耗时增大100倍，与面元细分倍数一致，而并行程序计算耗时增大约2.4倍，原因在于并行程序中，需要将计算参数从主机内存拷贝至显卡内存，GPU计算完成后还需将计算结果从显卡内存拷贝至主机内存，两次内存拷贝会耗费一定的时间，而在统计并行程序计算时间时，将内存拷贝时间统计在内，导致并行耗时的增大倍数与计算量不一致。
Subdivision counts Serial program runtime/s Parallel program runtime/s Ratio 500 2.48 0.18 13.78 2 000 9.98 0.19 52.53 5 000 24.82 0.21 118.19 10 000 49.51 0.22 225.05 50 000 248.12 0.44 563.91 Table 1. CUDA computation speed ratio of skin projection module
-
文中根据蒙皮和尾焰的成像特性和辐射计算方法，建立了空中目标红外图像计算方法，其中，采用正向光线追迹方法计算目标蒙皮的红外投影成像，采用反向光线追迹方法计算尾焰的红外辐射成像，并通过细分加密的方式提高了成像精细度。基于VS 2013+CUDA 8.0软件平台开发了目标红外成像计算软件，可以用于计算目标真实和表观光谱红外图像，对软件计算结果的准确性进行了验证，并分析了CUDA并行计算的加速能力，计算结果表明：软件可以正确生成设定观测角度下的目标红外图像，尾焰辐射计算结果与实验结果符合较好，CUDA并行可以显著提高计算速度，尤其当计算量较大时，并行程序的加速比可达500以上。
软件生成的红外图像可用于目标红外特性分析，也可为目标的红外探测、识别提供参考。
DownLoad: