CN112036037B - 一种倾斜地球同步轨道的长期演化快速分析方法 - Google Patents

Abstract

一种倾斜地球同步轨道的长期演化快速分析方法，首先推导各摄动项的一次或二次平均摄动势函数，包括地球非球形摄动带谐项

地球非球形摄动田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄的1:1共振部分、日月引力摄动勒让德展开截取到4阶项、太阳光压摄动；其次，通过拉格朗日型轨道摄动方程，结合平均摄动势函数，建立IGSO轨道的二次平均半解析轨道递推器；对比分析不同摄动源和阶数考虑下的IGSO轨道长期演化情况，简化模型，进一步提高轨道递推效率；最后，借助高效的轨道递推，绘制一系列完整轨道根数以及初始历元时刻组合的动力学网格图，其中可划分的轨道根数为(e,i,Ω,ω)，它们的全部二元组合为e‑i,ω‑Ω；e‑ω,i‑Ω；e‑Ω,i‑ω，有关初始历元时刻InitialEpoch的全部二元组合为InitialEpoch‑e,InitialEpoch‑i,InitialEpoch‑Ω,InitialEpoch‑ω，根据动力学网格图完成对IGSO轨道长期演化的快速且全面的分析。

Description

一种倾斜地球同步轨道的长期演化快速分析方法

技术领域

本发明涉及一种倾斜地球同步轨道的长期演化快速分析方法，尤其涉及航天器轨道半解析递推和动力学网格图分析方法，属于航天器轨道动力学领域。

背景技术

倾斜地球同步轨道(IGSO)在卫星通讯和卫星导航方面具有重要的应用价值。对IGSO轨道的长期演化分析为相应卫星的任务设计和任务后处理提供了重要信息，而长期演化分析的快速性有助于提供更多信息并有利于任务设计的快速迭代。在地球卫星轨道的长期演化分析中，采用半解析轨道递推方法绘制相应的动力学网格图是一种有效的手段。文献(Gkolias Ioannis,Colombo Camilla.“Towards a sustainable exploitation of thegeosynchronous orbital region”,Celestial Mechanics and Dynamical Astronomy(2019)131:19)针对地球同步轨道(GEO)区域建立了一次平均的半解析轨道递推模型，并绘制了部分轨道根数组合下的动力学网格图，分析了轨道的长期演化。方法理论研究无误，但对于需要递推大量轨道进行绘制的动力学网格图，使用该文献中的方法，往往需要耗费数十个小时才能完成绘制一张动力学网格图所需的计算，并占用大量的计算资源，甚至需要借助超级计算机。因此在实际研究工作中，该文献中的方法不能方便快速地绘制动力学网格图，从而无法对IGSO轨道的长期演化进行快速、全面的分析，不利于任务设计的快速迭代。为解决这一问题，针对IGSO轨道，必须建立计算速度更快的半解析轨道递推器，从而能够方便快速地绘制更多组合的动力学网格图，实现对IGSO轨道长期演化快速、全面的分析。此外，在大部分已发表的相关文献中，针对初始历元时刻的动力学网格图鲜少提出并研究。因此针对IGSO轨道，建立计算高效的半解析轨道递推器，快速绘制一系列完整轨道根数以及初始历元时刻组合的动力学网格图，对于IGSO轨道上卫星的任务设计和任务后处理具有重要意义。

发明内容

本发明的目的是为了解决现有技术对IGSO轨道长期演化的分析速度慢、分析不全面的问题，提供一种倾斜地球同步轨道的长期演化快速分析方法。该方法针对IGSO轨道，建立二次平均的半解析轨道递推器，借助该计算高效的轨道递推器，实现对多种组合的动力学网格图的快速绘制，从而完成对轨道长期演化的快速、全面的分析。此方法以其动力学建模简单易操作，对IGSO轨道长期演化分析的快速性和充分性的优点，为IGSO轨道上卫星的任务设计和任务后处理提供参考。

本发明的目的是通过以下技术方案实现的。

一种倾斜地球同步轨道的长期演化快速分析方法，首先推导各摄动项的一次或二次平均摄动势函数，包括地球非球形摄动带谐项

地球非球形摄动田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄的1:1共振部分、日月引力摄动勒让德展开截取到4阶项、太阳光压摄动；其次，通过拉格朗日型轨道摄动方程，结合平均摄动势函数，建立IGSO轨道的二次平均半解析轨道递推器；对比分析不同摄动源和阶数考虑下的IGSO轨道长期演化情况，简化模型，进一步提高轨道递推效率；最后，借助高效的轨道递推，绘制一系列完整轨道根数以及初始历元时刻组合的动力学网格图，其中可划分的轨道根数为(e,i,Ω,ω)，它们的全部二元组合为e-i,ω-Ω；e-ω,i-Ω；e-Ω,i-ω，有关初始历元时刻InitialEpoch的全部二元组合为InitialEpoch-e,InitialEpoch-i,InitialEpoch-Ω,InitialEpoch-ω，根据动力学网格图完成对IGSO轨道长期演化的快速且全面的分析。

一种倾斜地球同步轨道的长期演化快速分析方法，包括如下步骤：

步骤一：推导平均摄动势函数；

