一种地球轨道中航天器的轨迹预测方法及系统与流程

本发明涉及航天动力学技术领域：，特别涉及一种地球轨道中航天器的轨迹预测方法及系统。背景技术：：在未来空间任务中，航天器自主程度更高，协同和交互操作更加频繁，研究新型高性能空间轨道计算方法对于航天工程具有重大应用价值。对于轨道递推，研究人员关注的焦点在于计算精度和计算效率。过去的几十年中，已经有大量研究人员对轨道递推计算方法进行研究。传统的数值方法中，为了得到高精度的计算结果，需要选择非常小的计算步长，并构造大规模的代数矩阵方程，这在实际操作中往往难以实现。技术实现要素：本发明的目的是提供一种地球轨道中航天器的轨迹预测方法及系统，以提高轨迹预测的精度和效果。为实现上述目的，本发明提供了如下方案：一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。可选的，所述从预测时间步长内选取多个配点，具体包括：利用公式ξm＝cos[(m-1)π/(m-1)],m＝1,...,m，在[-1,1]区间中选取m个节点；分别令ξ＝ξ1,…,ξm，利用公式确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξm表示分别表示第1个、第m个和第m个节点，m表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。可选的，所述根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息，具体包括：判断预测时间步长是否小于一个轨道周期，获得第一判断结果；若所述第一判断结果表示是，则，分别令t＝t1,t2,…,tm，根据航天器的初始位置和初始速度，利用公式确定航天器在每个配点时的初始预测状态信息；其中，x(t)表示实际时间t的初始预测状态信息，r(t)和v(t)分别表示实际时间t的初始预测位置矢量和初始预测速度矢量，r0和v0为航天器的初始位置和初始速度，t1、t2和tm分别表示第1个、第2个和第m个配点；若所述第一判断结果表示否，则分别令t＝t1,t2,…,tm，根据航天器的初始位置r0和初始速度v0，利用公式确定航天器在每个配点时的初始预测状态信息；其中，a表示轨道长半轴，a＝-μ/(||v0||2-2μ/||r0||)，e为偏心率，e为偏近点角，通过求解关系式获得，为偏近点角的时间导数，μ为地球引力常数。可选的，所述采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；其中，un和un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息；t0、t1和g分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，t0＝diag[-1,-1,...]，g＝eun-f(un,t)，j为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tm]，t1、t2和tm分别表示第1个、第2个和第m个配点，f为一阶形式的轨道动力学方程的右端函数，e和分别为第一常值矩阵和第二常值矩阵，e＝diag[(lb)b-1,(lb)b-1,...,(lb)b-1,...,(lb)b-1]，l为微分算子，l-1为积分算子，diag[·]表示对角矩阵，b表示变换矩阵，b＝[φ(t1)t,φ(t2)t,...,φ(tm)t]t，φ(t)为横向量，φ(t)＝[t1(t),t2(t)...tm(t)]，t为实际时间，t＝t1,t2,…,tm,，t1(t)、t2(t)和tm(t)分别为横向量φ(t)中的第1个、第2个和第m个元素，n为插值函数阶数，m为配点个数；判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的的状态信息作为航天器在每个配点时的修正后的状态信息。可选的，所述利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：利用公式tm(t)＝cos[marccos(ξ)]，计算横向量φ(t)中的每个元素；其中，tm(t)表示横向量φ(t)中的第m个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，ti和ti+1为第i个预测时间步长的左端点和右端点。可选的，所述利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：利用公式t0(t)＝1计算横向量φ(t)中的第0个元素t0(t)；利用公式t1(t)＝ξ计算横向量φ(t)中的第1个元素t1(t)；根据横向量φ(t)中的第0个元素和第1个元素，利用公式tm+1(t)＝2ξtm(t)-tm-1(t)，采用迭代的方式计算横向量φ(t)中的第3个到第m个元素；其中，tm-1(t)、tm(t)和tm+1(t)分别表示横向量φ(t)中的第m-1个、第m个和第m+1个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，ti和ti+1为第i个预测时间步长的左端点和右端点。可选的，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，具体包括：采用chebyshev多项式作为差值函数，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。可选的，所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，之后还包括：将当前的预测时间步长输出的航天器在每个时刻的修正后的状态信息，作为下一个预测时间步长的初始状态信息，对下个预测时间步长的航天器的轨迹进行预测；所述初始状态信息为航天器的初始位置和初始速度。一种地球轨道中航天器的轨迹预测系统，所述预测系统包括：配点选取模块，用于从预测时间步长内选取多个配点；初始预测模块，用于根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；修正模块，用于采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；插值模块，用于对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。可选的，所述配点选取模块，具体包括：节点选取子模块，用于利用公式ξm＝cos[(m-1)π/(m-1)],m＝1,...,m，在[-1,1]区间中选取m个节点；配点选取子模块，用于分别令ξ＝ξ1,…,ξm，利用公式确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξm表示分别表示第1个、第m个和第m个节点，m表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明公开了一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。本发明采用高精度的反馈加速picard迭代方法对初始预测状态信息进行修正，无需对矩阵求逆，在提高轨迹预测精度的同时，能够降低轨道计算中的复杂加速度项的计算耗时，并且结合配点法的思想，进一步减少了计算量，本发明为航天器提供了精确、实时的轨道预测方法。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。图1为本发明提供的一种地球轨道中航天器的轨迹预测方法的流程图；图2为本发明提供的一种地球轨道中航天器的轨迹预测方法的原理图；图3为本发明提供的一个预测时间步长的轨道递推结果；图4为本发明提供的本发明获取的计算结果的hamilton(哈密顿)量的相对误差；图5为本发明提供的反馈加速picard迭代方法与rk12(10)计算结果的位置误差对比图。具体实施方式本发明的目的是提供一种地球轨道中航天器的轨迹预测方法及系统，以提高轨迹预测的精度和效果。为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。如图1和2所示，本发明提供一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：步骤101，从预测时间步长内选取多个配点。配点法也是求解轨道递推问题的一个重要的工具。在配点法中，问题的解是通过多项式、谐波函数、或者其他的基函数估计得到的。相比于其他纯粹的数值方法，配点法的优势在于能够在相对更大的时间尺度上得到半解析的估计解，同时只需要存储较少的基函数系数。newton-kantorovich/chebyshev(牛顿·坎托罗维奇/切比雪夫)伪谱法和rbfs(radialbasisfunctions，径向基插值方法)方法分别使用chebyshev(切比雪夫)多项式和径向基函数作为基函数对配点插值，都能够精确高效的求解轨道动力学问题。1-1选择初始化轨道方案：根据用户或任务需求，选择单次轨道预测时间(预测时间步长)δt＝ti+1-ti，其中ti为单次轨道预测初始时刻，ti+1为单次轨道预测末端时刻。当单次轨道预测时间δt小于一个轨道周期时，采用冷启动方案，以降低初始化过程计算量，跳转至步骤1-2；当单次轨道预测时间δt大于一个轨道周期时，采用热启动方案，以保证后期迭代修正过程能够收敛，跳转至步骤1-3。1-2冷启动方案：冷启动方案是利用初始位置r0及初始速度v0得到直线型估计轨道及其整体状态信息x(t)，其中包括轨道位置信息r(t)及速度信息v(t)。实现方法如下：在已知航天器初始位置r0和速度v0的条件下，假设航天器从起点出发做匀速直线运动，则航天器在t时刻的状态信息x(t)为：x(t)由t时刻的位置矢量r(t)和速度矢量v(t)组成，r0和v0为初始位置和速度，t0为初始时刻。1-3热启动方案：热启动方案是利用初始位置r0及初始速度v0得到开普勒模型的估计轨道及其整体状态信息x(t)，其中包括轨道位置信息r(t)及速度信息v(t)。实现方法如下：基于线性化中心引力场的开普勒轨道为周期性椭圆轨道，具有解析解如下：||r||＝a(1-ecose)(2)e-esine＝2π(t/t)(3)联立以上二式可得轨道位置信息r(t)和速度信息v(t)，从而得到线性化估计轨道整体状态信息x(t)。以上两式中e为偏近点角，a为轨道半长轴，e为偏心率，这些参数均可由经典轨道动力学理论得到。2-1在估计轨道上选取chebyshev-gauss-lobatto(切比雪夫·高斯·洛巴托，cgl)配点：基于轨道状态信息x(t)，通过在时间域上[ti,ti+1]选取离散时间点，得到m个离散化的cgl配点状态信息x(tm)，其中m＝1,2,…,m。实现方法如下：将局部时间区间ti≤t≤ti+1缩放至[-1,1]区间，缩放规则为：其中ti和ti+1为第i个步长的左右端点，t为实际时间，ξ为缩放后的时间。在[-1,1]区间内通过如下公式选取cgl节点ξm＝cos[(m-1)π/(m-1)],m＝1,...,m(5)并通过缩放公式反推出时间节点tm，从而得到cgl配点状态信息x(tm)。步骤101，所述从预测时间步长内选取多个配点，具体包括：利用公式ξm＝cos[(m-1)π/(m-1)],m＝1,...,m，在[-1,1]区间中选取m个节点；分别令ξ＝ξ1,…,ξm，利用公式确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξm表示分别表示第1个、第m个和第m个节点，m表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。步骤102，根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息。步骤102所述根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息，具体包括：判断预测时间步长是否小于一个轨道周期，获得第一判断结果；若所述第一判断结果表示是，则，分别令t＝t1,t2,…,tm，根据航天器的初始位置和初始速度，利用公式确定航天器在每个配点时的初始预测状态信息；其中，x(t)表示实际时间t的初始预测状态信息，r(t)和v(t)分别表示实际时间t的初始预测位置矢量和初始预测速度矢量，r0和v0为航天器的初始位置和初始速度，t1、t2和tm分别表示第1个、第2个和第m个配点；若所述第一判断结果表示否，则分别令t＝t1,t2,…,tm，根据航天器的初始位置r0和初始速度v0，利用公式确定航天器在每个配点时的初始预测状态信息；其中，a表示轨道长半轴，a＝-μ/(||v0||2-2μ/||r0||)，e为偏心率，e为偏近点角，通过求解关系式获得，为偏近点角的时间导数，μ为地球引力常数步骤103，采用反馈加速picard(比卡)迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息。基于步骤102中航天器在每个配点时的初始预测状态信息x(tm)，通过反馈picard迭代方程对其进行修正，得到精确的c航天器在每个配点时的初始预测状态信息实现方法如下：选取第一类chebyshev(切比雪夫)正交多项式作为插值基函数，cgl节点作为时间域节点。通过试算确定基函数个数和配点个数，一般情况下二者相同，从而为基函数系数的确定提供足够多的约束，同时也便于迭代计算并得到高精度解。记基函数个数和配点个数均为m，cgl配点的确定规则为：ξm＝cos[(m-1)π/(m-1)],m＝1,...,m(6)(i)构造反馈picard迭代方程中的常值矩阵e和常值矩阵e的构造表达式为e＝diag[(lb)b-1,(lb)b-1,...,(lb)b-1,...,(lb)b-1]，其中包括m个(lb)b-1；常值矩阵的构造表达式为其中包括m个(l-1b)b-1。上式中b＝[φ(t1)t,φ(t2)t,...,φ(tm)t]t，l为微分算子，l-1为积分算子。diag[·]表示以括号中元素构成的对角矩阵。φ(t)为横向量，表示第一类chebyshev正交基函数系在t时刻的值，φ(t)＝[t1(t),t2(t)...tm(t)]，由函数tn(ξ)＝cos[narccos(ξ)]或者如下迭代关系确定：t0(ξ)＝1，t1(ξ)＝ξ，tn+1(ξ)＝2ξtn(ξ)-tn-1(ξ)(7)ξ为缩放后的时间尺度。由于第一类chebyshev多项式和cgl配点的定义域均为[-1,1],在实际计算中，将局部时间区间ti≤t≤ti+1缩放至[-1,1]区间进行迭代，缩放规则为：其中ti和ti+1为第i个步长的左右端点，t为实际时间，ξ为缩放后的时间，ξ＝t1,t2,...,tm∈[-1,1]为cgl配点。(ii)构造反馈picard迭代方程中的变量矩阵t0,t1,g：变量矩阵t0的表达式为t0＝diag[-1,-1,...]，为mn×mn维单位阵，变量矩阵t1的表达式为其中j为雅可比阵，变量矩阵g的表达式为g＝eun-f(un,t)，其中t＝[t1,t2,...,tm]。以上各式中n为插值函数阶数，m为配点个数，一般使m＝n。另外t＝[t1,t2,...,tm]，f为如下一阶形式的轨道动力学方程右端函数，un为第n次迭代修正的状态矢量，其初始估计u0由步骤2-1中cgl配点状态信息x(tm)，m＝1,2,…,m，合并得到，其中u0＝[x(t1)t,x(t2)t,...x(tm)t]t，上标t代表矩阵转置。(iii)反馈picard修正：将常值矩阵e、变量矩阵t0,t1,g代入下式，对状态矢量un修正，得到修正状态矢量un+1上式的推导过程如下：对于一阶形式的轨道动力学方程组其中x＝x(t),x(t)＝[x1(t),x2(t),...,xd(t),...xn(t)]。根据变分迭代法可得到如下迭代公式其中λ(τ)为待求的拉格朗日乘子。令上式右端项变分为零,得到关于λ(τ)的约束条件：接下来，根据(13)式给出变分迭代法的变型。将λ(τ)估计成taylor级数形式，其估计式如下：λ(τ)≈t0[λ]+t1[λ](τ-t)+...+tk[λ](τ-t)k(14)tk[λ]是λ(τ)的k阶微分变换，即根据约束条件，tk[λ]可以由如下迭代过程得到其中令将λ(τ)的估计式代入原迭代方程，可以得到截取一阶估计式得到通过cgl配点及第一类chebyshev多项式对其进行离散化后即可得到反馈picard修正公式。(iv)判断迭代修正是否收敛：计算当前修正状态矢量un+1与上一步状态矢量un之间的范数误差εn+1＝||un+1-un||。设定计算误差限为δ，若εn+1＜δ，结束迭代修正计算，进入步骤(v)；若εn+1≥δ，则转到步骤(ii)继续计算。(v)提取精确的cgl配点：从迭代修正收敛的状态矢量u提取精确的cgl配点信息其提取规则为构造规则的逆过程，如下式所示：步骤104所述采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息，具体包括：利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息；其中，un和un+1分别为航天器在每个配点时的第n次和第n+1次修正后的状态信息；t0、t1和g分别为第一变量矩阵、第二变量矩阵和第三变量矩阵，t0＝diag[-1,-1,...]，g＝eun-f(un,t)，j为雅可比矩阵，t为配点向量，t＝[t1,t2,...,tm]，t1、t2和tm分别表示第1个、第2个和第m个配点，f为一阶形式的轨道动力学方程的右端函数，e和分别为第一常值矩阵和第二常值矩阵，e＝diag[(lb)b-1,(lb)b-1,...,(lb)b-1,...,(lb)b-1]，l为微分算子，l-1为积分算子，diag[·]表示对角矩阵，b表示变换矩阵，b＝[φ(t1)t,φ(t2)t,...,φ(tm)t]t，φ(t)为横向量，φ(t)＝[t1(t),t2(t)...tm(t)]，t为实际时间，t＝t1,t2,…,tm,，t1(t)、t2(t)和tm(t)分别为横向量φ(t)中的第1个、第2个和第m个元素，n为插值函数阶数，m为配点个数；判断航天器在每个配点时的第n+1次修正后的状态信息与第n次修正后的状态信息的差值是否小于误差限值，获得第二判断结果；若所述第二判断结果表示否，则令n的数值增加1，返回步骤“利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息”；若所述第二判断结果表示是，则输出航天器在每个配点时的第n+1次修正后的的状态信息作为航天器在每个配点时的修正后的状态信息。其中，所述利用公式对航天器在每个配点时的初始预测状态信息进行第n+1次修正，获得航天器在每个配点时的第n+1次修正后的状态信息，之前还包括：利用公式tm(t)＝cos[marccos(ξ)]，计算横向量φ(t)中的每个元素；其中，tm(t)表示横向量φ(t)中的第m个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，ti和ti+1为第i个预测时间步长的左端点和右端点。或利用公式t0(t)＝1计算横向量φ(t)中的第0个元素t0(t)；利用公式t1(t)＝ξ计算横向量φ(t)中的第1个元素t1(t)；根据横向量φ(t)中的第0个元素和第1个元素，利用公式tm+1(t)＝2ξtm(t)-tm-1(t)，采用迭代的方式计算横向量φ(t)中的第3个到第m个元素；其中，tm-1(t)、tm(t)和tm+1(t)分别表示横向量φ(t)中的第m-1个、第m个和第m+1个元素，ξ为实际时间t对应的[-1,1]区间内的缩放后的时间，ti和ti+1为第i个预测时间步长的左端点和右端点。步骤104，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。基于步骤103中精确的cgl配点的状态信息以chebyshev多项式为插值函数，生成精确的预测轨道及其整体状态信息实现方法如下：步骤104所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，具体包括：采用chebyshev多项式作为差值函数，对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。步骤104之后还包括：判断是否结束轨道预测：若当前轨道预测时间长度满足用户或任务需求，结束轨道预测。若需要做更长时间的轨道预测，基于步骤103中预测轨道整体状态信息获得末端时间节点的位置信息r(tf)和速度信息v(tf)。转到步骤101，并使用末端时间节点的位置信息r(tf)和速度信息v(tf)更新步骤1中的初始位置r0及初始速度v0，重复步骤101至步骤104。步骤104所述对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息，之后还包括：将当前的预测时间步长输出的航天器在每个时刻的修正后的状态信息，作为下一个预测时间步长的初始状态信息，对下个预测时间步长的航天器的轨迹进行预测；所述初始状态信息为航天器的初始位置和初始速度。下面结合附图对本发明方法的实施方式作详细说明。本发明针对低轨情形下的航天器轨道递推问题，提出了一种反馈加速picard迭代方法。其具体实施方式包括如下步骤：第一步：对于leo情形下的轨道递推问题，设计修正变分迭代法，得到轨道递推公式。首先推导航天器运动方程的变分迭代法递推公式，根据泛函取极值条件和泰勒展开，得到拉格朗日乘子的估计形式，将其代入原递推公式，得到修正变分迭代法的递推方程。计算e,t0,t1常值矩阵，代入迭代方程(11)。根据工程需要，配置计算步长和迭代终止条件参数。第二步：根据初值条件和修正迭代方程在第一个时间步长内进行轨道计算。由初始位置r0和初始速度得到u0，代入迭代方程(11)式进行轨道计算。第三步：计算状态变量修正量的模，和迭代终止条件对比，判定是否继续进行迭代。将修正量与迭代终止条件比对，若满足迭代终止条件ε≤δ，记录un+1，转入下一区间进行迭代。若ε＞δ，回到第二步继续迭代。第四步：根据上一步长终点信息，估计下一步长的初始状态信息，计算下一步长内的运行轨道。根据上一步长终点处位置、速度信息估计下一步长起点处状态信息，从而估计出u0，转入第二步继续计算，直至终点，从而得到航天器在给定时间区间的完整运行轨道。实例1：在leo情形下，不考虑空气阻力，利用球谐重力场40阶模型(emg2008)进行轨道递推，将常值矩阵e,t0,t1及表1中的计算参数代入公式(11)进行计算，得到一个周期内的轨道，如图3所示。为评估本方法的计算精度，对数值结果hamilton量的变化进行分析。数值结果中hamilton量的相对误差为εh＝δh/h0，其中h0是系统在初始时刻的hamilton量，δh为计算过程中hamilton量的偏差。εh的变化如图4所示。hamilton量的相对计算误差均被限制在10-13附近，这几乎达到了机器计算精度。表1实例1采用的参数值实例2：在实例1的重力场模型中加入大气阻力模型，并将计算结果与runge-kutta12(10)，即龙格库塔12(10)阶方法进行对比。此处采用的大气阻力模型为简单的球体阻力模型，相关参数如表2所示。其中大气密度通过msis-e-90大气模型得到。计算结果的位置误差如图5所示，可以看到在近地轨道，本发明所提方法的最大位置误差为10-6米。计算成本方面，本发明的方法相比于runge-kutta12(10)方法可以节约大约95％的计算时间。其中位置误差的计算方法如下：其中δri是采样点处的位置误差，ra是远地点相对于地心的距离，norbits为轨道个数。表2实例2采用的参数值参数数值rk12(10)的计算步长50s节点个数25计算步长(tp)(s)500迭代终止条件1e-7弹道系数0.01m2/kg本发明的方法不仅适用于低地球轨道，对其他的轨道类型也是适用的。本发明的前2阶形式，如表3所示。当用户允许的计算时间较长时，可采用0阶估计方法，反之，可采用1阶估计方法。表3低阶修正公式本发明还提供一种地球轨道中航天器的轨迹预测系统，所述预测系统包括：配点选取模块，用于从预测时间步长内选取多个配点。所述配点选取模块，具体包括：节点选取子模块，用于利用公式ξm＝cos[(m-1)π/(m-1)],m＝1,...,m，在[-1,1]区间中选取m个节点；配点选取子模块，用于分别令ξ＝ξ1,…,ξm，利用公式确定每个节点对应的第i个预测时间步长内的实际时间作为配点；其中，ξ1、ξm和ξm表示分别表示第1个、第m个和第m个节点，m表示节点的个数，ti和ti+1为第i个预测时间步长的左端点和右端点，t表示实际时间，ξ表示实际时间t对应的[-1,1]区间内的缩放后的时间。初始预测模块，用于根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息。修正模块，用于采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息。插值模块，用于对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明公开了一种地球轨道中航天器的轨迹预测方法，所述预测方法包括如下步骤：从预测时间步长内选取多个配点；根据航天器的初始位置和速度确定航天器在每个配点时的初始预测状态信息；采用反馈加速picard迭代方法对航天器在每个配点时的初始预测状态信息进行修正，获得航天器在每个配点时的修正后的状态信息；对航天器在每个配点时的修正后的状态信息进行插值，获得航天器在每个时刻的修正后的状态信息。本发明采用高精度的反馈加速picard迭代方法对初始预测状态信息进行修正，无需对矩阵求逆，在提高轨迹预测精度的同时，能够降低轨道计算中的复杂加速度项的计算耗时，并且结合配点法的思想，进一步减少了计算量，本发明为航天器提供了精确、实时的轨道预测方法。本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。当前第1页1 2 3