CN111854765B - 一种中轨道导航卫星轨道长期预报方法 - Google Patents

一种中轨道导航卫星轨道长期预报方法 Download PDF

Info

Publication number
CN111854765B
CN111854765B CN202010864123.5A CN202010864123A CN111854765B CN 111854765 B CN111854765 B CN 111854765B CN 202010864123 A CN202010864123 A CN 202010864123A CN 111854765 B CN111854765 B CN 111854765B
Authority
CN
China
Prior art keywords
orbit
average
perturbation
change rate
conservative
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010864123.5A
Other languages
English (en)
Other versions
CN111854765A (zh
Inventor
胡敏
徐家辉
张锐
王许煜
李玖阳
马聪慧
云朝明
潘升东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Publication of CN111854765A publication Critical patent/CN111854765A/zh
Application granted granted Critical
Publication of CN111854765B publication Critical patent/CN111854765B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/02—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/10—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration
    • G01C21/12—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by using measurements of speed or acceleration executed aboard the object being navigated; Dead reckoning

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种中轨道导航卫星轨道长期预报方法，包括获取保守摄动力下平均轨道根数的一阶变化率，保守摄动力包括地球非球形摄动力和日月三体引力摄动力；获取非保守摄动力下平均轨道根数的一阶变化率，非保守摄动力包括太阳光压摄动力；根据保守摄动力下平均轨道根数的一阶变化率和非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数。本发明通过对保守摄动力和非保守摄动力分别求解，更加能反映卫星轨道长期演化的特性；通过平均化的方法分离出短周期项，积分步长可以设置为数个周期，在兼顾长期预报精度的情况下，有效提高计算效率。

Description

一种中轨道导航卫星轨道长期预报方法
技术领域
本申请涉及一种中轨道导航卫星轨道长期预报方法，属于卫星轨道预报技术领域。
背景技术
轨道区域可按轨道高度划分如下：低地球轨道(Low Earth Orbit，LEO)区域，空间区域的高度范围为0到2000km；中地球轨道(Medium Earth Orbit，MEO)区域，空间区域的高度范围为2000km到35586km；地球同步轨道(Geosynchronous Orbit，GEO)区域，空间区域的高度范围为35586km到35986km(GEO高度35786km上下200km的区域)。中轨道区域即中地球轨道区域是一个很大的空间区域，是导航卫星的主要运行区域。
卫星轨道预报根据时间尺度可分为长期和短期预报，轨道长期预报模型经历了从早期公式推导获得解析解到小步长数值法积分外推，再到半分析法积分外推兼顾二者优点几个过程。刘林提出了一种卫星轨道预报的分析方法，用平均根数表示t时刻卫星的位置和速度，并直接用卫星直角坐标的位置和速度分量表示地球非球形摄动的周期项，可以有效避免小偏心率和轨道倾角而出现的奇点问题。周建华提出了一种基于神经网络混合建模的中长期轨道预报方法，在原动力学模型基础上引入神经元网络模型作为补偿，以提高预报精度。王大为对不同摄动力进行半分析法处理，并将完整摄动力模型的半分析方法应用于不同类型轨道，以研究在不同初始轨道根数半长轴和轨道倾角的组合下，半分析法在100年内轨道预报的精度。张斌斌以碎片所在的轨道高度和面质比对碎片进行分类，以碎片的划分单元为研究对象，通过平均法去除短周期项，在保证轨道预报进度的情况下显著提高计算效率。现有的轨道长期预报模型没有专门针对中轨道卫星的长期预报方法，同时考虑的摄动力太多，计算效率较低。
发明内容
本申请的目的在于，提供一种中轨道导航卫星轨道长期预报方法，以解决现有轨道长期预报方法存在的考虑的摄动力太多，计算效率较低技术问题。
本发明的中轨道导航卫星轨道长期预报方法，包括：
获取保守摄动力下平均轨道根数的一阶变化率，所述保守摄动力包括地球非球形摄动力和日月三体引力摄动力；
获取非保守摄动力下平均轨道根数的一阶变化率，所述非保守摄动力包括太阳光压摄动力；
根据所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数。
优选地，所述获取保守摄动力下平均轨道根数的一阶变化率，具体为：
获取轨道根数变化率随轨道根数变化的参数化运动方程；
获取所述参数化运动方程的一阶解，基于所述参数化运动方程的一阶解，计算保守摄动力下平均轨道根数的一阶变化率。
优选地，所述获取轨道根数变化率随轨道根数变化的参数化运动方程，具体为：
所述参数化运动方程根据第一公式确定；所述第一公式为：
Figure BDA0002649165020000021
式中，ai为第i个平均轨道根数，
Figure BDA0002649165020000022
为第i个平均轨道根数的导数、n为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)、
Figure BDA0002649165020000023
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、
Figure BDA0002649165020000024
r为卫星地心距矢量、
Figure BDA0002649165020000025
为由所述保守摄动力得到的摄动函数。
优选地，所述获取所述参数化运动方程的一阶解，具体为：
所述参数化运动方法的一阶解根据第二公式确定；所述第二公式为：
Figure BDA0002649165020000031
式中，ai为第i个平均轨道根数、t为时间、n(a)为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)、＜Fi(a,t)＞为摄动函数，＜·＞为平均算子，其定义为摄动函数对平经度λ一个周期内求平均，所述摄动函数是根据第三公式确定；所述第三公式为：
Figure BDA0002649165020000032
优选地，所述基于所述参数化运动方程的一阶解，计算保守摄动力下平均轨道根数的一阶变化率，具体为：
基于所述参数化运动方程的一阶解，计算地球非球形带谐项摄动力下平均轨道根数的一阶变化率；
基于所述参数化运动方程的一阶解，计算日月三体引力摄动力下平均轨道根数的一阶变化率；
基于所述参数化运动方程的一阶解，计算地球非球形田谐共振项摄动力下平均轨道根数的一阶变化率。
优选地，所述基于所述参数化运动方程的一阶解，计算地球非球形带谐项摄动力下平均轨道根数的一阶变化率，具体为：
根据第四公式确定地球非球形带谐项的平均势函数；所述第四公式为：
Figure BDA0002649165020000033
式中，μ为地球引力常数、a为导航卫星轨道半长轴、N表示带谐项的最高阶数，δ0s为(1,0,0,0,0,0)、R为地球平均赤道半径，Jn为重力场带谐项系数，Vns,
Figure BDA0002649165020000034
Qns,Gs为系数；
对所述地球非球形带谐项的平均势函数
Figure BDA0002649165020000035
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到地球非球形带谐项摄动力下平均轨道根数的一阶变化率。
优选地，所述基于所述参数化运动方程的一阶解，计算日月三体引力摄动力下平均轨道根数的一阶变化率，具体为：
根据第五公式确定日月三体引力摄动力下的一阶平均摄动势函数；所述第五公式为：
Figure BDA0002649165020000041
式中，μ3为第三体引力常数、R3为第三体到地球的距离、N表示第三体考虑的最高阶数、δ0s值为(1,0,0,0,0,0)、a为导航卫星轨道半长轴、Vns，
Figure BDA0002649165020000042
Qns，Gs为系数；
对所述日月三体引力摄动力下的一阶平均摄动势函数
Figure BDA0002649165020000043
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到第三体引力摄动力下平均轨道根数的一阶变化率。
优选地，所述基于所述参数化运动方程的一阶解，计算地球非球形田谐共振项摄动力下平均轨道根数的一阶变化率，具体为：
根据第六公式确定地球非球形田谐共振项摄动力下的平均摄动函数；所述第六公式为：
Figure BDA0002649165020000044
式中，i为虚数、Re(·)为选取函数的实数部分、a为平均轨道根数的6维向量、λ为平经度、θ为地球自转角、m,j为整数、U(a,θ,t)为关于a,θ,t的摄动函数表达式、
Figure BDA0002649165020000051
对所述地球非球形田谐共振项摄动力下的平均摄动函数
Figure BDA0002649165020000052
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到地球非球形田谐共振项摄动力下的平均轨道根数的一阶变化率。
优选地，所述获取非保守摄动力下平均轨道根数的一阶变化率，具体为：
根据第七公式确定
Figure BDA0002649165020000053
所述第七公式为：
Figure BDA0002649165020000054
式中，
Figure BDA0002649165020000055
为第i个平均轨道根数的导数、ai为第i个平均轨道根数、
Figure BDA0002649165020000056
r为卫星地心距、a为半长轴、
Figure BDA0002649165020000057
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、L为真经度。
优选地，所述根据所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数，具体为：
将所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率叠加；
基于叠加的结果，利用数值积分Adams-Cowell方法对导航卫星的轨道状态外推更新，得到中轨道导航卫星的平均轨道根数。
本发明的中轨道导航卫星轨道长期预报方法相较于现有技术，具有如下有益效果：
通过对保守力摄动力和非保守摄动力分别求解，更加能反映卫星轨道长期演化的特性。
通过平均化的方法分离出参数化运动方程中的短周期项，积分步长可以设置为数个周期，在兼顾长期预报精度的情况下，有效提高计算效率。
附图说明
图1为本发明实施例中轨道导航卫星轨道长期预报方法的流程图；
图2为本发明实施例中轨道导航卫星轨道长期预报方法中摄动加速度量级随轨道高度变化示意图；
图3为本发明实施例中轨道导航卫星轨道长期预报方法中与分点根数对应的分点坐标系。
具体实施方式
下面结合实施例详述本发明，但本发明并不局限于这些实施例。
图1为本发明的中轨道导航卫星轨道长期预报方法的流程图。
本发明的中轨道导航卫星轨道长期预报方法，包括：
步骤1、获取保守摄动力下平均轨道根数的一阶变化率，保守摄动力包括地球非球形摄动力和日月三体引力摄动力，具体为：
步骤1.1、获取轨道根数变化率随轨道根数变化的参数化运动方程，参数化运动方程根据第一公式确定；第一公式为：
Figure BDA0002649165020000061
式中，ai为第i个平均轨道根数、
Figure BDA0002649165020000062
为第i个平均轨道根数的导数、n为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)、
Figure BDA0002649165020000063
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、
Figure BDA0002649165020000064
r为卫星地心距矢量、
Figure BDA0002649165020000065
为由所述保守摄动力得到的摄动函数。
步骤1.2、获取参数化运动方程的一阶解，基于参数化运动方程的一阶解，计算保守摄动力下平均轨道根数的一阶变化率，具体为：
步骤1.2.1、参数化运动方法的一阶解根据第二公式确定；第二公式为：
Figure BDA0002649165020000071
式中，ai为第i个平均轨道根数、t为时间、n(a)为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)、＜Fi(a,t)＞为摄动函数，＜·＞为平均算子，其定义为摄动函数对平经度λ一个周期内求平均，所述摄动函数是根据第三公式确定；所述第三公式为：
Figure BDA0002649165020000072
Figure BDA0002649165020000073
步骤1.2.2、基于参数化运动方程的一阶解，计算地球非球形带谐项摄动力下平均轨道根数的一阶变化率，具体为：
根据第四公式确定地球非球形带谐项的平均势函数；第四公式为：
Figure BDA0002649165020000074
式中，μ为地球引力常数、a为导航卫星轨道半长轴、N表示带谐项的最高阶数，δ0s值为(1,0,0,0,0,0)、R为地球平均赤道半径，Jn为重力场带谐项系数，Vns,
Figure BDA0002649165020000075
Qns,Gs为系数；
对地球非球形带谐项的平均势函数
Figure BDA0002649165020000076
求导，并代入到拉格朗日形式的参数化运动方程中，得到地球非球形带谐项摄动力下平均轨道根数的一阶变化率。
步骤1.2.3、基于参数化运动方程的一阶解，计算日月三体引力摄动力下平均轨道根数的一阶变化率，具体为：
根据第五公式确定日月三体引力摄动力下的一阶平均摄动势函数；第五公式为：
Figure BDA0002649165020000081
式中，μ3为第三体引力常数、R3为第三体到地球的距离、N表示第三体考虑的最高阶数、δ0s值为(1,0,0,0,0,0)、a为导航卫星轨道半长轴、Vns，
Figure BDA0002649165020000082
Qns，Gs为系数；
对日月三体引力摄动力下的一阶平均摄动势函数
Figure BDA0002649165020000083
求导，并代入到拉格朗日形式的参数化运动方程中，得到第三体引力摄动力下平均轨道根数的一阶变化率。
步骤1.2.4、基于参数化运动方程的一阶解，计算地球非球形田谐共振项摄动力下平均轨道根数的一阶变化率，具体为：
根据第六公式确定地球非球形田谐共振项摄动力下的平均摄动函数；第六公式为：
Figure BDA0002649165020000084
式中，i为虚数、Re(·)为选取函数的实数部分、a为平均轨道根数的6维向量、λ为平经度、θ为地球自转角、m,j为整数、U(a,θ,t)为关于a,θ,t的摄动函数表达式、
Figure BDA0002649165020000085
对地球非球形田谐共振项摄动力下的平均摄动函数
Figure BDA0002649165020000086
求导，并代入到拉格朗日形式的参数化运动方程中，得到地球非球形田谐共振项摄动力下的平均轨道根数的一阶变化率。
步骤2、获取非保守摄动力下平均轨道根数的一阶变化率，非保守摄动力包括太阳光压摄动力，具体为：
根据第七公式确定
Figure BDA0002649165020000095
第七公式为：
Figure BDA0002649165020000091
式中，
Figure BDA0002649165020000092
为第i个平均轨道根数的导数、ai为第i个平均轨道根数、
Figure BDA0002649165020000093
r为卫星地心距、a为半长轴、
Figure BDA0002649165020000094
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、L为真经度。
步骤3、根据保守摄动力下平均轨道根数的一阶变化率和非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数，具体为：
将保守摄动力下平均轨道根数的一阶变化率和非保守摄动力下平均轨道根数的一阶变化率叠加；
基于叠加的结果，利用数值积分Adams-Cowell方法对导航卫星的轨道状态外推更新，得到中轨道导航卫星的平均轨道根数。
下面将以具体的实施例说明本发明的中轨道导航卫星轨道长期预报方法。
本发明实施例的中轨道导航卫星轨道长期预报的方法，包括：
步骤S1：确定中轨道区域的摄动力，摄动力包括保守摄动力和非保守摄动力。
中轨道区域的轨道动力学与低轨道区域的轨道动力学有很大不同，主要是由于不需要考虑大气阻力。在近地空间的这一区域，其它扰动(例如地球非球形、三体引力和太阳光压)主导并确定了轨道的演化。废弃卫星在轨运行中受到的作用力，决定了其长期演化的运动状态。图2给出了航天器受到的摄动力加速度量级(加速度取log)随轨道高度的变化。
从图2中可以看出，从低轨道到同步轨道区域，地球非球形摄动力都是最主要的摄动力；随着轨道高度的增加，J2项摄动加速度开始衰减，而日、月三体引力摄动加速度则增加，太阳光压变化不大，但对轨道演化的影响不可忽略。中轨道导航卫星轨道高度在19000km到24000km之间，因此，考虑的摄动力为地球非球形摄动、日月三体引力摄动和太阳光压摄动。导航卫星在轨过程还受到其它微小摄动力的作用，但其加速度的量级太小，对中轨道导航卫轨道长期预报时可以忽略；
asd＝a0+ans+aS+aM+asr (1)
式中，a0表示地球中心引力项；ans表示地球非球形摄动项；aS、aM分别表示太阳、月球引力摄动项；asr表示太阳光压摄动项。通过上述可以确定，对中轨道卫星产生影响的摄动力包括保守摄动力和非保守摄动力，其中保守摄动力为地球非球形摄动力和太阳、月球三体引力摄动力。
步骤S2：利用分点根数描述的参数运动方程；
采用开普勒根数描述轨道状态的运动方程时，在偏心率和轨道倾角趋于0时无法求解，而选取分点根数描述航天器的轨道状态则可以有效避免运动方程出现奇异的情况。分点根数记为(a,h,k,p,q,λ)，与开普勒轨道根数的转换关系如下：
Figure BDA0002649165020000101
与分点根数对应，定义分点坐标系
Figure BDA0002649165020000102
顺行轨道的分点坐标系三轴的方向矢量定义如图3所示，OE-XYZ为J2000坐标系。
定义方向余弦(α,β,γ)为zB轴和分点坐标系
Figure BDA0002649165020000103
轴夹角的余弦值：
Figure BDA0002649165020000104
Figure BDA0002649165020000105
Figure BDA0002649165020000106
在惯性坐标系下航天器的运动方程为：
Figure BDA0002649165020000111
式中，
Figure BDA0002649165020000112
是一种类似势能的函数，是由保守摄动力得到的摄动函数，fnon为航天器受到的非保守摄动力加速度。
运动方程(4)是由航天器的位置和速度描述的，将其转化为轨道根数变化率随轨道根数变化的参数化运动方程。用(a1,a2…a6)分别表示分点根数(a,h,k,p,q,λ)，参数化运动方程用分点根数描述的形式如下：
Figure BDA0002649165020000113
式中：n为卫星运行的平均角速度；δi6值为(0,0,0,0,0,1)。
(ai,aj)的运算规则为：
Figure BDA0002649165020000114
参数化运动方程的摄动函数
Figure BDA0002649165020000115
是关于(a,h,k,λ,α,β,γ)的函数，可直接求导得出，而对于p和q，需要通过链式求导获得，具体为将摄动函数
Figure BDA0002649165020000116
表示为关于轨道根数(a,h,k,λ)和方向余弦(α,β,λ)的函数。为了方便描述，引入中间变量A、B和C，分别定义为：
Figure BDA0002649165020000117
Figure BDA0002649165020000118
C＝1+p2+q2 (7)
引入交叉微分算子U,αβ，定义为：
Figure BDA0002649165020000119
然后，通过应用链式规则获得摄动函数
Figure BDA00026491650200001110
对于p和q的偏导数：
Figure BDA0002649165020000121
Figure BDA0002649165020000122
步骤S3：求解保守摄动力作用下平均根数的一阶变化率；
步骤S3.1：构造运动方程的一阶解；
为了通过平均法分离短周期项，将密切根数在平均根数的领域内展开：
Figure BDA0002649165020000123
式中，
Figure BDA0002649165020000124
表示第i个平均轨道密切根数，ai表示第i个平均轨道根数；
Figure BDA0002649165020000125
表示一个很小的短周期变化量，
Figure BDA0002649165020000126
表示第i个轨道根数的第j阶短周期变化项。
在参数运动方程(5)中引入参数ε，ε取值范围为[0,1]，参数运动方程可改写为：
Figure BDA0002649165020000127
εFi项给出了摄动作用下密切根数的变化率，对于平均根数的参数运动方程，我们假定以下形式：
Figure BDA0002649165020000128
结合密切根数运动方程(10)和平均根数运动方程(12)，并将公式(11)中密切根数函数在平均分点根数的邻域内展开：
Figure BDA0002649165020000131
其中：
Figure BDA0002649165020000132
Figure BDA0002649165020000133
类似的，我们将平均运动角速率函数
Figure BDA0002649165020000134
在平均分点根数的邻域内展开：
Figure BDA0002649165020000135
其中：
Figure BDA0002649165020000136
Figure BDA0002649165020000137
Figure BDA0002649165020000138
在获得上述展开形式的运动方程后，接下来就能求解平均运动方程。首先，将公式(10)对时间t微分，通过方程(12)获得密切根数变化率的表达式：
Figure BDA0002649165020000141
然后，将公式(13)到(16)代入公式(11)，获得密切根数变化率的另一个表达式：
Figure BDA0002649165020000142
任意参数ε属于[0,1]，密切根数变化率的两个表达式都相等，因此，对于j＝0,1,2,...，εj项的系数也应该相等，分别得到j阶的方程如下：
Figure BDA0002649165020000143
Figure BDA0002649165020000144
Figure BDA0002649165020000145
通过对公式(19)求平均可以得到平均分点根数的变化率
Figure BDA0002649165020000146
同时，因为短周期项
Figure BDA0002649165020000147
和N1(a)是关于快变量λ的周期性函数，对公式(19)的第一行求平均，得到一阶项的平均分点根数的变化率：
Figure BDA0002649165020000148
＜·＞是平均算子，其定义为摄动函数对平经度λ一个周期内求平均：
Figure BDA0002649165020000149
一阶平均运动方程的形式可写为：
Figure BDA00026491650200001410
步骤S3.2：地球非球形带谐项摄动下平均根数的一阶变化率
对于地球非球形带谐项，其平均势函数为：
Figure BDA0002649165020000151
式中：μ为地球引力常数、a为航天器轨道半长轴，N表示带谐项的最高阶数，R为地球平均赤道半径，Jn为重力场带谐项系数，δ0s值为(1,0,0,0,0,0)，Vns,
Figure BDA0002649165020000152
Qns,Gs为系数。
对于系数Vns，如果n-s是奇数，则Vns＝0；若n-s是偶数，则有：
Figure BDA0002649165020000153
迭代公式(24)的初始启动项为：
V0,0＝1
Figure BDA0002649165020000154
Figure BDA0002649165020000155
为Hansen系数，迭代公式为：
Figure BDA0002649165020000156
式中：
Figure BDA0002649165020000157
系数Qns＝Qns(γ)，可通过下述迭代公式进行计算：
Figure BDA0002649165020000161
系数Qns的初始项为Q0,0＝1，多项式Gs的迭代计算公式为：
Gs＝(kα+hβ)Gs-1-(hα-kβ)Hs-1,G0＝1
Hs＝(hα-kβ)Gs-1+(kα+hβ)Hs-1,H0＝0 (28)
将平均势函数
Figure BDA0002649165020000162
求导，并代入到拉格朗日形式的运动方程，得到地球非球形带谐项摄动作用下的一阶平均运动方程：
Figure BDA0002649165020000163
Figure BDA0002649165020000164
Figure BDA0002649165020000165
Figure BDA0002649165020000166
Figure BDA0002649165020000167
Figure BDA0002649165020000168
步骤S3.3：太阳/月球引力作用下平均根数的一阶变化率
对于三体引力摄动，重新定义方向余弦(α,β,γ)为第三体质心相对于地球质心的位置矢量R3和分点坐标系
Figure BDA0002649165020000169
轴夹角的余弦值：
Figure BDA0002649165020000171
Figure BDA0002649165020000172
Figure BDA0002649165020000173
采用平均算子对第三体引力摄动函数平均化，得到其一阶平均摄动势函数如下：
Figure BDA0002649165020000174
式中：μ3为第三体引力常数，R3为第三体到地球的距离，N表示第三体考虑的最高阶数、δ0s值为(1,0,0,0,0,0)、a为导航卫星轨道半长轴、Vns，
Figure BDA0002649165020000175
Qns，Gs为系数。
系数Vns、Qns和Gs利用公式(24)、(27)和(28)计算获得，系数
Figure BDA0002649165020000176
的迭代公式为：
Figure BDA0002649165020000177
初始项为：
Figure BDA0002649165020000178
对第三体引力的平均摄动势函数
Figure BDA0002649165020000179
求导，并代入到拉格朗日形式的运动方程，得到第三体引力摄动作用下的一阶平均运动方程：
Figure BDA0002649165020000181
Figure BDA0002649165020000182
Figure BDA0002649165020000183
Figure BDA0002649165020000184
Figure BDA0002649165020000185
Figure BDA0002649165020000186
步骤S3.4：地球非球形田谐共振项作用下平均根数的一阶变化率
定义地球自转角为θ，航天器赤经为αB，则航天器的地理经度为赤经与地球自转角的差值。
Figure BDA0002649165020000187
田谐项势函数可用地球转转角θ来描述，其平均摄动函数的形式如下：
Figure BDA0002649165020000188
式中：i为虚数，Re(·)为选取函数的实数部分，a为平均轨道根数的6维向量、λ为平经度、θ为地球自转角、m,j为整数、U(a,θ,t)为关于a,θ,t的摄动函数表达式、
Figure BDA0002649165020000189
Figure BDA00026491650200001810
为航天器运行角速度、
Figure BDA00026491650200001811
为地球自转角速度。
集合B为满足公式(36)的正整数集，当航天器运行的轨道周期与地球自转周期的比值为简单的整数比(或着偏差很小)，就会发生田谐项共振，此时田谐项对轨道的长期摄动影响不能忽略。
Figure BDA0002649165020000191
对田谐项平均摄动势函数
Figure BDA0002649165020000192
求导，并代入到拉格朗日形式的运动方程，得到其一阶平均运动方程：
Figure BDA0002649165020000193
Figure BDA0002649165020000194
Figure BDA0002649165020000195
Figure BDA0002649165020000196
Figure BDA0002649165020000197
Figure BDA0002649165020000198
步骤S4：求解非保守摄动力作用下平均根数的一阶变化率
太阳光压摄动力为非保守摄动力，不能通过平均算子获得解析的平均摄动函数，根据运动方程(5)，太阳光压摄动作用下平均根数变化率通过对平经度λ积分的方法求解。
Figure BDA0002649165020000199
分点根数对速度的偏导数和摄动力fnon是航天器位置矢量的函数，用真近点角f描述更加方便。引入真经度L和偏经度F两个辅助经度，它们与开普勒偏真近地点角f和近地点角E的关系如下：
L＝f+ω+IΩ
F＝E+ω+IΩ (39)
为了求解偏经度F，需要先求解关于偏经度F和平经度λ的开普勒方程(40)，再通过公式(41)求解真经度L。
λ＝F+h cos F-k sin F (40)
Figure BDA0002649165020000201
式中：b＝(1+B)-1。
根据公式(39)到(41)平经度λ和四个经度之间的关系，将公式(38)转化为关于真经度L的积分：
Figure BDA0002649165020000202
其中，地心距r也可以用真经度L表示：
Figure BDA0002649165020000203
因此，分点根数和摄动力fnon都是关于真经度L的函数，再利用数值积分方法对公式(42)的平均根数变化率进行求解。
步骤S5：对各摄动的一阶解叠加，采用数值积分方法外推；
对步骤S3、S4得到的不同摄动作用下平均根数的一阶变化率叠加，再利用数值积分Adams-Cowell方法对航天器轨道状态外推更新，可以计算出航天器预报时刻的平均轨道根数。
本发明提供的一种中轨道导航卫星轨道长期预报方法具有以下优点：
(1)通过平均化的方法分离出短周期项，积分步长可以设置为数个周期，在兼顾长期预报精度的情况下，有效提高计算效率。
(2)通过对保守力和非保守力分别求解，更加能反映卫星轨道长期演化的特性。
以上所述，仅是本申请的几个实施例，并非对本申请做任何形式的限制，虽然本申请以较佳实施例揭示如上，然而并非用以限制本申请，任何熟悉本专业的技术人员，在不脱离本申请技术方案的范围内，利用上述揭示的技术内容做出些许的变动或修饰均等同于等效实施案例，均属于技术方案范围内。