考虑到IGSO轨道是一类不易出现奇异的轨道，且从方便应用的角度出发，选取经典轨道根数σ＝(a,e,i,Ω,ω,M)作为变量，其中a为轨道半长轴，e为轨道偏心率，i为轨道倾角，Ω为轨道升交点赤经，ω为轨道近地点幅角，M为轨道平近点角。对任意函数F，记

分别为函数F的一次平均结果和二次平均结果。

针对IGSO轨道的二次平均摄动势函数

表示为：

其中，

为一次平均的地球非球形带谐项摄动势函数，

为地球非球形田谐项摄动1:1共振部分摄动势函数，

为二次平均的太阳第三体引力摄动势函数，

为二次平均的月球第三体引力摄动势函数，

为一次平均的太阳光压摄动势函数。

1)、推导

摄动势函数：

地球非球形带谐项摄动势函数R_Zonal表示为：

其中，μ为地球引力常数，a_e为地球赤道半径，r为地球质心到卫星的矢径的大小，

为卫星的地心纬度，l为整数系数，J_l为第l阶带谐项系数，

为

的第l阶勒让德多项式。

考虑主要带谐项J₂,J₃,J₄，J₂项是一阶摄动因素，是最主要的摄动源；J₃,J₄项各为奇次和偶次带谐项的代表，J₃,J₄项对卫星轨道影响的规律能够反映出其它奇次和偶次带谐项的影响规律。将

带入式(2)中，取l＝2,3,4，并将勒让德多项式展开，得主要带谐项J₂,J₃,J₄的摄动势函数：

根据椭圆运动关系式，时间t与真近点角f之间的微分关系为：

其中，

为轨道平均角速度。

利用式(4)得如下形式的函数，所述函数是关于一个卫星轨道周期的平均值：

其中，k₁,k₂为系数，k₁,k₂＝0,1,2,…。

将式(5)带入式(3)得主要带谐项J₂,J₃,J₄关于一个卫星轨道周期的一次平均摄动势函数：

在式(6)中已消去真近点角f。除了主要带谐项J₂,J₃,J₄的一阶项外，还需要考虑J₂的二阶项

所对应的摄动势函数，根据Brouwer的理论，同理得：

一次平均的地球非球形带谐项摄动势函数

表示为：

2)、推导

摄动势函数：

地球非球形田谐项摄动势函数R_Tesseral表示为：

其中，S_lmpq(ω,M,Ω,θ_G),θ_lmpq均为中间变量；F_lmp(i)为倾角函数，仅和轨道倾角i有关，是关于i的三角函数的有限形式；G_lpq(e)为汉森系数，仅和偏心率e有关，是关于e的无限形式，但有些是有限的；C_lm,S_lm为引力位球谐系数；θ_G为格林尼治恒星时；l,m,p,q均为整数系数。

由于IGSO轨道位于地球同步轨道区域，轨道周期等同于地球自转周期，因此地球非球形田谐项摄动对IGSO卫星轨道的主要影响为1:1共振项引起的长周期项。引入变量λ＝ω+Ω+M-θ_G，在式(9)中取J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄对应摄动势函数的1:1共振部分，分别为

和

则地球非球形田谐项摄动1:1共振部分摄动势函数

表示为：

3)、推导

摄动势函数：

太阳引力和月球引力对IGSO卫星轨道的摄动均为第三体引力摄动，其摄动势函数可以统一进行推导。以地心赤道惯性坐标系为基准坐标系，将太阳、月球、地球均视为质点，给出日月第三体引力摄动势函数为：

其中，μ′为第三体引力常数，r′为地球质心到第三体的矢径的大小，ψ为地球质心到卫星的矢径r和地球质心到第三体的矢径r′的夹角，P_l(cosψ)为cosψ的第l阶勒让德多项式。

第三体在地心赤道惯性坐标系下的轨道根数为(a′,e′,i′,Ω′,ω′,f′)，其中，a′为半长轴，e′为偏心率，i′为倾角，Ω′为升交点轨道赤经，ω′为近地点幅角，f′为真近点角。将cosψ用卫星的轨道根数和第三体的轨道根数表示为：

其中，

分别为卫星轨道近焦点坐标系x轴和y轴的方向向量：

A,B是和卫星真近点角无关的项，用c_x,s_x分别表示变量x的余弦函数cosx和正弦函数sinx，则A,B表示为：

其中，△Ω＝Ω-Ω′,u′＝ω′+f′。

将式(19)截取到勒让德多项式的4阶并展开得：

将式(25)表示的日月第三体引力摄动势函数在卫星一个运动周期内做平均，消除卫星真近点角：

其中，AF_2nd,AF_3rd,AF_4th的表达式涉及到如下形式的函数，所述函数为关于一个卫星轨道周期的平均值：

利用式(20)、(26)、(27)得AF_2nd,AF_3rd,AF_4th的表达式为：

将式(28)和式(23)、(24)代入式(26)即得到一次平均后的日月第三体引力摄动势函数

一次平均后的日月第三体引力摄动势函数

已经不再包含与f有关的短周期项，但仍然包含与第三体轨道真近点角f′有关的短周期项，因此将

再对第三体轨道的一个运动周期做一次平均，以消去第三体轨道真近点角f′。

根据式(23)、(24)，辅助变量A,B和卫星真近点角f无关，但包含第三体轨道的真近点角f′。将式(26)表示的一次平均后的日月第三体引力摄动势函数在第三体轨道的一个运动周期做一次平均：

其中，DAF_2nd,DAF_3rd,DAF_4th的表达式涉及到式(5)中所述的函数。

利用式(23)、(24)、(28)、(29)和式(5)得DAF_2nd,DAF_3rd,DAF_4th的表达式：

将式(30)带入式(29)即得二次平均后的日月第三体引力摄动势函数

二次平均后的日月第三体引力摄动势函数

中需要计算第三体在地心赤道惯性坐标系中的轨道根数。对太阳轨道和月球轨道做出如下简化：太阳轨道视为圆轨道，月球轨道视为长期进动椭圆模型并保留黄白交角。

在地心赤道惯性坐标系下，太阳轨道的平均轨道根数为：

其中，AU为天文单位，T为儒略世纪，d为平太阳日，其和儒略世纪满足关系：T＝36525d。

在地心黄道惯性坐标系下，月球轨道的平均轨道根数为：

其中，km为千米单位。

对于月球轨道，式(32)为其相对黄道面的轨道根数，其中(i_M,Ω_M,ω_M)需装换到地球赤道面内，令对应根数为(i′_M,Ω′_M,ω′_M)，则：

其中，ε为平黄赤交角，ζ为一辅助角。

对于太阳引力摄动，已假设太阳轨道为圆轨道，偏心率为0，将式(31)代入式(29)可得二次平均后的太阳引力摄动势函数

为：

对于月球引力摄动，将式(32)、(33)代入式(29)可得二次平均后的月球引力摄动势函数

为：

4)、推导

摄动势函数：

对于IGSO轨道，不考虑地影对太阳光压摄动的影响。

太阳光压摄动势函数R_SRP表示为：

其中，

表示

的高阶小量，μ_SRP为：

其中，C_R为与卫星表面反射系数相关的参数，取值为1到2；S为卫星受光压作用的有效面积；m为卫星质量；P_SRP为地球表面处的太阳光压常数；△₀为日地平均距离。

式(36)中所包含的cosψ含义同式(20)，将式(20)代入式(36)中并去掉高阶小量，得：

将式(37)代入式(38)，忽略式(38)中等式右端的第二大项，仅保留第一大项，也即取到勒让德展开式的一阶项，取r′为日地平均距离△₀，得太阳光压摄动势函数为：

将式(39)表示的太阳光压摄动势函数在一个卫星轨道周期内做一次平均，利用式(27)，得一次平均后的太阳光压摄动势函数：

步骤二：利用拉格朗日型轨道摄动方程建立IGSO轨道的二次平均半解析轨道递推器；

针对所选用的经典轨道根数σ＝(a,e,i,Ω,ω,M)的拉格朗日型轨道摄动方程为：

根据式(41)，需要求摄动势函数对轨道根数的偏导数。步骤一中已给出IGSO轨道上二次平均后的摄动势函数

摄动势函数

对轨道根数的偏导数表示为：

将式(42)中的偏导数代入式(41)，即得到IGSO轨道二次平均后的半解析递推模型。将所得的二次平均半解析递推模型结合数值积分器，得到同时考虑地球非球形摄动带谐项

地球非球形摄动田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄的1:1共振部分、日月引力摄动勒让德展开截取到4阶项、太阳光压摄动的IGSO轨道的二次平均半解析轨道递推器。

步骤三：IGSO轨道摄动模型对比分析；

步骤二得到的二次平均半解析轨道递推器中考虑的摄动源和对应的摄动阶数为：1)地球非球形摄动带谐项，可选阶数为

2)地球非球形摄动田谐项，可选阶数为J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄；3)太阳引力摄动项，可选阶数为2阶、4阶；4)月球引力摄动项，可选阶数为2阶、3阶、4阶；5)太阳光压摄动项。

针对任意IGSO轨道，在使用二次平均半解析轨道递推器进行长期轨道演化分析时，通过考虑不同的摄动源和对应的摄动阶数，对比分析不同摄动和摄动阶数对IGSO轨道长期演化的影响。根据式(1)、(8)、(18)、(34)、(35)、(40)，不同摄动源和摄动阶数的摄动势函数之间是以加和的形式组成总的摄动势函数

且式(41)表明，各摄动对轨道根数的影响也是加和的形式，因此在步骤二得到的二次平均半解析轨道递推器中，可以选择是否包含某个特定的摄动源或某一特定的摄动阶数对轨道长期演化的影响。对比方法为仅改变某一特定摄动项，其他条件不变，通过轨道长期演化结果对比分析该摄动项的影响是否较小，以此判断在绘制动力学网格图时是否可忽略，以进一步加快计算。

步骤四：绘制一系列完整轨道根数以及初始时刻组合的动力学网格图；

对于IGSO轨道，在绘制动力学网格图时，根据步骤三的分析结论，选取满足精度要求且摄动项考虑最少的摄动模型，以获得最快的计算速度。若卫星面质比较大，则必须考虑太阳光压摄动对轨道长期演化的影响。

在IGSO轨道的长期演化中，初始轨道根数对轨道的长期演化起到决定性的影响。对轨道根数(e,i,Ω,ω)进行划分，其全部二元组合为e-i,ω-Ω；e-ω,i-Ω；e-Ω,i-ω。此外，初始历元时刻也会影响IGSO轨道的长期演化，因此还需要绘制有关初始历元时刻InitialEpoch的动力学网格图，其全部二元组合为InitialEpoch-e,InitialEpoch-i,InitialEpoch-Ω,InitialEpoch-ω。对于每一个划分的二元组合，其余初始根数或初始历元时刻保持不变，绘制两张动力学网格图，指标分别为轨道长期演化中的最大偏心率和轨道寿命。当IGSO轨道在长期演化过程中的偏心率达到设置的临界偏心率时，此时轨道近地点已再入大气层，轨道将逐渐衰减，视为轨道寿命终止。对IGSO轨道长期演化起决定性作用的初始轨道根数(e,i,Ω,ω)以及有较大的影响的初始历元时刻InitialEpoch间一系列二元组合的动力学网格图充分反映了IGSO轨道的长期演化规律，而经过步骤三进一步简化后的二次平均半解析轨道递推器则保证了整个过程的快速性。

步骤五：利用步骤四得到的动力学网格图，完成对IGSO轨道长期演化的快速分析。

步骤四已采用简化后的二次平均半解析轨道递推器快速绘制了有关初始轨道根数和初始历元时刻的动力学网格图，这些动力学网格图充分反映了IGSO轨道长期演化的最大偏心率分布以及轨道寿命分布。对IGSO轨道长期演化的最大偏心率分布进行分析，为IGSO轨道卫星的任务设计提供参考；对IGSO轨道长期演化的轨道寿命分布进行分析，为IGSO轨道卫星的任务后处理提供参考。步骤四中动力学网格图的快速绘制保证了本步骤的快速性，有利于任务设计的快速迭代。

综上，针对倾斜地球同步轨道长期演化快速分析问题，推导了一次或二次平均后的各摄动源所对应的摄动势函数，结合拉格朗日型轨道摄动方程建立了针对IGSO轨道的二次平均半解析轨道递推器。通过轨道摄动模型对比分析，进一步简化了二次平均半解析轨道递推器以获得更快的计算速度。对初始轨道根数和初始历元时刻间的二元组合进行划分，利用上述的高效半解析轨道递推器绘制一系列动力学网格图，既充分反映了IGSO轨道的长期演化规律，也满足了快速性的需求，进而完成对IGSO轨道长期演化的快速分析。

有益效果

1、本发明公开的一种倾斜地球同步轨道的长期演化快速分析方法，提供了一种简单易操作的建立二次平均半解析轨道递推模型的方法。其中，将日月引力摄动势函数进行了两次平均，并使用简单的解析星历，所建立的二次平均半解析轨道递推器相比于已有的一次平均半解析轨道递推器，大幅提高了计算速度。

2、本发明公开的一种倾斜地球同步轨道的长期演化快速分析方法，通过对比分析不同摄动和摄动阶数对IGSO轨道长期演化的影响，忽略部分摄动项，简化了二次平均半解析轨道递推器，进一步提高了绘制动力学网格图的速度。

3、本发明公开的一种倾斜地球同步轨道的长期演化快速分析方法，使用简化后的二次平均半解析轨道递推器进行轨道递推，快速地给出了一系列动力学网格图。其中，给出了初始轨道根数完整二元组合的动力学网格图，相比于已有文献中部分组合的动力学网格图，更为充分得反映了IGSO轨道的长期演化规律，此外，还给出了已有文献中鲜少考虑的关于初始历元时刻的动力学网格图，使得对IGSO轨道长期演化规律的分析更加充分。由于简化后的二次平均半解析轨道递推器非常高效，使得上述过程具有快速性特点，进而保证了对IGSO轨道长期演化分析的快速性。

附图说明

图1为实施例中给定IGSO轨道在多种摄动模型下的偏心率演化情况；

图2为实施例中给定IGSO轨道在多种摄动模型下的轨道倾角演化情况；

图3为实施例中给定IGSO轨道在多种摄动模型下的升交点赤经演化情况；

图4为实施例中给定IGSO轨道在多种摄动模型下的近地点幅角演化情况；

图5为实施例中给定IGSO轨道在多种摄动模型下的e-ω相位图演化情况；

图6为实施例中给定IGSO轨道在多种摄动模型下的偏心率矢量演化情况；

图7为实施例中二元组合e-i的动力学网格图，其中：图7(a)为轨道长期演化的最大偏心率分布情况、图7(b)为轨道寿命分布情况；

图8为实施例中二元组合ω-的动力学网格图，其中：图8(a)为轨道长期演化的最大偏心率分布情况、图8(b)为轨道寿命分布情况；

图9为实施例中二元组合e-ω的动力学网格图，其中：图9(a)为轨道长期演化的最大偏心率分布情况、图9(b)为轨道寿命分布情况；

图10为实施例中二元组合i-的动力学网格图，其中：图10(a)为轨道长期演化的最大偏心率分布情况、图10(b)为轨道寿命分布情况；

图11为实施例中二元组合e-的动力学网格图，其中：图11(a)为轨道长期演化的最大偏心率分布情况、图11(b)为轨道寿命分布情况；

图12为实施例中二元组合i-ω的动力学网格图，其中：图12(a)为轨道长期演化的最大偏心率分布情况、图12(b)为轨道寿命分布情况；

图13为实施例中二元组合InitialEpoch-的动力学网格图，其中：图13(a)为轨道长期演化的最大偏心率分布情况、图13(b)为轨道寿命分布情况；

具体实施方式

为了更好地说明本发明的目的和优点，下面结合附图对本发明的具体实施方式作进一步的说明。

本实施例公开的一种倾斜地球同步轨道的长期演化快速分析方法，为验证该方法，首先，选取一条典型的IGSO轨道，为具有一定偏心率的苔原轨道，作为主要研究对象。初始轨道根数、卫星表面反射系数以及面质比等参数如下表所示。初始历元时刻InitialEpoch选为UTC时间的2020年5月15日11点17分04秒，对应的儒略日为2458984.97018519。

表1选取的IGSO轨道的初始参数

步骤一：推导平均摄动势函数；

由式(8)、(18)、(34)、(35)、(40)，可得如式(1)所示的二次平均摄动势函数。

步骤二：利用拉格朗日型轨道摄动方程建立IGSO轨道的二次平均半解析轨道递推器；

由式(41)、(42)，结合数值积分器即得二次平均半解析轨道递推器。在本实施例中，数值积分器选择龙格库塔7(8)变步长积分器，误差控制设置为10^-13。

步骤三：IGSO轨道摄动模型对比分析；

将表1中的参数和初始历元时刻InitialEpoch输入IGSO轨道的二次平均半解析轨道递推器，然后按如下摄动项设置分别对本实施例给定的IGSO轨道进行长期递推，为验证二次平均半解析轨道递推器的正确性，使用已有的一次平均半解析轨道递推器做相同的轨道递推，用作对比：

(a)SA：已有的一次平均半解析轨道递推模型，用于验证所得二次平均半解析轨道递推器的正确性。

(b)DA-44：带谐项

田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄；日月引力摄动均取到4阶项；太阳光压摄动。

(c)DA-22：带谐项

田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄；日月引力摄动均只取到2阶项；太阳光压摄动。

(d)DA-22-NoTes：带谐项

不考虑田谐项；日月引力摄动均只取到2阶项；太阳光压摄动。

(e)DA-22-NoTes-J2：带谐项J₂；不考虑田谐项；日月引力摄动均只取到2阶项；太阳光压摄动。

(f)DA-22-NoTes-J2-NoSRP：带谐项J₂；不考虑田谐项；日月引力摄动均只取到2阶项；不考虑太阳光压摄动。

在本实施例中，轨道长期演化递推时长为120年。

不同摄动项设置(a)-(f)下的该IGSO轨道长期演化中偏心率演化情况表示在图1中、轨道倾角演化情况表示在图2中、升交点赤经演化情况表示在图3中、近地点幅角演化情况表示在图4中、e-ω相位图表示在图5中、偏心率矢量表示在图6中。DA-44和SA的曲线基本重合，即本发明建立的二次平均半解析轨道递推器和已有的一次平均半解析轨道递推器的结果高度吻合，验证了二次平均半解析轨道递推器的正确性。DA-22和DA-44的曲线在80年内基本吻合，且它们的最大偏心率演化基本一致，表明在本发明所考虑的精度内，仅考虑太阳引力摄动和月球引力摄动二次平均摄动势函数中的2阶项。DA-22-NoTes和DA-22的曲线在整个120年内都基本吻合，表明在本发明所考虑的精度内，不考虑地球非球形田谐项摄动。DA-22-NoTes-J2和DA-22-NoTes的曲线基本重合，表明对于IGSO轨道的长期演化，地球非球形带谐项摄动只需要考虑到J2项。本实施例中，卫星的面质比较小，DA-22-NoTes-J2-NoSRP和DA-22-NoTes-J2的曲线基本重合，表明当IGSO轨道上卫星面质比较小时，对于轨道长期演化不考虑太阳光压摄动。

对该IGSO轨道在不同摄动设置下长期演化的对比分析表明，在一般情况下，使用DA-22-NoTes-J2设置的二次平均半解析轨道递推器对IGSO轨道递推，相比于DA-44，满足精度要求且计算速度更快。若卫星的面质比较小，不考虑太阳光压摄动，使用DA-22-NoTes-J2-NoSRP设置的二次平均半解析轨道递推器对IGSO轨道递推，进一步提高计算速度。

步骤四：绘制一系列完整轨道根数组合的动力学网格图；

对本实施例中的IGSO轨道，卫星面质比较小，使用DA-22-NoTes-J2-NoSRP设置下的二次平均半解析轨道递推器对IGSO轨道递推。考虑轨道根数(e,i,Ω,ω)的全部二元组合：e-i,ω-；e-ω,i-；e-,i-ω，对于有关初始历元时刻InitialEpoch的二元组合，本实施例仅考虑其中一种：InitialEpoch-Ω。每一组二元组合对应两张动力学网格图，指标分别为轨道长期演化中的最大偏心率和轨道寿命，除考虑的二元组合外，其余初始根数或初始历元时刻保持不变。为保证动力学网格图的精度，对每个划分的参数进行了足够多的划分，如下表所示：

