一种双重自适应切比雪夫皮卡迭代法的轨道预测方法

1.本发明涉及航天轨道设计技术领域，具体是一种双重自适应切比雪夫皮卡迭代法的轨道预测方法。


背景技术：

2.卫星轨道预测问题是航天领域一个基本且重要的实际问题，通常只能通过使用数值积分方法来解决。常用的数值积分器分为单步法、多步法。而单步法与多步法只有使用较小的步长才能保证计算不发散，在卫星轨道预测方面会面临精度与效率的问题。
3.随着计算机及数值理论的快速发展，皮卡迭代法由于其快速收敛性与可并行性得到航天轨道计算领域的学者的关注。该方法最早可追溯到1963年clenshaw及其合作者的工作，他们结合切比雪夫多项式和皮卡迭代发展了常微分方程数值求解的有效迭代算法。bai和junkins在2010年提出了修正切比雪夫皮卡迭代方法(mcpi)。基于mcpi方法，woollands和junkins针对摄动二体动力学问题，提出了一种求解初边值问题的自适应切比雪夫皮卡法。而这些切比雪夫皮卡方法总是选择足够大的节点数目n，使得力函数可以很好地用多项式逼近，然后进行数值picard迭代和积分误差反馈。足够大的节点数目n虽然可确保插值误差足够小，但同时也会带来计算成本增加、迭代效率降低等问题。


技术实现要素：

4.针对上述现有技术中的不足，本发明提供一种双重自适应切比雪夫皮卡迭代法的轨道预测方法，由于切比雪夫皮卡迭代法通常以粗略的近似解开始迭代，意味着在最初的几次迭代中迭代误差相对较大。在这种情况下，根据误差分析，不需要高精度的插值，这表明可以在开始时采用较小的n的点来节省计算量。当迭代误差减小时，再自适应地增加插值节点的数量，通过平衡迭代误差与插值误差来自适应地决定插值节点的数量以提高切比雪夫皮卡迭代法的效率，降低计算成本。
5.为实现上述目的，本发明提供一种双重自适应切比雪夫皮卡迭代法的轨道预测方法，包括如下步骤：
6.步骤1，获得航天器的初始位置与初始速度并选择航天器的初始预估轨道；
7.步骤2，自适应设置预测时间步长，并设置预测时间步长内的初始节点数目为n0；
8.步骤3，令n＝n0，其中，n为预测时间步长内的节点数目；
9.步骤4，基于状态初值与速度初值，采用切比雪夫皮卡迭代法对航天器在各时间步长内各节点上的状态信息、速度信息进行迭代预测；
10.每次迭代过程中，基于当前迭代步的预测结果与上一迭代步的预测结果得到迭代误差e1，并基于上一迭代步过程中的力函数得到插值误差e2；
11.判断e1《ε是否成立：
12.若是则输出当前迭代步预测结果中的状态信息与速度信息；
13.否则，继续判断e2≤ke1是否成立，若成立则进行下一次迭代，否则令n＝n+δn并更
新当前迭代步预测结果中的状态信息与速度信息后再进行下一次迭代；
14.其中，ε为容差，k为迭代因子，δn为节点增量；
15.步骤5，基于步骤4输出的各节点的状态信息与速度信息，通过插值技术获得航天器在任意时刻的状态信息与速度信息。
16.在其中一个实施例，步骤4中，每次迭代过程中的状态信息、速度信息进行迭代预测，具体为：
[0017][0018]
式中，vi为第i次迭代的速度信息，v0为速度初值，s为积分矩阵，[f
i-1
]为第i-1次迭代中力函数的值矩阵，xi为第i次迭代的状态信息，x0为速度初值。
[0019]
在其中一个实施例，所述积分矩阵具体为：
[0020][0021]
其中：
[0022][0023][0024]
[0025]
式中，t0表示预测时间步长的初始时间，tf表示预测时间步长的终止时间，矩阵t表示切比雪夫值矩阵，τ0～τn表示切比雪夫点，tn(
·
)～表示0～n次切比雪夫多项式，矩阵w表示对角矩阵，矩阵c表示转换矩阵，t
t
表示矩阵t的转置。
[0026]
在其中一个实施例，步骤4中，所述迭代误差为：
[0027]
e1＝||x
i-x
i-1
||+||v
i-v
i-1
||
[0028]
所述插值误差为：
[0029]
e2＝||e[f
i-1
]||
[0030]
式中，x
i-1
为第i-1次迭代的状态信息，v
i-1
为第i-1次迭代的速度信息，e为误差矩阵。
[0031]
在其中一个实施例，所述误差矩阵具体为：
[0032][0033]
其中：
[0034][0035][0036]
式中，i
n+1
表示n+1维单位矩阵，矩阵表示截断切比雪夫值矩阵，矩阵表示截断对角矩阵，表示矩阵的转置。
[0037]
在其中一个实施例，步骤4中，所述更新当前迭代步预测结果中的状态信息与速度信息，具体为：
[0038][0039]
式中，等式左边的xi、vi为第i次迭代更新后采样点数目增加后修正的状态信息与速度信息，等式右边的xi、vi为第i次迭代更新后采样点数目增加前的状态信息与速度信息，p为重采样矩阵，即将采样点数从n变更为n+δn。
[0040]
在其中一个实施例，步骤2中，所述自适应设置预测时间步长，具体为：
[0041]
步骤2.1，从近地点开始将每个轨道周期分为x段，并设置每个预测时间步长上节点数目的最大数量为n
max
，其中，x为奇数；
[0042]
步骤2.2，令x＝3；
[0043]
步骤2.3，根据指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度；
[0044]
步骤2.4，以近地点为起始点，沿着开普勒轨道向前传播轨道的1/x，利用n1个节点
的切比雪夫多项式逼近该1/x段开普勒轨道上的引力函数：
[0045]
步骤2.5，判断切比雪夫多项式的最后三个系数是否均小于给定误差的百分之一：
[0046]
若是，则采用当前1/x的轨道作为预测时间步长；
[0047]
否则，则令节点数目n1＝n1+δn
′
后重复步骤2.4-2.5直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或n1≥n
max
，其中，若n1≥n
max
则进行步骤2.6；
[0048]
步骤2.6，令x＝x+2后重复步骤2.3-2.5，直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一；
[0049]
其中，若初始和/或最终时间不在预先计算的分段点上，则缩短第一个预测时间步长和/或最后一个预测时间步长。
[0050]
与现有技术相比，本发明具有如下有益技术效果：
[0051]
本发明通过在迭代前期采用较小的节点数目n来节省计算量，随着迭代过程中迭代误差减小时，再自适应地增加插值节点n的数量，并通过平衡迭代误差与插值误差来自适应地决定插值节点的数量以提高切比雪夫皮卡迭代法的效率，进而可以在轨道预测的过程中，不仅能保持高精度预测，还能有效地降低计算成本。
附图说明
[0052]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
[0053]
图1为本发明实施例中预测方法的流程图；
[0054]
图2为本发明实施例中自适应更新节点数的流程图；
[0055]
图3为本发明实施例中mcpi和rmcpi方法的迭代和插值误差的示意图，其中，(a)为一阶系统的误差示意图，(b)为二阶系统的误差示意图；
[0056]
图4为本发明实施例中leo(10个轨道)的收敛速度比较示意图，其中，(a)为apc-g与rapc-g的收敛速度比较示意图，(b)为apc-qg与rapc-qg的收敛速度比较示意图；
[0057]
图5为本发明实施例中heo(10个轨道)的收敛速度比较示意图，其中，(a)为apc-g与rapc-g的收敛速度比较示意图，(b)为apc-qg与rapc-qg的收敛速度比较示意图；
[0058]
图6为本发明实施例中geo(10个轨道)的收敛速度比较示意图，其中，(a)为apc-g与rapc-g的收敛速度比较示意图，(b)为apc-qg与rapc-qg的收敛速度比较示意图；
[0059]
图7为本发明实施例中将leo和heo传播10个开普勒轨道周期所需的等效力函数计算次数示意图，其中，(a)为leo的等效力函数计算次数示意图，(b)为heo的等效力函数计算次数示意图；
[0060]
图8为本发明实施例中leo(10轨道)的哈密顿误差比较示意图，其中，(a)为apc-g与rapc-g的哈密顿误差比较示意图，(b)为apc-qg与rapc-qg的哈密顿误差比较示意图；
[0061]
图9为本发明实施例中heo(10轨道)的哈密顿误差比较示意图，其中，(a)为apc-g与rapc-g的哈密顿误差比较示意图，(b)为apc-qg与rapc-qg的哈密顿误差比较示意图；
[0062]
图10为本发明实施例中geo(10轨道)的哈密顿误差比较示意图，其中，(a)为apc-g与rapc-g的哈密顿误差比较示意图，(b)为apc-qg与rapc-qg的哈密顿误差比较示意图；
[0063]
图11为本发明实施例中在七周(100个轨道)的哈密顿量中实现13位精度具有70度和阶次球谐重力场的molniya轨道的传播示意图。
[0064]
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
[0065]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0066]
在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0067]
在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接，还可以是物理连接或无线通信连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
[0068]
另外，本发明各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
[0069]
如图1所示为本实施例公开的一种双重自适应切比雪夫皮卡迭代法的轨道预测方法，具体包括如下步骤：
[0070]
步骤1，获得航天器的初始位置与初始速度并选择航天器的初始预估轨道，本实施例中，采用开普勒轨道作为预估轨道。
[0071]
步骤2，自适应设置预测时间步长，并设置预测时间步长内的初始节点数目为n0。其中，自适应设置预测时间步长的具体实施例过程为：
[0072]
步骤2.1，预测时间步长根据真近点角的均匀划分来确定，从近地点开始将每个轨道周期分为x段，并设置每个预测时间步长上节点数目的最大数量为n
max
，其中，x为奇数；
[0073]
步骤2.2，令x＝3；
[0074]
步骤2.3，根据指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度；
[0075]
步骤2.4，以近地点为起始点，沿着开普勒轨道向前传播轨道的1/x，利用n1个节点的切比雪夫多项式逼近该1/x段开普勒轨道上的引力函数：
[0076]
步骤2.5，判断切比雪夫多项式的最后三个系数是否均小于给定误差的百分之一：
[0077]
若是，则采用当前1/x的轨道作为预测时间步长；
[0078]
否则，则令节点数目n1＝n1+δn
′
后重复步骤2.4-2.5直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或n1≥n
max
，其中，若n1≥n
max
则进行步骤2.6；
[0079]
步骤2.6，令x＝x+2后重复步骤2.3-2.5，直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一；
[0080]
其中，δn
′
为预测时间步长自适应设置过程中的节点增量，若初始和/或最终时间
不在预先计算的分段点上，则缩短第一个预测时间步长和/或最后一个预测时间步长。
[0081]
下面结合具体的示例对预测时间步长以及节点数目的设置进行进一步的说明。
[0082]
例如：从近地点开始将每个轨道周期分为奇数段，从三段开始，最开始将轨道周期划分为三段。每个预测时间步长的最大节点数n
max
＝40。根据用户指定的初始条件，首先转换为经典轨道元素，然后用于计算开普勒轨道近地点的位置和速度。以近地点为起始点，沿着开普勒轨道向前传播轨道的三分之一，利用10个切比雪夫点的切比雪夫多项式逼近该段轨道上的引力函数，如果切比雪夫多项式的最后三个系数均小于给定误差的百分之一，则采用三分之一的轨道作为预测时间步长，并将预测时间步长上节点的数目设为10。如果达不到误差要求，则将节点数目加倍。例如20、40个节点，直到达到误差要求或节点数目达到最大值40。若节点数目达到最大值还未满足误差要求，则减少预测时间步长为五分之一的轨道，节点数目再从10个节点开始计算，预测时间步长按下一个奇数分之一的轨道进行减少，依次类推，直至达到误差要求。一旦确定了预测时间步长和每个预测时间步长上的节点数，就可以开始逐段利用切比雪夫皮卡迭代法进行轨道传播计算。
[0083]
步骤3，令n＝n0，其中，n为预测时间步长内的节点数目；
[0084]
步骤4，基于状态初值与速度初值，采用切比雪夫皮卡迭代法对航天器在各时间步长内各节点上的状态信息、速度信息进行迭代预测；
[0085]
参考图2，每次迭代过程中，基于当前迭代步的预测结果与上一迭代步的预测结果得到迭代误差e1，并基于上一迭代步过程中的力函数得到插值误差e2；
[0086]
判断e1《ε是否成立：
[0087]
若是则输出当前迭代步预测结果中各节点的状态信息与速度信息；
[0088]
否则，继续判断e2≤ke1是否成立，若成立则进行下一次迭代，否则令n＝n+δn并更新当前迭代步预测结果中的状态信息与速度信息后再进行下一次迭代；
[0089]
其中，ε为容差，k为迭代因子，δn为节点增量；
[0090]
步骤5，基于步骤4输出的各节点的状态信息与速度信息，过插值技术获得航天器在任意时刻的状态信息与速度信息。
[0091]
在具体实施过程中，将轨道预测问题转换为二阶微分方程组的求解的问题，采用切比雪夫皮卡迭代法对二阶微分方程组进行求解，得到的最终的近似解即为轨道预测。其具体实施过程为：
[0092]
首先，将轨道预测问题转换为二阶微分方程组的求解的问题，为：
[0093]
x
″
(t)＝f(t,x(t),x
′
(t)),t0≤t≤tfꢀꢀꢀ
(1)
[0094]
式中，x
″
(t)表示位置状态函数的二阶导数即加速度函数，x(t)表示位置状态函数，x
′
(t)表示位置状态函数的一阶导数即速度函数，t表示时间节点，t0表示预测时间步长的初始时间，tf表示预测时间步长的终止时间，表示一个四元向量值力函数；
[0095]
初始状态信息x0和初始速度信息v0的获取方式可以直接获取或通过包含六要素的tle两行轨道数据转换得到，在初始条件x(t0)＝x0，x
′
(t0)＝v0下，可以将式(1)改写成为一阶状态空间方程组，为：
[0096]
[0097]
将其转换为积分方程，为：
[0098][0099]
其中，s表示积分变量，给定初始的状态信息x0,速度信息v0和以开普勒轨道作为初始轨道x0(t)、v0(t)，可通过皮卡迭代产生一列近似解xi(t)、vi(t)，i＝1,2，
···
，其中，二阶微分方程组的切比雪夫皮卡迭代公式为：
[0100][0101]
当f是光滑、可积的单值函数且区间[t0,tf]的长度足够小时，这一列近似解xi(t)、vi(t)，i＝1,2，
···
可以收敛到精确解。
[0102]
切比雪夫皮卡迭代法是在皮卡迭代的基础上，通过基于切比雪夫节点和第一类切比雪夫多项式逼近未知轨迹和非线性函数而发展起来的。
[0103]
假设给定一个初始轨迹x0，未知轨迹xi和沿轨迹x
i-1
的力函数f由正交切比雪夫多项式的有限项逼近，为：
[0104][0105]
式中，是由最小二乘法确定的系数，tk是第一类k阶切比雪夫多项式，其中，如果节点数目n足够大，则可以获得力函数的近机器精度近似。由于切比雪夫多项式的积分可以解析计算，力函数的积分也可以近似为具有接近机器精度的截断切比雪夫级数，从而获得更精确的近似轨迹xi。为了确定系数mcpi方法利用切比雪夫节点，它们是tn的极值，为：
[0106][0107]
式中，τj表示切比雪夫节点，j表示节点次序。
[0108]
切比雪夫节点具有非均匀余弦分布，可以最大限度地减少由龙格现象引起的误差。由于正交性，系数可以通过矩阵向量乘积而不是离散样本点的矩阵求逆来获得:
[0109]
[f
i-1
]＝[f(sj,x
i-1
(sj),v
i-1
(sj))]
(n+1)
×m[0110]
sj＝s(τj)
[0111]
式中，[f
i-1
]为第i-1次迭代中力函数的值矩阵，x
i-1
(sj)表示第i-1次迭代位置函数在sj处的值，v
i-1
(sj)表示第i-1次迭代速度函数在sj处的值，sj表示第j个修正的切比雪夫点，s(τj)表示转换函数在第j个切比雪夫点处的值。
[0112]
借助切比雪夫积分关系可以得到解轨迹，其中，切比雪夫积分关系为：
[0113]
∫t0(s)ds＝t1(s)
[0114][0115][0116]
式中，tk(s)为第一类k次切比雪夫多项式。
[0117]
解轨迹可以通过截断切比雪夫级数形式给出，为：
[0118][0119]
其中，系数可以通过矩阵向量的形式获得。为了呈现切比雪夫皮卡迭代法的向量矩阵形式，定义矩阵如下：
[0120][0121][0122][0123]
式中，矩阵t表示切比雪夫值矩阵，τ0～τn表示切比雪夫点，t0(
·
)～tn(
·
)表示0～n次切比雪夫多项式，矩阵w表示对角矩阵，矩阵c表示转换矩阵。
[0124]
令其中，t
t
表示矩阵t的转置，可以得到轨迹解的矩阵形式为：
[0125][0126]
其中，为速度初始值，为位置初始值，βi为第i次迭代速度函数的系数，αi为第i次迭代位置函数的系数。
[0127]
具体到本实施例中，区别于上述现有的切比雪夫皮卡迭代法，提出一种新的向量矩阵形式，该方法不使用系数，而是使用离散状态值。为此，本实施例引入了一个矩阵，为：
[0128][0129]
该矩阵s为积分矩阵，其功能等同于积分符号。本实施例中跳过了计算系数环节，直接获得节点处的函数值，可以降低矩阵向量乘积运算。
[0130]
在定义积分矩阵s后，令xi和vi分别表示xi和vi在修正切比雪夫点{sj:j＝0,1,
···
,n},i＝1,2,
···
上的值，以及：
[0131][0132]
可以得到每次迭代过程中的状态信息、速度信息进行迭代预测的矩阵形式，具体为：
[0133][0134]
式中，vi为第i次迭代的速度信息，v0为速度初值，s为积分矩阵，[f
i-1
]为第i-1次迭代中力函数的值矩阵，xi为第i次迭代的状态信息，x0为速度初值。
[0135]
在迭代过程中，假设已知轨迹x
i-1
。令表示非线性函数的插值，式(4)可以更精确地写为：
[0136][0137]
用e(t)＝v(t)-vi(t)定义第i个速度信息近似解的误差。可以得到：
[0138][0139]
从式(6)可知，e1(t)与精确的皮卡迭代有关，而e2(t)与力函数插值有关，由此可见，切比雪夫皮卡迭代法的误差可以分解为两部分，其一是皮卡迭代过程中的迭代误差e1(t)，另一个则为与力函数插值相关的插值误差e2(t)。
[0140]
基于上述误差分解，本实施例提出了一种自适应地确定每次迭代中的节点数目的方法。一旦将经典的切比雪夫皮卡迭代法与自适应方案相结合，获得新的迭代法，其具体实施原理为：
[0141]
切比雪夫皮卡迭代法(例如mcpi方法)总是选择足够大的节点数目n，使得力函数可以很好地用多项式逼近，然后进行皮卡迭代和积分误差反馈修正。足够大的n可确保插值误差足够小，但也会增加函数计算的计算成本。由于切比雪夫皮卡迭代法通常以粗略的近似解开始迭代，这意味着在最初的几次迭代中迭代误差相对较大。在这种情况下，根据误差分析，不需要高精度的插值，这表明可以在开始时采用较小的n的点来节省计算量。当迭代误差减小时，再自适应地增加插值节点的数量。通过平衡迭代误差与插值误差来自适应地决定插值节点的数量以提高切比雪夫皮卡迭代法的效率，其具体的实施过程为；
[0142]
从一个相对较小的节点数目n和一个初始的近似轨迹x0(t)开始，执行一次皮卡迭代以获得一个新的近似轨迹x1(t)，以n次多项式表示；
[0143]
然后计算迭代误差e1(t)和插值误差e2(t)，如果插值误差大于迭代误差，即e2(t)》e1(t)，则将n值增加δn并将近似轨迹x1(t)表示为(n+δn)次多项式的形式，即将更高阶的系数设置为0；
[0144]
然后继续进行皮卡迭代以产生新的轨迹x2(t)，如果e2(t)≤e1(t)，不做任何更改，只需进行皮卡迭代以产生新的轨迹x2(t)；
[0145]
依次类推，直到达到给定的误差限。如果e2(t)达到机器精度、给定的误差限或n达到给定的最大值，则n的值不增加。需要注意的是，当选择的节点数目n使得e2(t)与e1(t)大小相当时，可能会影响迭代收敛，这意味着需要更多的迭代次数。为了避免这种情况，在实际运行中，可以要求插值误差e2(t)应该小于迭代误差e1(t)乘上一个非常小的因子k《1，即e2(t)≤ke1(t)，使得最终不会影响收敛速度。
[0146]
为了自适应地确定节点数目n，得到迭代误差与插值误差的数值估计至关重要，对于迭代误差的近似，可以直接利用当前和先前的状态值，并通过它们的差异来近似迭代误差，即为：
[0147]
e1＝||x
i-x
i-1
||+||v
i-v
i-1
||
ꢀꢀꢀ
(7)
[0148]
式中，xi为第i次迭代的状态信息，x
i-1
为第i-1次迭代的状态信息，vi为第i次迭代的速度信息，v
i-1
为第i-1次迭代的速度信息。
[0149]
而对于插值误差的估计，切比雪夫迭代法的插值误差主要由以下项主导，为：
[0150][0151]
因为其它项中都包含一个较小乘积因子，例如因此，只需要考虑力函数的插值误差即可，从而提高计算效率。
[0152]
估计插值误差的详细过程为：首先计算力函数拟合的最小二乘系数，然后计算删除最后一个系数后的多项式与力函数的精确值在插值点上的误差。本实施例中提出了一个用于估计插值误差的向量矩阵格式。为此，首先令：
[0153][0154][0155][0156]
式中，矩阵表示截断切比雪夫值矩阵，τ0～τn表示切比雪夫点，t0(
·
)～t
n-1
(
·
)表示0～(n-1)次切比雪夫多项式，矩阵表示截断对角矩阵，矩阵w表示对角矩阵；
[0157]
定义误差矩阵e，为：
[0158][0159]
式中，i
n+1
表示n+1维单位矩阵，表示矩阵的转置。本实施例中定义该误差矩阵来计算插值误差，具有操作简单且计算高效的效果。
[0160]
力函数在修正的切比雪夫点sj，j＝0,1，
···
，n上的插值误差则为：
[0161]
e2＝||e[f
i-1
]||
[0162]
该插值误差e2的计算原理为通过系数截断，用低一阶多项式来近似原多项式的插值误差。
[0163]
在步骤4中，更新当前迭代步预测结果中的状态信息与速度信息，具体为：在重新确定节点后通过重采样矩阵来更新当前迭代步的状态信息与速度信息。本实施例中的重采样矩阵为：
[0164]
p＝[p
i,j
]
(n+1)
×
(n+1)
[0165][0166]
式中，p为重采样矩阵，p
i,j
表示重采样矩阵的第i行j列元素，(n+1)
×
(n+1)表示重采样矩阵的规模，ω
j,n
表示第j个重心插值系数，ti表示第i个需要插值的节点，sj表示第j个修正的切比雪夫点，ω
l,n
表示第l个重心插值系数，s
l
表示第l个修正的切比雪夫点；
[0167]
因此更新当前迭代步预测结果中的状态信息与速度信息的过程为：
[0168][0169]
式中，等式左边的xi、vi为第i次迭代更新后采样点数目增加后修正的状态信息与速度信息，等式右边的xi、vi为第i次迭代更新后采样点数目增加前的状态信息与速度信息，p为重采样矩阵，即将采样点数从n变更为n+δn。
[0170]
下面结合具体的算例对本发明作出进一步的说明。
[0171]
首先，通过考虑一阶和二阶非线性系统来考虑新方法的效率：
[0172]
一阶非线性系统为：
[0173][0174]
一阶非线性系统为：
[0175][0176]
其中，t0＝1，tf＝256π，ε＝0.001。
[0177]
这两个问题常被作为一个已知解析解的基准来进行算法准确性研究。对于mcpi方法，方法使用阶n＝130的chebyshev多项式来逼近区间内的解16π的长度。对与本发明方法，设置一阶系统的n0＝700、δn＝200、n
max
＝1500、k＝10-8
；二阶系统设置为n0＝50、δn＝20、n
max
＝130、k＝10-8
。对于所有情况，停止迭代的容差设置为ε＝10-11
。
[0178]
mcpi和本发明方法(rmcpi)的迭代误差和插值误差如图3所示。对于一阶和二阶系统，可以观察到mcpi方法的收敛性几乎不受早期迭代中使用较小点的影响。这说明点数调整思想在本发明方法中行之有效。表1中给出了力函数计算次数的比较。发现对于一阶和二阶系统，在不影响迭代次数和精度的情况下，力函数计算次数分别减少了22.2％和26.2％。
[0179]
表1：力函数计算次数的比较
[0180][0181]
另外，将本发明方法应用于扰动轨道传播，由于求解椭圆引力扰动轨道的二阶系统，本算例主要考虑两种二阶apc方法，apc-iig和apc-qiig，以及相应的修正版本。为简化符号，apc-iig和apc-qiig方法分别缩写为apc-g和apc-qg，基于本发明的方法修正apc-iig
和apc-qiig方法分别表示为rapc-g和rapc-qg。
[0182]
为了全面展示本发明方法的效率，算例测试了三种不同的场景，分别为：低地球轨道(leo)、高偏心轨道(heo)和地球同步轨道(geo)。采用三种轨道状态的常用参数和初始条件，列于表2中。
[0183]
表2：轨道参数与初始条件
[0184][0185]
为了获得更好的性能，本算例使用了变增量δn，并且针对不同场景的各预测方法方法的参数设置列在表3中，rapc表示，rapc-q表示。其中，还可以通过数值实验进一步优化设置。
[0186]
表3：参数设置
[0187][0188]
首先检查新方法在每个部分上的收敛性。为了方便比较，插值误差和迭代误差均归一化为无量纲误差。在rapc-g方法的实现中迭代误差与插值误差由下式给出：
[0189]
迭代误差为：
[0190][0191]
插值误差为：
[0192][0193]
式中，du、vu和au分别是距离单位、速度单位和加速度单位，||
·
||表示最大范数。rapc-qg方法的误差可以类似地定义。图4-6显示了本发明在一个时间步长中的收敛结果，三个轨道测试均采用70阶球谐重力模型(地球重力模型2008)。结果表明，在大多数情况下，本发明的收敛性与经典切比雪夫皮卡迭代法具有很好的一致性。特别地，在leo和geo情况
下，rapc-qg方法比apc-qg方法收敛得更快。结果验证了每次迭代的节点自适应方案是成功的。
[0194]
图7显示了leo和heo传播10个开普勒轨道周期所需的等效力函数计算的数量。在每幅图中，第五条(rapc)表示在基于apc方法的每次迭代中仅采用节点自适应方案时的力函数计算次数。与第三条(apc-q)相比，表明所提出方法的贡献与准线性化技术的贡献相当。比较第四条(apcqg)和最后一条(rapc-qg)发现，本发明提出的的节点自适应方案可以进一步减少当前最先进的自适应切比雪夫皮卡迭代法的力函数计算次数。
[0195]
本算例还通过哈密顿误差、函数计算的总次数和计算时间进一步证明了本发明的性能。哈密顿量的相对误差，定义为其中δh是哈密顿量在每个采样点上的差异，h0是初始时刻的哈密顿量。计算时间是十次程序运行的平均时间。三个轨道测试用例的结果显示在表4-6和图8-10中。从图8-10清楚地表明，本发明所提出的方法在所有三种情况下具有相当的哈密顿误差。从表4-6中观察到，在达到相同精度下，与经典apc方法相比，所有三种情况下力函数计算数量显著减少。对于leo和heo情况，计算时间也相应减少。对于meo案例，计算时间的减少不明显，因为此时个方法中函数计算次数都相对较少，在时间消耗上不占主导地位。综上数值结果显示了本发明方法的效率优势。
[0196]
表4:leo(10个轨道)性能比较
[0197] apc-grapc-gapc-qgrapc-qg最大哈密顿误差1.00
×
10-14
9.15
×
10-15
3.12
×
10-14
3.27
×
10-14
力函数计算次数229887640168206308计算时间(s)1.63241.13311.43581.0387
[0198]
表5:heo(10个轨道)性能比较
[0199] apc-grapc-gapc-qgrapc-qg最大哈密顿误差2.43
×
10-11
4.14
×
10-12
2.43
×
10-11
4.12
×
10-12
力函数计算次数13054215888901659计算时间(s)0.74980.60590.72270.5613
[0200]
表6:geo(10个轨道)性能比较
[0201] apc-grapc-gapc-qgrapc-qg最大哈密顿误差5.32
×
10-15
1.36
×
10-14
4.93
×
10-15
1.10
×
10-14
力函数计算次数8317550640计算时间(s)0.14890.14610.12370.1257
[0202]
为了测试长时稳定性的性能，本算例计算传播100个开普勒轨道周期(长于7周)的molniya轨道(a＝26,554km，e＝0.72，i＝63deg，ω＝0deg，ω＝0deg，m＝0deg)，且考虑到70阶球谐重力场。rapc-qg方法的哈密顿误差，如图11所示，精度达到了13位，这与apc-qg方法的一致。因此，在这种情况下，每次迭代中的节点自适应技术不会影响切比雪夫皮卡迭代法的良好稳定性。并且rapc-qg的力函数计算次数和计算时间都小于apc-qg方法，如表6所示。
[0203]
表6:molniya轨道性能比较
[0204] apc-qgrapc-qg最大哈密顿误差8.12
×
10-14
7.13
×
10-14
力函数计算次数11094130855计算时间(s)17.995115.2593
[0205]
综上所述，本发明提出了一类双重自适应切比雪夫皮卡迭代法来传播扰动轨道。该方法包括一个自适应方案来决定每次皮卡迭代中的节点数。本方法还能结合了切比雪夫皮卡迭代法的现有增强技术，例如准线性化、局部重力近似等。在达到相同精度的情况下，双重自适应切比雪夫皮卡迭代法比经典自适应切比雪夫皮卡迭代法更有效，且具有长时稳定性。
[0206]
以上所述仅为本发明的优选实施例，并非因此限制本发明的专利范围，凡是在本发明的发明构思下，利用本发明说明书及附图内容所作的等效结构变换，或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。