Claims (8)

1.一种中轨道导航卫星轨道长期预报方法，其特征在于，包括：
获取保守摄动力下平均轨道根数的一阶变化率，所述保守摄动力包括地球非球形摄动力和日月三体引力摄动力；
获取非保守摄动力下平均轨道根数的一阶变化率，所述非保守摄动力包括太阳光压摄动力；
根据所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数；所述获取保守摄动力下平均轨道根数的一阶变化率，具体为：
获取轨道根数变化率随轨道根数变化的参数化运动方程；
获取所述参数化运动方程的一阶解，基于所述参数化运动方程的一阶解，计算保守摄动力下平均轨道根数的一阶变化率；
所述获取非保守摄动力下平均轨道根数的一阶变化率，具体为：
根据第七公式确定
Figure FDA0003548294310000011
所述第七公式为：
Figure FDA0003548294310000012
式中，
Figure FDA0003548294310000013
为第i个平均轨道根数的导数、ai为第i个平均轨道根数、
Figure FDA0003548294310000014
r为卫星地心距、a为导航卫星轨道 半长轴、
Figure FDA0003548294310000015
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、L为真经度。
2.根据权利要求1所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述获取轨道根数变化率随轨道根数变化的参数化运动方程，具体为：
所述参数化运动方程根据第一公式确定；所述第一公式为：
Figure FDA0003548294310000021
式中，ai为第i个平均轨道根数、
Figure FDA0003548294310000022
为第i个平均轨道根数的导数、
Figure FDA0003548294310000023
为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)、
Figure FDA0003548294310000028
为卫星地心距矢量的导数、fnon为导航卫星受到的所述非保守摄动力的加速度、
Figure FDA0003548294310000024
r为卫星地心距矢量、
Figure FDA0003548294310000025
为由所述保守摄动力得到的摄动函数。
3.根据权利要求2所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述获取所述参数化运动方程的一阶解，具体为：
所述参数化运动方程的一阶解根据第二公式确定；所述第二公式为：
Figure FDA0003548294310000026
式中，ai为第i个平均轨道根数、t为时间、n(a)为导航卫星运行的平均角速度、δi6值为(0,0,0,0,0,1)，<Fi(a,t)>为摄动函数，<·>为平均算子，其定义为摄动函数对平经度λ一个周期内求平均，所述摄动函数是根据第三公式确定；所述第三公式为：
Figure FDA0003548294310000027
4.根据权利要求3所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述基于所述参数化运动方程的一阶解，计算保守摄动力下平均轨道根数的一阶变化率，具体为：
基于所述参数化运动方程的一阶解，计算地球非球形带谐项摄动力下平均轨道根数的一阶变化率；
基于所述参数化运动方程的一阶解，计算日月三体引力摄动力下平均轨道根数的一阶变化率；
基于所述参数化运动方程的一阶解，计算地球非球形田谐共振项摄动力下平均轨道根数的一阶变化率。
5.根据权利要求4所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述基于所述参数化运动方程的一阶解，计算地球非球形带谐项摄动力下平均轨道根数的一阶变化率，具体为：
根据第四公式确定地球非球形带谐项的平均势函数；所述第四公式为：
Figure FDA0003548294310000031
式中，μ为地球引力常数、a为导航卫星轨道半长轴、N表示带谐项的最高阶数，δ0s值为(1,0,0,0,0,0)、R为地球平均赤道半径，Jn为重力场带谐项系数，
Figure FDA0003548294310000032
Qns,Gs为系数；
对所述地球非球形带谐项的平均势函数
Figure FDA0003548294310000033
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到地球非球形带谐项摄动力下平均轨道根数的一阶变化率。
6.根据权利要求4所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述基于所述参数化运动方程的一阶解，计算日月三体引力摄动力下平均轨道根数的一阶变化率，具体为：
根据第五公式确定日月三体引力摄动力下的一阶平均摄动势函数；所述第五公式为：
Figure FDA0003548294310000041
式中，μ3为第三体引力常数、R3为第三体到地球的距离、N1表示第三体考虑的最高阶数、δ0s值为(1,0,0,0,0,0)、a为导航卫星轨道半长轴、Vns，
Figure FDA0003548294310000042
Qns，Gs为系数；
对所述日月三体引力摄动力下的一阶平均摄动势函数
Figure FDA0003548294310000043
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到第三体引力摄动力下平均轨道根数的一阶变化率。
7.根据权利要求4所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述基于所述参数化运动方程的一阶解，计算地球非球形田谐共振项摄动力下平均轨道根数的一阶变化率，具体为：
根据第六公式确定地球非球形田谐共振项摄动力下的平均摄动函数；所述第六公式为：
Figure FDA0003548294310000044
式中，i为虚数、Re(·)为选取函数的实数部分、a为平均轨道根数的6维向量、λ为平经度、θ为地球自转角、m,j为整数、U(a,θ,t)为关于a,θ,t的摄动函数表达式、
Figure FDA0003548294310000051
对所述地球非球形田谐共振项摄动力下的平均摄动函数
Figure FDA0003548294310000052
求导，并代入到拉格朗日形式的所述参数化运动方程中，得到地球非球形田谐共振项摄动力下的平均轨道根数的一阶变化率。
8.根据权利要求1～7任一项所述的中轨道导航卫星轨道长期预报方法，其特征在于，所述根据所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率计算导航卫星的平均轨道根数，具体为：
将所述保守摄动力下平均轨道根数的一阶变化率和所述非保守摄动力下平均轨道根数的一阶变化率叠加；
基于叠加的结果，利用数值积分Adams-Cowell方法对导航卫星的轨道状态外推更新，得到中轨道导航卫星的平均轨道根数。
CN202010864123.5A 2020-06-08 2020-08-25 一种中轨道导航卫星轨道长期预报方法 Active CN111854765B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN2020105151107 2020-06-08
CN202010515110 2020-06-08

Publications (2)

Publication Number Publication Date
CN111854765A CN111854765A (zh) 2020-10-30
CN111854765B true CN111854765B (zh) 2022-04-26

Family

ID=72968173

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010864123.5A Active CN111854765B (zh) 2020-06-08 2020-08-25 一种中轨道导航卫星轨道长期预报方法

Country Status (1)

Country Link
CN (1) CN111854765B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114357788B (zh) * 2022-01-10 2023-08-01 中国空间技术研究院 低轨巨型星座偏差演化分析方法及装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010078476A (ja) * 2008-09-26 2010-04-08 Seiko Epson Corp 長期予測軌道データの信頼性判定方法、測位方法、測位装置及び測位システム
CN103499349A (zh) * 2013-09-29 2014-01-08 桂林电子科技大学 基于广播星历参数外推的卫星轨道中长期预报方法和系统
CN108959665A (zh) * 2017-05-17 2018-12-07 上海微小卫星工程中心 适用于低轨卫星的轨道预报误差经验模型生成方法及系统
CN109063380A (zh) * 2018-09-12 2018-12-21 北京理工大学 一种静止轨道电推进卫星故障检测方法及位置保持方法
CN109613574A (zh) * 2018-11-13 2019-04-12 中国人民解放军战略支援部队航天工程大学 计算北斗中轨卫星坟墓轨道最早穿越其他全球卫星导航系统轨道时间的方法
CN110132261A (zh) * 2018-11-16 2019-08-16 中国西安卫星测控中心 一种基于数值拟合的高精度星上轨道预报方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109255096B (zh) * 2018-07-25 2022-10-04 西北工业大学 一种基于微分代数的地球同步卫星轨道不确定演化方法
CN111125874B (zh) * 2019-11-18 2023-08-22 中国人民解放军63686部队 一种动平台高精度测轨预报方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010078476A (ja) * 2008-09-26 2010-04-08 Seiko Epson Corp 長期予測軌道データの信頼性判定方法、測位方法、測位装置及び測位システム
CN103499349A (zh) * 2013-09-29 2014-01-08 桂林电子科技大学 基于广播星历参数外推的卫星轨道中长期预报方法和系统
CN108959665A (zh) * 2017-05-17 2018-12-07 上海微小卫星工程中心 适用于低轨卫星的轨道预报误差经验模型生成方法及系统
CN109063380A (zh) * 2018-09-12 2018-12-21 北京理工大学 一种静止轨道电推进卫星故障检测方法及位置保持方法
CN109613574A (zh) * 2018-11-13 2019-04-12 中国人民解放军战略支援部队航天工程大学 计算北斗中轨卫星坟墓轨道最早穿越其他全球卫星导航系统轨道时间的方法
CN110132261A (zh) * 2018-11-16 2019-08-16 中国西安卫星测控中心 一种基于数值拟合的高精度星上轨道预报方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
"Long-term evolution of orbits about a precessing oblate planet. 2. The case of variable precession";Michael Efroimsky 等;《Celestial Mechanics and Dynamical Astronomy》;20061231;第96卷(第3期);正文第259-286页 *
"The effect of secular resonances on the long-term orbital evolution of uncontrollable objects on satellite radio navigation systems in the MEO region";Bordovitsyna 等;《Solar System Research》;20121231;第46卷(第5期);正文第329-339页 *
地球非球形对卫星轨道的长期影响及补偿研究;项军华等;《飞行力学》;20070615(第02期);正文第85-87页 *
平均轨道根数与密切轨道根数的互换;胡敏等;《飞行器测控学报》;20120415(第02期);正文第77-80页 *

Also Published As

Publication number Publication date
CN111854765A (zh) 2020-10-30

Similar Documents

Publication Publication Date Title
CN109255096B (zh) 一种基于微分代数的地球同步卫星轨道不确定演化方法
Van Patten et al. A possible experiment with two counter-orbiting drag-free satellites to obtain a new test of einstein's general theory of relativity and improved measurements in geodesy
CN109032176B (zh) 一种基于微分代数的地球同步轨道确定和参数确定方法
CN110595485A (zh) 基于两行根数的低轨卫星长期轨道预报方法
CN109738919B (zh) 一种用于gps接收机自主预测星历的方法
CN107421550B (zh) 一种基于星间测距的地球-Lagrange联合星座自主定轨方法
CN106679674B (zh) 基于星历模型的地月L2点Halo轨道阴影分析方法
Ciufolini et al. Measurement of dragging of inertial frames and gravitomagnetic field using laser-ranged satellites
Lorell et al. Mariner 9 celestial mechanics experiment: Gravity field and pole direction of Mars
Ning et al. A novel differential Doppler measurement-aided autonomous celestial navigation method for spacecraft during approach phase
CN110816896B (zh) 一种卫星星上简易轨道外推方法
CN102230969A (zh) 一种卫星星座星间链路的长时间自主维持方法
CN113589832B (zh) 对地表固定区域目标稳定观测覆盖的星座快速设计方法
CN111854765B (zh) 一种中轨道导航卫星轨道长期预报方法
Sengupta et al. Satellite orbit design and maintenance for terrestrial coverage
Flores et al. A method for accurate and efficient propagation of satellite orbits: A case study for a Molniya orbit
Patel et al. Controlled short-period orbits around Earth-Moon equilateral libration points for Lunar Occultations
Golikov THEONA—a numerical-analytical theory of motion of artificial satellites of celestial bodies
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
Cappuccio et al. BepiColombo gravity and rotation experiment in a pseudo drag-free system
CN112393835B (zh) 一种基于扩展卡尔曼滤波的小卫星在轨推力标定方法
CN114386282A (zh) 半分析法的低轨巨型星座轨道动力学分析方法及装置
Arbinger et al. Impact of orbit prediction accuracy on low earth remote sensing flight dynamics operations
CN111547274A (zh) 一种航天器高精度自主目标预报方法
CN112036037A (zh) 一种倾斜地球同步轨道的长期演化快速分析方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant