CN115392540A - 一种用于月球轨道交会制导的快速预报方法 - Google Patents

一种用于月球轨道交会制导的快速预报方法 Download PDF

Info

Publication number
CN115392540A
CN115392540A CN202210897104.1A CN202210897104A CN115392540A CN 115392540 A CN115392540 A CN 115392540A CN 202210897104 A CN202210897104 A CN 202210897104A CN 115392540 A CN115392540 A CN 115392540A
Authority
CN
China
Prior art keywords
acceleration
center
gravitational
pseudo
calculating
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.)
Pending
Application number
CN202210897104.1A
Other languages
English (en)
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.)
Beijing Institute of Spacecraft System Engineering
Original Assignee
Beijing Institute of Spacecraft System Engineering
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 Beijing Institute of Spacecraft System Engineering filed Critical Beijing Institute of Spacecraft System Engineering
Priority to CN202210897104.1A priority Critical patent/CN115392540A/zh
Publication of CN115392540A publication Critical patent/CN115392540A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0631—Resource planning, allocation, distributing or scheduling for enterprises or organisations

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • General Physics & Mathematics (AREA)
  • Economics (AREA)
  • Strategic Management (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Data Mining & Analysis (AREA)
  • Development Economics (AREA)
  • Quality & Reliability (AREA)
  • General Business, Economics & Management (AREA)
  • Operations Research (AREA)
  • Marketing (AREA)
  • Tourism & Hospitality (AREA)
  • Computational Mathematics (AREA)
  • Remote Sensing (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Game Theory and Decision Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Educational Administration (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computing Systems (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供了一种用于月球轨道交会制导的快速预报方法，能够实现高精度、快速预报，满足任务的要求。本发明的预报方法，将用于地球轨道卫星GPS自主定轨的引力场加速度计算方法应用到月球轨道交会领域，针对该算法无法直接应用到月球轨道交会制导策略设计所存在的问题，逐一进行改进和升级，适用于月球轨道交会飞控任务实施需要的高精度、快速算法以满足任务的要求。具体地，本发明首次将该算法移植到月球轨道应用领域，这是前人没有涉及的；针对月球引力场的特殊性，各引力摄动加速度量级相当，无法单独用J2项进行近似描述的问题，在算法中增加了J22项摄动加速度的考虑，从而使得算法的计算精度和数值稳定性得到了大幅度提升。

Description

一种用于月球轨道交会制导的快速预报方法
技术领域
本发明涉及深空探测轨道设计技术领域，具体涉及一种用于月球轨道交会制导的快速预报方法。
背景技术
100km到200km左右的近圆环月轨道是探月工程任务中最常用的轨道类型，其受到的最主要的摄动力是月球的非球形引力场的摄动。引力摄动加速度计算模型通常采用球谐函数表示，每次加速度计算都要计算一组Legendre多项式递推公式。递推的计算量是随引力场阶数增加呈几何级数增长的。由于月球引力场的不均匀性，没有像地球或者火星引力场类似J2主项，为了确保计算的精度，通常需要引入更高阶的引力场模型进行计算。
工程经验表明，对于轨道设计而言，通常需要保留32x32的阶数才能保证轨道设计结果的精度；而对于轨道确定和飞控变轨控制策略实施而言，至少需要100x100阶的引力场阶数才能满足任务要求。而引力场阶数的增加会造成计算量呈几何级数增长，造成计算时长也相应大幅度增加。在月球轨道交会对接远程导引变轨策略飞控实施阶段，需要在短时间内根据实时测定轨结果完成轨控策略的在线设计、生成和注入工作。由于控制变量多、积分时间长、引力场模型精度要求高，传统数值积分的方法由于计算速度慢，很难适应应急条件下高精度、快速策略重构设计的要求。
发明内容
有鉴于此，本发明提供了一种用于月球轨道交会制导的快速预报方法，能够实现高精度、快速预报，满足任务的要求。
为实现上述目的，本发明技术方案如下：
本发明的一种用于月球轨道交会制导的快速预报方法，首先，将引力场加速度刨除中心引力场和J2项、J22项残余引力场引力加速度后的残余加速度折算成等效的“伪中心”参数，从而将残余的引力场加速度通过非线性映射成“伪中心”；然后，在给定的经度和维度点和指定的高度范围内，通过多项式插值，计算在该高度范围内“伪中心”关于高度的函数的插值多项式系数；再将月球表面按经度和维度分割成若干网格点，计算并保存每个经纬度点在给定高度范围内的“伪中心”函数的拟合多项式系数；最后，对于任意空间位置点，通过计算其附近6个经纬度网格点在该位置点瞬时高度对应的“伪中心”通过二维插值获得该点对应的伪中心，并通过非线性映射折算出该点对应的真实引力场加速度；计算加速度矢量对位置矢量的偏导数矩阵；基于真实引力场加速度以及加速度矢量对位置矢量的偏导数矩阵，计算得到探测器的位置和速度。
其中，包括如下步骤：
步骤1，给定探测器初始位置速度；计算月球中心引力场二体问题的加速度及其偏导数；
步骤2，计算相应的J2项和J22项摄动加速度和偏导数；
步骤3，计算去掉J2项和J22项的伪中心；
步骤4，将三维空间分解为经度、纬度和高度三个方向，对伪中心进行离散化，获得高度方向的伪中心多项式拟合；
步骤5，根据给定的探测器位置，得到高度h、经度和纬度的离散网格点；对高度h、经度和纬度的二维离散网格点，进行二维六点插值，得到对应网格点的伪中心插值结果；
步骤6，根据插值结果，计算得到对应的引力加速度，并加入J2项和J22项的影响，最终得到真实加速度；
步骤8，根据计算得到的状态转移矩阵偏导数，获得任意时刻的位置、速度和状态转移矩阵，至此完成轨道预报问题的求解。
其中，所述步骤1中，中心引力场二体问题的加速度a2B表示为：
Figure BDA0003769338000000031
其中，μ为地球引力场常数，其中x、y、z分别为探测器位置矢量的三方向分量r＝(xy z)，r为位置矢量模值，
Figure BDA0003769338000000032
a2B对位置矢量的偏导数为：
Figure BDA0003769338000000033
其中，所述步骤2中，J2项的摄动加速度为：
Figure BDA0003769338000000034
其中，J2为二阶带谐项系数，Re为地球赤道半径；
J2项的加速度偏导数为：
Figure BDA0003769338000000041
Figure BDA0003769338000000042
Figure BDA0003769338000000043
Figure BDA0003769338000000044
J22项的摄动加速度为
Figure BDA0003769338000000045
其中，J22为二阶二次田谐系数，λ22为二阶田谐项主轴的地理经度，即赤道椭圆常州的方位角，均为常数；Λ22为中间计算量，Λ22＝xcosλ22+ysinλ22；
J22摄动加速度对位置矢量的偏导数为
Figure BDA0003769338000000051
Figure BDA0003769338000000052
Figure BDA0003769338000000053
Figure BDA0003769338000000054
其中，所述步骤3中，月心固连坐标系下三维空间任意一点位置r和引力加速度a0与伪中心c0的映射为
Figure BDA0003769338000000061
其中r为探测器实际的月心距，表示伪月心距ρ加上伪中心c；μ为地球引力场常数；a0为去掉
Figure BDA0003769338000000062
和
Figure BDA0003769338000000063
后的摄动引力加速度，满足下式：
Figure BDA0003769338000000064
其中，aCBF是采用非球形引力场加速度函数计算的探测器相对月心固联坐标系的引力加速度，包含中心引力场加速度和摄动加速度，CBF表示中心天体固联坐标系；所述伪中心c与摄动加速度aCBF的关系为
Figure BDA0003769338000000065
有益效果：
1、本发明的预报方法，将用于地球轨道卫星GPS自主定轨的引力场加速度计算方法应用到月球轨道交会领域，针对该算法无法直接应用到月球轨道交会制导策略设计所存在的问题，逐一进行改进和升级，适用于月球轨道交会飞控任务实施需要的高精度、快速算法以满足任务的要求。具体地，本发明首次将该算法移植到月球轨道应用领域，这是前人没有涉及的；针对月球引力场的特殊性，各引力摄动加速度量级相当，无法单独用J2项进行近似描述的问题，在算法中增加了J22项摄动加速度的考虑，从而使得算法的计算精度和数值稳定性得到了大幅度提升。
2、本发明针对原算法无法计算敏度矩阵从而无法用于制导策略迭代修正的问题，在原始算法基础上增加了加速度偏导数的插值计算方法，从而可以满足月球轨道交会制导律设计的需要。
3、本发明方法针对月球轨道交会对接远程导引任务的需求，将在传统地球引力场计算引力加速度的GAAF算法首次引入月球引力场计算，并针对月球引力场与地球引力场的差异，以及月球轨道导引策略计算的需求对算法进行了改进和升级，在不损失精度的情况下可以有效提升环月轨道预报速度达2个数量，满足在轨实时任务规划的要求，特别适用于应急条件下轨控策略快速重构的要求。
附图说明
图1为本发明实施例方法流程图。
图2为本发明二维六点插值方法示意图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
本发明提供了一种用于月球轨道交会制导的快速预报方法，针对月球轨道交会对接远程导引任务的需求，将在传统地球引力场计算引力加速度的GAAF算法首次引入月球引力场计算，并针对月球引力场与地球引力场的差异，以及月球轨道导引策略计算的需求对算法进行了改进和升级，通过将J22项引力场系数引入GAAF算法以及推导了GAAF引力加速度的偏导数差值公式用于快速计算状态转移矩阵，实现探测器位置和速度的高精度快速预报。本发明具体流程是：
首先，将引力场加速度刨除中心引力场和J2项、J22项残余引力场引力加速度后的残余加速度折算成等效的“伪中心”参数，从而将残余的引力场加速度通过非线性映射成“伪中心”；
然后，在给定的经度和维度点和指定的高度范围内，通过多项式插值，计算在该高度范围内“伪中心”关于高度的函数的插值多项式系数；
再将月球表面按经度和维度分割成若干网格点，计算并保存每个经纬度点在给定高度范围内的“伪中心”函数的拟合多项式系数；
最后，对于任意空间位置点，通过计算其附近6个经纬度网格点在该位置点瞬时高度对应的“伪中心”通过二维插值获得该点对应的伪中心，并通过非线性映射折算出该点对应的真实引力场加速度；计算加速度矢量对位置矢量的偏导数矩阵；基于真实引力场加速度以及加速度矢量对位置矢量的偏导数矩阵，计算得到探测器的位置和速度。
具体地，本发明方法原理为：
给定航天器在t0时刻的位置和速度(r0,v0)以及航天器加速度函数a3×1(r,v,t)后，航天器环月轨道预报问题定义为如下42维一阶常微分方程组的初值问题：
Figure BDA0003769338000000081
其中(r,v)为航天器在t时刻的位置和速度，Φ6×6(t,t0)为t0时刻至t时刻的状态转移矩阵。
方程组的初值为：
Figure BDA0003769338000000082
环月轨道航天器受到的加速度主要包括月球引力场加速度aCBF，地球三体摄动加速度aEarth3B，太阳三体摄动加速度aSun3B以及太阳光压加速度aSRP，即：
a3×1(r,v,t)＝aCBF+aEarth3B+aSun3B+aSRP (3)
对于低轨环月近圆轨道航天器，月球引力场加速度aCBF及其偏导数
Figure BDA0003769338000000091
(引力场加速度仅与位置r矢量相关，因此
Figure BDA0003769338000000092
涉及计算高阶勒让德多项式，是影响计算速度的主要因素。本发明给出引力场加速度aCBF以及相应的偏导数矩阵
Figure BDA0003769338000000093
快速插值计算方法。针对月球引力场的特殊性，本文将其摄动加速度的两个主项：J2项(二阶带谐项)引力摄动加速度
Figure BDA0003769338000000094
以及J22项(二阶二次田谐项)引力加速度
Figure BDA0003769338000000095
以及相应的加速度偏导数
Figure BDA0003769338000000096
和
Figure BDA0003769338000000097
采用解析表达方法给出，剩余的其它加速度a0：
Figure BDA0003769338000000098
将其映射为伪中心后，在高度方向和经纬度三个方向进行离散化处理。其中，高度方向离散化采用多项式拟合系数存储；经纬度方向采用等间隔经纬度离散点存储。在计算真实加速度的时候，首先通过解析方法计算出两个主项加速度
Figure BDA0003769338000000099
和
Figure BDA00037693380000000910
然后通过插值离散数据反演方法获得剩余加速度a0，将两者相加后，就可以计算出总的引力摄动加速度aCBF。相应的加速度偏导数
Figure BDA00037693380000000911
也采用相同的思路进行处理，这里不再赘述。将计算出的aCBF和
Figure BDA00037693380000000912
带入式(1)中，在初值(2)条件下，就可以进行轨道预报工作了。轨道预报的计算速度主要是由右函数中加速度的计算速度决定的。由于a0部分采用了数值插值方法进行处理，大幅地提高了计算速度，从而大幅地提高了轨道预报的速度。而且，由于插值计算保留了所有的高阶摄动加速度的特性，因此算法不会因为插值带来的速度提升而降低计算精度。实际测试表明，与100x100阶引力场的传统轨道预报模型相比，轨道预报2天，本算法轨道预报位置误差小于100m，速度误差小于0.01m/s，状态转移矩阵误差小于0.1％，计算速度提高近两个数量级。
本实施例的轨道预报方法流程图如图1所示，包括如下步骤：
步骤1，给定探测器初始位置速度
Figure BDA00037693380000000913
计算月球中心引力场二体问题的加速度及其偏导数。其中，中心引力场二体问题的加速度a2B可以表示为：
Figure BDA0003769338000000101
其中，μ为地球引力场常数，其中x、y、z分别为探测器位置矢量的三方向分量r＝(xy z)，r为位置矢量模值，
Figure BDA0003769338000000102
a2B对位置矢量的偏导数为：
Figure BDA0003769338000000103
步骤2，计算相应的J2项和J22项摄动加速度和偏导数；其中，J2项的摄动加速度为：
Figure BDA0003769338000000104
其中，J2为二阶带谐项系数，Re为地球赤道半径。
J2项的加速度偏导数为：
Figure BDA0003769338000000105
Figure BDA0003769338000000111
Figure BDA0003769338000000112
Figure BDA0003769338000000113
J22项的摄动加速度为
Figure BDA0003769338000000114
其中，J22为二阶二次田谐系数，λ22为二阶田谐项主轴的地理经度，即赤道椭圆常州的方位角，均为常数；Λ22为中间计算量，Λ22＝xcosλ22+ysinλ22。
J22摄动加速度对位置矢量的偏导数为
Figure BDA0003769338000000121
Figure BDA0003769338000000122
Figure BDA0003769338000000123
Figure BDA0003769338000000124
步骤3，计算去掉J2和J22项的伪中心；
其中，引力加速度可以表达为：
Figure BDA0003769338000000131
其中，aCBF是采用非球形引力场加速度函数计算的探测器相对月心固联坐标系的引力加速度，包含中心引力场加速度和摄动加速度。CBF表示中心天体固联坐标系。ρ定义为探测器相对引力加速度aCBF的伪月心距离(Pseudo－Radius)，ρ＝|ρ|。
探测器实际的月心距r定义为伪月心距加上伪中心(Pseudo－Center)c：
c＝r-ρ (18)
伪中心c与摄动加速度aCBF的关系为：
Figure BDA0003769338000000132
计算去掉
Figure BDA0003769338000000133
和
Figure BDA0003769338000000134
后的摄动引力加速度a0为
Figure BDA0003769338000000135
进一步可得到月心固连坐标系下三维空间任意一点位置r和引力加速度a0与伪中心c0(c0是与a0对应的伪中心)的映射为
Figure BDA0003769338000000136
步骤4，将三维空间分解为经度、纬度和高度三个方向，对伪中心进行离散化，获得高度方向的伪中心多项式拟合，具体过程为：
在经度和纬度方向，将球面离散为若干网格点，对于给定的网格点坐标
Figure BDA0003769338000000138
其高度方向的伪中心c0＝[cx cy cz]采用有理多项式进行拟合：
Figure BDA0003769338000000137
i＝x,y,z代表伪中心矢量的x,y,z三个方向，Dk为分母多项式，Nk为分子多项式，a1…an-1、b1…bn-1均为多项式系数。拟合方法采用通用的Pade近似函数逼近。其中，
Figure BDA0003769338000000141
Hmax和Hmin分别为拟合的高度区间，对于低轨近圆环月轨道，通常取Hmin＝10km，Hmax＝350km可以满足大多数任务的要求。不同m和n的取值，会影响拟合精度和计算速度。通过大量的计算和测试，取n＝7和m＝6可以在存储量、计算精度与计算速度之间取得较好的权衡。
步骤5，根据给定的探测器位置，得到高度h、经度和纬度的离散网格点；对高度h、经度和纬度的二维离散网格点，进行二维六点插值，得到对应网格点的伪中心插值结果。采用双变量六点插值方法对任意点
Figure BDA0003769338000000142
的插值算法如下：
Figure BDA0003769338000000143
其中，p和q为经度和纬度的归一化参数：
Figure BDA0003769338000000144
λ0为任意点
Figure BDA0003769338000000145
对应的最近左下角离散网格点的经度，
Figure BDA0003769338000000146
为任意点
Figure BDA0003769338000000147
对应的最近左下角离散网格点的纬度，Δλ和
Figure BDA0003769338000000148
为经度和纬度离散网格点宽度，通常取为Δλ＝0.5°，
Figure BDA0003769338000000149
c(h)参见公式(22)。
c0,0、c0,1、c1,0、c1,1、c0,-1和c-1,0为c0同高度上的经纬度方向六个离散点，详见图2，具体表达式为：
Figure BDA0003769338000000151
步骤6，根据插值结果，计算得到对应的引力加速度，并加入J2和J22项的影响，最终得到真实加速度，具体过程为：
根据插值结果c0，计算对应的伪中心距ρ0：
ρ0＝r-c0 (27)
再计算对应的剩余加速度a0：
Figure BDA0003769338000000152
加入J2和J22项的影响，最终得到引力加速度
Figure BDA0003769338000000153
步骤7，计算加速度矢量对位置矢量的偏导数矩阵，具体过程如下：
首先计算伪中心c对h的偏导数如下：
Figure BDA0003769338000000154
由于不同插值点所在的高度区间的拟合函数都不同，因此，对于6个二维插值点，均可以计算相应的伪中心对x的偏导数
Figure BDA0003769338000000155
Figure BDA0003769338000000156
其中，k＝x,y,z代表伪中心矢量的x,y,z三个方向。
可以计算得到
Figure BDA0003769338000000157
Figure BDA0003769338000000161
对于经纬度方向的6个插值点，可以分别插值得到
Figure BDA0003769338000000162
总加速度的偏导数的计算公式如下：
Figure BDA0003769338000000163
其中，
Figure BDA0003769338000000164
上式中，
Figure BDA0003769338000000165
Figure BDA0003769338000000166
Figure BDA0003769338000000171
Figure BDA0003769338000000172
参照式(6)计算；
Figure BDA0003769338000000173
按式(34)计算；
Figure BDA0003769338000000174
按式计算；
Figure BDA0003769338000000175
按式(13)计算。
步骤8，根据计算得到的状态转移矩阵偏导数
Figure BDA0003769338000000176
带入公式(1)，以公式(2)为初值，求解相应的常微分方程组初值问题，就可以获得任意t时刻的位置r(t)、速度v(t)和状态转移矩阵Φ(t,t0)，至此完成轨道预报问题的求解。
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种用于月球轨道交会制导的快速预报方法，其特征在于，首先，将引力场加速度刨除中心引力场和J2项、J22项残余引力场引力加速度后的残余加速度折算成等效的“伪中心”参数，从而将残余的引力场加速度通过非线性映射成“伪中心”；然后，在给定的经度和维度点和指定的高度范围内，通过多项式插值，计算在该高度范围内“伪中心”关于高度的函数的插值多项式系数；再将月球表面按经度和维度分割成若干网格点，计算并保存每个经纬度点在给定高度范围内的“伪中心”函数的拟合多项式系数；最后，对于任意空间位置点，通过计算其附近6个经纬度网格点在该位置点瞬时高度对应的“伪中心”通过二维插值获得该点对应的伪中心，并通过非线性映射折算出该点对应的真实引力场加速度；计算加速度矢量对位置矢量的偏导数矩阵；基于真实引力场加速度以及加速度矢量对位置矢量的偏导数矩阵，计算得到探测器的位置和速度。
2.如权利要求1所述的预报方法，其特征在于，包括如下步骤：
步骤1，给定探测器初始位置速度；计算月球中心引力场二体问题的加速度及其偏导数；
步骤2，计算相应的J2项和J22项摄动加速度和偏导数；
步骤3，计算去掉J2项和J22项的伪中心；
步骤4，将三维空间分解为经度、纬度和高度三个方向，对伪中心进行离散化，获得高度方向的伪中心多项式拟合；
步骤5，根据给定的探测器位置，得到高度h、经度和纬度的离散网格点；对高度h、经度和纬度的二维离散网格点，进行二维六点插值，得到对应网格点的伪中心插值结果；
步骤6，根据插值结果，计算得到对应的引力加速度，并加入J2项和J22项的影响，最终得到真实加速度；
步骤8，根据计算得到的状态转移矩阵偏导数，获得任意时刻的位置、速度和状态转移矩阵，至此完成轨道预报问题的求解。
3.如权利要求2所述的预报方法，其特征在于，所述步骤1中，中心引力场二体问题的加速度a2B表示为：
Figure FDA0003769337990000021
其中，μ为地球引力场常数，其中x、y、z分别为探测器位置矢量的三方向分量r＝(x yz)，r为位置矢量模值，
Figure FDA0003769337990000022
a2B对位置矢量的偏导数为：
Figure FDA0003769337990000023
4.如权利要求2或3所述的预报方法，其特征在于，所述步骤2中，J2项的摄动加速度为：
Figure FDA0003769337990000024
其中，J2为二阶带谐项系数，Re为地球赤道半径；
J2项的加速度偏导数为：
Figure FDA0003769337990000031
Figure FDA0003769337990000032
Figure FDA0003769337990000033
Figure FDA0003769337990000034
J22项的摄动加速度为
Figure FDA0003769337990000035
其中，J22为二阶二次田谐系数，λ22为二阶田谐项主轴的地理经度，即赤道椭圆常州的方位角，均为常数；Λ22为中间计算量，Λ22＝xcosλ22+ysinλ22；
J22摄动加速度对位置矢量的偏导数为
Figure FDA0003769337990000041
Figure FDA0003769337990000042
Figure FDA0003769337990000043
Figure FDA0003769337990000044
5.如权利要求4所述的预报方法，其特征在于，所述步骤3中，月心固连坐标系下三维空间任意一点位置r和引力加速度a0与伪中心c0的映射为
Figure FDA0003769337990000051
其中r为探测器实际的月心距，表示伪月心距ρ加上伪中心c；μ为地球引力场常数；a0为去掉
Figure FDA0003769337990000052
和
Figure FDA0003769337990000053
后的摄动引力加速度，满足下式：
Figure FDA0003769337990000054
其中，aCBF是采用非球形引力场加速度函数计算的探测器相对月心固联坐标系的引力加速度，包含中心引力场加速度和摄动加速度，CBF表示中心天体固联坐标系；所述伪中心c与摄动加速度aCBF的关系为
Figure FDA0003769337990000055
CN202210897104.1A 2022-07-28 2022-07-28 一种用于月球轨道交会制导的快速预报方法 Pending CN115392540A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210897104.1A CN115392540A (zh) 2022-07-28 2022-07-28 一种用于月球轨道交会制导的快速预报方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210897104.1A CN115392540A (zh) 2022-07-28 2022-07-28 一种用于月球轨道交会制导的快速预报方法

Publications (1)

Publication Number Publication Date
CN115392540A true CN115392540A (zh) 2022-11-25

Family

ID=84116129

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210897104.1A Pending CN115392540A (zh) 2022-07-28 2022-07-28 一种用于月球轨道交会制导的快速预报方法

Country Status (1)

Country Link
CN (1) CN115392540A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116701823A (zh) * 2023-08-07 2023-09-05 长沙翔宇信息科技有限公司 交会点空间范围估算方法、装置、终端设备及存储介质

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116701823A (zh) * 2023-08-07 2023-09-05 长沙翔宇信息科技有限公司 交会点空间范围估算方法、装置、终端设备及存储介质
CN116701823B (zh) * 2023-08-07 2023-10-27 长沙翔宇信息科技有限公司 交会点空间范围估算方法、装置、终端设备及存储介质

Similar Documents

Publication Publication Date Title
CN112257343B (zh) 一种高精度地面轨迹重复轨道优化方法及系统
CN106697333B (zh) 一种航天器轨道控制策略的鲁棒性分析方法
CN102878995B (zh) 一种静止轨道卫星自主导航方法
CN107421550A (zh) 一种基于星间测距的地球‑Lagrange联合星座自主定轨方法
CN104792340A (zh) 一种星敏感器安装误差矩阵与导航系统星地联合标定与校正的方法
CN109656133B (zh) 一种针对空间走廊跟踪观测的分布式卫星群优化设计方法
Wolf et al. Performance trades for Mars pinpoint landing
CN105737858A (zh) 一种机载惯导系统姿态参数校准方法与装置
CN112161632B (zh) 一种基于相对位置矢量测量的卫星编队初始定位方法
CN110096726B (zh) 基于月球借力的geo卫星应急转移轨道快速优化设计方法
CN111301715A (zh) 基于霍曼变轨的同轨道特定相位分布的星座布局与轨道调整方法、装置及计算机存储介质
CN109708663B (zh) 基于空天飞机sins辅助的星敏感器在线标定方法
CN103853887A (zh) 一种冻结轨道的偏心率的卫星轨道确定方法
CN110779532A (zh) 一种应用于近地轨道卫星的地磁导航系统及方法
CN103968834A (zh) 一种近地停泊轨道上深空探测器的自主天文导航方法
Qin et al. An innovative navigation scheme of powered descent phase for Mars pinpoint landing
CN115392540A (zh) 一种用于月球轨道交会制导的快速预报方法
CN108959665B (zh) 适用于低轨卫星的轨道预报误差经验模型生成方法及系统
CN111814313B (zh) 一种高精度引力场中回归轨道设计方法
Nordkvist et al. Attitude feedback tracking with optimal attitude state estimation
CN116125503A (zh) 一种高精度卫星轨道确定及预报算法
D'Souza et al. Orion cislunar guidance and navigation
CN114200491A (zh) 一种基于导航数据的应急航天器星历确定方法及系统
Li et al. Navigation and guidance of orbital transfer vehicle
CN110579784B (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