高轨卫星的月球引力摄动效应研究

赵东方，张捍卫

（1. 河南理工大学 测绘与国土信息学院，河南 焦作 454003；2. 河南理工大学 资源与环境学院，河南 焦作 454003）

0 引言

近年来，随着我国北斗卫星导航系统（BeiDou navigation satellite system, BDS）的逐步建立，BDS越来越多地应用到测绘、导航制导、气象等领域[1]。而各项应用对定位的要求越来越高，使得对卫星定轨的要求也越来越高。而卫星在围绕地球运行过程中，受到的主要摄动包括地球非球形引力摄动、第三体引力摄动、辐射压力摄动、潮汐摄动、大气阻力摄动等[2]。各种摄动因素对不同轨道的影响是不同的。

除去地球和空间探测器外，第三个天体引力对空间探测器产生的摄动，称为第三体引力摄动。而对于围绕地球的卫星而言，第三体引力摄动主要是由日月引力产生的，称为日月引力摄动。绕地运行卫星轨道高度越高，第三体引力摄动的影响就越大。对于高轨绕地运行的卫星，第三体引力摄动的影响十分可观。第三体引力摄动对高轨卫星的影响将仅次于地球中心引力和地球非球形引力摄动[3]。

文献[4]指出对于围绕地球运行的卫星而言，第三体引力摄动影响的大小，主要取决于卫星轨道的高度、形状、轨道面位置和拱线相对于月地、日地连线的位置。文献[5]指出，第三体引力摄动对偏心率和倾角有影响，而地球扁率项J2对其没有影响。

因此，对第三体引力摄动的研究，可以探究第三体引力摄动对绕地运行卫星轨道的影响规律和大小，对于卫星轨道保持、自主定轨等都有积极的意义。

由于第三体引力摄动对于高轨卫星的影响更为显著，而BDS 中有多颗地球静止轨道（geostationary Earth orbit, GEO）卫星和倾斜地球同步轨道（inclined geosynchronous orbits, IGSO）卫星的高度都在35 000 km 左右[6]，属于高轨卫星，受到第三体引力的影响较大。因此，对第三体引力摄动的研究，也对BDS 也有积极的意义。

本文主要研究的是在不同卫星轨道面与地月连线向量夹角下，月球引力摄动对高轨卫星的影响。

1 基本原理

卫星运动的微分方程都很复杂，除最简单的二体等少数问题外，到目前为止都不能给出严格解。而卫星的实际运动受到各种摄动因素的影响，其实际微分方程十分复杂，很难给出准确的解析解。而运用数值解法，可以给出满足一定精度的离散解，可以定性与定量地分析卫星运动的各种规律[7]。

文献[8]指出，利用龙格-库特塔-费尔布尔格（Runge-Kutta-Fehlbrg, RKF）单步法，先求出12 步解，然后采用12 阶亚当斯（Adams）预估-校正公式往后积分，步长采用75 s，每隔15 min输出一组结果。结果表明，二体意义下对地球静止轨道卫星和倾斜地球同步轨道卫星的数值积分结果与真值相差甚小，均在亚毫米级或以内。文献[9]采用2005-06-07（年积日第158 天）的全球定位系统（global positioning system, GPS）PRN1 卫星星历，利用Adams 数值积分对PRN1卫星的运动方程进行求解。PRN1 卫星星历由国际全球卫星导航系统（global navigation satellite system, GNSS）服务组织（International GNSS Service, IGS）提供。积分初值采用IGS 网站下载的G 文件中，年积日第158 天12 时的数据作为初值，积分步长采用75 s，每隔15 min 输出一组结果，积分时间为20 h，将结果与IGS 精密星历比较。结果显示，该数值积分结果与精密星历差值在2 cm 左右。说明Adams 数值积分方法对于确定卫星轨道是可靠的。

1.1 RKF7、RKF8 单步积分法

RKF 方法是一种嵌套技术的龙格-库特塔（Runge-Kutta, RK）方法，利用差分格式中右函数系数可有不同选择的特点，即同时给出n阶和n+1阶两组RK 计算公式，用两组公式计算结果之差来估计截断误差，根据截断误差的大小来控制步长。本文采用的是RKF7、RKF8 来进行起步，RKF 的7 阶（RKF7）和8 阶（RKF8）公式[10]为

式中：T为截断误差；f10、f11、f12分别为式（2）中k=10、k=11、k=12 时fk的值。

1.2 Adams 预估校正法

在卫星运动方程的数值解法中，单步法的优点是起步和变步长比较容易，但是计算量较大。所以单步法通常用作多步法计算的起算数据。当采用单步法计算出足够的起算数据后，就可采用效率更高的多步法进行积分[11]。本文使用的是12 阶Admas 预估-校正法，预估用显式的亚当斯-巴什福斯（Admas-Bashforth）公式为

校正用隐式的亚当斯-莫尔顿（Adams-Moulton）公式为

1.3 滑动拉格朗日插值法

利用式（7）在X、Y、Z方向上分别进行插值处理，即可得到任意时刻的卫星坐标。为了提高插值精度，常利用高阶拉格朗日插值多项式，但是这又会出现插值区间两端不稳定，易发生震荡或跳跃现象，这种情况被称为“龙格”现象。

为了解决“龙格”现象，可以采用滑动拉格朗日插值法解决插值区间两端不稳定的问题。该方法使待内插点始终处于插值区间的中央，将插值区间作为一个固定“窗口”，该“窗口”长度固定。“窗口”每向后移动一个单位，内插点也随之向后移动一个单位[13]。本文采用8 阶滑动拉格朗日插值。文献[14]指出利用8 阶滑动式拉格朗日插方法对卫星插值可以达到毫米级。

2 数值计算

2.1 处理流程

卫星初始状态经过RKF7、RKF8 单步法积分，求得足够的多步法积分的起步数据；再由Adams预估-校正法求得卫星在一定时间内，一定时间间隔的卫星位置矢量和卫星速度矢量，转化为卫星轨道根数；最后用8 阶滑动拉格朗日插值法求得卫星加速度。

2.2 数据准备

在利用数值积分进行卫星轨道积分过程中，除了考虑月球引力外，还需考虑的力学模型包括表1中所示的力学模型。表2 是一组卫星初始时刻的轨道参数，其中a为卫星轨道的长半轴，e为卫星轨道的偏心率，M为卫星轨道的平近点角，ω为卫星轨道的近地点幅角，Ω为卫星轨道的升交点赤经，i为卫星轨道的倾角，α为地月连线向量与卫星轨道面的夹角，夹角α范围为 0 ≤α≤ 180° 。

使用表1 中的力学模型，分别对表2 中各个卫星进行轨道数值积分，卫星质量都为1 500 kg，积分时间长度为48 h，历元间隔为75 s，输出数据的间隔为15 min。

表1 力学模型

表2 卫星初始时刻轨道根数

3 结果分析

通过计算，得到了月球引力摄动加速度和两个周期轨道根数随夹角α和时间的变化情况，如图1至图12 所示。两个周期时间可根据开普勒第三定律计算。

开普勒第三定律为

式中：T为卫星运行周期；a为卫星轨道的长半轴，G为引力常数，M0为地球质量。

根据开普勒第三定律可以计算出，在不同夹角α下，卫星的2 个周期时间大约为2 874 min。

3.1 卫星月球引力摄动加速度的变化规律

图1 为卫星月球引力摄动加速度随夹角α变化的情况，图2 为卫星月球引力摄动加速度随时间变化的情况。从图2 中可以看出，月球引力摄动加速度的大小为3.0×10-6～8.4×10-6m/s2；当α从0 和180°趋向90°时，月球引力摄动加速度逐渐减小。在夹角α一定的情况下，月球引力摄动加速度随时间呈一定的周期性变化，震动幅度随着夹角α趋向90°而减小。