表2动力学网格图参数划分

在本实施例中，临界偏心率设置为0.846。当IGSO轨道演化过程中的偏心率达到临界偏心率时，则认为轨道寿命终止。动力学网格图中的每一个点都是对一条对应轨道长期演化递推的结果，因此对于每一组二元组合需要递推40000条左右的轨道。使用DA-22-NoTes-J2-NoSRP设置下的二次平均半解析轨道递推器进行轨道递推，在本实施例中，计算环境为Intel(R)Core(TM)i7-8700 CPU@3.20GHz，仅需要约30分钟即完成40000条左右的时长为120年的轨道递推，可快速绘制IGSO轨道的动力学网格图。所绘制的动力学网格图为：图7为二元组合e-i的动力学网格图，图8为二元组合ω-Ω的动力学网格图，图9为二元组合e-ω的动力学网格图，图10为二元组合i-Ω的动力学网格图，图11为二元组合e-Ω的动力学网格图，图12为二元组合i-ω的动力学网格图，图13为二元组合InitialEpoch-Ω的动力学网格图。

步骤五：利用步骤四得到的动力学网格图，完成对IGSO轨道长期演化的快速分析。

步骤四已得到一系列动力学网格图，其中，图7(b)、图8(b)、图9(b)、图10(b)、图11(b)、图12(b)表示的轨道寿命分布情况反映出IGSO轨道的轨道寿命随不同初始轨道根数而剧烈变化的现象，表明IGSO轨道的轨道寿命对初始轨道根数十分敏感，此外，部分轨道的轨道寿命较短，甚至低于20年，这为IGSO轨道上卫星的任务后处理提供了新思路，而IGSO轨道上的任务设计则需要避开这些偏心率变化较大的轨道区域。图7(a)、图8(a)、图9(a)、图10(a)、图11(a)、图12(a)也为IGSO轨道上的任务设计提供了参考。本实施例中给定的IGSO轨道的轨道寿命接近100年，但根据图7(b)、图8(b)、图9(b)、图10(b)、图11(b)、图12(b)，在考虑任务后处理时，适当改变初始轨道根数以减少其轨道寿命。

图13为二元组合InitialEpoch-Ω的动力学网格图，其中，图13(a)和图13(b)均表明初始历元时刻对IGSO轨道的长期演化有不可忽视的影响，且影响呈现出周期性，约18.6年，即月球轨道周期。

图7、图8、图9、图10、图11、图12和图13的绘制总共需要递推约280000条IGSO轨道，借助简化后的二次平均半解析轨道递推器，在本实施例中总共耗时约6.5小时，相比于已有一次平均半解析轨道递推器递推40000条IGSO轨道所需的数十小时，表明本发明所述方法具有显著的快速性。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种倾斜地球同步轨道的长期演化快速分析方法，其特征在于：包括如下步骤：

步骤一：推导平均摄动势函数；

考虑到IGSO轨道是一类不易出现奇异的轨道，且从方便应用的角度出发，选取经典轨道根数σ＝(a,e,i,Ω,ω,M)作为变量，其中a为轨道半长轴，e为轨道偏心率，i为轨道倾角，Ω为轨道升交点赤经，ω为轨道近地点幅角，M为轨道平近点角；对任意函数F，记

分别为函数F的一次平均结果和二次平均结果；

针对IGSO轨道的二次平均摄动势函数

表示为：

其中，

为一次平均的地球非球形带谐项摄动势函数，

为地球非球形田谐项摄动1:1共振部分摄动势函数，

为二次平均的太阳第三体引力摄动势函数，

为二次平均的月球第三体引力摄动势函数，

为一次平均的太阳光压摄动势函数；

1)、推导

摄动势函数：

地球非球形带谐项摄动势函数R_Zonal表示为：

其中，μ为地球引力常数，a_e为地球赤道半径，r为地球质心到卫星的矢径的大小，

为卫星的地心纬度，l为整数系数，J_l为第l阶带谐项系数，

为

的第l阶勒让德多项式；

考虑主要带谐项J₂,J₃,J₄，J₂项是一阶摄动因素，是最主要的摄动源；J₃,J₄项各为奇次和偶次带谐项的代表，J₃,J₄项对卫星轨道影响的规律能够反映出其它奇次和偶次带谐项的影响规律；将

带入式(2)中，取l＝2,3,4，并将勒让德多项式展开，得主要带谐项J₂,J₃,J₄的摄动势函数：

根据椭圆运动关系式，时间t与真近点角f之间的微分关系为：

其中，

为轨道平均角速度；

利用式(4)得如下形式的函数，所述函数是关于一个卫星轨道周期的平均值：

其中，k₁,k₂为系数，k₁,k₂＝0,1,2,；

将式(5)带入式(3)得主要带谐项J₂,J₃,J₄关于一个卫星轨道周期的一次平均摄动势函数：

在式(6)中已消去真近点角f；除了主要带谐项J₂,J₃,J₄的一阶项外，还需要考虑J₂的二阶项

所对应的摄动势函数，根据Brouwer的理论，同理得：

一次平均的地球非球形带谐项摄动势函数

表示为：

2)、推导

摄动势函数：

地球非球形田谐项摄动势函数R_Tesseral表示为：

其中，S_lmpq(ω,M,Ω,θ_G),θ_lmpq均为中间变量；F_lmp(i)为倾角函数，仅和轨道倾角i有关，是关于i的三角函数的有限形式；G_lpq(e)为汉森系数，仅和偏心率e有关，是关于e的无限形式，但有些是有限的；C_lm,S_lm为引力位球谐系数；θ_G为格林尼治恒星时；l,m,p,q均为整数系数；

由于IGSO轨道位于地球同步轨道区域，轨道周期等同于地球自转周期，因此地球非球形田谐项摄动对IGSO卫星轨道的主要影响为1:1共振项引起的长周期项；引入变量λ＝ω+Ω+M-θ_G，在式(9)中取J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄对应摄动势函数的1:1共振部分，分别为

和

则地球非球形田谐项摄动1:1共振部分摄动势函数

表示为：

3)、推导

摄动势函数：

太阳引力和月球引力对IGSO卫星轨道的摄动均为第三体引力摄动，其摄动势函数统一进行推导；以地心赤道惯性坐标系为基准坐标系，将太阳、月球、地球均视为质点，给出日月第三体引力摄动势函数为：

其中，μ′为第三体引力常数，r′为地球质心到第三体的矢径的大小，ψ为地球质心到卫星的矢径r和地球质心到第三体的矢径r′的夹角，P_l(cosψ)为cosψ的第l阶勒让德多项式；

第三体在地心赤道惯性坐标系下的轨道根数为(a′,e′,i′,Ω′,ω′,f′)，其中，a′为半长轴，e′为偏心率，i′为倾角，Ω′为升交点轨道赤经，ω′为近地点幅角，f′为真近点角；将cosψ用卫星的轨道根数和第三体的轨道根数表示为：

其中，

分别为卫星轨道近焦点坐标系x轴和y轴的方向向量：

A,B是和卫星真近点角无关的项，用c_x,s_x分别表示变量x的余弦函数cosx和正弦函数sinx，则A,B表示为：

其中，ΔΩ＝Ω-Ω′,u′＝ω′+f′；

将式(19)截取到勒让德多项式的4阶并展开得：

将式(25)表示的日月第三体引力摄动势函数在卫星一个运动周期内做平均，消除卫星真近点角：

其中，AF_2nd,AF_3rd,AF_4th的表达式涉及到如下形式的函数，所述函数为关于一个卫星轨道周期的平均值：

利用式(20)、(26)、(27)得AF_2nd,AF_3rd,AF_4th的表达式为：

将式(28)和式(23)、(24)代入式(26)即得到一次平均后的日月第三体引力摄动势函数

一次平均后的日月第三体引力摄动势函数

已经不再包含与f有关的短周期项，但仍然包含与第三体轨道真近点角f′有关的短周期项，因此将

再对第三体轨道的一个运动周期做一次平均，以消去第三体轨道真近点角f′；

根据式(23)、(24)，辅助变量A,B和卫星真近点角f无关，但包含第三体轨道的真近点角f′；将式(26)表示的一次平均后的日月第三体引力摄动势函数在第三体轨道的一个运动周期做一次平均：

其中，DAF_2nd,DAF_3rd,DAF_4th的表达式涉及到式(5)中所述的函数；

利用式(23)、(24)、(28)、(29)和式(5)得DAF_2nd,DAF_3rd,DAF_4th的表达式：

将式(30)带入式(29)即得二次平均后的日月第三体引力摄动势函数

二次平均后的日月第三体引力摄动势函数

中需要计算第三体在地心赤道惯性坐标系中的轨道根数；对太阳轨道和月球轨道做出如下简化：太阳轨道视为圆轨道，月球轨道视为长期进动椭圆模型并保留黄白交角；

在地心赤道惯性坐标系下，太阳轨道的平均轨道根数为：

其中，AU为天文单位，T为儒略世纪，d为平太阳日，其和儒略世纪满足关系：T＝36525d；

在地心黄道惯性坐标系下，月球轨道的平均轨道根数为：

其中，km为千米单位；

对于月球轨道，式(32)为其相对黄道面的轨道根数，其中(i_M,Ω_M,ω_M)需装换到地球赤道面内，令对应根数为(i′_M,Ω′_M,ω′_M)，则：

其中，ε为平黄赤交角，ζ为一辅助角；

对于太阳引力摄动，已假设太阳轨道为圆轨道，偏心率为0，将式(31)代入式(29)可得二次平均后的太阳引力摄动势函数

为：

对于月球引力摄动，将式(32)、(33)代入式(29)可得二次平均后的月球引力摄动势函数

为：

4)、推导

摄动势函数：

对于IGSO轨道，不考虑地影对太阳光压摄动的影响；

太阳光压摄动势函数R_SRP表示为：

其中，

表示

的高阶小量，μ_SRP为：

其中，C_R为与卫星表面反射系数相关的参数，取值为1到2；S为卫星受光压作用的有效面积；m为卫星质量；P_SRP为地球表面处的太阳光压常数；Δ₀为日地平均距离；

式(36)中所包含的cosψ含义同式(20)，将式(20)代入式(36)中并去掉高阶小量，得：

将式(37)代入式(38)，忽略式(38)中等式右端的第二大项，仅保留第一大项，也即取到勒让德展开式的一阶项，取r′为日地平均距离Δ₀，得太阳光压摄动势函数为：

将式(39)表示的太阳光压摄动势函数在一个卫星轨道周期内做一次平均，利用式(27)，得一次平均后的太阳光压摄动势函数：

步骤二：利用拉格朗日型轨道摄动方程建立IGSO轨道的二次平均半解析轨道递推器；

针对所选用的经典轨道根数σ＝(a,e,i,Ω,ω,M)的拉格朗日型轨道摄动方程为：

根据式(41)，需要求摄动势函数对轨道根数的偏导数；步骤一中已给出IGSO轨道上二次平均后的摄动势函数

摄动势函数

对轨道根数的偏导数表示为：

将式(42)中的偏导数代入式(41)，即得到IGSO轨道二次平均后的半解析递推模型；将所得的二次平均半解析递推模型结合数值积分器，得到同时考虑地球非球形摄动带谐项J₂,J₃,J₄,

地球非球形摄动田谐项J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄的1:1共振部分、日月引力摄动勒让德展开截取到4阶项、太阳光压摄动的IGSO轨道的二次平均半解析轨道递推器；

步骤三：IGSO轨道摄动模型对比分析；

步骤二得到的二次平均半解析轨道递推器中考虑的摄动源和对应的摄动阶数为：1)地球非球形摄动带谐项，可选阶数为J₂,J₃,J₄,

2)地球非球形摄动田谐项，可选阶数为J₂₂,J₃₁,J₃₂,J₃₃,J₄₁,J₄₂,J₄₃,J₄₄；3)太阳引力摄动项，可选阶数为2阶、4阶；4)月球引力摄动项，可选阶数为2阶、3阶、4阶；5)太阳光压摄动项；

针对任意IGSO轨道，在使用二次平均半解析轨道递推器进行长期轨道演化分析时，通过考虑不同的摄动源和对应的摄动阶数，对比分析不同摄动和摄动阶数对IGSO轨道长期演化的影响；根据式(1)、(8)、(18)、(34)、(35)、(40)，不同摄动源和摄动阶数的摄动势函数之间是以加和的形式组成总的摄动势函数

且式(41)表明，各摄动对轨道根数的影响也是加和的形式，因此在步骤二得到的二次平均半解析轨道递推器中，选择是否包含某个特定的摄动源或某一特定的摄动阶数对轨道长期演化的影响；对比方法为仅改变某一特定摄动项，其他条件不变，通过轨道长期演化结果对比分析该摄动项的影响，以此判断在绘制动力学网格图时是否可忽略，以进一步加快计算；

步骤四：绘制一系列完整轨道根数以及初始时刻组合的动力学网格图；

对于IGSO轨道，在绘制动力学网格图时，根据步骤三的分析结论，选取满足精度要求且摄动项考虑最少的摄动模型，以获得最快的计算速度；若卫星面质比大于一数值，则必须考虑太阳光压摄动对轨道长期演化的影响；

在IGSO轨道的长期演化中，初始轨道根数对轨道的长期演化起到决定性的影响；对轨道根数(e,i,Ω,ω)进行划分，其全部二元组合为e-i,ω-Ω；e-ω,i-Ω；e-Ω,i-ω；此外，初始历元时刻也会影响IGSO轨道的长期演化，因此还需要绘制有关初始历元时刻InitialEpoch的动力学网格图，其全部二元组合为InitialEpoch-e,InitialEpoch-i,InitialEpoch-Ω,InitialEpoch-ω；对于每一个划分的二元组合，其余初始根数或初始历元时刻保持不变，绘制两张动力学网格图，指标分别为轨道长期演化中的最大偏心率和轨道寿命；当IGSO轨道在长期演化过程中的偏心率达到设置的临界偏心率时，此时轨道近地点已再入大气层，轨道将逐渐衰减，视为轨道寿命终止；对IGSO轨道长期演化起决定性作用的初始轨道根数(e,i,Ω,ω)以及有较大的影响的初始历元时刻InitialEpoch间一系列二元组合的动力学网格图充分反映了IGSO轨道的长期演化规律，而经过步骤三进一步简化后的二次平均半解析轨道递推器则保证了整个过程的快速性；

步骤五：利用步骤四得到的动力学网格图，完成对IGSO轨道长期演化的快速分析；

步骤四已采用简化后的二次平均半解析轨道递推器快速绘制了有关初始轨道根数和初始历元时刻的动力学网格图，这些动力学网格图充分反映了IGSO轨道长期演化的最大偏心率分布以及轨道寿命分布；对IGSO轨道长期演化的最大偏心率分布进行分析，为IGSO轨道卫星的任务设计提供参考；对IGSO轨道长期演化的轨道寿命分布进行分析，为IGSO轨道卫星的任务后处理提供参考；步骤四中动力学网格图的快速绘制保证了本步骤的快速性，有利于任务设计的快速迭代。

Images