飞机的红外图像仿真

王 彪，童中翔，王超哲，马 榜

(空军工程大学航空航天工程学院，陕西 西安710038)

1 引言

随着红外制导技术的不断发展，红外制导导弹的抗干扰能力不断提高，目前最新一代的红外制导导弹普遍采用了红外成像制导技术［1］。红外成像制导探测器能够利用目标的形状特征、灰度分布规律、运动特点等诸多物理特征形成对目标的识别，进而形成制导信号，其抗干扰能力得到了极大提高，因此，对目标的红外图像特征进行研究已变得十分重要。

目前红外图像的获取方式有三种，一是在真实条件下进行实际拍摄，二是在地面进行模拟测试，三是理论建模与仿真。理论建模与仿真法与前两种方法相比，其所需经费较少，且不易受场地、设备、天气等因素的限制，可以在一定的精度下模拟出飞机在任意飞行条件下的红外图像，在研究目标红外图像的过程中具有重要意义。

针对飞行器的红外图像仿真方法，广大专家学者进行了广泛的研究。在机体红外图像方面，文献［2］借助Fluent软件生成了蒙皮的红外图像，文献［3］、［4］通过建立热平衡方程生成了机体红外图像，但建立的模型不能反映飞机外形对气动加热的影响。在排气系统红外图像方面，文献［5］对发动机热喷流进行了红外成像仿真，文献［6］采用粒子系统对尾焰图像进行了仿真研究。然而由于飞机在空中受到多种因素的影响，其红外辐射模型的建立比较复杂，模型的精确性与实时性往往不能同时满足，因此在建立仿真模型的过程中，如何快速地生成满足精度要求的图像是模型的关键。本文将整机辐射分解为机体辐射与排气系统辐射两部分，建立了一种可用于任意视线方向的飞机红外图像的生成模型，并利用OpenGL图形库生成了可视化的红外图像，通过验证模型计算速度较快，且满足一定的精确性。

2 机体红外辐射计算

2．1 蒙皮温度场建模

飞机在空中飞行的过程中要受到高速气流的气动加热、内热源加热、环境辐射热量与蒙皮向外辐射热量的综合作用，如图1所示，使得机体的各部位呈现出不均匀的温度分布，因此，求解蒙皮的温度场是生成机体红外图像的关键。

图1 飞机热交换示意图

本文中首先将机体表面分成n个面元，认为每个面元具有均匀的温度与辐射量。综合考虑各因素对面元温度变化的影响，建立式(1)的热平衡方程。

式中，mi为面元的质量;ci为面元的比热容;Ti为面元温度;qcv，qin，qenv，qcon，qrad分别为气动加热、内热源加热、环境辐射热、相邻面元之间的热传递与蒙皮向外辐射热量的热流密度。

由于机体蒙皮很薄，面元之间的热传导面积很小，且蒙皮温度场是连续的，相邻面元之间的温差不大，因此忽略面元之间的热传递，即认为:

2．1．1 气动加热

高速气流流过飞机表面时，会因气流压缩与摩擦产生大量的热量，热量大小与面元所在位置的机体形状相关。根据机体各部位形状不同，将机体分为机身、机头、翼前缘、翼面四个部分分别求解其气动加热大小。

a)机身及翼面气动加热计算

机身可近似等价为圆柱形，气流的流动方向沿圆柱轴向，其热流密度可按照平板气动加热计算。翼面的热流密度也可近似按照平板气动加热公式计算，其热流密度公式为:

式中，Ti为面元i的温度;Tr与αox分别为恢复温度和热交换系数，其计算式为:

式中，Te，ue，ρe分别为附面层边缘气体的温度、速度和密度;pr为普朗特数;cp为空气的定压比热;Rex为特征长度x处的雷诺数;ρ*是参考温度T*对应的密度，可由气体状态方程求解:

其中，R为气体常数;T*可取Eckert参考温度:

b)机头气动加热计算

机头可近似等价为圆锥形，其热流密度公式和平板热流密度公式的区别在于换热系数不同，圆锥的换热系数公式为［7］:

c)翼前缘气动加热计算

翼前缘可近似等价为一个无限后掠圆柱，其驻点线上的热流密度qor与半球的驻点热流密度qob有如下近似关系:

半球的驻点热流密度qob可按照Kemp－Riddle经验公式计算:

式中，R0为半球半径;ρsl为海平面大气密度;Ts为驻点温度。

对于后掠机翼，当后掠角在一定范围内时(0°≤Λ≤60°)，对后掠角的修正因子如下［8］:

2．1．2 内热源加热

