空间目标轨道外热流计算及辐射特性研究

郑鸿儒，马 岩 ，张 帅，王建超，曲友阳

（1.北京跟踪与通信技术研究所,北京 100094；2.长光卫星技术股份有限公司,吉林 长春 130000）

1 引言

目前各国对空间资源的争夺愈演愈烈，空间态势感知技术对国防的战略意义日益凸显，如何对空间目标进行有效监视成为亟需解决的问题之一。空间目标由于距离太远，往往在探测器上成像只有几个像素，可提取的目标信息十分有限[1]。为了维持载荷正常工作，空间目标不可避免地会向外辐射多余热量，因此，通过目标红外信息获得目标的工作状态是空间态势感知系统的重要组成部分[2]。

然而，空间目标红外探测试验成本较高，通过建立理论模型，开展仿真分析已成为研究目标红外特性的重要手段。其中，目标轨道外热流计算是重要一环，主要包括太阳辐射、地球红外辐射和地球反照辐射。太阳辐射及遮挡的影响相对容易计算，而地球红外辐射和反照辐射则需要关联目标姿态和受照射情况，计算过程十分复杂。国内外的相关学者开展了大量研究，主要分为针对简单结构的解析法或积分法，以及处理复杂结构的蒙特卡洛（Monte Carlo）法。前者在处理常见六面体结构时比较方便，如易桦[3]等提出了一种针对圆轨道航天器外热流的计算方法，对偏航姿态瞬态外热流进行了分析。GARZÓN[4]等人针对3U 卫星建立了简化的热流和温度模型，分析了贝塔角和导热率的影响。李志松[5]等人针对微纳卫星开展了轨道外热流计算及在轨温度分析。李世俊[6]、吴愉华[7]等人将相机简化为六面体结构，对太阳同步轨道和地球同步轨道外热流分别进行了分析，得到变姿态条件下的热流值。

积分方法虽然能够快速获得初步结果，但对于处理复杂结构及遮挡问题时则比较困难，且精度不够。蒙特卡洛方法是商业仿真软件常用的计算方法，在处理复杂结构航天器及多次反射等情况具有明显优势。Atra 等人[8]对卫星热环境的建模和分析进行了综述，并对比分析了几款商业软件的计算效果。韩玉阁等人[9]利用蒙特卡洛方法计算卫星红外辐射特征，结果表明卫星的散热面是判断卫星是否失效的依据。潘晴[10]等使用反向蒙特卡洛方法对带天线结构的立方体卫星进行了轨道外热流的计算，并与商业软件进行了对比分析，以验证其精度。刘巨[11]使用STK 软件获得太阳矢量关系，使用IDEAS/TMG 模块获得了轨道周期内的热流变化。到目前为止，国内采用自主编程对空间目标外热流的精细化计算研究的较少。

本文采用蒙特卡洛方法，基于非结构化四面体网格，编写空间目标轨道外热流仿真软件，采用OpenMP 并行加速光线计算，获得目标在任意时间、任意轨道、任意姿态下的轨道外热流，并进行了对比验证。进一步地，对存在遮挡情况下的目标表面热流情况进行了仿真分析，并结合表面材料属性，对表面温度和辐射特性开展了研究。

2 计算方法

本文卫星目标轨道外热流的计算软件主要采用矢量坐标变换法，计算顺序如下：首先读取卫星轨道参数，根据卫星轨道信息，计算得到太阳和卫星在J2000 坐标系下的位置坐标，获得转换矩阵。然后，读取网格文件，遍历每一个表面网格，计算每个表面受到的太阳辐射、地球辐射和地球反照辐射。判断卫星所在位置是否处于地影区，如果在地影区则将太阳辐射和地球反照太阳辐射值设置为0，不在地影区则总辐射为3 种辐射产生的热流之和。最后判断表面是否被遮挡，如果被遮挡，则热流值设置为0。

软件运行时卫星轨道信息输入为6 个（半长轴、偏心率、倾角、升交点赤经、近地点幅角、真近点角），并在内部完成递推。在不进行姿态调整时，本体系与轨道系重合。轨道系定义为，+X指向飞行方向，+Z指向地心，+Y遵循右手定则。

卫星轨道外热流计算中，暂不考虑目标自身温度产生的辐射。下面将逐一介绍太阳直接辐射、地球红外辐射和地球反照辐射的计算方法。

2.1 太阳辐射计算

太阳辐射及其遮挡情况计算与地球辐射和反照辐射计算不同，不需要对大量光线进行随机计算。太阳在J2000 系的位置计算过程参见文献[7]。在获得太阳矢量Ssun后，目标表面接收到的太阳辐射的热流密度可表示为

式中，Sc为太阳辐射照度，φ为表面网格法向与太阳光向量的夹角，ρ为面元反射率。太阳辐射照度随季节变化，通常取平均太阳辐射照度=1 367W/m2进行计算：

式中Rs为日地距离，为平均日地距离。目标在光照区受到太阳辐射，在地影区则设置Sc为0。地影区通过计算当前时刻卫星视角下的地球张角βe判断，当卫星与地球连线和卫星与太阳连线的夹角大于 βe时，卫星受到光照。βe的计算方法如下：

式中Re为地球半径，Rsat为卫星与地心的距离。

2.2 地球红外辐射计算

地球红外辐射及反照外热流示意图如图1 所示。计算地球红外辐射时认为地球处于热平衡状态，任意表面的辐射强度均匀且相等，则地球辐射热流的表达式[3,12]为：

图1 地球红外辐射及反照外热流示意图Fig.1 Schematic diagram of earth’s infrared radiation and albedo external heat flow

其中，ρe为地球反射率，本文设定为0.3，ε为面元红外发射率；α1和 α2分别为地球表面微元 ds与目标表面微元 dA的连线与二者法向之间的夹角，L为目标表面微元与地球表面之间的距离。令，为地球红外辐射角系数，是主要求解对象。

通过积分方法计算时，主要判断面元法向与目标-地球连线的夹角 β 与 βe的余角之间的关系，令k=Re/Rsat，则 ϕe可由下式得到：

当(π-arccosk)≤β≤ 时，ϕe=0。

地球红外角系数公式的推导详见文献[13]。在使用蒙特卡洛方法时，则通过统计光线实现地球红外角系数的计算，此时地球辐射热流表达式变换为：

式中，N是面元发出的总光线数量，Nabs是被吸收的光线数量，θi是光线发射位置的天顶角。

2.3 地球反照辐射计算

地球反照辐射计算中将地球设定为漫反射，依据兰贝特余弦定理，则目标表面接收到的地球反照辐射热流[3,12]可表示为：

令 ϕer为地球反照角系数，(cosα2,0)·max(cosη,0)/πL2ds。由于涉及光线与地球表面交点，以及和太阳的位置关系，计算较为复杂。在工程上，地球反照辐射热流可以由地球辐射热流得到:

其中，ϕ为卫星-地球连线与太阳光的夹角。在使用蒙特卡洛方法计算时，地球反照辐射热流变换为：

2.4 温度场计算

详细计算空间目标温度场的分布十分复杂，本文对计算条件进行了一定简化。目标在大气层外时，忽略热对流；除帆板外，忽略星体各表面之间的导热和辐射；除太阳和地球外，忽略其他天体辐射的影响。在这种条件下，空间目标接收到的外热流即为太阳辐射、地球辐射和地球反照辐射。星体表面热平衡方程如下：

对于太阳能帆板电池片表面，平衡方程如下：

对于太阳能帆板背板表面，平衡方程如下：

其中，αs、αe、αa分别为目标表面对太阳辐射、地球辐射、地球反照的吸收率，太阳辐射和地球反照辐射吸收率取值为材料吸收率，地球辐射吸收系数为表面发射率 ε；qs,i、qa,i、qe,i分别为第i个面元接收到的太阳辐射、地球反照、地球辐射的热流；Ai是第i个面元的面积，Ns是该表面上面元的总数目，A是该表面所有面元的总面积，σ为斯忒芬－玻耳兹曼常数，σ=5.67×10-8W/(m2·K4)；T为目标表面温度；T1为帆板电池的表面温度；T2为帆板背板的表面温度；ηs为太阳能帆板的光电转换率，设为0.2；qin为目标表面的内热源，此处等效为面热源，对于本体表面，其值为20 W/m2。太阳能帆板蜂窝材料厚度h设为20 mm，热导率ζ设为1.7 W/(m·K)[14]，表面材料属性如表1 所示。

表1 表面材料热参数[15]Tab.1 Thermal parameters of surface material

2.5 计算结果验证

为了验证仿真软件计算结果的准确性，针对800 km 太阳同步轨道和地球同步轨道开展轨道外热流计算，并与文献[5]中给出的结果进行对比分析。仿真中采用的光线数为10 000，轨道递推间隔为90 s，仿真日期选择为春分日，对比结果如图2～图3（彩图见期刊电子版）所示。其中大写字母（线条）代表文献中给出的各表面结果，小写字母（符号）代表本文仿真结果。其中q为包括太阳辐射、地球辐射和地球反照辐射的轨道总热流。由图2～图3 可以看出，两种轨道高度下轨道外热流的计算值与文献中的计算值误差小于5%，具有较高的一致性。

图2 800 km 太阳同步轨道外热流对比Fig.2 Comparison of external heat flow of the 800 km sun-synchronous orbit

图3 地球同步轨道外热流对比Fig.3 Comparison of external heat flow of the geosynchronous orbit

此外，对地球辐射角系数也进行了验证，如表2 所示，给出了与文献[16]的对比结果。可以看出，两者一致性较好，可以认为本软件计算精度较高。

表2 典型位置地球红外角系数随平板俯仰角变化对比Tab.2 Comparison of earth infrared angular coefficients varying with plate pitch angle at typical positions

2.6 计算域

进一步地，建立带帆板的模拟卫星模型网格，计算域如图4 所示。卫星本体系与轨道系重合（本体系原点为卫星质心）。帆板为固定式帆板，分布在+Y和-Y侧，卫星本体尺寸为30 cm×20 cm×30 cm，帆板尺寸为30 cm×40 cm×2 cm，帆板边缘与星体±Y侧平面4 cm，网格数量为28 万。

图4 计算域Fig.4 Computational domain

2.7 遮挡计算

在计算目标本体结构对太阳辐射的遮挡时，由网格中心点沿着卫星指向太阳的矢量方向发射一条光线。通过给定足够大的步长，保证矢量终点在计算域外。然后遍历目标所有面元，如果光线矢量穿插了其他本体表面，则将该光线所在面元太阳辐射热流值设置为0。

在计算对地球红外辐射和地球反照辐射的遮挡时，由表面面元向 2π 空间内随机发射N条光线，跟踪每条光线的运动。在计算光线运动路径时，采用临近网格检索方法，通过计算光线路径与四面体网格的4 个三角面元相交情况，当光线穿插本体其他表面时删除光线，以此循环，直到光线到达计算域边界。此时，根据光线、目标-地球连线、地球-太阳连线等的关系，统计地球红外辐射和地球反照辐射量。

3 分析与讨论

3.1 不同模式下轨道外热流计算

空间目标的辐射特性除与所处位置、结构形状、表面属性因素相关外，还受姿态变化的影响。对于天基目标观测尤为明显。通常，三轴稳定卫星在轨长期模式可分为三轴对地、三轴对日等。本节对固定帆板的卫星在两种模式下的春分日轨道外热流进行了计算。对日坐标系定义为，-Z轴指向太阳，+Y轴指向黄北极，X轴遵循右手定则。轨道参数为535 km 太阳同步轨道，降交点地方时为10:30，计算起始时刻为春分日UTC 时间12:30:00，迭代时间间隔为120 s，轨道周期为5 720 s，其中2 838 s 至4 901 s 为地影区。

图5～图10（彩图见期刊电子版）给出了卫星本体各表面在一个轨道周期内的3 种热流变化情况，其中mode 1 代表三轴对地模式，图中以实线表示，mode 2 代表三轴对日模式，图中以点划线表示。可以看出，在对地模式下，除+Y面以外，本体其余各表面均受到太阳辐射，由于太阳辐射在总辐射中占比较大，因此，在卫星设计上一般选择+Y面为散热面。在各面中，-Y面受到太阳帆板遮挡，太阳辐射在0～550 W/m2范围内波动，相较于其他受晒表面幅值较小。+Z面长期对地，只有进出地影区时短时间受晒，峰值在510 W/m2左右。除-Z面以外，其他表面均受到地球辐射及地球反照辐射影响。对地模式下，各表面受到的地球辐射恒定，反照辐射随时间变化。

图5 +X 面外热流曲线Fig.5 External heat flow of +X surface

图6 -X 面外热流曲线Fig.6 External heat flow of -X surface

图7 +Y 面外热流曲线Fig.7 External heat flow of +Y surface

图8 -Y 面外热流曲线Fig.8 External heat flow of -Y surface

图10 -Z 面外热流曲线Fig.10 External heat flow of -Z surface

在对日模式下，各表面热流情况与对地模式有显著不同。太阳辐射仅存在于-Z面，在光照区恒定，约为1 378 W/m2，其他表面只受地球辐射和反照辐射，与太阳辐射相比量级较小，因此受热相对稳定，在卫星热设计上，±X、±Y、+Z面均可设计为散热面。对于地球辐射和反照辐射来说，对日模式和对地模式相比，主要不同是由于姿态变化产生的在时间和幅值上的差别。对于载荷在+Z面的卫星来说，从图9 可以看出，卫星在两种模式下受到的3 种辐射值均变化不大，热环境比较温和。

3.2 遮挡影响分析

图11（彩图见期刊电子版）给出了计算起始时刻对地模式下星体各表面受到的总外热流云图。从图中可以看出，由于卫星降交点地方时选择为10:30 am，太阳从-Y侧照射星体，星体和帆板的-Z、+X，-Y侧受光照，+X面热流值在1 000 W/m2以上，其他面在600 W/m2左右。受帆板遮挡影响，-Y侧部分位置不被太阳直接照射，+Y侧帆板靠近星体边缘被遮挡，说明本文编写的软件可以很好地处理复杂结构体遮挡情况。

图11 外热流分布及遮挡影响示意图Fig.11 Schematic diagram of external heat flow distribution and occluding effects

表3 给出了一个轨道周期内各表面的平均热流情况。对于三轴对地模式（mode 1），除+Y面平均热流在90 W/m2左右外，各表面热流分布相对均匀，在340～400 W/m2左右变化。而对于对日模式，-Z侧的星体表面和帆板表面热流值在972 W/m2左右，其他表面则小于150 W/m2。针对结构遮挡的影响开展了仿真分析，无遮挡情况下不考虑光线被星体结构或帆板遮挡。计算结果表明，遮挡的影响主要表现在对地模式下，由于帆板的遮挡，-Y表面遮挡后热流值降低了53.79 W/m2，由于星本体的遮挡，+Y-Z侧帆板表面遮挡后热流值降低了32.05 W/m2。对本体其他表面影响微弱，对于对日模式则无影响。

表3 各平面一个轨道周期内的平均外热流Tab.3 Average external heat fluxes of each surface in one orbital period W/m2

3.3 表面温度

结合各表面热流变化及表1 中给出的表面材料热参数，根据式（9）～式（11）计算得到各表面温度。图12～图13（彩图见期刊电子版）分别给出了对地和对日模式下的各表面温度。由于各表面的质量和比热容未知，本节计算结果仅考虑稳态解，相较于实际情况波动较大。可以看出，对地模式下各表面温度在不同时刻变化较大，其中，星本体-Z面温度变化最大，约为180 K，+Y侧温度变化最小，约为22 K。帆板由于存在导热，电池片侧（-Z）表面温度在光照区和地影区的差值相较于本体-Z面较小，约为120 K；随着时间变化，+X和-X侧表面受到阳光照射，温度变化较大，幅值约在140 K 波动。-Y面在光照区和地影区温度都比较恒定，分别为283 K 和200 K，+Z表面温度由于在进出地影区时受到太阳照射，存在两个峰值，其余时间波动较小。

图12 对地模式下的各表面温度Fig.12 Temperature of each surface in the earth-pointing mode

图13 对日模式下的各表面温度Fig.13 Temperature of each surface in the sun-pointing mode

在对日模式下，除帆板和本体-Z面在光照区和地影区变化较大外，其余表面变化幅度较小，其中波动最大的是本体+Z面，波动区间约为132 K。对于帆板和本体-Z面，光照区温度在340 K 左右波动，地影区在220 K 左右波动，整个区间相对恒定。

3.4 表面温度对比验证

为了验证轨道外热流计算的准确性，对“吉林一号”某卫星帆板温度进行仿真与在轨实测值验证分析。卫星运行在535 km 太阳同步轨道，降交点地方时为12:00，帆板位于-Y侧，帆板运行模式为对日模式。在实际应用中帆板处于非平衡状态，因此表面温度使用以下简化公式进行求解：

式中，c为比热容，取值为350 J/(kg·K)，m是帆板质量，取值为1.5 kg，As为帆板面积，此处取值为0.3，Q为帆板表面总的吸收热流，表面吸收率设为0.9。采用蒙特卡洛方法获得表面外热流后，代入公式迭代求解帆板温度，并与卫星遥测数据做对比，结果如图14 所示。

图14 太阳帆板温度对比Fig.14 Temperature comparison of solar panel

可以看出，计算结果在各个轨道周期内的变化幅值和规律与遥测数据符合较好。最高温度约为95 °C，最低值约为-76 °C，在光照区升高，在地影区下降，不断循环。本文计算结果可以为太阳帆板温度场分析提供可信参考。同时也可以看出，计算值和遥测值在部分区域存在一定偏差，可能是参数设置与实际情况的差异影响了计算结果，需要开展进一步研究。

3.5 红外辐射强度分布

目标自身红外辐射由表面温度和表面材料的发射率决定，在获得温度变化曲线后，目标自身红外光谱辐射出射度可由式（13）表示：

其中，c1为第一辐射常量，值为3.742×10-16W ·m2，c2为第二辐射常量，其值为1.438 8×10-2m·K，λ为波长，单位为μm，Ti是面元表面温度。εi(λ)是面元光谱发射率。

结合图12 和图13 给出的两种模式下的各表面温度，认为各表面为漫反射，根据朗伯余弦定律计算得到两种模式下一个轨道周期内的各方向上的辐射强度，如图15～图18（彩图见期刊电子版）所示。

图15 对地模式下的各方向辐射强度（3～5 μm）Fig.15 Radiation intensity in each direction in the earthpointing mode (3～5 μm)

图16 对日模式下的各方向辐射强度（3～5 μm）Fig.16 Radiation intensity in each direction in the sunpointing mode (3～5 μm)

图17 对地模式下的各方向辐射强度（8～14 μm）Fig.17 Radiation intensity in each direction in the earthpointing mode (8～14 μm)

图18 对日模式下的各方向辐射强度（8～14 μm）Fig.18 Radiation intensity in each direction in the sunpointing mode (8～14 μm)

可以看出，在常用的3～5 μm 和8～14 μm 两个探测波段中，两种模式下+Z和-Z方向的辐射强度较高，3～5 μm 波段最大值在1.5 W·sr-1左右。8～14 μm波段最大值在22 W·sr-1左右。主要原因是帆板温度较高且面积较大，是红外信号的主要来源。8～14 μm 波段辐射强度明显强于3～5 μm 波段。同一谱段内，对地模式下，不同时间内，+X、-X和-Y方向也存在红外信号较强时刻。但对日模式下，除+Z和-Z方向的辐射强度较高外，其他方向红外辐射值较小，目标探测存在困难。

4 结论

本文针对绕地卫星轨道外热流，采用蒙特卡洛法、基于非结构四面体网格和OpenMP 并行算法编写了仿真软件。针对文献中的典型工况进行了对比分析，验证了软件计算结果的准确性。进一步，针对535km太阳同步轨道，考虑卫星结构遮挡的影响，对不同姿态控制策略下的卫星轨道外热流进行了仿真分析，并将太阳帆板表面温度与在轨遥测数据进行了对比分析。得出如下结论：

（1）不同姿态模式下目标的轨道外热流区别较大。对于535km太阳同步，10:30地方时轨道来说，对地模式下除+Y面外，各表面热流值随时间变化波动较大，而对日模式下除-Z面本体及帆板表面波动较大外，其他表面变化较小。

（2）蒙特卡洛算法对空间目标复杂结构及结构遮挡具有很好的适应性，对地模式下，考虑遮挡后，-Y表面热流值降低了53.79 W/m2，+Y-Z侧帆板表面热流值降低了32.05W/m2。

（3）不同模式下目标各表面温度特性不同。对地模式下各表面温度随时间波动较大，使红外观测窗口规划提高了难度。在两种模式下，帆板在光照区温度较高，具有明显的红外特征，便于开展空间目标红外观测。