一种用于月球轨道交会制导的快速预报方法与流程

1.本发明涉及深空探测轨道设计技术领域，具体涉及一种用于月球轨道交会制导的快速预报方法。


背景技术：

2.100km到200km左右的近圆环月轨道是探月工程任务中最常用的轨道类型，其受到的最主要的摄动力是月球的非球形引力场的摄动。引力摄动加速度计算模型通常采用球谐函数表示，每次加速度计算都要计算一组legendre多项式递推公式。递推的计算量是随引力场阶数增加呈几何级数增长的。由于月球引力场的不均匀性，没有像地球或者火星引力场类似j2主项，为了确保计算的精度，通常需要引入更高阶的引力场模型进行计算。
3.工程经验表明，对于轨道设计而言，通常需要保留32x32的阶数才能保证轨道设计结果的精度；而对于轨道确定和飞控变轨控制策略实施而言，至少需要100x100阶的引力场阶数才能满足任务要求。而引力场阶数的增加会造成计算量呈几何级数增长，造成计算时长也相应大幅度增加。在月球轨道交会对接远程导引变轨策略飞控实施阶段，需要在短时间内根据实时测定轨结果完成轨控策略的在线设计、生成和注入工作。由于控制变量多、积分时间长、引力场模型精度要求高，传统数值积分的方法由于计算速度慢，很难适应应急条件下高精度、快速策略重构设计的要求。


技术实现要素：

