基于影锥特征角的近圆轨道卫星进出地影时间的快速判定法

谢文杰，邵圣祥，王 霞，宋易敏

(1. 中国人民解放军63611 部队，库尔勒841001； 2. 中国人民解放军63610 部队，库尔勒841001)

1 引言

当太阳光线受到地球遮挡后会在地球背后形成阴影区域，称为地影。 许多绕地运行的卫星都要周期性地进出地影。 卫星位于地影区时，由于星体受到的太阳辐射大幅减少或消失，不仅影响卫星服务功能的提供[1-3]，而且限制航天发射和测控任务的实施[4-7]。 因此，准确、快速计算卫星进出地影的时刻及阴影时间对制定卫星温控策略、确定太阳辐射压大小、规划航天任务等至关重要[8-11]。

在工程应用中主要有2 种地影模型：柱形地影模型和锥形地影模型[12]。 柱形地影模型比较简单，只能描述卫星被完全照亮或遮挡2 种状态，适用于低轨卫星的阴影预报[1]，高轨卫星建议谨慎采用[2]。 锥形地影模型考虑到太阳被遮掩的程度，构造了相对复杂和真实的阴影结构，导致阴影判断条件增多、计算量加大[13-14]，但是却成为各种高精度应用的必然选择[3]。 目前关于锥形地影模型阴影判定依据，一种是利用卫星的角度得出太阳、地球的视半径、角距和视面积来判断[2]；另一种是根据卫星到阴影轴的距离与影锥在基本面上的投影半径进行判断[4]。 由于卫星时刻处于运动之中，根据卫星位置确定的天体视半径、投影基本面也随之不断改变，所以在迭代判定卫星阴影状态的过程中需要更新每一时刻的阴影条件相关参数。 现有文献没有对锥形地影模型不同阴影计算方法的预报精度、计算效率进行定量分析。

本文借鉴贾向华[1]论述的阴影特征角和判定角的概念，把变化率很小的影锥半锥角作为比较的参考基准，提出了一种新的锥形阴影判定条件。 只对背向太阳侧的卫星轨道进行地影计算，并给出了基于时间窗口数据筛选用于近圆轨道卫星的地影快速计算方法；使用STK 和MATLAB 工具对3 种轨道高度的卫星进行地影仿真计算，将不同算法计算结果与STK 结果进行比较和分析，以验证本文算法的精度和效率。

2 锥形地影模型及其判定条件

将太阳视为有一定大小的发光球体，忽略地球扁率和大气的影响，卫星为一个不计尺寸的点目标，太阳光以直线形式传播，这样因地球遮挡在空间形成的圆锥形阴影如图1 所示。 在地球的背后形成了太阳光完全照射不到的区域和受到部分太阳光照射的区域，前者称为本影，后者称为半影。 以日心Os和地心Oe连线为轴，整个地影由本影圆锥和半影圆锥同轴互相嵌套构成。V1为半影圆锥的顶点，α为半影锥的半锥角；V2为本影圆锥的顶点，β为本影锥的半锥角。 本文结合圆锥形阴影几何特征定义阴影特征角和判定角，根据空间星体几何关系推导出它们的求解方法和卫星进出地影的判定条件。

图1 锥形地影模型Fig.1 Conical Earth shadow model

定义rs、re、rt分别为太阳、地球和卫星在惯性坐标系(如J2000.0 坐标系)中的位置矢量，则太阳、卫星相对于地球的坐标s、st分别为式(1)、式(2)：

Rs、Re分别为太阳和地球的半径，由相似三角形和比例的性质可以计算得到式(3)、式(4)：

半影判定角是影锥顶点V1和卫星位置Ot的连线与阴影轴的夹角，角的开口方向与α一致，记作θ；本影判定角是影锥顶点V2和卫星位置Ot的连线与阴影轴所夹的非钝角，记作φ。s1、s2分别为地心Oe指向影锥顶点V1、V2的矢量，φ为辅助变量。 以上夹角和矢量由式(5)～式(9)确定：

在地球一个公转周期中，地球与太阳之间的距离在1.47×108km 到1.52×108km 之间变化[8]，则α的变化区间为0.264 8°～0.273 8°，β的变化区间为0.260 0°～0.268 8°，α的变化率为5.707 8×10-10°/s、β的变化率为5.580 9×10-10°/s，可见α、β随时间变化很小，在短期内可以视为常值。 由图1可以看出，地球的本影区实际上内嵌在半影圆锥空间之中，所以确定卫星是否位于地球的阴影内，首先要判断卫星是否在外层圆锥范围内，然后再判断卫星是否在内层圆锥区域内。下面给出地影的判定条件，如式(10)所示：

