王彩云,赵焕玥,李晓飞,刘纯胜,王立国,王佳宁
(1.南京航空航天大学航天学院,江苏 南京 210001;2.北京电子工程总体研究所,北京 100854)
由于航天技术的发展,各军事强国发射进入太空的空间目标数量剧增。空间目标的体积愈加小型化,结构更加复杂化,形状越来越多样化,空间目标探测被赋予新的挑战。红外探测作为一种重要的光学探测手段,因其对目标背景温差的高敏感性,探测的高隐蔽性以及可以进行多光谱探测的优点成为了理想的空间探测方式[1]。与外场试验相比较,红外成像仿真可以获得特定位置、时间和天气下的红外图像,不仅能减少试验成本,还能提高试验效率,具有广泛的应用前景与研究价值。
国内在空间目标的红外成像起步较晚,但是也取得了很大进步。例如文献[2]~[4]中阐述了空间目标的红外辐射特性分析建模方法,对其红外辐射强度进行了数值求解,并将计算结果转化成灰度值,最后采用计算机图形学算法生成红外图像,为空间目标红外探测提供了有力的工具。但是在红外成像方面,已有的研究大都基于理论数值计算阶段,尚未开发一套从目标表面温度分布求解到目标红外辐射特性分析求解,再转入目标红外成像的系统仿真软件,而国外成熟的红外成像商业软件由于版权限制,价格昂贵,也不是最佳选择[5]。
所以基于经济因素和可行性的综合考虑,本文提出了一种基于三角剖分的空间目标红外成像仿真方法,建立了一套完整的空间目标红外成像仿真方法。此方法首先根据空间目标热平衡方程[6],利用节点网络法求解出目标的各节点温度,然后利用STK得出探测器、目标、地球与太阳的位置关系,随后由目标节点温度求得目标红外辐射强度,最后利用改进Delaunay三角剖分的方法重建目标得到红外图像。以某卫星[7]为例进行了成像仿真实验,成像结果验证了本文方法的有效性,为进一步的空间目标的探测与识别提供了基础。
空间目标的红外辐射特性主要取决于目标的表面温度分布。要对空间目标进行红外成像,首先应该建立目标表面温度求解模型,然后根据温度分布求解其红外辐射特性。
为了简化计算,忽略空间目标自身工作的发热,目标表面温度的分布主要由太阳直接辐射[8]、地球反射太阳辐射[9]、地球直接辐射和目标表面向宇宙空间的辐射四个因素决定。
(1)太阳直接辐射
太阳的半径为6.960×105km,可以被看作是温度为5900 K的黑体。把太阳与地球的距离当作其平均距离(1.4961×108km),地球大气层外法线方向的太阳总辐照度Es为:
(1)
式中,Eλ是太阳的光谱辐射照度。
另外,因为地球半径远远小于太阳与地球距离,可以将太阳光视为辐射度均匀的平行光,则在目标表面的第i个单元吸收的太阳直接辐射qi1为:
qi1=εsEsSiFs,i
(2)
式中,εs表示材料对太阳辐射的吸收率;Si为表面面积;Fs,i是外表面对于太阳直接辐射的角系数。
(2)地球反射太阳辐射
除了太阳直接辐射,目标表面还会接收地球反射太阳辐射,被视为地球表面的漫反射,则该表面单元所吸收的地球反射太阳辐射qi2为:
qi2=εrEsρFse,i
(3)
式中,εr为吸收率;ρ为地球反射率;Fse,i为角系数。
(3)地球直接辐射
地球可以被视为一个温度280 K的黑体,其空间某处的总辐射照度Ee为:
(4)
式中,Re为 6450 km,是地球及其大气系统总半径;h为目标和大气层的距离,则该表面单元吸收的地球直接辐射qi3为:
qi3=εeEeSiFe,i
(5)
式中,εe为吸收率;Fe,i为角系数。
(4)目标向宇宙空间的辐射
若目标的温度不低于热力学零度,则目标表面一定会向外辐射能量。第i个单元表面向宇宙空间的辐射能qi4为:
qi4=εσTi4(t)Si
(6)
式中,ε为发射率;σ为斯蒂芬-玻尔兹曼常数;Ti(t) 为该部分在时间t的温度。
(5)内能的增加
空间目标的内能会根据温度的变化而发生变化,在任一瞬间,第i个单元内能的增加量为:
(7)
式中,c表示表面的比热容;mi为质量。
根据上文的描述,热平衡方程如下:
qi1+qi2+qi3-qi4=qi5
(8)
(9)
(2)将高次项变为一次项:
Ti4(t)=4Ti3(t-Δt)Ti(t)-3Ti4(t-Δt)
(10)
如果给定太阳、目标、地球的位置关系,就能计算出式(8)的角系数Fs,i、Fse,i和Fe,i的值,然后由式(8)~(10),可知对目标任意表面单元,只要知道其在t-Δt时刻的温度值Ti(t-Δt),其在t时刻的温度值Ti(t)就可以推算得到。所以,若知道了目标表面的初始温度分布以及目标与太阳、地球的方位,目标的表面温度分布就能求出。
本文采用节点网络法对空间目标的表面温度分布进行计算。节点网络法是把空间表面划分为许多热力学特性基本一致的单元,把单元作为节点,对节点和背景的热交换进行分析,建立各节点的瞬时热平衡方程。在给出的初始条件下,递推计算这些方程组,便可以计算出目标表面节点的温度值和其变化规律[9]。
空间目标对空间某处的红外辐射包含两部分:目标自身红外辐射和目标反射背景辐射。
(1)目标自身红外辐射
在空间目标表面温度分布被计算出以后,可以采用普朗克公式,在其红外波段范围的积分得到:
(11)
式中,λ1,λ2为红外波段范围的上下限;C1为第一辐射常数,其值为3.742×108W·μm4/m2;T为单元表面温度;C2为第二辐射常数,其值为1.439×104μm·K。
目标自身红外辐射强度采用漫辐射源的辐射特性得出:
(12)
式中,ΔA为目标表面单元面积;θi为目标表面单元法线与辐射方向的夹角。
(2)目标反射辐射
目标对背景的反射辐射亮度为:
(13)
式中,ρe为等效反射系数;E为外界对目标辐射量的总和。目标反射强度可采用漫反射源的辐射特性得到:
Ir=∑LrΔAcosθi
(14)
式中,ΔA为目标表面单元面积;θi为目标表面单元面法线与反射方向的夹角。
因此,空间目标的总辐射强度为:
I=Ifa+Ir
(15)
前一节计算出的辐射强度并不是红外成像仿真需要的最终结果,需要把经过成像系统效应作用后的辐射强度转化成灰度值。
现在通常使用的是均匀量化,辐射强度和灰度值之间是一个线性关系。一幅图像中像面微元对应的最高和最低辐射强度分别为Imin和Imax,设图像灰度等级为0~255级,则对应于辐射强度为I(j)的像面微元的灰度值为:
(16)
前面已经计算出红外图像中节点的灰度值,下一步是将这些离散的节点导出,重建为空间目标,并且选取探测器所在角度进行成像,这一步要采用Delaunay三角剖分[11]方法。
Delaunay三角剖分算法常用的分为两种:基于投影的三角剖分算法[12],此算法容易出现拼接错误并且不利于复杂曲面的重构;另一种为基于Delaunay的生长算法,即在三维离散点中构造一个种子三角形,随后根据Delaunay准则寻找最优点扩展形成一个新的三角形,以此三角形为基础再次向外搜索最优点,直到所有的散射点都加入到三角网中构成三角形[13]。此算法的优点是思想简单,能对复杂三维曲面进行重构,但是此方法在搜索满足条件的第三点时计算复杂,耗费大量时间,且易出现曲面重叠、空洞等现象[14]。所以本文将采用三维栅格法[15]和自适应空间外接球方法[16]对其进行改进,提出一种改进Delaunay三角剖分算法。
三维栅格法首先读取三维点云坐标,求得点云在三个坐标轴X、Y和Z轴上的最大值和最小值,确定能将所有数据包含在内的最小长方体,然后选区合适的边长值,将长方体分割成若干个相同大小的立方体。最后依据各点的三位坐标值将点存储在相对应的小立方体编号数组中,完成了点云数据的空间划分。该方法可以根据点云的疏密程度自由调整子区域的大小,在领域点搜索的过程中,可以减小搜索范围,提高效率。
自适应空间外接球方法的思想是找到点集中距离生长边的中心点最近的点组成新三角形,确定此三角形外接球并找到球内部的点云数据,最后通过迭代找到最优扩展点。
综上,本文的改进Delaunay三角剖分算法的步骤为:
(1)利用上述三维栅格法完成某个点的k邻域寻找;
(2)构造种子三角形,利用种子三角形来生长三角网;
(3)采用自适应空间外接球方法寻找三角形生长的最优点;
(4)对未使用的点集进行迭代生长。
本文采用这种方法先对空间目标进行重建,将红外辐射强度转化的灰度值分布于空间目标的表面,随后采用STK计算出探测器方位,根据其方位矢量得出二维图像。
本文以某卫星为例,进行红外辐射成像仿真实验,首先进行目标三维建模并计算卫星表面温度分布,随后利用第二节中的公式求出目标红外辐射特性大小,最后采用改进的Delaunay三角剖分成像方法得到目标的红外图像。
首先表1中列出了探测目标的物性参数。表2列出了目标轨道与探测轨道参数。
表1 目标物性参数表
表2 目标与探测轨道参数
目标和探测器轨道均为太阳同步轨道,其中ad为轨道半长轴;ed为轨道偏心率;id为轨道倾角;Ωd为升交点赤经;ωd为近地点辐角;τd为升交点地方时。
为了验证本文计算空间目标红外辐射强度的模型与成像方法是否合理有效。图1、图2首先比较了在3~5μm与8~12μm两个波段下目标本体面元与太阳帆板正面面元总红外辐射强度变化,从2009 年 6 月 5 日 15∶30开始,横坐标为距开始时间的秒数,纵坐标为红外辐射强度大小。由于红外辐射强度随目标周期性运动也呈现周期性变化,图1展现了红外辐射强度的两个周期变化曲线。
图1 两种波段下目标本体面元总辐射强度随时间变化曲线
由图1(a)与图1(b)对比可知,8~12 μm波段探测的目标本体面元辐射强度比3~5 μm更强,采用8~12 μm波段更易于探测目标。在某一时间段,目标本体面元辐射强度达到最小保持不变,这时面元处于阴影区。
图2 两种波段下目标太阳帆板面元总辐射强度随时间变化曲线
对比图1、图2,可以看出在光照区,目标太阳帆板正面面元总辐射强度在8~12 μm波段大于目标本体面元的总辐射强度,在3~5 μm波段并不明显,这是因为太阳辐射能量主要集中在中短波,目标本体反射率较大,所以在3~5 μm波段反射了更多的太阳能量。
图3为目标红外辐射成像仿真结果,仿真从5400 s(17∶00)开始,直到114000 s(18∶40)结束,每隔1200 s成一次像。
图3 不同时间目标双波段成像结果
由图3可以看出,同一时间8~12 μm波段辐射强度大于3~5 μm波段。目标的成像角度随时间不断变化,在光照区,目标的太阳帆板正面的亮度大于目标本体,满足理论实际。综上所述,本文提出的基于改进Delaunay三角剖分的空间目标红外辐射成像方法有效。
为了比较传统的Delaunay三角剖分算法与本文改进后的Delaunay三角剖分算法的重建效率,进行了两种算法对同一目标的重建时间对比实验,实验平台均为Matlab R2017a,且均在安装 64 位Windows 10 操作系统的同一台计算机上完成,表3是运行时间比较结果。
表3 算法运行时间比较
由表3可以看出,本文改进后的算法重建时间比传统Delaunay三角剖分算法重建时间短,效率更高。
本文提出了一种基于改进Delaunay三角剖分的空间目标红外辐射成像方法,建立了一套完整的空间目标红外辐射成像方法。此方法首先根据热平衡方程,利用节点网络法求解出空间目标的各节点温度,然后利用STK得出探测器、目标、地球与太阳的位置关系,随后由目标节点温度求得目标红外辐射强度,最后利用改进的Delaunay三角剖分的方法重建目标得到红外图像。本文以某卫星为例进行了成像仿真实验,并对实验结果进行了分析,验证了本文方法的有效性,此方法对空间目标的红外探测以及识别具有重要意义。