图1 月球引力摄动加速度随夹角α 变化情况

图2 月球引力摄动加速度随时间变化情况

3.2 卫星轨道长半轴变化规律

图 3 为两个周期长半轴增量随夹角α变化的情况，图4 为长半轴改变量随时间变化的情况。从图3、图4 可以看出，月球引力使得长半轴长度减小；两个周期长半轴长度增量绝对值最大位于α=180°处，为0.59 km；长半轴长度增量绝对值最小位于α=90°处，为0.23 km；长半轴增量绝对值随着夹角α趋向90°而减小。长半轴改变量随时间呈现周期性变化，震动幅度随着夹角α趋向90°而减小。

图3 轨道长半轴两个周期增量随夹角α 变化情况

图4 月球引力引起长半轴改变量随时间变化情况

3.3 卫星轨道偏心率变化规律

图5 为两个周期卫星轨道偏心率增量随夹角α变化的情况，图6 为偏心率改变量随时间变化的情况。从图5、图6 可以看出，当0≤α＜90°时，两个周期偏心率增量逐渐变大；当90°＜α≤180°时，两个周期偏心率增量逐渐变小。偏心率改变量随时间震荡变化，震荡幅度在α=0 和α=180°附近最大，随着α趋向90°震荡幅度逐渐减小。

图5 轨道偏心率两个周期增量随夹角α 变化情况

图6 月球引力引起偏心率改变量随时间变化情况

3.4 卫星轨道平近点角变化规律

图7 为两个周期平近点角增量随夹角α变化的情况，图8 为平近点角改变量随时间变化的情况。从图7、图8 可以看出，两个周期平近点角增量先增加后减小，再增加；增量绝对值在α=0处取得最大值为0.046°。平近点角改变量随时间呈一定的周期性波动，震荡幅度在α=0 和α=180°附近最大，随着夹角α趋向90°震荡幅度逐渐减小。

图8 月球引力引起平近点角改变量随时间变化情况

3.5 卫星轨道升交点赤经变化规律

图9 为两个周期升交点赤经增量随夹角α变化的情况，图10 为升交点赤经改变量随时间变化的情况。从图9、图10 可以看出，当α=120°时，升交点赤经增量绝对值最大为 0.000 25°。当0≤α＜90°时，升交点赤经改变量随时间减小；当 90°＜α≤180°时，升交点赤经改变量随时间增加。

图9 轨道升交点赤经两个周期增量随夹角α 变化情况

图10 月球引力引起升交点赤经改变量随时间变化情况

3.6 卫星轨道倾角变化规律

图11 为两个周期轨道倾角增量随夹角α变化的情况，图12 为轨道倾角改变量随时间变化的情况。从图11、图12 可以看出，当α趋向90°时，两个周期轨道倾角增量逐渐变大；在夹角α=90°时，两个周期轨道倾角增量绝对值最大为0.000 09°；当夹角α=0 或α=180°时，两个周期轨道倾角增量绝对值最小为0.000 014°。轨道倾角改变量随时间都逐渐增加。

图11 轨道倾角两个周期增量随夹角α 变化情况

图12 月球引力引起倾角改变量随时间变化情况

4 结束语

本文以轨道长半轴为42 167.26 km 的一组卫星为研究对象，通过数值方法对其进行积分。计算结果表明，经过两个周期，当卫星轨道面与地月连线共面时，月球引力对加速度、轨道长半轴、轨道偏心率和轨道平近点角的影响最大；当卫星轨道面与地月连线垂直时，月球引力对轨道倾角的影响最大。其中，月球引力对卫星轨道的长半轴和平近点角影响较大，分别可达0.59 km 和0.046°。这将有助于认识和了解短时期月球引力对高轨航天器造成的影响。