空间目标地基观测红外辐射特性研究

郑鸿儒，马 岩 ，张 帅，陈亚涛

（1.北京跟踪与通信技术研究所,北京 100094；2.长光卫星技术股份有限公司,吉林 长春 130000；3.北京航空航天大学 宇航学院,北京 100191）

1 引言

近年来，随着卫星技术的不断发展，各国对空间资源的争夺越来越激烈。空间目标的探测、识别和监视对于国防的意义日益凸显[1,2]。开展空间目标的红外辐射特性研究对于敌方目标探测和我方红外隐身都具有重要意义。

国内外空间目标红外辐射特性研究主要包括目标表面温度场计算、目标红外辐射计算，以及通过地基或天基平台的目标观测等。目标表面温度计算主要是通过积分法或蒙特卡洛法等获得表面热流，然后结合材料属性迭代求解表面温度。获得表面温度后，结合表面发射和反射特性即可计算出目标红外辐射值。进而可以得到地基或天基平台的可观测性以及目标光谱辐射特性。其可为分析目标工作状态，反演目标表面材料属性提供参考。

刘建斌等人[3]将目标分别视为朗伯表面和随机粗糙表面，计算得到了风云一号卫星在地面观测位置处的光照度。申文涛等人[4]采用随机起伏表面算法模拟材料表面，并结合双向反射分布函数（BRDF）计算了卫星在可见光和红外谱段的辐射特性，并与地面卫星模型试验进行了对比。杨帆等人[5]计算了探测器接收到的红外辐射光谱，采用最小二乘算法反演了目标各组件的有效面积和表面温度。李文豪等人[6]分析了散热面在不同功率下在距离5 km 处入瞳的中波和长波辐射照度。谷牧等人[7]考虑了地基探测器的可见范围及大气衰减的影响，分析了地基探测器入瞳处影响红外光谱辐射照度的主要因素。孙成明等人[8]对风云一号卫星通过有限元法进行了建模，并与缩比模型散射辐照度的实验结果进行了对比分析，结果基本一致。Skinner 等人[9-10]对地球同步轨道卫星的红外光谱进行了长期观测和分析，获得了等效面积和等效温度等信息。

然而，以上研究要么对目标的建模过程介绍较少，要么采用结构化面网格进行计算，当空间目标结构比较复杂时建模及处理遮挡的难度较大。其次，在计算轨道外热流时，通常将太阳光线与表面夹角的余弦值作为太阳辐射角系数，但该方法仅适用于凸表面。最后，以往研究中对于目标红外辐射研究较多，对考虑大气衰减的地基探测场景研究较少。

基于此，本文基于非结构化四面体网格开展有限元建模，编写空间目标红外辐射特性仿真程序，采用蒙特卡洛算法追踪光线，计算轨道外热流。结合表面材料属性计算得到目标表面温度和辐射特性。进一步，考虑大气衰减效应，对目标在地基探测器入瞳处接收到的红外光谱辐射特性进行研究，给出升轨和降轨过程中典型位置的红外辐射光谱特性。

2 计算方法

本文采用有限元方法开展了空间目标红外光谱特性研究，探测器接收到的红外辐射主要有目标自身辐射、目标反射太阳辐射、目标反射地球辐射和目标反射地球反照辐射[11]。

2.1 仿真算法流程

首先，需要对目标模型进行有限元建模，采用非结构四面体网格对空间目标进行离散，光线的步入和追踪均以网格为单位进行，网格结构也是信息存储结构的基础。程序中网格元、面元和节点分别定义为相应的结构体，并通过指针或序号相互映射。软件通过读取CGNS（CFD General Notation System）格式网格文件，自主建立网格-面元的映射关系，判断相邻网格和边界网格等，并将空间目标表面温度和材料属性等参数离散到各表面中心点。

然后，通过矢量坐标变换得到各面元中心点坐标在目标本体坐标系下与地球、太阳的位置关系。通过蒙特卡洛方法，由目标各表面面元发出给定数目的光线，并通过计算光线是否被遮挡，以及光线方向和太阳矢量、地球矢量的关系，汇总所有光线计算结果，得到目标接收到的轨道外热流。此处的外热流包括太阳辐射、地球辐射及地球反照辐射，忽略其他天体的影响。

