考虑太阳辐射影响的飞机实时红外探测成像仿真

2019-12-13张昊春曲博岩

应用光学 2019年6期

张昊春，曲博岩，金亮

(哈尔滨工业大学能源科学与工程学院，黑龙江哈尔滨 150001)

引言

随着光电子技术的发展，红外成像、光电制导及跟踪、激光测距、激光雷达、声呐等技术开始广泛应用于陆海空多种作战平台，从而大大改变了现代战场的攻防形态。在很长一段时间里，针对战斗机等飞行目标的探测手段主要分为雷达探测和非雷达探测两种。而鉴于目前所开发的新一代飞行器多采用隐形材料，使得雷达的探测有效性明显降低，而红外探测技术以其不受无线电干扰，机动性较好，探测精确性高等优势，显著地弥补了传统雷达对飞行器目标探测能力的不足[1]，进而得到了发展。

飞行器在高速飞行时机身与周围空气摩擦将产生大量的热，加上机体表面接收到的来自太阳、地表和大气背景辐射的影响，因此机身对外会发出强烈的红外辐射。哈尔滨工业大学的夏新林等[2]通过研究飞机蒙皮温度场中的传热过程，建立起飞机内部热传导、飞机蒙皮表面气动对流及辐射的瞬态传热模型。北京航空航天大学的吕建伟等[3]通过反向的蒙特卡洛法(RMC)，引入表面多重遮挡的算法，计算具有复杂几何外形的F22机身红外辐射特性。国防科技大学的卢建等[4]通过3DMAX软件建立了飞行器目标的三维坐标模型以及传感器光学成像系统仿真模型。王科伟等[5]介绍了红外成像波门形心跟踪算法的原理，通过建立跟踪误差模型，讨论分析了红外成像波门形心跟踪算法的误差。成刚等[6]对多光谱成像探测技术进行了详细研究，针对迷彩伪装目标侦察设计了无人机载多光谱摄像机，重点分析了多光谱侦察探测、识别的计算方法。韩军等[7]提出一种提高现有光电观瞄系统成像空间分辨率的新方法。在不改变阵列探测器件像元尺寸和不移动探测器的前提下，利用双光楔的较大移动使像平面发生微小位移，实现对探测器各相邻像元和不感光间隔目标进行微位移采样，提高目标信息的采样率，通过超分辨重建技术来提高系统的分辨能力，达到提高光电观瞄系统空间分辨率的目的。Nicolas Douchin[8]论述了如何利用SE-Workbench-EO软件模拟飞机红外辐射特征以及获得实时渲染画面，并采用宽带、窄带模型及Curtis-Godson逼近法计算出飞机羽流区内的辐射亮度。Haynes[9]介绍了一种紫外至红外波段的场景生成软件CAMEO-SIM的各组件及其相应的功能，说明了该软件中加入红外痕迹及阴影、尾焰、传感器效应等因素的方法。Gretzmacher Floris M等[10]采用一种改进的欧几里得距离测量的视觉图像来进行伪装评估，在图像中形状相似的选择区域可以可视化。Bastiaansen M C等[11]采用丹麦国防研究机构(DDRE)开发的计算机伪装评估方法CAMEVA，输入由高分辨率目标和适当数量的背景组成的单个图像。Jacobs P.A.M等[12]为了确定伪装材料跟踪背景温度变化的能力，在各种气象条件下进行伪装和背景温度的测量，以确定伪装材料有效降低目标特征的情况。

截至目前，在红外成像领域，较多学者针对飞行器自身的红外辐射特性进行了成像仿真研究，而较少考虑了太阳辐射与天空背景对成像的影响，导致成像结果有一定误差。而当探测时刻、地点或季节改变时，太阳相对飞机的位置变化将对飞机红外成像产生较大的影响，综合考虑太阳辐射等环境因素，针对不同的飞行姿态及速度，不同的探测时刻、地点对飞机红外成像的影响进行研究。

1 飞机红外辐射特性分析

1.1 飞机模型

采用COMSOL软件建立J-16飞机简化后的三维几何模型，并对机体外表面划分非结构网格。机体表面网格沿飞机中轴线呈左右对称分布，并将机头、发动机两处温度变化梯度较大的位置上的网格进行加密处理，几何模型与划分结果分别如图1和图2所示。

图1 飞机几何模型示意图Fig.1 Diagram of aircraft geometric model

图2 网格划分结果示意图Fig.2 Diagram of grating division result

1.2 红外辐射特性

气动加热的形成条件为：飞机处于高速运动过程中时，机体周围的气流将有一部分动能不可逆地转化为热能，从而在蒙皮表面产生一层热层，与机身蒙皮进行换热[4]。目前对于飞机蒙皮温度分布的计算主要有数值仿真与经验公式两种方法。第一种是通过ANSYS等软件进行仿真模拟，进而得出蒙皮的温度场分布。在工程应用中，为了简化计算，常常采用经验公式((1)式)计算得到驻点温度Tr，以之代替蒙皮温度。

(1)

式中：Ta为处于飞机飞行高度上的大气温度，单位K，可以查表获得；k即气体绝热指数，文中取为1.4；γ即恢复系数，其中层流一般取0.85，湍流取为0.83；Ma即飞机飞行速度马赫数。由普朗克定律可知，蒙皮区域单个面元上的光谱辐射出射度等于：

(2)

式中：λ为波长，单位μm；T为绝对温度，单位K；c1、c2分别为第一、第二辐射常数。另外，为计算飞机蒙皮的定向红外辐射亮度，在t时刻将沿视线方向能够观察到的若干面元的辐射强度作积分处理，即得蒙皮区域的总辐射强度：

(3)

式中：ε即面元的发射率，文中均取为0.9；Ai为面元i的面积；ωi为面元i的外法线相对于视线的交角。因为目标与其所处背景的红外辐射强度分布不同，随着探测器成像单元接收到的辐射强度增大，对应的像素点灰度值增加，根据这一点，在捕捉到的红外图像中应用灰度等级表示目标区域的辐射能量分布。灰度等级映射关系由飞机蒙皮区域红外辐射亮度的计算结果建立，计算公式如下[15]：

GL=[r+(1-r)·(L-Lmin)/

(Lmax-Lmin)]×255

(4)

式中：GL表示某一像素点的灰度值，取值范围为0～255；r表示环境泛光的红外等效值，取值范围为0～1；Lmin和Lmax表示成像域辐射亮度的下限与上限，L即为该像素点处红外辐射亮度的初始值。最后根据像素点辐射温度和灰度值的转化关系得出飞机蒙皮红外仿真图像，由红外测温原理知，飞机蒙皮任一面元辐射温度Ts为[16]

(5)

式中：τa为探测器到达飞机蒙皮位置处的大气透过率；ε(χ)为λ1～λ2波段内飞机蒙皮面元法线方向的平均发射率；χ为入射角，即探测方向与面元的夹角；Tr为蒙皮温度，根据(1)式进行计算；n的取值与波段相关，在2 μm～5 μm波段，n=8.68，在8 μm～13 μm波段，n=4.09。

2 太阳辐射模型

2.1 太阳高度角及方位角

太阳高度角Hs及方位角As的计算公式如下：

sinHs=sinφsinδ+cosφcosδcost

(6)

cosAs=(sinHs·sinφ-sinδ)/(cosHs·cosφ)

(7)

式中：φ为飞机所处的纬度；δ为太阳赤纬；t为太阳时角。太阳赤纬的计算公式如下：

δ=0.006 918-0.399 912cos(b)+

0.070 257sin(b)-0.007 758cos(2b)+

0.000 907sin(2b)-0.002 697cos(3b)+0.001 48sin(3b)

(8)

式中：b=2π(N-1)/365，N从每年的1月1号算起，直到仿真时所取日期为止的天数(1月1日则N=1，1月2日N=2，以此类推)。时角的计算公式如下：

t=ts·15°

(9)

式中：ts为视太阳时，其值为平太阳时(即地方时)减去12。太阳时角在正午时为 0 (即位于天空的正中间)，上午为负，下午为正，日出时为-90°，日落时为+90°，平均每小时时角变换为15°。以中国东北某城市(E125°，N44°，海拔170 m)为例，零方位角为S方向，计算6:00～18:00期间内太阳高度角与方位角的变化曲线，结果如图3和图4。

图3 某市太阳高度角Fig.3 Solar height angle in a certain city

图4 某市太阳方位角Fig.4 Solar azimuth angle in a certain city

2.2 太阳辐射面积比例因子

由于所取飞机模型的结构比较复杂，各个面元之间彼此存在着遮挡的情况，故需定义一个与面元遮挡关系相关的参量，即太阳辐射的面积比例因子。假设面元a由于b的遮挡，只有总面积中的一部分能够接收到太阳的辐射能量，于是面积比例因子D定义为

(10)

式中：Ai即某面元的总面积大小；Ai,sun为面元能够接受到太阳辐射的面积大小。遮挡可分为完全遮挡和不完全遮挡两种情况，如果此面元完全没有被其他面元遮挡，则D=1；而当此面元完全被其他面元遮挡时，D=0，由此可知D的取值范围在0到1之间。

基于反向蒙特卡洛法确定每个面元的面积比例因子。首先，在一个面元中取总数为N的随机分布的点，然后根据这个面元所接收的太阳光的入射方向，令所取的每个点沿着其反方向射出一束光线，再通过空间几何关系推断出每束光线与其余的面元是否相交。只要有一束光线与其余面元相交，即说明该随机点被遮挡(所有的光束都相交即为完全遮挡)；如每一束都不相交，说明该点可以接受到全部的太阳辐射。对所有光线的判断循环完成后，统计面元被遮挡的点的数目，从而得到可被照射到的点的数目，最终求出每个面元的D值。根据以上方法，假定太阳入射方向向量为(-1,-1,-1)时，D值的计算结果如图5所示。

图5 D值分布情况示意图Fig.5 Distribution schematic diagram of solar radiation area scale factor

2.3 太阳辐射值

太阳常数定义为太阳在大气层外产生的总辐照度Eall=1 353 W/m2，太阳的平均半径约6.363 8×105km，日地平均距离约为1.5×108km[14]。当计算到达飞机所在位置的太阳辐照度时，需要考虑的影响因素诸如太阳的高度角、观察者的海平面高度、天空中云雾、尘埃等，情况比较复杂。在工程计算中，太阳辐射值Es(W/m2)一般使用如下经验公式来进行估算：

Es=[1-A(U*,β)]·(0.349Eall)sinβ+

(11)

表1 不同地表状态反射率Table 1 Reflectivity of different surface states

3 红外成像仿真

3.1 坐标系转换

在整个成像仿真过程中，涉及到各坐标系之间的转换，其中包括用于标定飞机三维坐标的大地坐标系，描述飞机尺寸和飞行姿态的飞机坐标系，表示探测器三维坐标的探测器坐标系，确定成像位置的成像面坐标系以及位于焦平面的像元坐标系。

为表明各个坐标系之间的相对位置关系，画出成像仿真坐标系统示意图，如图6所示。其中主要涉及的坐标系包括大地坐标系Xd-Yd-Zd，飞机模型坐标系X-Y-Z，成像面坐标系Xp-Yp和探测系统坐标系Xs-Ys-Zs。

图6 系统坐标示意图Fig.6 System coordinate diagram

进行三维空间的坐标旋转变换要求已知旋转角度和旋转轴，假设在右手系中，全部坐标点绕x轴、y轴、z轴旋转的角度分别为α、β和γ，旋转变换的变换矩阵R如下式[9]：

(12)

假设坐标系中的全部坐标点相对于x轴、y轴、z轴的平移量分别为tx、ty、tz，平移变换的变换矩阵T如下式：

(13)

由(10)式和(11)式，得出飞机模型坐标系到大地坐标系的转换关系式：

(14)

探测器坐标系到成像坐标系的转换关系式如下：

(15)

式中f为焦距。接下来在成像坐标系中确定像元的位置坐标，假设像元坐标系下某个像元在x和y方向上的宽度分别为Δm和Δn，于是根据下式可将成像坐标系下的投影点映射到像元坐标系。

(16)

式中[x]代表对x向下取整。飞机机身蒙皮离散为很多个三角面元，如果三角面元至少有2个点位于某成像单元的有效成像范围内，则该面元在这个成像单元内成像[15]。将机身蒙皮上各个面元逐一投影到像元坐标系上，即可得到飞机蒙皮面元在该坐标系下的3个顶点坐标，由这3个顶点确定的三角面之和即为目标的成像区域，而探测器阵列上其余位置则为背景区域。

3.2 仿真结果及分析

假定在大地坐标系下，飞机的三维坐标为(0,0,1 000)，探测器坐标为(0,0,0)，飞行速度1马赫，大气模式为中纬度夏季模式，探测地点为E120°，N30°，地表状态为干燥裸地，探测时间为北京时间2019年6月1日星期六，8:00AM，探测波段取8 μm～12 μm，仿真结果如图7所示，其中图7(a)为灰度等级由(5)式转化而得的红外辐射温度(单位K)，图7(b)为处理后相应的RGB图像。

图7 飞行速度1 Ma时仿真结果Fig.7 Simulation results at flight speed of 1 Ma

根据计算结果，结合图像进行分析知，此时飞机处于高速运行状态，气动加热效果显著，机身的发射辐射中来自太阳辐射的份额占总辐射强度的比例较小，所以蒙皮的辐射亮度分布主要由气动加热层定向辐射强度决定；另一方面由于太阳光倾斜照射，飞机机头被直接照射的区域亮度明显高于其他区域。经计算，图7中框区域内的辐射亮度比其他区域辐射亮度高18.2%。

飞机在飞行过程中速度的变化将会影响气动加热效果，进而影响仿真结果，另外，飞行高度的变化将会使得探测器光学路径长度改变，也会对结果产生影响。现保持上述计算条件不变的情况下，改变飞机的运行速度和高度，分别研究二者带来的影响，其中飞行速度分别取为0.5 Ma、1 Ma(图7)、3 Ma和6 Ma，将得到的仿真结果和图7统一转化为辐射温度，结果如图8至图10所示。飞行高度分别取1 000 m(图7)、4 000 m和8 000 m，最终得到的结果如图11所示。

图8 飞行速度0.5 Ma时的仿真结果Fig.8 Simulation results at flight speed of 0.5 Ma

图9 飞行速度3 Ma时的仿真结果Fig.9 Simulation results at flight speed of 3 Ma

图10 飞行速度6 Ma时的仿真结果Fig.10 Simulation results at flight speed of 6 Ma

图11 改变飞行高度时的仿真结果Fig.11 Simulation results when flying altitude is changed

从图中可以看出，随着飞行速度的提高，飞机蒙皮辐射亮度逐渐增大，这一现象在翼面上更加明显(图8至图10)，这是由于当速度提高时，飞机与表面气流摩擦生成的热量更高，换热效果更加明显，从而使蒙皮区域温度和辐射亮度提高。图8至图10中框区域内的辐射亮度比其他区域辐射亮度分别高21.4%，17.9%和16.1%。另一方面，随着飞行高度的增加，目标在像平面中所占面积逐渐减小，当飞机从1 000 m高度上升到4 000 m时，蒙皮区域红外辐射亮度整体上变化很小，分析认为，这是由两方面因素导致的，首先随着飞行高度的增加，目标与探测器之间大气散射路径增大，大气介质中辐射能量的吸收和散射衰减增强，探测器接收到的辐射能量减小，所以仿真得到的目标红外辐射亮度减小；同时，随着飞行高度增加，对应于单个成像单元的飞机表面划分的面元数量增加，从而使像平面接收的每个成像单元红外辐射能量增加。所以，由于投影面元叠加和大气衰减的双重作用，从1 000 m到4 000 m高空的红外仿真图像整体变化很小。而当飞机从4 000 m爬升至8 000 m时，蒙皮区域整体的辐射亮度明显下降，此时大气衰减成为主导因素。

图12 飞机偏航，俯仰及旋转示意图Fig.12 Diagram of aircraft yaw, pitch and rotation

接下来研究探测时间以及飞行姿态对成像效果的影响。假定飞行速度为0.5马赫，飞行高度为1 000 m，大气模式为中纬度夏季模式，探测地点为E120°，N30°，地表状态为干燥裸地。以上条件均控制不变，现设定日期为6月1日，探测波段仍取8 μm～12 μm，探测时间及飞行姿态控制参数见表2，得到的仿真结果如图13至图16所示。规定飞机左偏航时偏航角为正，抬头时俯仰角为正，右滚转时旋转角为正，如下图所示。

图13 12:00AM的仿真结果Fig.13 Simulation result at 12:00AM

图14 14:00PM的仿真结果Fig.14 Simulation result at 14:00PM

图15 16:00PM的仿真结果Fig.15 Simulation result at 16:00PM

图16 18:00PM的仿真结果Fig.16 Simulation result at 18:00PM

表2 仿真计算参数设置Table 2

由计算结果分析知，因飞行速度较低，气动加热效果不明显，机身的发射辐射中来自太阳辐射的份额占总辐射强度的比例较大，所以蒙皮的辐射亮度分布主要由接收的太阳辐射条件(高度角，面积比例因子等)决定。当探测时间为中午12点时，机身被太阳照射的面积上辐射亮度呈对称分布，同时因机头处D值较小，所以辐射亮度较小，经计算得图8中框区域的辐射亮度比其他区域辐射亮度低7.9%。