内热源对蒙皮面元的加热以两种方式进行:热辐射与通过金属结构的热传导。将内热源看作是温度为T0、发射率为ε0、有一定形状的辐射灰体，通过第一种方式辐射到面元i的热流密度为:

式中，αi为蒙皮面元的吸收率;F0－i为内热源到面元i的辐射角系数。

内热源通过第二种方式传导到面元i的热流密度为:

式中，λc为综合导热系数;l0为内热源到面元的热传导距离。

则内热源对蒙皮加热的总热流密度为:

2．1．3 环境辐射热

对于高空飞行的飞机，天空散射太阳辐射、地面辐射及地面反射太阳辐射影响不大，本文中予以忽略，环境辐射对蒙皮温度的影响只考虑太阳直射辐射与大气辐射，即:式中，Esun为太阳常数(Esun=1353 W/m2);τh为从大气层外边界到飞行高度h的路径上大气对太阳辐射的透过率，计算方法可参考文献［9］;θs为太阳直射光线与面元法线方向的夹角，当0≤θs＜π/2时，qsun按式(17)计算，π/2＜θs≤π时，qsun=0。

2．1．4 蒙皮自身辐射热量Ti

将蒙皮看作发射率为εi，温度为Ti的辐射灰体，则蒙皮面元自身辐射热量为:

2．1．5 热平衡方程求解

根据以上描述，可完成式(1)中热平衡方程的建立，方程的右侧实际上是Ti的函数，其中含有T4

i、Ti、T－1i项及与Ti无关的常数项，记作:

其中，C1，C2，C3，C4为与Ti无关的系数。

因为本文所求的是稳态时的温度场，因此令dTi/dt=0，得到以下方程:

方程(20)虽然是一个5次方程，但根据其物理意义，方程有唯一解，且Ti＞T∞。因此在求解时，设定合适的门限值σ与步长B，以T∞为初始值，当满足式(21)时，认为T∞+kB即为方程(20)的解。

2．2 机体红外辐射亮度计算

得到蒙皮温度场后，将蒙皮看作灰体，应用普朗克定律和基尔霍夫定律可得蒙皮自身辐射的红外辐射亮度，其计算公式为:

式中，ε为蒙皮发射率，与蒙皮材料有关，本文取ε=0．8;c1与c2为辐射常数。

蒙皮的辐射亮度不仅与其自身辐射有关，还与反射的环境辐射有关，本文主要考虑蒙皮面元对太阳辐射的反射，其反射辐射亮度为:

式中，θD为面元法向与视线方向的夹角;其余各参数与式(16)一致;F(λ1，λ2，Tsun)为黑体辐射函数，其计算式为:

式中，M(λ，Tsun)为波长为λ，温度为Tsun的黑体光谱辐射出射度，可由普朗克公式计算得出，Tsun=5900 K。

蒙皮面元的红外辐射亮度等于自身辐射亮度与反射太阳辐射亮度之和，即:

3 排气系统红外辐射计算

3．1 尾焰的温度和组分分布计算

本文采用半经验法计算尾焰流场，将尾焰分成初始段与主体段两部分，其中初始段可分为核心区与混合区，如图2所示。设核心区长度为LC，外边界锥角为α，尾喷管半径为R0，初始段内边界半径为rc，尾焰外边界半径为rm，飞机飞行速度为Vplane，尾焰流速为Vnozzle，当地声速为Vsonic。

图2 尾焰流场模型

各几何参数的计算公式如下:核心区长度:

边界锥角(α0=8．3°):

本文将大气简化成由CO2、H2O、N2、O2四种气体构成，根据涡流的传播规律，对于坐标为(x，r)的点，其温度计算式为:

其中，T0为发动机喷口的气体温度;Ta为飞机所在环境的大气温度;Tw为x轴上温度，其计算式为:

气体的重量分数计算式为:

式中，gi为尾焰气体第i种组分的重量成分;g0i为发动机喷口处第i种气体的重量成分;gai为飞机所在环境的大气中第i种气体的重量成分。

3．2 排气系统红外辐射亮度计算

由于红外成像设备在某一时刻拍摄到的尾焰红外图像实际上是一幅二维图像，图像中每个像素点的灰度值取决于尾焰经由包络面发射向成像设备的红外辐射亮度大小，因此只需对尾焰包络面划分网格，如图3所示，网格的数量取决于计算的精度。

图3 尾焰包络面网格

计算前首先需要对每个结点进行可见性判断，判断方法如下:设结点P的坐标为(xp，yp，zp)，外法向为 珗np=(npx，npy，npz);探测器所在位置的坐标为(xD，yD，zD)，法向为 珗nD=(nDx，nDy，nDz)，探测器接收角为φ，点P到探测器的光束方向为珗ndir，如图4所示。若 珗np与 珗ndir的内积(珗np，珗ndir)＜0且 珗ndir与 珗nD的夹角满足:

则点P可见，反之点P不可见。

图4 可见性判断

若点P可见，求取珗ndir所在直线与包络面的另一个交点Q，交点Q有两种情况，一是Q点在尾焰扩散边界面上，二是Q点在尾喷口截面上。第一种情况下，Q点的辐射亮度为周围环境大气的光谱辐射亮度。第二种情况下，将尾喷管看作灰色辐射体，Q点的辐射亮度按照下式计算:

式中，εQ为尾喷管热空腔的发射率;TQ为热空腔温度。计算得到LQ后，将QP路径划分为n段，如图5所示。近似认为每一小段路径Rk，k+1中气体的温度和组分相同，对QP路径应用Malkmus统计窄谱带模型与C－G近似［10］，得到尾焰通过点P的红外辐射亮度大小。

图5 QP路径分段

4 探测器成像仿真

4．1 灰度值量化

分别计算出机体和尾焰的红外辐射亮度后，对所得的结果进行比较，得出辐射亮度的最大值Lmax与最小值Lmin，设任意一点i的辐射亮度为Li，i点的灰度可设置为:

Gmax与Gmin分别为灰度值上限与下限，取值范围在［0，255］。

4．2 噪声模型

噪声的灰度值变化近似符合高斯分布，噪声仿真方法如下:

假设成像系统的感光元件由m×n个像素组成，需要首先生成m×n个均值为μ、方差为σ的随机高斯数，作为噪声灰度值矩阵，然后将图像分为m×n个单元，在每个单元原有灰度值的基础上叠加噪声灰度值矩阵，即可得到叠加噪声的红外图像。若噪声类型是不随帧改变的空间噪声，在生成下一帧图像时，不重新生成噪声灰度值矩阵。若噪声类型是随机三维噪声，则在生成下一帧图像时，重新生成噪声灰度值矩阵。

4．3 三维图像投影

真实成像系统所成的图像是二维图像，因此需要对三维图像进行投影，本文采取透视投影法对三维图像进行投影，投影时定义图所示的平截头体ABCD－A'B'C'D'，将三维图像置于平截头体中，然后将平截头体投影到探测器所在位置。图中w/h为图像宽高比，θ1与θ2为探测器的视场角。

图6 透视投影示意图

5 仿真算例及分析

本文模拟了F16飞机在不同条件下的红外图像，飞行条件设置如下:飞行高度3 km，飞行速度1．2 Ma，太阳位于目标体轴系的y轴正方向，发动机尾喷口气体流速600 m/s，探测器与目标之间的距离为900 m，像素数为640×480，探测器视场角为3°×2．25°。图7(a)与图7(b)分别为探测器位于目标正上方与目标侧后方时所成的8～12μm波段的红外图像;图7(c)与图7(d)分别为探测器位于目标正上方与目标侧后方时所成的3～5μm波段的红外图像;图7(e)与图7(f)分别为探测器位于目标正上方，但不考虑太阳辐射时所成的3～5μm与8～12μm波段的红外图像。

保持探测器的视场角为3°×2．25°不变，改变探测器与目标之间的距离，图8(a)为探测器与目标距离1800 m时所成的8～12μm波段的红外图像;图8(b)为探测器与目标距离3200 m时所成的8～12μm波段的红外图像，

按照4．2节中的噪声模型生成640×480像素的随机噪声，图9(a)与图9(b)分别为叠加噪声前后飞机的侧向红外图像。

分别对比图7(a)与图7(c)、图7(b)与图7(d)可知，在8～12μm波段，目标轮廓较明显，机体辐射占主要部分，而在3～5μm波段，机体辐射几乎不可见，尾焰辐射占主要部分。

分别对比图7(a)与图7(b)、图7(c)与图7(d)可知，探测器与目标之间的距离不变时，探测器角度不同，目标红外图像所占像素点的数量及目标轮廓存在很大差别。

分别对比图7(c)与图7(e)、图7(a)与图7(f)可知，太阳辐射在3～5μm波段对于机体辐射的影响较大，而在8～12μm波段对于机体辐射的几乎没有影响。而对比图7(c)与图7(d)，图7(d)中的机体辐射亮度相对于图7(c)中较低，说明机体反射太阳辐射的能量与机体角度相关。

对比图8(a)与图8(b)可知，视场角不变的情况下，探测器距目标越远，目标图像所占的像素点越少，目标的轮廓特征及灰度变化越不明显。

对比图9(a)与图9(b)可知，叠加噪声后，目标与背景的对比变弱，不利于探测器对目标的识别。

图7 探测器与目标距离400 m时的红外图像

图8 改变探测器与飞机距离时的红外图像

图9 叠加噪声前后图像对比

6 结语

本文研究了飞机红外图像的仿真方法，然后以F16飞机为例，生成了飞机的红外图像，并对不同条件下的红外图像进行了对比分析。在计算蒙皮表面温度场时，采用求解热平衡方程法，并基于传热学原理提出了蒙皮气动加热的计算方法;计算排气系统红外辐射时，采用基于包络面的红外辐射计算方法，结合非均匀热气体的计算模型，有效减小了计算量。与传统方法相比，该方法能在准确性的基础上以较快的速度生成飞机的红外图像，对于分析飞机的红外图像及进行红外成像制导导弹的仿真特征具有较强的应用价值。

本文只研究了飞机目标本身的红外图像，并未对背景的红外图像以及红外辐射的大气传输过程进行研究，在下一步的研究中，应将飞机目标本身的红外图像、背景的红外图像与红外辐射的大气传输过程相结合，从而模拟出更加真实的红外场景。

［1］ GE Wei，CAO Dongjie，HAO Hongxu．Application of IR control and guidance technology in precise attack weapons［J］．Acta Armamentarll，2010，36(S2):35－38．(in Chinese)葛炜，曹东杰，郝宏旭．红外制导技术在精确打击武器中的应用［J］．兵工学报，2010，36(S2):35－38．

［2］ ZHANG Ke，LI Ning，LIU Fumei，et al．Method of aircraft skin infrared radiation image generating［J］．Laser＆Infrared，2011，41(3):272－277．(in Chinese)张可，黎宁，刘福美，等．一种飞行器蒙皮红外辐射图像生成方法［J］．激光与红外，2011，41(3):272－277．

［3］ WANG Fei，HE Jing，WANG Xin－sai，et al．Simulation model of IR imaging for the aeroplane［J］．Infrared and Laser Engineering，2007，36(3):352－356．(in Chinese)王飞，贺菁，王新赛，等．空中飞机红外成像仿真模型研究［J］．红外与激光工程，2007，36(3):352－356．

［4］ ZHANG Ke．Research on technology of creating aircraft’s skin infrared images［D］．Nanjing:Nanjing University of Aeronautics and Astronautics，2009．(in Chinese)张可．飞行器蒙皮红外辐射图像生成技术研究［D］．南京:南京航空航天大学，2010．

［5］ MEI Fei，JIANG Yong，CHEN Shiguo，et al．Infrared imaging prediction model for aero－engine exhaust plume［J］．Laser＆Infrared，2012，42(8):909－913．(in Chinese)梅飞，江勇，陈世国，等．一种新的三维实时红外尾焰仿真方法［J］．红外技术，2012，42(8):909－913．

［6］ YU Yang，TANG Xinyi，LIU Peng，et al．A new way of real－time 3D simulation of infrared plume［J］．Infrared Technology，2009，31(10):577－580．(in Chinese)于洋，汤心溢，刘鹏，等．一种航空发动机喷流红外成像仿 真 模 型［J］．激 光 与 红 外，2009，31(10):577－580．

［7］ 姜贵庆，刘连元．高速气流传热与烧蚀热防护［M］．北京:国防工业出版社，2003．

［8］ REN Qingmei，YANG Zhibin，CHENG Zhu，et al．Develop－ment of the platform for analysis coupling aero heating and structural temperature field［J］．Structure＆Environment Engineering，2009，36(5):33－38．(in Chinese)任青梅，杨志斌，成竹，等．气动加热与结构温度场耦合分析平台研发技术［J］．强度与环境，2009，36(5):33－38．

［9］ MAO Xia，HU Haiyong，HUANG Kang．Calculation method for airplane IR radiation and atmospheric transmit－ttance［J］．Journal of Beijing University of Aeronautics and Astronautics，2009，35(10):1228－1231．(in Chinese)毛峡，胡海勇，黄康．飞机红外辐射及大气透过率计算方法［J］．北京航空航天大学学报，2009，35(10):1228－1231．

［10］LI Jianxun，TONG Zhongxiang，WANG Chaozhe，et al．Calculation and simulation on infrared radiation of hot jet from engine［J］．Spectroscopy，2013，33(1):7－13．(in Chinese)李建勋，童中翔，王超哲，等．发动机热喷流红外辐射计算与仿真［J］．光谱学与光谱分析，2013，33(1):7－13．