详细计算表面温度特性十分复杂，本文在计算时忽略了星体各表面面元之间的导热和辐射换热，仅考虑太阳帆板正反面的导热。使用牛顿迭代法求解表面热平衡方程，并根据温度特性计算红外辐射量。在计算轨道外热流、各表面对给定探测器入瞳方向的辐射强度的过程中使用Open-MP 并行加速算法进行加速。

最后，汇总各表面面元的参数计算结果，得到目标温度及红外辐射特性。

2.2 自身辐射

表面自身辐射由表面温度和表面材料发射率决定，在获得各表面温度后，目标自身红外光谱辐射出射度可由下式表示：

其中：c1为第一辐射常量，其值为3.742×10-16W·m2；c2为第二辐射常量，其值为1.438 8×10-2m·K；λ为波长，单位为 μm；Ti是面元表面温度；εi(λ)是面元光谱发射率，根据文献[5]知，在红外谱段发射率基本不变，本文设置为常值。

2.3 反射太阳辐射

计算目标表面反射太阳辐射时，通常将太阳视为一个温度为5 900 K 的黑体，由普朗克公式得到的光谱辐射出射度为：

其中，Ts是太阳温度。则太阳在面元上的光谱辐射强度为：

其中：Rs是太阳半径；Ds是太阳-地球之间的距离；ϕs为太阳辐射角系数。面元反射太阳辐射光谱辐射出射度为：

其中，ρs为表面材料对太阳辐射的反射率。

2.4 反射地球辐射

面元反射地球辐射包括地球红外辐射和地球反照太阳辐射。地球红外辐射可以将地球视为一个温度为255 K 的黑体，由普朗克公式得到光谱辐射出射度为：

其中，Te是地球温度。则地球自身辐射在面元上的光谱辐射强度为：

式中Rsat为卫星与地心的距离，Re是地球半径，ϕe是地球红外辐射角系数。地球反射辐射在面元上的光谱辐射强度为：

式中 ρa为地球对太阳辐射的平均反射率，ϕa是地球反照辐射角系数。

面元反射地球辐射光谱辐射出射度为：

其中，ρe为表面材料对地球红外辐射的反射率。

3 表面属性及模型参数

3.1 表面反射特性

对于卫星表面反射参数设置，本文将太阳帆板背阳面设置为漫反射，根据朗伯余弦定律计算。对于镜面属性较强的太阳帆板电池片表面和多层包覆的卫星本体表面采用双向反射分布函数（BRDF）进行计算。由于BRDF 的模型较多，计算过程也比较复杂，本文选用文献[12]中介绍的改进Phong 模型进行计算，选用GaSa 材料的模型参数计算帆板正面电池片反射特性，选用聚酰亚胺薄膜的模型参数计算卫星本体反射。计算方法如下式所示：

其中，fr(θin;θr,φ) 为表面的BRDF，θin、θr、φ分别为入射天顶角、观测天顶角、观测方位角；ρd为材料的漫反射系数；ρsp为镜面反射系数；α为镜向指数；β为观测方向与镜反射方向的夹角，具体定义可见参考文献[12]。目标表面与探测器入瞳的位置关系示意如图1 所示。

图1 目标表面与探测器入瞳的位置关系示意图Fig.1 Schematic diagram of the position relationship between the target surface and the probe's entrance pupil

对于地球红外辐射和地球反照辐射，由于每个方向基本都有能量贡献，因此本文只对太阳直射采用BRDF 计算，其他辐射采用漫反射计算面元的反射红外辐射强度[13]。则目标在探测器方向上的红外光谱辐射强度对于漫反射表面表示为：

对于采用BRDF 计算的表面表示为：

其中，i是目标表面离散的面元序号，N是面元总个数。

3.2 模型及参数设置

本文计算空间目标辐射特性采用的模型主要参考风云一号C 卫星（FY-1C），卫星本体尺寸为：1.42 m×1.42 m×1.2 m，太阳帆板位于两侧，尺寸为4.5 m×0.2 m×1.2 m。卫星在太阳同步轨道运行。太阳帆板沿飞行方向前后两侧固定安装。帆板所在平面与轨道面平行，以保证在整个轨道周期内太阳光照角度恒定[14]。轨道坐标系定义为，+X指向飞行方向，+Z指向地心，+Y符合右手定则。仿真软件运行时需要输入卫星轨道信息（半长轴、偏心率、倾角、升交点赤经、近地点幅角、真近点角）。软件内部完成递推后。认为本体系与轨道系重合。仿真中使用的6 个轨道参数如表1 所示。表面参数设置如表2 所示。卫星本体设置为多层材料。

表1 轨道参数Tab.1 Orbit parameters

表2 表面材料热参数Tab.2 Thermal parameters of surface material

卫星运行轨道降交点地方时约为8:34，太阳帆板电池片朝向为-Y平面。计算使用的模型如图2 所示，使用的网格如图3 所示，计算域为包括卫星本体在内的13 m×3 m×3m 的区域，总四面体网格数约为13.5 万，卫星目标表面面元网格数约1.4 万。

图2 仿真算例中使用的简化模型Fig.2 The simplified model used in simulation

图3 仿真算例中的网格划分Fig.3 Grid division in simulation cases

4 分析与讨论

4.1 表面温度计算结果

首先，为了验证算法的准确性，根据前述的仿真条件，计算得到不考虑遮挡影响的正六面体目标在一个轨道周期内各表面受到的轨道外热流，并将其与商业软件的计算结果进行对比，如图4（彩图见期刊电子版）所示。其中，线条为本文软件计算结果，符号为对应表面的商业软件计算结果，计算起始时刻为2000-03-21 04:58:00。可以看出，本软件的轨道外热流计算结果与商业软件相近，误差较小。

图4 不考虑遮挡效应时本文计算得到的轨道外热流与商用软件的结果对比Fig.4 Comparison of calculation results of orbit external heat flow by proposed method and commercial software without taking into account the occulusion effect

采用图3 给出的网格模型计算考虑遮挡效应的轨道外热流分布，如图5（彩图见期刊电子版）所示。可以看出，由于目标三轴稳定对地，且帆板位置固定不变，故本体-Y面与帆板电池片受到的热流最大，为1 200 W/m2左右，且在阳照区和地影区都比较稳定。其中，+X和-X侧帆板表面在阳照区相对于图4 存在热流波动，这是由于卫星本体对阳光的遮挡造成的。-Z表面在阳照区、+Z表面在卫星星下点进出地影区时存在被阳光直射区间，因此，热流随时间变化较大，其他表面受到的热流变化幅度较小。

图5 考虑遮挡效应时轨道外热流计算结果Fig.5 Calculation results of orbit external heat flow with taking into account the occlusion effect

结合表面材料参数，计算得到的表面温度，如图6（彩图见期刊电子版）所示。此处的计算结果仅考虑稳态解。假设本体各表面间相互绝热，内热源等效为各本体表面受到的热流，设为50 W/m2。帆板导热系数设置为 1.7 W/(m·K)。可以看出，除+Y表面外，其余各表面温度在一个轨道周期内变化较大。温度最高的位置在帆板上，温度在311～330 K 范围内变化，本体-Y表面阳照区温度约为316 K。帆板在地影区和阳照区温差最大，变化幅度约为150 K。其他表面温度变化趋势与热流变化基本相同。

图6 各表面温度Fig.6 Temperature of each surface

4.2 红外辐射强度计算结果

进一步地，计算目标表面在3～5 μm 和8～14 μm的红外辐射强度，如图7～8（彩图见期刊电子版）所示。此处选择的位置点为阳照区过赤道时刻，方位角定义为绕Z轴逆时针旋转方向，0 点位于+X轴，天顶角定义为绕Y轴逆时针旋转方向，0 点位于+Z轴，探测器入瞳法线指向目标，计算步长为10°。

图7 3～5 μm 波段红外辐射强度Fig.7 Infrared radiation intensity in 3～5 μm band

图8 8～14 μm 波段红外辐射强度Fig.8 Infrared radiation intensity in 8～14 μm band

从图中可以看出，在3～5 μm 中波和8～14 μm长波波段的辐射强度分布都存在天顶角90°、方位角90°和270°两个峰值，与文献[4]分布趋势一致。这主要是因为帆板反、正两面的温度均较高，而空间目标反射的太阳辐射主要在可见光波段，地球辐射相对较小，因此，目标自身温度产生的红外辐射占比较高。

对于中波波段，辐射强度最大值在45 W/sr左右。对于长波波段，辐射强度最大值在770 W/sr 左右。两个波段在天顶角90°、方位角90°，即本体+Y轴方向强度稍高。这是因为帆板正反面温度相近的情况下，太阳帆板背板的发射率设置稍高于正面电池片。

4.3 地基探测红外计算结果