式中，θ＝α或者φ＝β时，表示卫星正好处于进/出半影或者本影位置。

3 卫星地影快速计算方法

太阳、地球和卫星在空间沿各自的轨道运动，只有在满足特定的位置关系时，卫星才可能进入到地影中。 因此，实际上仅有一部分卫星轨道穿越地球阴影区，显然对所有轨道数据进行计算不是合理的做法。 本文提出一种数据筛选方法，先剔除卫星轨道中不可能位于阴影区的数据，然后根据第2 节提出的锥形阴影判定条件对剩余的轨道数据进行计算，从而大幅度减少数据量，以获得较高的计算效率。

太阳、地球和卫星的空间几何关系可以用太阳地心矢量s和卫星地心矢量st之间的夹角σ来描述，如图2 所示。 卫星从地球向阳的一侧朝着地球背阳一侧运动，当卫星位于位置O′t时(图1)，σ为90°，卫星与半影圆锥母线交于点M，M到地心的距离为d，其计算公式如式(11)所示：

图2 太阳、地球和卫星的空间位置关系Fig.2 Relative position of the Sun，the Earth and the satellite

图3 夹角σ 与卫星地影的关系Fig.3 Relationship between angle σ and Earth shadow of the satellite

把上述星体几何约束关系转换为对应的卫星轨道数据时间窗口，可以剔除不可能发生地影的时段。 如图2 所示，以地球为中心，太阳在黄道面上运动，卫星在其轨道面上运动，太阳地心矢量s在卫星轨道面上的投影与卫星地心矢量st所成的夹角η和σ的变化趋势一致。 根据几何中的三垂线定理，当σ＝90°时，η＝90°。 这样可以用夹角η在同一个轨道平面中等价表征太阳、地球和卫星的位置关系，并且方便运用天体运动规律计算卫星从某一时刻初始位置到达特定位置所需经历的时间。

为方便计算，太阳、卫星轨道数据时间采用约简儒略日表示。 地心卫星轨道坐标系定义为：坐标原点位于地心，x轴在轨道面内且指向近地点，y轴在轨道面内且与x轴垂直，指向由x轴逆时针旋转90°确定，z轴与x轴、y轴构成右手螺旋系。 地心赤道坐标系到地心卫星轨道坐标系的转换如式(12)所示：

式中Xe为地心赤道坐标系矢量，Xo为地心卫星轨道坐标系矢量，RX、RZ分别为绕x、z轴旋转的算子，T 表示转置矩阵[8]，Ω是升交点赤经，i是倾角，ω是近地点辐角。

目前大多数卫星都采用近圆轨道[8，15]。 这样可以近似认为卫星作匀速圆周运动，同时把太阳也作匀速运动处理。 太阳和卫星在圆形轨道上的角速度分别如式(13)、(14)所示：

式中365.25 代表地球公转周期，地球引力系数GM＝398 600.441 5 km3/s2，a为卫星轨道半长轴。 在卫星轨道坐标系中，使用式(15)计算太阳投影和卫星的真近点角ν，以确定某一时刻太阳、卫星在轨道面的相对位置及其夹角。

式中x、y分别为太阳或卫星在地心卫星轨道平面中的位置坐标。 假设卫星从某一位置运行到另一位置绕地心旋转了Δν弧度，则所需时间t为式(16)：

当卫星使用逆行轨道，即与太阳运动方向相反时，式分母取＋号。

综上，依据初始时刻太阳、卫星的相对位置关系可以计算出卫星在之后周期运动中可能会穿越地影的时间窗口，并对这些时间段内的轨道数据进行阴影判定，此时阴影判定条件中的第2 个不等式可以忽略。

4 仿真计算与分析

在适当的近似假设下，基于锥形地影模型，利用卫星处的太阳、地球的视半径和视面积判定卫星阴影状态的方法为天体视半径判定法[2，4，8](算法1)；利用半影锥、本影锥在基本面上的投影半径及卫星到阴影轴的距离进行卫星地影解算的方法为影锥基本面投影半径判定法[4，8](算法2)。 这2 种算法简单实用，输入高精度太阳和卫星星历，能够计算得到满足工程应用要求的解。 定义本文提出的基于影锥特征角的阴影判定条件和快速计算方法为影锥特征角快速判定法(算法3)。 使用MATLAB R2015b软件编写计算机程序对卫星进出地影的时间进行预报和分析。

选择低轨、中轨和高轨3 颗近圆轨道卫星进行地影时间计算和分析，卫星数据如表1 所示。首先使用STK11.0.1 计算2019 年7 月30 日0 时(UTC)至2019 年8 月2 日0 时(UTC)期间3 颗卫星的地影状态，然后利用其高精度轨道外推模型获取该时间范围内间隔为1 s 的太阳和卫星J2000.0 坐标系位置矢量数据。

STK(Satellite Tool Kit)是一个经过航天任务验证的软件，在仿真分析方面具有极高的准确性[16]。 表2 为3 种算法计算得出的卫星绕地球飞行一圈，经历一次半影、本影的平均持续时间，并以STK 地影计算结果为基准，计算每种算法的相对误差。 数据表明3 种算法地影时间计算精度基本相当，半影时间最大误差5.447%，本影时间最大误差0.972%，阴影时间越短、计算的相对误差也越大。

表1 卫星数据Table 1 Satellite data

表2 卫星穿越地影的平均持续时间及相对误差Table 2 Average duration of satellite through penumbra/umbra and relative errors s

由图1 可知，卫星在周期运动中按照半影-本影-半影的顺序进出地影。 以卫星第一次进入半影的时刻和最后一次离开半影的时刻分别作为卫星进地影的时刻和出地影的时刻，将3 种算法的计算结果和STK 预报的地影时刻进行比较得到卫星进出地影时间差，图4～图6 分别是3 颗卫星在仿真计算期间进出地影时间差的变化。 从表2和图4～图6 可以看出，算法1 和算法2 对3 种卫星的计算结果完全相同；对于低轨、中轨、高轨卫星，算法3与算法1、算法2 预报的卫星进/出地影时刻的最大时间差分别为1 s、2 s、4 s，且算法3 预报的地影时刻准确度稍高于算法1、算法2。

在处理相同的数据量下，比较3 种算法的计算效率。 计算机的软硬件配置为Intel Core i7-9700 处理器、32 GB 内存、Windows 7 Service Pack 1 64 位操作系统。 表3 记录了输入3 组时间间隔为1 s、共计259 201 点太阳和卫星位置矢量数据时，MATLAB R2015b 软件测量得到的每一种算法的平均完成时间。 由表3 可知，3 种算法的计算效率高低依次为：算法3、算法1 和算法2，算法3 比算法1 效率提高25.305%，比算法2 提高52.172%，算法3 具有非常明显的时间性能优势。

图4 32060 卫星进出地影时间差值变化Fig.4 Changes of time difference of satellite 32060 entry/exit the Earth shadow

图5 29061 卫星进出地影时间差值变化Fig.5 Changes of time difference of satellite 29061 entry/exit the Earth shadow

图6 26880 卫星进出地影时间差值变化Fig.6 Changes of time difference of satellite 26880 entry/exit the Earth shadow

表3 3 种算法计算时间的比较Table 3 Comparison of time calculated by three methods s

由仿真结果可以看出，算法3 地影判定依据正确，预报时间窗口设计合理，在略优于现有其他算法预报精度的条件下，具有较高的计算效率，适用于航天工程中大数据量、高强度地影时间计算。

5 结论

1)新的锥形地影条件以影锥半锥角为特征角和判定依据，在迭代过程中按时间段计算特征角，可大幅减少条件中参数更新的次数；根据发生地影时太阳、地球和卫星存在的几何边界条件，确定与这些特定空间位置对应的轨道数据的时间窗口，可大大减少参与地影计算的数据量。

2)仿真计算结果表明，本文建立的地影快速计算方法预报精度略优于天体视半径判定法和影锥基本面投影半径判定法，并且计算效率分别提高25.305%和52.172%。

3)尽管影锥特征角快速判定法以地球作为遮挡物进行分析和计算，但该算法同样适用于其他天体(如月球)遮挡的情况。 下一步应研究椭圆轨道几何特征和运动规律，进一步拓宽地影快速计算方法应用范围。