CN115204449B - 一种基于自适应勒让德皮卡迭代法的轨道预测方法 - Google Patents
一种基于自适应勒让德皮卡迭代法的轨道预测方法 Download PDFInfo
- Publication number
- CN115204449B CN115204449B CN202210582728.4A CN202210582728A CN115204449B CN 115204449 B CN115204449 B CN 115204449B CN 202210582728 A CN202210582728 A CN 202210582728A CN 115204449 B CN115204449 B CN 115204449B
- Authority
- CN
- China
- Prior art keywords
- legend
- iteration
- pick
- approximate solution
- information
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 72
- 230000003044 adaptive effect Effects 0.000 claims abstract description 13
- 238000005516 engineering process Methods 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 42
- 238000004364 calculation method Methods 0.000 claims description 16
- 230000008569 process Effects 0.000 claims description 15
- 230000010354 integration Effects 0.000 claims description 13
- 230000009466 transformation Effects 0.000 claims description 6
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 3
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000012804 iterative process Methods 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 230000011218 segmentation Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000008094 contradictory effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000000725 suspension Substances 0.000 description 1
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"
-
- 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
- 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/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
-
- 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Human Resources & Organizations (AREA)
- Operations Research (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Development Economics (AREA)
- Game Theory and Decision Science (AREA)
- Entrepreneurship & Innovation (AREA)
- Marketing (AREA)
- Computing Systems (AREA)
- Quality & Reliability (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Probability & Statistics with Applications (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于自适应勒让德皮卡迭代法的轨道预测方法,包括如下步骤:步骤1,根据航天器的初始位置与初始速度选择航天器的初始预估轨道;步骤2,设置预测时间步长以及预测时间步长内的节点数目;步骤3,基于带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息,并通过插值技术获得航天器在任意时刻的状态信息与速度信息。本发明应用于航天轨道设计领域,通过基于多次迭代后的近似轨迹获得理想的轨迹预测,不仅能提供足够的精度使预测迭代收敛加速,同时还大大减少了计算量。
Description
技术领域
本发明涉及航天轨道设计技术领域,具体是一种基于自适应勒让德皮卡迭代法的轨道预测方法。
背景技术
卫星轨道预测问题是航天领域一个基本且重要的实际问题,通常只能通过使用数值积分方法来解决。从历史上看,解决此问题的数值积分器可分为单步法或多步法。而这些只有使用较小的步长才能保证计算不发散,且必须采用高阶方法才能保证一定的计算精度,这导致轨迹预测过程需要大量的计算节点,且各个节点的计算难以并行处理。因此在应对长时间轨道计算和复杂动力学模型时,单步法、多步法需要消耗大量计算资源和在线计算时间,在卫星轨道预测方面面临精度差与效率低的问题。
发明内容
针对上述现有技术中的不足,本发明提供一种基于自适应勒让德皮卡迭代法的轨道预测方法,该方法不仅能提供足够的精度使预测迭代收敛加速,同时还大大减少了计算量。
为实现上述目的,本发明提供一种基于自适应勒让德皮卡迭代法的轨道预测方法,包括如下步骤:
步骤1,获得航天器的初始位置与初始速度并选择航天器的初始预估轨道;
步骤2,自适应设置预测时间步长以及预测时间步长内的节点数目;
步骤3,基于带积分反馈项的勒让德皮卡迭代法得到航天器在各时间步长内各节点上的状态信息与速度信息,并通过插值技术获得航天器在任意时刻的状态信息与速度信息。
在其中一个实施例,步骤2中,设置预测时间步长以及预测时间步长内的节点数目具体包括:
步骤2.1,预测时间步长根据真近点角的均匀划分来确定,从近地点开始将每个轨道周期分为X段,并设置每个预测时间步长上节点数目的最大数量为Y0,其中,X为奇数;
步骤2.2,令X=3;
步骤2.3,根据指定的初始条件,首先转换为经典轨道元素,然后用于计算开普勒轨道近地点的位置和速度;
步骤2.4,以近地点为起始点,沿着开普勒轨道向前传播轨道的1/X,利用Y(Y<Y0)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数:
若切比雪夫多项式的最后三个系数均小于给定误差的百分之一,则采用当前1/X的轨道作为预测时间步长,每个预测时间步长上节点的数目设为Y;
若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一,则令节点数目Y=Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y0,其中,若Y≥Y0则进行步骤2.5;
步骤2.5,令X=X+2后重复步骤2.3-2.4,直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一;
其中,若初始和/或最终时间不在预先计算的分段点上,则缩短第一个预测时间步长和/或最后一个预测时间步长。
在其中一个实施例,步骤3中,所述基于带积分反馈项的勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息,具体包括:
步骤3.1,获取航天器的初始状态信息x0和初始速度信息v0;
步骤3.2,将轨道预测问题转换为二阶微分方程组的求解问题,采用勒让德皮卡迭代法得到各节点的速度信息近似解以及状态信息近似解;
步骤3.3,采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值,再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及状态信息近似解,得到带积分误差反馈的速度信息近似解与状态信息近似解,反复迭代该过程直到获得满足给定精度的速度信息近似解与状态信息近似解,并将其作为航天器在各个离散的节点上的速度信息与状态信息。
在其中一个实施例,步骤3.2的具体包括:
将轨道预测问题转换为二阶微分方程组的求解的问题,即:
x″(t)=f(t,x(t),x′(t)),t0≤t≤tf (1)
式中,x″(t)表示位置状态函数的二阶导数即加速度函数,x(t)表示位置状态函数,x′(t)表示位置状态函数的一阶导数即速度函数,t表示时间节点,t0表示预测时间步长的初始时间,tf表示预测时间步长的终止时间,表示一个四元向量值函数;
在初始条件x(t0)=x0,x′(t0)=v0下,可以将式(1)改写成为一阶状态空间方程组,为:
将其转换为积分方程,为:
给定初始的状态信息x0、速度信息v0和以开普勒轨道作为初始轨道可通过勒让德皮卡迭代产生一列近似解xi(t)、vi(t),i=1,2,···,其中,二阶微分方程组的勒让德皮卡迭代公式为:
式中,vi(t)表示第i次勒让德皮卡迭代的速度信息近似解,xi-1(τ)表示第i-1次勒让德皮卡迭代的状态信息近似解,vi-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解,xi(t)表示第i次勒让德皮卡迭代的状态信息近似解,vi(τ)表示第i次勒让德皮卡迭代的速度信息近似解;
式中,s表示积分变量,g(·)为力函数,g(τ,x(τ))=f(t(τ),x(t(τ))以及x(τ)=x(t(τ);
采用勒让德多项式来逼近未知轨迹xi和沿轨迹xi-1的力函数g并使用勒让德多项式LN+1在(-1,1)上的零点来表示状态轨迹,其中,将k次勒让德多项式记为Lk,零点记为τj,j=0,1,···,N,通过勒让德多项式插值来近似力函数,即:
式中,Lk(τ)表示k次勒让德多项式;
式中,L′N+1(τj)表示(N+1)次勒让德多项式在点τj处的导数值;
通过结合g(τ,xi-1(τ))的插值及式(5),得到xi(t)、vi(t)的近似值,即式(5)可以重写为:
因此xi(τj)、vi(τj)的值可以通过后向离散勒让德变换来得到,为:
获取积分矩阵S,将式(10)转换为简化的矩阵向量形式,为:
基于勒让德皮卡迭代,得到状态信息与速度信息第i次迭代的近似解,为:
勒让德多项式的导数递推关系为:
(2k+1)Lk(s)=L′k+1(s)-L′k-1(s),k≥1
则有:
式中,k表示自然数,Lk(s)表示k次勒让德多项式,L′k+1(s)表示(k+1)次勒让德多项式的导数,L′k-1(s)表示(k-1)次勒让德多项式的导数,L0(s)表示零次勒让德多项式,L1(τ)表示一次勒让德多项式,L0(τ)表示零次勒让德多项式,Lk+1(τ)表示(k+1)次勒让德多项式,Lk-1(τ)表示(k-1)次勒让德多项式;
在其中一个实施例,积分矩阵S的确定过程具体为:
首先,令:
根据勒让德多项式和点{τj}定义两个矩阵,为:
再定义矩阵T,为:
用[vk]表示一个由向量vk组成的矩阵,因此,可以将式(7)转换为矩阵形式,为:
因此可以得到式(10)的矩阵形式,为:
最终得到积分矩阵,为:
其中,积分矩阵S仅取决于N。
在其中一个实施例,在步骤3.3中,得到带积分误差反馈的速度信息近似解与状态信息近似解的具体程为:
状态信息近似解与速度信息近似解的误差为:
又由多元函数的泰勒公式有:
式中,O(·)表示同阶无穷小;
通过速度信息的误差近似值修正勒让德皮卡迭代的近似解,从而得到带积分误差反馈的近似解,为:
在迭代过程中判断是否有max{||xi-xi-1||,||vi-vi-1||}<ε成立,ε为容差,若成立则退出迭代循环,否则继续迭代;
退出迭代时输出状态信息与速度信息的近似解xJ、vJ,即得到航天器在各个离散的节点上的状态信息与速度信息,最终通过插值技术获得航天器在任意时刻的状态信息与速度信息即能完成轨道预测;其中,J为退出迭代循环时的迭代次数。
在其中一个实施例,步骤1中,采用开普勒轨道作为初始预估轨道。
本发明提供的一种基于自适应勒让德皮卡迭代法的轨道预测方法,采用自适应策略确定时间步长及节点数目、基于勒让德高斯积分点利用勒让德正交多项式逼近力函数、结合积分误差反馈项等实现轨道预测的高精度,同时该方法易于实现并行计算,可大大降低轨道预测计算中的计算耗时,提高计算效率。本发明为卫星轨道预测提供了高精度、高效的轨道预测方法。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图示出的结构获得其他的附图。
图1为本发明实施例中轨道预测方法的流程图;
图2为本发明实施例中置预测时间步长以及预测时间步长内的节点数目的设置流程图;
图3为本发明实施例中状态信息与速度信息近似解的获取流程图;
图4为本发明实施例中Molniya轨道模拟示意图;
图5为本发明实施例中Hamilton量误差随递推时间的变化示意图。
本发明目的的实现、功能特点及优点将结合实施例,参照附图做进一步说明。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明的一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
另外,本发明各个实施例之间的技术方案可以相互结合,但是必须是以本领域普通技术人员能够实现为基础,当技术方案的结合出现相互矛盾或无法实现时应当认为这种技术方案的结合不存在,也不在本发明要求的保护范围之内。
如图1所示为本实施例公开的一种基于自适应勒让德皮卡迭代法的轨道预测方法,具体包括如下步骤:
步骤1,根据航天器的初始位置与初始速度选择航天器的初始预估轨道,本实施例中,采用开普勒轨道作为预估轨道。
步骤2,设置预测时间步长以及预测时间步长内的节点数目;
步骤3,基于带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息与速度信息,并通过插值技术获得航天器在任意时刻的状态信息与速度信息。
参考图2,在步骤2中,本实施例中采用一种自适应技术确定预测时间步长和每个预测时间步长上的节点数目,其具体实施过程为:
步骤2.1,预测时间步长根据真近点角的均匀划分来确定,从近地点开始将每个轨道周期分为X段,并设置每个预测时间步长上节点数目的最大数量为Y0,其中,X为奇数;
步骤2.2,令X=3;
步骤2.3,根据指定的初始条件,首先转换为经典轨道元素,然后用于计算开普勒轨道近地点的位置和速度;
步骤2.4,以近地点为起始点,沿着开普勒轨道向前传播轨道的1/X,利用Y(Y<Y0)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数:
若切比雪夫多项式的最后三个系数均小于给定误差的百分之一,则采用当前1/X的轨道作为预测时间步长,每个预测时间步长上节点的数目设为Y;
若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一,则令节点数目Y=Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y0,其中,若Y≥Y0则进行步骤2.5;
步骤2.5,令X=X+2后重复步骤2.3-2.4,直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一;
其中,若初始和/或最终时间不在预先计算的分段点上,则缩短第一个预测时间步长和/或最后一个预测时间步长。
下面结合具体的示例对预测时间步长以及节点数目的设置进行进一步的说明。
例如:从近地点开始将每个轨道周期分为奇数段,从三段开始,最开始将轨道周期划分为三段。每个预测时间步长的最大节点数40。根据用户指定的初始条件,首先转换为经典轨道元素,然后用于计算开普勒轨道近地点的位置和速度。以近地点为起始点,沿着开普勒轨道向前传播轨道的三分之一,利用10个切比雪夫点的切比雪夫多项式逼近该段轨道上的引力函数,如果切比雪夫多项式的最后三个系数均小于给定误差的百分之一,则采用三分之一的轨道作为预测时间步长,并将预测时间步长上节点的数目设为10。如果达不到误差要求,则将节点数目加倍。例如20、40个节点,直到达到误差要求或节点数目达到最大值40。若节点数目达到最大值还未满足误差要求,则减少预测时间步长为五分之一的轨道,节点数目再从10个节点开始计算,预测时间步长按下一个奇数分之一的轨道进行减少,依次类推,直至达到误差要求。一旦确定了预测时间步长和每个预测时间步长上的节点数,就可以开始逐段利用勒让德皮卡迭代进行轨道传播计算。
在步骤3,首先采用带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的速度信息与状态信息,在通过插值技术获得航天器在任意时刻的速度信息与状态信息,即完成了航天器的轨道预测。参考图3,其中采用带积分反馈项的自适应勒让德皮卡迭代法得到航天器在各个节点上的状态信息的具有实施过程为:
步骤3.1,获取航天器的初始状态信息x0和初始速度信息v0,可以直接获取或通过包含六要素的TLE两行轨道数据转换得到;
步骤3.2,将轨道预测问题转换为二阶微分方程组的求解的问题,采用勒让德皮卡迭代法得到各节点的速度信息近似解以及状态信息近似解;
步骤3.3,采用速度信息的近似解近似真实解,得到速度信息的误差近似值,再通过速度信息的误差近似值修正勒让德皮卡迭代的信息近似解以及状态信息近似解,得到带积分误差反馈的速度信息近似解与状态信息近似解,并将其作为航天器在各个离散的节点上的速度信息与状态信息。
本实施例中,步骤3.2的具体实施过程包括:
将轨道预测问题转换为二阶微分方程组的求解的问题,即:
x″(t)=f(t,x(t),x′(t)),t0≤t≤tf (1)
式中,x″(t)表示位置状态函数的二阶导数即加速度函数,x(t)表示位置状态函数,x′(t)表示位置状态函数的一阶导数即速度函数,t表示时间节点,t0表示预测时间步长的初始时间,tf表示预测时间步长的终止时间,表示一个四元向量值函数;
在初始条件x(t0)=x0,x′(t0)=v0下,可以将式(1)改写成为一阶状态空间方程组,为:
将其转换为积分方程,为:
给定初始的状态信息x0,速度信息v0和以开普勒轨道作为初始轨道,可通过勒让德皮卡迭代产生一列近似解xi(t)、vi(t),i=1,2,···,其中,二阶微分方程组的勒让德皮卡迭代公式为:
式中,vi(t)表示第i次勒让德皮卡迭代的速度信息近似解,xi-1(τ)表示第i-1次勒让德皮卡迭代的状态信息近似解,vi-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解,xi(t)表示第i次勒让德皮卡迭代的状态信息近似解,vi(τ)表示第i次勒让德皮卡迭代的速度信息近似解;
式中,s表示积分变量,g(·)为力函数,g(τ,x(τ))=f(t(τ),x(t(τ))以及x(τ)=x(t(τ);
采用勒让德多项式来逼近未知轨迹xi和沿轨迹xi-1的力函数g并使用勒让德多项式LN+1在(-1,1)上的零点来表示状态轨迹,其中,将k次勒让德多项式记为Lk,零点记为τj,j=0,1,···,N,通过勒让德多项式插值来近似力函数,即:
式中,Lk(τ)表示k次勒让德多项式;
式中,L′N+1(xj)表示(N+1)次勒让德多项式在点τj处的导数值;
通过结合g(τ,xi-1(τ))的插值及式(5),得到xi(t)、vi(t)的近似值,即式(5)可以重写为:
因此xi(τj)、vi(τj)的值可以通过后向离散勒让德变换来得到,为:
勒让德多项式的导数递推关系为:
(2k+1)Lk(s)=L′k+1(s)-L′k-1(s),k≥1
则有:
式中,k表示自然数,Lk(s)表示k次勒让德多项式,L′k+1(s)表示(k+1)次勒让德多项式的导数,L′k-1(s)表示(k-1)次勒让德多项式的导数,L0(s)表示零次勒让德多项式,L1(τ)表示一次勒让德多项式,L0(τ)表示零次勒让德多项式,Lk+1(τ)表示(k+1)次勒让德多项式,Lk-1(τ)表示(k-1)次勒让德多项式;
现有的勒让德皮卡迭代中都是采用两个常数矩阵来形成矩阵向量形式,其具有计算量大的弊端,基于此,本实施例提出了仅采用一个常数矩阵来形成矩阵向量形式。其具体实施过程为:
通过获取积分矩阵S,将式(10)转换为简化的矩阵向量形式,为:
基于勒让德皮卡迭代,得到状态信息与速度信息第i次迭代的近似解,为:
在具体实施过程中,积分矩阵S的确定过程具体为:
首先,令:
根据勒让德多项式和点{τj}定义两个矩阵,为:
再定义矩阵T,为:
用[vk]表示一个由向量vk组成的矩阵,因此,可以将式(7)转换为矩阵形式,为:
因此可以得到式(10)的矩阵形式,为:
最终得到积分矩阵,为:
其中,积分矩阵S仅取决于N,即本实施例中可以仅通过一个常值矩阵来实现,使得LPI方法更容易理解,并且节省了一半矩阵向量乘积的计算量。
在步骤3.3中,采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值,再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及状态信息近似解,得到带积分误差反馈的速度信息近似解与状态信息近似解的具体实施过程为:
在步骤3.2中得到状态信息与速度信息第i次迭代的近似解的前提下,此时状态信息近似解与速度信息近似解的误差为:
又由多元函数的泰勒公式有:
式中,O(·)表示同阶无穷小;
通过速度信息的误差近似值修正勒让德皮卡迭代的近似解,从而得到带积分误差反馈的近似解,为:
迭代过程中判断是否有max{||xi-xi-1||,||vi-vi-1||}<ε成立,ε为容差,若成立则退出迭代循环,否则继续迭代;
退出迭代时输出状态信息与速度信息的近似解xJ、vJ,即得到航天器在各个离散的节点上的状态信息与速度信息,最终通过插值技术获得航天器在任意时刻的状态信息与速度信息即能完成轨道预测。其中,J为退出迭代循环时的迭代次数。
对于轨道问题,当使用复杂的引力模型时,雅可比矩阵的精确计算是非常昂贵的。通常,不需要超过4位精度的J,因为积分误差项本身就是随着迭代而衰减得很快。在近地球轨道的情况下,我们发现使用基于前次迭代的近似轨迹获得的理想二体问题的雅可比矩阵就可提供足够的精度使得迭代解收敛加速,使得计算大大减少。
下面结合具体的实例对本发明作出进一步的说明。
参考图4,模拟一个Molniya轨道。在实例中不考虑空气阻力,考虑球谐重力场70阶模型(EMG2008)进行轨道递推,向前递推预测10个开普勒周期,约5.9×105s。为评估本方法的计算精度,对数值结果Hamilton量的变化进行分析。数值结果中Hamilton量的相对误差为其中H0是系统在初始时刻的Hamilton量,H(t)为t时刻Hamilton量,e(t)随递推时间的变化如图5所示。Hamilton量的相对计算误差控制在3×10-14内,达到了机器计算精度。该实例如表1所示。
表1实例1采用的轨道参数及自适应确定的时间步长与节点数目参数
参数 | 数值 |
<![CDATA[初始位置r<sub>0</sub>(km)]]> | <![CDATA[[9.2210×10<sup>3</sup>,-2.5413×10<sup>-12</sup>,-4.9875×10<sup>-12</sup>]<sup>T</sup>]]> |
<![CDATA[初始速度v<sub>0</sub>(km/s)]]> | <![CDATA[[3.0433×10<sup>-15</sup>,3.9146,7.6829]<sup>T</sup>]]> |
偏心率 | 0.72 |
倾斜角 | 63° |
一个周期分段数 | 5 |
节点数目 | 40 |
迭代中止误差 | <![CDATA[10<sup>-15</sup>]]> |
以上所述仅为本发明的优选实施例,并非因此限制本发明的专利范围,凡是在本发明的发明构思下,利用本发明说明书及附图内容所作的等效结构变换,或直接/间接运用在其他相关的技术领域均包括在本发明的专利保护范围内。
Claims (6)
1.一种基于自适应勒让德皮卡迭代法的轨道预测方法,其特征在于,包括如下步骤:
步骤1,获得航天器的初始位置与初始速度并选择航天器的初始预估轨道;
步骤2,自适应设置预测时间步长以及预测时间步长内的节点数目;
步骤3,基于带积分反馈项的勒让德皮卡迭代法得到航天器在各时间步长内各节点上的位置信息与速度信息,并通过插值技术获得航天器在任意时刻的位置信息与速度信息;
步骤3中,所述基于带积分反馈项的勒让德皮卡迭代法得到航天器在各个节点上的位置信息与速度信息,具体包括:
步骤3.1,获取航天器的初始位置信息x0和初始速度信息v0;
步骤3.2,将轨道预测问题转换为二阶微分方程组的求解问题,采用勒让德皮卡迭代法得到各节点的速度信息近似解以及位置信息近似解;
步骤3.3,采用速度信息的近似解、位置信息的近似解得到速度信息的误差近似值,再通过速度信息的误差近似值修正勒让德皮卡迭代的速度信息近似解以及位置信息近似解,得到带积分误差反馈的速度信息近似解与位置信息近似解,反复迭代该过程直到获得满足给定精度的速度信息近似解与位置信息近似解,并将其作为航天器在各个离散的节点上的速度信息与位置信息;
步骤3.2的具体包括:
将轨道预测问题转换为二阶微分方程组的求解的问题,即:
x″(t)=f(t,x(t),x′(t)),t0≤t≤tf (1)
式中,x″(t)表示位置函数的二阶导数即加速度函数,x(t)表示位置函数,x′(t)表示位置函数的一阶导数即速度函数,t表示时间节点,t0表示预测时间步长的初始时间,tf表示预测时间步长的终止时间,表示一个四元向量值函数;
在初始条件x(t0)=x0,x′(t0)=v0下,可以将式(1)改写成为一阶状态空间方程组,为:
将其转换为积分方程,为:
给定初始的位置信息x0、速度信息v0和以开普勒轨道作为初始轨道可通过勒让德皮卡迭代产生一列近似解xi(t)、vi(t),i=1,2,···,其中,二阶微分方程组的勒让德皮卡迭代公式为:
式中,vi(t)表示第i次勒让德皮卡迭代的速度信息近似解,xi-1(τ)表示第i-1次勒让德皮卡迭代的位置信息近似解,vi-1(τ)表示第i-1次勒让德皮卡迭代的速度信息近似解,xi(t)表示第i次勒让德皮卡迭代的位置信息近似解,vi(τ)表示第i次勒让德皮卡迭代的速度信息近似解;
式中,s表示积分变量,g(·)为力函数,g(τ,x(τ))=f(t(τ),x(t(τ)))以及x(τ)=x(t(τ));
采用勒让德多项式来逼近未知轨迹xi和沿轨迹xi-1的力函数g并使用勒让德多项式LN+1在(-1,1)上的零点来表示状态轨迹,其中,将k次勒让德多项式记为Lk,零点记为τj,j=0,1,···,N,通过勒让德多项式插值来近似力函数,即:
式中,Lk(τ)表示k次勒让德多项式;
式中,L′N+1(τj)表示(N+1)次勒让德多项式在点τj处的导数值;
通过结合g(τ,xi-1(τ))的插值及式(5),得到xi(t)、vi(t)的近似值,即式(5)可以重写为:
因此xi(τj)、vi(τj)的值可以通过后向离散勒让德变换来得到,为:
获取积分矩阵S,将式(10)转换为简化的矩阵向量形式,为:
基于勒让德皮卡迭代,得到位置信息与速度信息第i次迭代的近似解,为:
2.根据权利要求1所述基于自适应勒让德皮卡迭代法的轨道预测方法,其特征在于,步骤2中,设置预测时间步长以及预测时间步长内的节点数目具体包括:
步骤2.1,预测时间步长根据真近点角的均匀划分来确定,从近地点开始将每个轨道周期分为X段,并设置每个预测时间步长上节点数目的最大数量为Y0,其中,X为奇数;
步骤2.2,令X=3;
步骤2.3,根据指定的初始条件,首先转换为经典轨道元素,然后用于计算开普勒轨道近地点的位置和速度;
步骤2.4,以近地点为起始点,沿着开普勒轨道向前传播轨道的1/X,利用Y(Y<Y0)个节点的切比雪夫多项式逼近该1/X段开普勒轨道上的引力函数:
若切比雪夫多项式的最后三个系数均小于给定误差的百分之一,则采用当前1/X的轨道作为预测时间步长,每个预测时间步长上节点的数目设为Y;
若切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一,则令节点数目Y=Y+c后重复步骤2.4直至切比雪夫多项式的最后三个系数均小于给定误差的百分之一或Y≥Y0,其中,若Y≥Y0则进行步骤2.5;
步骤2.5,令X=X+2后重复步骤2.3-2.4,直至切比雪夫多项式的最后三个系数大于或等于给定误差的百分之一;
其中,若初始和/或最终时间不在预先计算的分段点上,则缩短第一个预测时间步长和/或最后一个预测时间步长。
勒让德多项式的导数递推关系为:
(2k+1)Lk(s)=L′k+1(s)-L′k-1(s),k≥1
则有:
式中,k表示自然数,Lk(s)表示k次勒让德多项式,L′k+1(s)表示(k+1)次勒让德多项式的导数,L′k-1(s)表示(k-1)次勒让德多项式的导数,L0(s)表示零次勒让德多项式,L1(τ)表示一次勒让德多项式,L0(τ)表示零次勒让德多项式,Lk+1(τ)表示(k+1)次勒让德多项式,Lk-1(τ)表示(k-1)次勒让德多项式;
5.根据权利要求1至4任一项所述基于自适应勒让德皮卡迭代法的轨道预测方法,其特征在于,在步骤3.3中,得到带积分误差反馈的速度信息近似解与位置信息近似解的具体程为:
位置信息近似解与速度信息近似解的误差为:
又由多元函数的泰勒公式有:
式中,O(·)表示同阶无穷小;
通过速度信息的误差近似值修正勒让德皮卡迭代的近似解,从而得到带积分误差反馈的近似解,为:
在迭代过程中判断是否有max{||xi-xi-1||,||vi-vi-1||}<ε成立,ε为容差,若成立则退出迭代循环,否则继续迭代;
退出迭代时输出位置信息与速度信息的近似解xJ、vJ,即得到航天器在各个离散的节点上的位置信息与速度信息,最终通过插值技术获得航天器在任意时刻的位置信息与速度信息即能完成轨道预测;其中,J为退出迭代循环时的迭代次数。
6.根据权利要求1至4任一项所述基于自适应勒让德皮卡迭代法的轨道预测方法,其特征在于,步骤1中,采用开普勒轨道作为初始预估轨道。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210582728.4A CN115204449B (zh) | 2022-05-26 | 2022-05-26 | 一种基于自适应勒让德皮卡迭代法的轨道预测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210582728.4A CN115204449B (zh) | 2022-05-26 | 2022-05-26 | 一种基于自适应勒让德皮卡迭代法的轨道预测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115204449A CN115204449A (zh) | 2022-10-18 |
CN115204449B true CN115204449B (zh) | 2023-04-25 |
Family
ID=83576222
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210582728.4A Active CN115204449B (zh) | 2022-05-26 | 2022-05-26 | 一种基于自适应勒让德皮卡迭代法的轨道预测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115204449B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116502774B (zh) * | 2023-06-26 | 2023-09-12 | 南京信息工程大学 | 一种基于时间序列分解和勒让德投影的时间序列预测方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107194039A (zh) * | 2017-04-26 | 2017-09-22 | 西北工业大学 | 一种基于改进高斯伪谱法的空间柔性系统展开控制方法 |
CN107977486A (zh) * | 2017-11-06 | 2018-05-01 | 北京宇航系统工程研究所 | 一种地球扰动引力场球冠谐模型阶次扩展方法及系统 |
Family Cites Families (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103033827B (zh) * | 2012-12-13 | 2015-09-09 | 中国航天科工信息技术研究院 | 一种卫星位置解算方法 |
CN112141366B (zh) * | 2020-09-23 | 2022-03-25 | 西北工业大学 | 一种地球轨道中航天器的轨迹预测方法及系统 |
-
2022
- 2022-05-26 CN CN202210582728.4A patent/CN115204449B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107194039A (zh) * | 2017-04-26 | 2017-09-22 | 西北工业大学 | 一种基于改进高斯伪谱法的空间柔性系统展开控制方法 |
CN107977486A (zh) * | 2017-11-06 | 2018-05-01 | 北京宇航系统工程研究所 | 一种地球扰动引力场球冠谐模型阶次扩展方法及系统 |
Also Published As
Publication number | Publication date |
---|---|
CN115204449A (zh) | 2022-10-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Ozaki et al. | Tube stochastic optimal control for nonlinear constrained trajectory optimization problems | |
US8332451B2 (en) | Programmable CORDIC Processor | |
Tanner et al. | From EM to data augmentation: the emergence of MCMC Bayesian computation in the 1980s | |
US5479571A (en) | Neural node network and model, and method of teaching same | |
CN115204449B (zh) | 一种基于自适应勒让德皮卡迭代法的轨道预测方法 | |
Haykin et al. | Optimum nonlinear filtering | |
US20180013439A1 (en) | Method to perform convolutions between arbitrary vectors using clusters of weakly coupled oscillators | |
US8732223B2 (en) | Deriving a function that represents data points | |
US8447799B2 (en) | Process for QR transformation using a CORDIC processor | |
US6223194B1 (en) | Adaptive filter, step size control method thereof, and record medium therefor | |
Sveier et al. | Applied Runge–Kutta–Munthe-Kaas integration for the quaternion kinematics | |
Roa et al. | Reduced nonlinear model for orbit uncertainty propagation and estimation | |
Probe et al. | Terminal Convergence Approximation Modified Chebyshev Picard Iteration for efficient numerical integration of orbital trajectories | |
US8452830B2 (en) | Programmable CORDIC processor with stage re-use | |
Galrinho et al. | On estimating initial conditions in unstructured models | |
CN114897254A (zh) | 一种双重自适应切比雪夫皮卡迭代法的轨道预测方法 | |
CN113872567A (zh) | 一种基于核函数的复数仿射投影自适应信号处理方法 | |
Zucchelli et al. | A Bayesian Approach to Low-Thrust Maneuvering Spacecraft Tracking | |
Im et al. | A least squares-type density estimator using a polynomial function | |
Weisman | Nonlinear transformations and filtering theory for space operations | |
CN116776474A (zh) | 一种航天器轨道的并行自适应计算方法及系统 | |
Anzalone et al. | Numerical Integration Techniques in Orbital Mechanics Applications | |
Lee | Stochastic attitude smoothing and interpolation on SO (3) with matrix Fisher distributions | |
Tahavori et al. | Time-weighted balanced stochastic model reduction | |
CN116661021A (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 |