地基探测器位置选择为北纬43.8°，东经125.3°，高232 m 位置，设定最小俯仰角为20°后得到升轨过站弧段（UTC 时间）：2000-03-21 11:26:16 至11:34:13，降轨过站弧段：2000-03-21 23:48:26 至23:54:58。分别选取弧段的起始、中点和结束位置为观测点，两种情况下的各位置点方位、俯仰角以及距离如表3 所示。表3 中，序号1、2、3 分别代表升轨弧段的进站、中点和出站；序号4、5、6 代表降轨弧段的进站、中点及出站。地面站采用北东地坐标系。

表3 不同位置点处的方位角、俯仰角和距离Tab.3 Azimuth,elevation,and range at different positions

计算得到在3～5 μm 和8～14 μm 波段的红外光谱辐射强度分布，如图9～图10（彩图见期刊电子版）所示，光谱分辨率为0.02 μm。可以看出两个弧段都是中间位置处红外光谱辐射强度最大，升轨进站时的红外光谱辐射强度略低于出站，降轨进站时的红外光谱辐射强度高于出站。在升轨弧段对于观测贡献较大的+Z表面温度逐渐升高，其他可探测表面温度变化较小。而降轨弧段仅-Z表面温度变化较大，但其对于地基观测没有贡献，进、出站时刻的红外光谱辐射强度差别主要由地基探测位置可见的有效面积决定。

图9 升轨观测红外光谱辐射强度Fig.9 Infrared radiation spectral intensity of ascending detection

图10 降轨观测红外光谱辐射强度Fig.10 Infrared radiation spectral intensity of descending detection

从表3 可以看出，升轨弧段俯仰角较大，接近过顶，对于温度较高的-Y面帆板和卫星本体在探测器入瞳处的有效面积较小。光谱辐射强度最大值在中点位置处的9～10 μm 波段内产生，约为23 W/(sr·μm)，进、出站时刻最大值分别约为18 W/(sr·μm)和19 W/(sr·μm)。降轨弧段俯仰角较小，且探测器方向在卫星本体的-Y侧，探测器入瞳处的有效面积较大，因此光谱辐射强度较大。中间位置处最大值约为113 W/(sr·μm)，进、出站时刻最大值分别约为85 W/(sr·μm)和76 W/(sr·μm)。

进一步地，考虑到大气吸收及观测路径上大气的自身辐射，使用Modtran 软件计算大气透过率和沿程辐射值。计算条件设置为中纬度地区冬季大气，得到探测器入瞳处接收到的叠加大气衰减效应和沿程辐射的空间目标红外光谱辐射强度如图11～图12（彩图见期刊电子版）所示。

图11 考虑大气影响的升轨观测红外光谱辐射强度Fig.11 Infrared radiation spectral intensity of ascending detection with the atmosphere impact

图12 考虑大气影响的降轨观测红外光谱辐射强度Fig.12 Infrared radiation spectral intensity of descending detection with the atmosphere impact

可以看出，两种观测情况下大气影响均较大。对3～5 μm 中波波段影响比8～14 μm 长波波段大，光谱强度降低比例较多。对于两种观测条件，由于进站和出站时的俯仰角较小，沿程损失比中间时刻大。在4.3 μm 和14.7 μm 的CO2吸收带附近，以及9.5 μm 的O3吸收带附近损失较大，在进行地基红外探测时应避免使用附近波段。同时也看到，地基观测无法避免地会受到大气衰减和背景辐射干扰，导致很难准确反演目标信息，因此，联合天基平台开展空间目标联合探测很有必要。

5 结论

本文针对空间目标的红外辐射特性，基于有限元方法编写了仿真计算软件，计算得到了目标在太阳同步轨道上接受的外热流，结合表面材料属性，计算得到了轨道周期内的表面温度变化。进一步地，结合改进的Phong 模型和漫反射模型，计算了目标在中波波段和长波波段的红外辐射强度空间分布。结果表明，对地三轴稳定太阳同步卫星采用沿飞行方向前后布置固定太阳帆板时，阳照区帆板表面温度较高，约在311～330 K 范围内变化，是红外观测的主要辐射来源。采用8～14 μm 波段进行探测时，红外辐射强度最大值在770 W/sr 左右，远大于3～5 μm 波段。最后，基于地基探测场景，使用Modtran 软件计算得到了在不用观测路径上的大气衰减和沿程辐射，计算得到空间目标在阳照区和地影区观测场景中的红外光谱辐射特性。本文研究结果可为空间目标信息反演及探测波段优选提供助力。