4.有鉴于此，本发明提供了一种用于月球轨道交会制导的快速预报方法，能够实现高精度、快速预报，满足任务的要求。
5.为实现上述目的，本发明技术方案如下：
6.本发明的一种用于月球轨道交会制导的快速预报方法，首先，将引力场加速度刨除中心引力场和j2项、j
22
项残余引力场引力加速度后的残余加速度折算成等效的“伪中心”参数，从而将残余的引力场加速度通过非线性映射成“伪中心”；然后，在给定的经度和维度点和指定的高度范围内，通过多项式插值，计算在该高度范围内“伪中心”关于高度的函数的插值多项式系数；再将月球表面按经度和维度分割成若干网格点，计算并保存每个经纬度点在给定高度范围内的“伪中心”函数的拟合多项式系数；最后，对于任意空间位置点，通过计算其附近6个经纬度网格点在该位置点瞬时高度对应的“伪中心”通过二维插值获得该点对应的伪中心，并通过非线性映射折算出该点对应的真实引力场加速度；计算加速度矢量对位置矢量的偏导数矩阵；基于真实引力场加速度以及加速度矢量对位置矢量的偏导数矩阵，计算得到探测器的位置和速度。
7.其中，包括如下步骤：
8.步骤1，给定探测器初始位置速度；计算月球中心引力场二体问题的加速度及其偏导数；
9.步骤2，计算相应的j2项和j
22
项摄动加速度和偏导数；
10.步骤3，计算去掉j2项和j
22
项的伪中心；
11.步骤4，将三维空间分解为经度、纬度和高度三个方向，对伪中心进行离散化，获得高度方向的伪中心多项式拟合；
12.步骤5，根据给定的探测器位置，得到高度h、经度和纬度的离散网格点；对高度h、经度和纬度的二维离散网格点，进行二维六点插值，得到对应网格点的伪中心插值结果；
13.步骤6，根据插值结果，计算得到对应的引力加速度，并加入j2项和j
22
项的影响，最终得到真实加速度；
14.步骤8，根据计算得到的状态转移矩阵偏导数，获得任意时刻的位置、速度和状态转移矩阵，至此完成轨道预报问题的求解。
15.其中，所述步骤1中，中心引力场二体问题的加速度a
2b
表示为：
[0016][0017]
其中，μ为地球引力场常数，其中x、y、z分别为探测器位置矢量的三方向分量r＝(x y z)，r为位置矢量模值，
[0018]a2b
对位置矢量的偏导数为：
[0019][0020]
其中，所述步骤2中，j2项的摄动加速度为：
[0021][0022]
其中，j2为二阶带谐项系数，re为地球赤道半径；
[0023]
j2项的加速度偏导数为：
[0024][0025][0026][0027][0028]j22
项的摄动加速度为
[0029][0030]
其中，j
22
为二阶二次田谐系数，λ
22
为二阶田谐项主轴的地理经度，即赤道椭圆常州的方位角，均为常数；λ
22
为中间计算量，λ
22
＝xcosλ
22
+ysinλ
22
；
[0031]j22
摄动加速度对位置矢量的偏导数为
[0032][0033][0034][0035][0036]
其中，所述步骤3中，月心固连坐标系下三维空间任意一点位置r和引力加速度a0与伪中心c0的映射为其中r为探测器实际的月心距，表示伪月心距ρ加上伪中心c；μ为地球引力场常数；a0为去掉和后的摄动引力加速度，满足下式：
[0037][0038]
其中，a
cbf
是采用非球形引力场加速度函数计算的探测器相对月心固联坐标系的引力加速度，包含中心引力场加速度和摄动加速度，cbf表示中心天体固联坐标系；所述伪中心c与摄动加速度a
cbf
的关系为
[0039]
有益效果：
[0040]
1、本发明的预报方法，将用于地球轨道卫星gps自主定轨的引力场加速度计算方法应用到月球轨道交会领域，针对该算法无法直接应用到月球轨道交会制导策略设计所存在的问题，逐一进行改进和升级，适用于月球轨道交会飞控任务实施需要的高精度、快速算法以满足任务的要求。具体地，本发明首次将该算法移植到月球轨道应用领域，这是前人没有涉及的；针对月球引力场的特殊性，各引力摄动加速度量级相当，无法单独用j2项进行近似描述的问题，在算法中增加了j22项摄动加速度的考虑，从而使得算法的计算精度和数值稳定性得到了大幅度提升。
[0041]
2、本发明针对原算法无法计算敏度矩阵从而无法用于制导策略迭代修正的问题，在原始算法基础上增加了加速度偏导数的插值计算方法，从而可以满足月球轨道交会制导律设计的需要。
[0042]
3、本发明方法针对月球轨道交会对接远程导引任务的需求，将在传统地球引力场计算引力加速度的gaaf算法首次引入月球引力场计算，并针对月球引力场与地球引力场的差异，以及月球轨道导引策略计算的需求对算法进行了改进和升级，在不损失精度的情况下可以有效提升环月轨道预报速度达2个数量，满足在轨实时任务规划的要求，特别适用于应急条件下轨控策略快速重构的要求。
附图说明
[0043]
图1为本发明实施例方法流程图。
[0044]
图2为本发明二维六点插值方法示意图。
具体实施方式
[0045]
下面结合附图并举实施例，对本发明进行详细描述。
[0046]
本发明提供了一种用于月球轨道交会制导的快速预报方法，针对月球轨道交会对接远程导引任务的需求，将在传统地球引力场计算引力加速度的gaaf算法首次引入月球引力场计算，并针对月球引力场与地球引力场的差异，以及月球轨道导引策略计算的需求对算法进行了改进和升级，通过将j
22
项引力场系数引入gaaf算法以及推导了gaaf引力加速度的偏导数差值公式用于快速计算状态转移矩阵，实现探测器位置和速度的高精度快速预报。本发明具体流程是：
[0047]
首先，将引力场加速度刨除中心引力场和j2项、j
22
项残余引力场引力加速度后的残余加速度折算成等效的“伪中心”参数，从而将残余的引力场加速度通过非线性映射成“伪中心”；
[0048]
然后，在给定的经度和维度点和指定的高度范围内，通过多项式插值，计算在该高度范围内“伪中心”关于高度的函数的插值多项式系数；
[0049]
再将月球表面按经度和维度分割成若干网格点，计算并保存每个经纬度点在给定高度范围内的“伪中心”函数的拟合多项式系数；
[0050]
最后，对于任意空间位置点，通过计算其附近6个经纬度网格点在该位置点瞬时高度对应的“伪中心”通过二维插值获得该点对应的伪中心，并通过非线性映射折算出该点对应的真实引力场加速度；计算加速度矢量对位置矢量的偏导数矩阵；基于真实引力场加速度以及加速度矢量对位置矢量的偏导数矩阵，计算得到探测器的位置和速度。
[0051]
具体地，本发明方法原理为：
[0052]
给定航天器在t0时刻的位置和速度(r0,v0)以及航天器加速度函数a3×1(r,v,t)后，航天器环月轨道预报问题定义为如下42维一阶常微分方程组的初值问题：
[0053][0054]
其中(r,v)为航天器在t时刻的位置和速度，φ6×6(t,t0)为t0时刻至t时刻的状态转移矩阵。
[0055]
方程组的初值为：
[0056][0057]
环月轨道航天器受到的加速度主要包括月球引力场加速度a
cbf
，地球三体摄动加速度a
earth3b
，太阳三体摄动加速度a
sun3b
以及太阳光压加速度a
srp
，即：
[0058]
a3×1(r,v,t)＝a
cbf
+a
earth3b
+a
sun3b
+a
srp
ꢀꢀꢀ
(3)
[0059]
对于低轨环月近圆轨道航天器，月球引力场加速度a
cbf
及其偏导数(引力场加速度仅与位置r矢量相关，因此涉及计算高阶勒让德多项式，是影响计算速度的主要因素。本发明给出引力场加速度a
cbf
以及相应的偏导数矩阵快速插值计算方法。针对月球引力场的特殊性，本文将其摄动加速度的两个主项：j2项(二阶带谐项)引力摄动加速度以及j
22
项(二阶二次田谐项)引力加速度以及相应的加速度偏导数和采用解析表达方法给出，剩余的其它加速度a0：
[0060][0061]
将其映射为伪中心后，在高度方向和经纬度三个方向进行离散化处理。其中，高度方向离散化采用多项式拟合系数存储；经纬度方向采用等间隔经纬度离散点存储。在计算真实加速度的时候，首先通过解析方法计算出两个主项加速度和然后通过插值离
散数据反演方法获得剩余加速度a0，将两者相加后，就可以计算出总的引力摄动加速度a
cbf
。相应的加速度偏导数也采用相同的思路进行处理，这里不再赘述。将计算出的a
cbf
和带入式(1)中，在初值(2)条件下，就可以进行轨道预报工作了。轨道预报的计算速度主要是由右函数中加速度的计算速度决定的。由于a0部分采用了数值插值方法进行处理，大幅地提高了计算速度，从而大幅地提高了轨道预报的速度。而且，由于插值计算保留了所有的高阶摄动加速度的特性，因此算法不会因为插值带来的速度提升而降低计算精度。实际测试表明，与100x100阶引力场的传统轨道预报模型相比，轨道预报2天，本算法轨道预报位置误差小于100m，速度误差小于0.01m/s，状态转移矩阵误差小于0.1％，计算速度提高近两个数量级。
[0062]
本实施例的轨道预报方法流程图如图1所示，包括如下步骤：
[0063]
步骤1，给定探测器初始位置速度计算月球中心引力场二体问题的加速度及其偏导数。其中，中心引力场二体问题的加速度a
2b
可以表示为：
[0064][0065]
其中，μ为地球引力场常数，其中x、y、z分别为探测器位置矢量的三方向分量r＝(x y z)，r为位置矢量模值，
[0066]a2b
对位置矢量的偏导数为：
[0067][0068]
步骤2，计算相应的j2项和j
22
项摄动加速度和偏导数；其中，j2项的摄动加速度为：
[0069][0070]
其中，j2为二阶带谐项系数，re为地球赤道半径。
[0071]
j2项的加速度偏导数为：
[0072][0073][0074][0075][0076]j22
项的摄动加速度为
[0077][0078]
其中，j
22
为二阶二次田谐系数，λ
22
为二阶田谐项主轴的地理经度，即赤道椭圆常州的方位角，均为常数；λ
22
为中间计算量，λ
22
＝xcosλ
22
+ysinλ
22
。
[0079]j22
摄动加速度对位置矢量的偏导数为
[0080][0081][0082][0083][0084]
步骤3，计算去掉j2和j
22
项的伪中心；
[0085]
其中，引力加速度可以表达为：
[0086][0087]
其中，a
cbf
是采用非球形引力场加速度函数计算的探测器相对月心固联坐标系的
引力加速度，包含中心引力场加速度和摄动加速度。cbf表示中心天体固联坐标系。ρ定义为探测器相对引力加速度a
cbf
的伪月心距离(pseudo－radius)，ρ＝|ρ|。
[0088]
探测器实际的月心距r定义为伪月心距加上伪中心(pseudo－center)c：
[0089]
c＝r-ρ
ꢀꢀꢀ
(18)
[0090]
伪中心c与摄动加速度a
cbf
的关系为：
[0091][0092]
计算去掉和后的摄动引力加速度a0为
[0093][0094]
进一步可得到月心固连坐标系下三维空间任意一点位置r和引力加速度a0与伪中心c0(c0是与a0对应的伪中心)的映射为
[0095][0096]
步骤4，将三维空间分解为经度、纬度和高度三个方向，对伪中心进行离散化，获得高度方向的伪中心多项式拟合，具体过程为：
[0097]
在经度和纬度方向，将球面离散为若干网格点，对于给定的网格点坐标其高度方向的伪中心c0＝[c
x c
y cz]采用有理多项式进行拟合：
[0098][0099]
i＝x,y,z代表伪中心矢量的x,y,z三个方向，dk为分母多项式，nk为分子多项式，a1…an-1
、b1…bn-1
均为多项式系数。拟合方法采用通用的pade近似函数逼近。其中，
[0100][0101]hmax
和h
min
分别为拟合的高度区间，对于低轨近圆环月轨道，通常取h
min
＝10km，h
max
＝350km可以满足大多数任务的要求。不同m和n的取值，会影响拟合精度和计算速度。通过大量的计算和测试，取n＝7和m＝6可以在存储量、计算精度与计算速度之间取得较好的权衡。
[0102]
步骤5，根据给定的探测器位置，得到高度h、经度和纬度的离散网格点；对高度h、经度和纬度的二维离散网格点，进行二维六点插值，得到对应网格点的伪中心插值结果。采用双变量六点插值方法对任意点的插值算法如下：
[0103][0104]
其中，p和q为经度和纬度的归一化参数：
[0105][0106]
λ0为任意点对应的最近左下角离散网格点的经度，为任意点对应的最近左下角离散网格点的纬度，δλ和为经度和纬度离散网格点宽度，通常取为δλ＝0.5
°
，c(h)参见公式(22)。
[0107]c0,0
、c
0,1
、c
1,0
、c
1,1
、c
0,-1
和c-1,0
为c0同高度上的经纬度方向六个离散点，详见图2，具体表达式为：
[0108][0109]
步骤6，根据插值结果，计算得到对应的引力加速度，并加入j2和j
22
项的影响，最终得到真实加速度，具体过程为：
[0110]
根据插值结果c0，计算对应的伪中心距ρ0：
[0111]
ρ0＝r-c0ꢀꢀꢀ
(27)
[0112]
再计算对应的剩余加速度a0：
[0113][0114]
加入j2和j
22
项的影响，最终得到引力加速度
[0115][0116]
步骤7，计算加速度矢量对位置矢量的偏导数矩阵，具体过程如下：
[0117]
首先计算伪中心c对h的偏导数如下：
[0118][0119]
由于不同插值点所在的高度区间的拟合函数都不同，因此，对于6个二维插值点，均可以计算相应的伪中心对x的偏导数
[0120][0121]
其中，k＝x,y,z代表伪中心矢量的x,y,z三个方向。
[0122]
可以计算得到
[0123][0124]
对于经纬度方向的6个插值点，可以分别插值得到
[0125]
总加速度的偏导数的计算公式如下：
[0126][0127]
其中，
[0128][0129]
上式中，
[0130][0131]
[0132][0133]
参照式(6)计算；按式(34)计算；按式计算；按式(13)计算。
[0134]
步骤8，根据计算得到的状态转移矩阵偏导数带入公式(1)，以公式(2)为初值，求解相应的常微分方程组初值问题，就可以获得任意t时刻的位置r(t)、速度v(t)和状态转移矩阵φ(t,t0)，至此完成轨道预报问题的求解。
[0135]
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。