HTML
-
中波红外4.3 μm附近波段由于大气中CO2的强烈吸收作用，地表辐射几乎难以到达星载传感器，同时火箭尾喷焰燃烧产物中的高温气体和粒子发出强烈的辐射，使其成为SBIRS-GEO等天基红外系统对发出高温喷焰的飞行器进行检测、跟踪和识别的重要波段之一。美国“空间中段监视”试验卫星(MSX)上搭载的空间红外成像望远镜SPIRIT III，其中的B1(4.22~4.36 μm)和B2 (4.24~4.46 μm)波段可用于跟踪空间飞行器等目标[1]。国内学者对具体的波段选择方法进行了研究，叶庆[2]通过分析目标背景对比度，提出了4.18~4.50 μm波段作为最佳探测波段的结论。
卫星观测火箭发射的场景中主要包括火箭尾喷焰、地表、大气和云。火箭发动机喷口喷出的高温高压气体、固体粒子燃烧产物，使得尾喷焰成为目标红外辐射的主要来源，其红外辐射图像仿真需经过流场组分特性计算、辐射传输计算和大气传输计算三个过程。刘尊洋[3]利用FLUENT软件计算不同飞行高度和速度的尾焰流场，利用HITEMP数据库和有限体积法（FVM）计算了尾焰在2.0~5.0 μm的辐射特性，并分析了不同飞行参数的影响。包醒东等[4]研究了喷焰中Al2O3颗粒与燃气的相互作用程度在不同高度下的变化规律，得出了不同飞行高度对气固两相喷焰的流动及红外辐射特性的影响。
在红外场景仿真技术方面，国外开发了SE Workbench和Vega Prime在内的一系列成熟的商业仿真软件，可涵盖材料数据库、场景几何建模、红外辐射建模、大气传输以及传感器建模的全流程，覆盖0.3~25 μm的红外波段，通用性强[5]。卫星观测火箭飞行时，视场范围较大，地物背景复杂多变，利用软件进行几何建模并嵌入其仿真流程，难度高且效率低。因此，文中选择利用大规模的卫星遥感数据集，通过波段迁移方法实现吸收波段地物背景的仿真。搭载在Terra和Aqua卫星上的中分辨率成像光谱仪(MODIS)，覆盖了0.4~14.4 μm的36个离散波段，其中部分波段与吸收波段相邻。MODIS数据因为多光谱、大宽幅、高重访、免费的特点被广泛应用[6]。Yan Chen等[7]利用MODIS对同一区域白天和夜晚的大气顶层辐亮度数据，通过联立方程反演了3.7、8.5、10.8、12.0 μm波段的地表发射率。Yao Liu等[8]利用光谱相关性，通过MODIS中的相邻波段数据仿真了SPIRIT III波段的背景图像。李文杰[9]使用单通道法对FY-3/MERSI数据集进行地表温度反演，计算了分辨率为5 km的地球背景红外辐射图。目前，关于尾喷焰和地物的仿真研究比较孤立，在目标、背景和传感器的关系上缺乏系统性，没有形成完整的动态探测场景。
文中建立了4.3 μm波段卫星观测火箭的场景生成模型。其中，背景辐射特性仿真方面建立了从吸收波段高分辨率地表发射率图像生成至大气传输计算的仿真流程；目标红外辐射特性仿真方面建立了从火箭轨迹生成、尾喷焰辐射传输计算到大气传输计算的仿真流程；通过构建完整的全场景内目标、背景、传感器的几何关系进行投影成像，形成了卫星观测助推段火箭飞行的动态场景，并对典型观测条件下的场景仿真结果进行了分析。
-
星载传感器入瞳处接收的辐亮度，在中红外波段主要包括地表辐射和路径辐射，地表某一点对于入瞳光谱辐亮度的贡献可以分解为[10]：
式中：
$L_{\lambda ,T}^{}$ 为总辐亮度；$L_{\lambda ,T}^E$ 和$L_\lambda ^P$ 分别为地表热辐射和路径辐射；${\tau _\lambda }$ 为光谱透过率；${\varepsilon _\lambda }$ 为光谱发射率；h为普朗克常数；c为光速；λ为波长(μm)；k为玻耳兹曼常数；T为热力学温度。太阳辐射集中在3.5 μm以下，波长越长，对总辐亮度的贡献越小。 文中的数据源是MODIS第22波段(Band22，3.929~3.989 μm)、第23波段(Band23，4.020~4.080 μm)以及需要仿真的4.18~4.5 μm波段(下称BandS)，光谱响应函数如图1所示。BandS的光谱响应函数使用高斯函数进行模拟，由于存在带外响应，实际计算时BandS展宽至4.1~4.6 μm。
Figure 1. Spectral response functions of Band22, Band23 and BandS
由公式(1)可知，背景的入瞳辐亮度取决于地表温度、发射率、大气透过率及路径辐射，其中发射率是物质的本质属性，大气透过率和路径辐射与大气模式以及地表点至传感器的几何路径等参数有关。文中使用的MODIS数据产品如表1所示。MOD11B1中Band22/23的地表发射率是通过万正明[11]提出的算法反演的，由于同时需要白天和夜晚的数据，需要从不同传感器匹配同一地表点，反演的空间分辨率降低至6 km左右。在星载红外相机的分辨率高于此分辨率的情况下，为了保证投影成像的精度，地表温度和发射率的空间分辨率应当进行相应的提高，甚至要达到百米量级。
Product name Data content applied to this paper Spatial resolution/km Projection method MOD021KM Calibrated and geolocated radiances in Band22 and Band23 1 Geographic Lat/Lon MOD11_L2 Surface temperature 1 Geographic Lat/Lon MOD03 Geolocation fields for 1 km IFOV 1 Geographic Lat/Lon MOD11B1 Surface emissivity in Band22 and Band23 6 Sinusoidal Table 1. Introduction to MODIS data
通过分析，在同一个场景中其他因素变化不大的情况下，地表点的辐亮度正向计算过程可用公式(2)简要描述：
式中：H为地表点的海拔高度；
${\theta _1}$ 为太阳到地表点的天顶角；${\theta _2}$ 为地表点到传感器的天顶角；${f_\lambda }$ 为从上述5个要素至辐亮度的映射，其中包含普朗克函数，是一个复杂的非线性函数。而对于发射率反演问题，则可表示为：式中：函数
${g_\lambda }$ 表征了从影响地表发射率反演结果的各要素到地表发射率之间的映射关系。基于上述分析，提出利用拟合能力极强的人工神经网络表征此映射关系的方法，将地表发射率分辨率从km级提升至100 m左右。基本原理是将视场内已知地表温度、海拔高度、太阳天顶角、传感器天顶角、表观辐亮度以及地表发射率的点构成训练数据集，训练神经网络的参数得到近似的映射${g_\lambda }$ 。当已知输入场景中邻近点的5个参数时，前向计算输出其地表发射率。网络结构如图2所示。对于MOD11B1中的地表发射率矩阵
${\varepsilon _{M \times N}}$ ，通过MOD11_L2、MOD03、MOD021 km数据向其经纬度网格插值，分别得到对应的地表温度矩阵${T_{M \times N}}$ 、海拔高度矩阵$ {H_{M \times N}} $ 、太阳天顶角矩阵${\theta _{1,M \times N}}$ 、传感器天顶角矩阵${\theta _{2,M \times N}}$ 以及表观辐亮度矩阵${L_{M \times N}}$ 。神经网络的输入层节点数为5 （$ T,H,{\theta _1},{\theta _2},L $ ）；两个隐层节点数均为10，激活函数为sigmoid；输出层节点为1（$\varepsilon $ ），激活函数为linear。
Figure 2. Schematic of neural networks for retrieving surface emissivity
MODIS数据中存在的无效值造成
${T_{M \times N}}$ 和${L_{M \times N}}$ 中也存在部分无效值，所以实际输入数据尺寸为$[ {M \times }$ ${ {N -N_{invalid}},5}] $ ，输出数据尺寸为$\left[ {M \times N - {N_{invalid}},1} \right]$ ，即样本数量为$M \times N - {N_{invalid}}$ ，其中${N_{invalid}}$ 为无效点数量。为验证上述反演算法的有效性，将已知发射率的点构成的数据集按70%、15%、15%的比例随机划分为训练集、验证集和测试集。均方误差统计结果如表2所示。
由表2可见，由于邻近点之间相对于传感器的观测几何关系相近，所提出的方法可以很好地表征潜在的大气传输物理模型，快速准确地生成高分辨率的地表发射率图像。图3为生成的100 m左右分辨率Band22/23地表发射率图像。
Number of samples MSE of Band22 MSE of Band23 Training set 4 074 3.06×10−4 2.95×10−4 Validation set 873 3.34×10−4 3.20×10−4 Test set 873 2.94×10−4 3.06×10−4 Table 2. MSE of training set, validation set and test set
Figure 3. Surface emissivity images with a resolution of 100 meters
-
20世纪80年代以来，美国喷气推进实验室JPL以及NASA等多家研究机构相继建立了地物波谱数据库。文中从MODIS UCSB发射率光谱库[12]中选取了4大类地物，包括土壤、植被、水体和人造材料共158条光谱，图4为其中41种地物在3.8~4.6 μm波段的光谱发射率。
Figure 4. Spectral emissivity in the 3.8-4.6 μm band of 4 categories including soil (a), vegetation (b), water (c) and manmade materials (d)
从图4可以得出以下结论：(1) 在3.8~4.6 μm，绝大部分土壤、植被和水体的发射率都很高，在0.7以上；(2) 部分人造材料发射率相对较低；(3) 大部分材料的发射率曲线平稳，可用直线拟合。基于此，采用参考文献[8]提出的迁移方法，通过Band22/23发射率拟合BandS发射率。参考文献[8]将所求波段的平均发射率作为其中所有微波段(microband)的发射率，文中为了提高精度并且便于匹配后续的光谱辐亮度计算，直接对BandS中每个10 nm带宽的微波段发射率进行拟合。
多元线性回归的公式为：
式中：X为已知的Band22/23发射率构成的矩阵，最后一列对应常数项；y为BandS中所有微波段发射率构成的矩阵；b为系数矩阵；
${N_{spec}}$ 是光谱数量；${N_{modisband}}$ 为2；微波段数量${N_{microband}}$ 为501。线性回归获得系数矩阵b后，已知任意地表点的Band22和Band23发射率，可计算BandS中每个微波段的发射率：使用统计量中的判定系数R2检验回归模型的有效性，R2越接近1表示线性拟合效果越好。如图5所示，R2基本保持在0.9以上。距离Band22/23越近的微波段，通过谱段关联法拟合得到的发射率置信度越高。
Figure 5. Coefficient of determination R2 of spectrum correlation algorithm
-
使用MODTRAN对高分辨率地表网格进行大气传输计算，如果每个点都计算一次会耗费大量时间，采用查询表（LUT）方法，将主要影响因素分别等间隔地设置参数表，对参数表的每一项进行计算，最后每个地表点根据插值确定在LUT中的位置，得到大气透过率和路径辐射。文中选取了海拔高度和观测天顶角两个参数建立LUT，步长分别为0.5 km和1°。此外，将最主要的吸收气体CO2的体积混合比从330 ppmV更新为380 ppmV。
-
火箭飞行轨迹是场景的关键，对尾喷焰辐射特性的动态变化有重要的影响。采用控制力加速度实现转弯的方法[13]，飞行加速度方面考虑气动阻力、重力、控制力和推力加速度。使用椭球地球模型计算重力，根据所在海拔的空气密度计算气动阻力大小，根据开始施加控制力的时间和发动机关机时间计算控制力加速度。
推力是发动机燃烧室燃料燃烧产物高速向后喷出而作用在火箭上的力。假设攻角为零，推力方向同速度方向，推力加速度大小为[12]：
式中：
${F_T}(t)$ 和$M(t)$ 分别为t时刻的推力大小和火箭剩余质量；${\rm d}M_f^{}$ 和$I_{sp}^{}$ 分别为推进剂秒消耗量和比冲；${M_{pl}}$ 和$M_b^{}$ 分别为有效载荷质量和助推器总质量；g为重力加速度。设置仿真时间间隔
${T_h}$ ，采用变步长的四阶龙格库塔法进行弹道外推，为火箭沿Y轴即竖直向上的速度赋微小的初始值：式中：u为火箭的状态向量。
-
通过8组分11化学反应模型计算尾喷焰流场的组分和温度、压强等的分布，采用视在光线法(LOS)计算尾焰红外辐射传输[13]，如图6所示。首先根据轨迹点的位置、速度方向、传感器位置等，确定视线与火箭轴线的夹角θ，建立垂直于视线的假想像平面，并对每个像元生成平行路径；然后对路径进行分层，每层的参数通过原始流场分布插值计算；最后沿路径逐层计算辐亮度与透过率[14]，最终到达假想像平面以进行后续投影成像。
Figure 6. Schematic of the LOS method
-
文中整合了第1节和第2节中提出的算法，开发了场景生成软件，其总体流程如图7所示。软件采用模块化设计，主体以及轨迹生成和背景仿真的各个模块基于Matlab平台编写，尾喷焰流场和辐射计算模块基于C++语言实现，大气传输通过调用MODTRAN程序进行计算。
Figure 7. Schematic of infrared dynamic scene generation of rocket plume observed by satellite
背景部分，首先通过MODIS原始数据生成较高分辨率的Band22/23地表发射率，然后拟合吸收波段BandS光谱发射率，最后结合观测几何关系利用MODTRAN计算所有地表点的辐亮度；目标部分，通过发射参数和火箭参数生成飞行轨迹，在轨迹上的各个时刻对尾喷焰图像进行仿真；最终，地表点和尾喷焰假想像平面中的点均通过投影成像的方式到达星载传感器入瞳处。
将地心坐标系作为世界坐标系，投影成像的过程即求解相机内参和外参矩阵：
式中：XO、YO和ZO分别为地表点或尾焰的地心系坐标；R、T分别为从地心坐标系到相机坐标系的旋转和偏移矩阵，组成外参矩阵；焦距fx、fy分别为从相机系到图像系的变换；u0、v0分别为从图像系到像素系的偏移量，两者共同组成内参矩阵；u、v为像素系的坐标；ZC为相机系Z轴分量。
星载传感器对尾喷焰和地面背景的观测几何关系和投影成像模型如图8所示。
Figure 8. Schematic of projection imaging of rocket plume and background from satellite
通过公式(8)计算视场内所有地表点和尾喷焰网格点在像素坐标系下的坐标，也即确定了它们属于图像中哪一个像元，其表观辐亮度计算公式为：
式中：
${L_{(m,n)}}$ 和${S_{(m,n)}}$ 为图像上(m, n)处像元的表观辐亮度和投影到地表所占的面积；$L_{(x,y)}^{BG}$ 和$S_{(x,y)}^{BG}$ 为属于像元(m, n)的地表点表观辐亮度和所占面积；$L_{(x,y)}^{Plume}$ 和$S_{(x,y)}^{Plume}$ 为LOS法输出的假想像平面中属于像元(m, n)的网格的表观辐亮度和网格面积。 -
此节对尾喷焰和背景的仿真结果可信度进行验证，通过对一个典型场景进行仿真，分析辐亮度图像序列中尾喷焰相对背景的变化规律，设置两个火箭发射点以及多条不同射向的飞行轨迹，分析不同场景和不同轨迹下尾喷焰总辐射强度的差异性。
-
为了验证所设计软件中火箭尾喷焰辐射的仿真结果，选择与参考文献[14]中相同的参数进行仿真，计算了1.5~5.5 μm波段的光谱辐射强度，与实测结果的比对如图9所示。可知光谱辐射强度曲线与实测结果规律一致，吻合度较高，4.3 μm附近峰值非常接近，验证了尾喷焰辐射仿真结果的准确性。
Figure 9. Comparison of the plume spectral radiation intensity between calculated and measured
对发射率计算方法进行验证，从JHU波谱库选取了包括土壤、植被、水体和人造材料在内的共470条发射率光谱作为测试样本。通过谱段关联方法计算BandS的光谱发射率，与波谱库中实测数据进行比较，分别统计了每类地物在不同波长下的平均相对误差。如图10所示，相对误差均保持在2%以内，表明通过Band22/23平均发射率拟合计算邻近的BandS光谱发射率的可行性。
Figure 10. Average relative error between the results of spectrum correlation algorithm and JHU spectrum library
最后，分析不同大气模式和观测条件下软件对于地表点入瞳辐亮度的计算误差。分别设定中纬度夏季和中纬度冬季两种大气模式，以及观测天顶角0°、30°和60°，地表温度均为300 K，生成上述470种不同地物类型的地表点到达卫星红外传感器入瞳处的辐亮度，并与直接使用JHU波谱库中BandS实测发射率计算得到的辐亮度结果进行比较，相对误差如图11所示。
Figure 11. Relative error of radiance under different atmospheric models and viewing zenith angles
从图11中可得到以下结论：(1) 由于夏季大气水汽含量高，路径辐射所占比重高于冬季，发射率误差引起的地面热辐射误差对于辐亮度误差的贡献小于冬季；(2) 随着观测天顶角的增大，大气透过率降低，地面热辐射误差对于辐亮度误差的贡献逐渐减小。(3) 吸收波段地面热辐射占总辐亮度的比例较小，使得辐亮度计算误差相比发射率误差更小，符合场景生成的精度要求，验证了仿真结果的可信度。
-
仿真场景情形为：卫星运行在地球同步轨道，其传感器视场中心提前调整至火箭发射点；视场内的地表海拔最低点13.5 m，最高点2.16 km，平均海拔1.12 km。采用BEM发动机[15]的喷口参数进行尾喷焰图像的仿真。火箭的发射参数、卫星及传感器参数如表3所示。
Satellite orbit Geographic Lat/Lon (0°, 100°E) Orbit altitude/km 35793 Infrared sensor Geographic Lat/Lon of line-of-sight center (42°N, 116°E) Number of pixels 256×256 Launch of the rocket Geographic Lat/Lon of launch site (42°N, 116°E) Altitude of launch site/km 1.43 Azimuth angle/(°) 45 Steering angle/(°) 60 Table 3. Scene simulation parameters
火箭助推段共飞行50 s，此时飞行高度为36 km。图12是轨迹上的部分时刻的辐亮度图像，目标及其放大后的形态用红线框出。可以看出，除部分像元由于MOD11_L2中存在温度无效值无法计算辐亮度以外，整个视场范围内地表背景均匀，辐亮度以路径辐射为主，地表热辐射在传感器响应函数范围内绝大部分被大气吸收。
Figure 12. Part of radiance images on the rocket trajectory
火箭在大气层以内飞行时，尾喷焰高温区域相对较小和大气吸收的双重作用，使得火箭尾喷焰所在像元在30 s (H=13.6 km)时的辐亮度仅比背景高出11%；在40 s (H=23.7 km)时比背景高出64%，达到可分辨状态；50 s (H=36.0 km)时助推器关机前的时刻，尾喷焰所在像元的辐亮度达到最高，在背景的4倍以上。
目标在图像中所占像元数理论上由面积决定，实际取决于尾喷焰的高温区域分布和目标检测跟踪算法的阈值等参数。如图12和图13所示，尾喷焰在图像中并不总是单像素，在20 s时虽然尾喷焰流域仅有25 m左右，但是高温区域恰好处在像元边界，所以占到2个像元；在50 s (H=36.0 km)也即一级分离前的时刻，尾喷焰膨胀至轴向长度接近200 m，占2个像元。如图14所示，尾喷焰在高空膨胀使得总辐射强度也急剧增大，两者变化规律几乎完全一致。
Figure 13. Plume length and number of pixels on the rocket trajectory
Figure 14. Total radiation intensity trend on the rocket trajectory
-
设置的第二个场景发射点位于更低纬度(32°N, 116°E)。每个发射点计算4条轨迹，发射转向角相同，发射方位角间隔为15°。两个场景尾喷焰总辐射强度变化的仿真结果分别如图15和图16所示。
Figure 15. Total radiation intensity trend on different trajectories of scene 1
Figure 16. Total radiation intensity trend on different trajectories of scene 2
同一场景内不同轨迹尾喷焰总辐射强度随飞行时间逐渐增大，但总体来说非常相近，变化规律一致。经过分析，(1) 地球扁率对重力的影响，使得沿不同方位发射的火箭飞行高度和速度差异随飞行时间逐渐增大，从而影响尾喷焰的温度分布和面积；(2) 卫星传感器对于尾喷焰的观测视向差异也随飞行时间逐渐增大，关机前时刻不同轨迹下的观测视向角度最大差异在6°左右。而不同场景之间传感器对于尾喷焰的观测视向差异则相对明显，尾喷焰总辐射强度差异较大。
-
文中构建了完整的包含卫星传感器、火箭尾喷焰目标和地面背景在内的动态场景，提出的地表发射率计算方法利用神经网络充分挖掘了大气传输过程中的潜在信息，提高了发射率图像的空间分辨率，一定程度上解决了吸收波段地面背景红外辐射特性研究的数据需求。对场景仿真结果的分析表明，大气吸收使得4.3 μm波段下地面背景辐亮度很低，但同时也造成了目标被卫星首次发现的时刻延后的问题。此外，尾喷焰在图像中的像元数也有所变化。
对于同一型火箭，不同场景和发射参数对尾喷焰辐射强度的变化都有影响。实际场景中的轨迹是通过发射点和落点位置根据实际需求计算的，在一定范围内可以调整。因此对不同的发射试验，飞行高度和速度随时间变化的差异更为复杂，这些因素共同作用使得同一型火箭的尾喷焰总辐射强度序列存在类内差异。
文中软件仿真生成的图像序列有助于卫星下视场景下目标和背景红外辐射特性的研究，也为基于卫星图像序列的目标检测、跟踪和识别研究提供了参考。
DownLoad: