CN112319859B - 一种基于自主滤波阶切换的非线性卫星轨道确定方法 - Google Patents

一种基于自主滤波阶切换的非线性卫星轨道确定方法 Download PDF

Info

Publication number
CN112319859B
CN112319859B CN202011166646.9A CN202011166646A CN112319859B CN 112319859 B CN112319859 B CN 112319859B CN 202011166646 A CN202011166646 A CN 202011166646A CN 112319859 B CN112319859 B CN 112319859B
Authority
CN
China
Prior art keywords
order
state
filtering
time
value
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
Application number
CN202011166646.9A
Other languages
English (en)
Other versions
CN112319859A (zh
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.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
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 Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202011166646.9A priority Critical patent/CN112319859B/zh
Publication of CN112319859A publication Critical patent/CN112319859A/zh
Application granted granted Critical
Publication of CN112319859B publication Critical patent/CN112319859B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Remote Sensing (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Automation & Control Theory (AREA)
  • Filters That Use Time-Delay Elements (AREA)

Abstract

本发明公开了一种基于自主滤波阶切换的非线性卫星轨道确定方法，属于卫星轨道确定领域。利用高阶方法适合于非线性较强时的精确状态估计，及低阶方法适合于滤波器进入平稳阶段的高效率滤波，通过设计自适应滤波阶切换策略，在保证精度的前提下，在一个估计过程中根据需要动态改变滤波阶。在滤波过程中，通过计算测量新息并检验滤波一致性特性，当系统非线性较强时，采用高阶算法保证精度；当滤波过程进入稳定状态时，采用低阶算法保证效率。该方法不仅保持原有非线性扩展卡尔曼滤波采用高阶时的高精度特点，同时继承了低阶运行时的高效率，普遍适用于卫星、空间碎片等的轨道确定，可加速高精度卫星轨道确定的执行速度。

Description

一种基于自主滤波阶切换的非线性卫星轨道确定方法
技术领域
本发明属于卫星轨道确定领域，涉及一种基于自主滤波阶切换的非线性卫星轨道确定方法。
背景技术
近年来，随着人类空间活动的精细化和多样化，对在轨卫星轨道确定提出了更高的要求，包括高精度和高效率，低测量需求，低初值依赖度等。然而，对于目前大多数的常规滤波器，例如卡尔曼滤波器、扩展卡尔曼滤波器、粒子滤波器等传统的滤波器等，总是无法满足所有这些要求。潜在的能够同时满足上述性能要求的滤波器只有无味卡尔曼滤波器和基于泰勒多项式展开技术的非线性扩展卡尔曼滤波器，两者在非线性较强时展现出了特别优良的性能，但是当系统的非线性过于强时，两种滤波器都将因为过于承重的计算成本而失效。然而，失效的原因却截然不同，一般来说，无味卡尔曼滤波器是因为过多的采样点导致的重复性计算而失效，而作为一种半解析方法，非线性扩展卡尔曼滤波器则是由于滤波阶过高导致用于逼近的近似多项式过于复杂。另外，无味卡尔曼滤波器只能达到二阶精度，而理论上来说非线性扩展卡尔曼滤波器可以达到任意高阶。
另一方面，非线性扩展卡尔曼滤波器具有良好的推广性，即其低阶算法和高阶算法的算法执行过程完全相同，唯一的区别在于算法执行过程的展开阶选择。当选择滤波阶为1时，非线性扩展卡尔曼滤波器将退化为常规的扩展卡尔曼滤波器，显然，低阶算法可以获得良好的计算效率却有可能遭遇精度损失，相反，当选择滤波阶较高时，非线性扩展卡尔曼滤波器可以获得较高的精度却面临着巨大的计算压力。因此，如何克服非线性扩展卡尔曼滤波器在不同滤波阶时遭遇的问题而同时获得低级算法和高阶算法的优势，具有非常重要的实际意义。
发明内容
本发明的目的在于克服上述现有技术中，非线性扩展卡尔曼滤波器不能同时获得低级算法和高阶算法的优势的缺点，提供一种基于自主滤波阶切换的非线性卫星轨道确定方法
为了达到上述目的，本发明采用以下技术方案予以实现：
一种基于自主滤波阶切换的非线性卫星轨道确定方法，包括如下步骤：
步骤一，建立方程：分别建立航天器轨道动力学方程和观测方程；
步骤二，非线性状态和测量预测：分别对步骤一的轨道动力学方程和观测方程进行泰勒展开逼近，得到轨道状态的高阶多项式映射和观测方程的高阶多项式映射，并根据上一时刻状态统计信息计算得到状态预测值；
步骤三，计算测量新息，并进行卡方检测，根据结果对步骤二的两个高阶多项式映射过程进行自主阶切换；当系统非线性状态较强时，自主阶切换算法采用高阶算法进行滤波；当系统非线性状态较弱或滤波过程进入稳定状态时，自主阶切换算法采用低阶算法进行滤波；
步骤四，状态更新：设置时间参数、初始状态和当前时刻测量值，采用步骤三选择的自主阶切换算法，对步骤二中得到的状态预测值进行更新，得到目标时刻航天器的轨道状态。
优选地，步骤一建立的轨道动力学方程为：
Figure GDA0003339982680000021
式(1)中，x表示n维状态变量，t表示时间；假设t0时刻航天器的状态为x0，该常微分方程的解表示为x(t)＝Φ(t；t0,x0)；
航天器的观测方程为：
z＝h(x,t)+u (2)
式(2)中，z表示t时刻的观测量，x表示t时刻的状态值，u表示t时刻的观测噪声。
优选地，步骤二所述的非线性状态预测包括两个，分别为
1)对轨道运动进行非线性状态预测，具体为：
假设tk时刻航天器的标称状态为
Figure GDA0003339982680000031
初始偏差为δxk，此时多项式形式的状态变量表示为
Figure GDA0003339982680000032
将该多项式变量在区间[tk,tk+1]上进行积分，得到微分方程在tk+1时刻的解为
Figure GDA0003339982680000033
对该解在标称状态
Figure GDA0003339982680000034
处做泰勒展开，得到高阶展开式
Figure GDA0003339982680000035
其中，
Figure GDA0003339982680000036
表示状态邻域[xk+1]对初始邻域大小δxk的非线性依赖关系；
2)对航天器的观测进行非线性预测，具体为：
时间为tk+1时的测量方程为：
zk+1＝h(xk+1,tk+1)+uk+1 (3)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，uk+1表示tk+1时刻的观测噪声，在tk+1时刻，将方程(3)在标称状态
Figure GDA0003339982680000037
处做泰勒展开，得到高阶展开式，
Figure GDA0003339982680000038
式(4)中，p表示m维测量矢量的索引值，
Figure GDA0003339982680000039
表示tk+1时刻的标称观测矢量，系数
Figure GDA00033399826800000310
表示泰勒展开系数，
Figure GDA00033399826800000311
表示tk时刻第i个状态偏差分量，其阶数为γi；γ＝[γ1,…,γn]T，表示各状态偏差分量的阶数。
优选地，高阶展开式
Figure GDA0003339982680000041
的γ阶多项式的近似解表示为：
Figure GDA0003339982680000042
式(5)中，i表示n维状态矢量的索引值；
Figure GDA0003339982680000043
表示标称状态；δxk＝[δxk,1，…,δxk,n]T，表示初始偏差；
Figure GDA0003339982680000044
表示相应的泰勒展开式系数，γ＝[γ1,…,γn]T表示各状态偏差分量的指数。
优选地，步骤二所述的非线性状态预测过程中，轨道运动的预测状态均值
Figure GDA0003339982680000045
和协方差矩阵
Figure GDA0003339982680000046
为：
Figure GDA0003339982680000047
Figure GDA0003339982680000048
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure GDA0003339982680000049
表示过程噪声的协方差矩阵；式(7)中，
Figure GDA00033399826800000410
和
Figure GDA00033399826800000411
其余的
Figure GDA00033399826800000412
和
Figure GDA00033399826800000413
与对应的
Figure GDA00033399826800000414
均相等。
优选地，预测观测值的均值
Figure GDA00033399826800000415
和协方差矩阵Pk+1如下：
Figure GDA00033399826800000416
Figure GDA00033399826800000417
Figure GDA00033399826800000418
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure GDA00033399826800000419
表示观测噪声的协方差矩阵；
Figure GDA00033399826800000420
表示状态矢量和测量矢量的协态矩阵，
Figure GDA00033399826800000421
表示测量值协方差矩阵；E[]表示期望值；式(9)中，
Figure GDA00033399826800000422
其余的
Figure GDA00033399826800000423
Figure GDA00033399826800000424
均与对应的
Figure GDA00033399826800000425
相等。
优选地，步骤三所述的滤波阶切换具体为：
定义测量新息矢量vk+1和相应的协方差矩阵Sk+1为
Figure GDA0003339982680000051
Figure GDA0003339982680000052
Figure GDA0003339982680000053
定义tk+1时刻的卡方指标归一化值为
Figure GDA0003339982680000054
定义tk+1时刻的卡方指标的时间平均值
Figure GDA0003339982680000055
为
Figure GDA0003339982680000056
其中，uk+1表示tk+1时刻的观测噪声；Rk+1表示tk+1时刻观测噪声的协方差矩阵，E[]表示期望值，L表示数据采集的窗口宽度，L等于3。
优选地，步骤三所述的滤波阶切换的具体过程包括两种状态，分别为：
1)滤波器需要保持滤波阶γ，此时
Figure GDA0003339982680000057
滤波器不能正常工作，滤波器自动增加滤波阶，其中
Figure GDA0003339982680000058
可根据m和α值查卡方值表得到；
2)滤波器在滤波阶为γ时正常工作，此时
Figure GDA0003339982680000059
滤波器正常工作，之后进行降低滤波阶存在的可能验证过程。
优选地，降低滤波阶存在的可能性验证过程，包括：当
Figure GDA00033399826800000510
滤波器将尝试降低滤波阶γ；否则将保持滤波阶不变；其中，
Figure GDA00033399826800000511
和
Figure GDA00033399826800000512
通过卡方分布的概率表得到。
优选地，步骤四所述的状态更新的具体过程为：
通过观测设备得到tk+1时刻新的观测值
Figure GDA00033399826800000513
时，将其融合到预测值中，得到卡尔曼滤波器的增益Kk+1、状态的估计值
Figure GDA00033399826800000514
和协方差矩阵
Figure GDA00033399826800000515
Figure GDA0003339982680000061
Figure GDA0003339982680000062
Figure GDA0003339982680000063
与现有技术相比，本发明具有以下有益效果：
本发明公开了一种基于自主滤波阶切换的非线性卫星轨道确定方法，利用高阶方法适合于非线性较强时的精确状态估计，及低阶方法适合于滤波器进入平稳阶段的高效率滤波，通过设计自适应滤波阶切换策略，在保证精度的前提下，在一个估计过程中根据需要动态改变滤波阶。在滤波过程中，通过计算测量新息并检验滤波一致性特性，动态地设计了一个自适应自主滤波阶切换策略，当系统非线性较强时，采用高阶算法保证精度；当滤波过程进入稳定状态时，采用低阶算法保证效率。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法是一种混合算法，同时获取了低阶算法和高阶算法的优点。确保非线性扩展卡尔曼滤波器在保证精度的前提下，加快滤波过程的执行速度。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法，在保证精度的前提下，可以降低滤波器在滤波平稳阶段的滤波阶，显著减少计算成本，提高计算效率。克服了传统的非线性扩展卡尔曼滤波器在一个估计过程中无法改变滤波阶的困境，不仅保持原有非线性扩展卡尔曼滤波采用高阶时的高精度特点，同时继承了低阶运行时的高效率。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法普遍适用于卫星、空间碎片等的轨道确定，可加速高精度卫星轨道确定的执行速度。
进一步地，通过对滤波过程中测量新息的监测在保证滤波一致性的要求下，动态地调整滤波平稳阶段的滤波阶，大大降低了滤波器的计算成本，进一步优化了卫星轨道确定的执行速度。
附图说明
图1为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的位置误差示意图；
图2为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的速度误差示意图；
图3为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的位置分量标准差估计值示意图；
图4为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的速度分量标准差估计值示意图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
下面结合附图对本发明做进一步详细描述：
实施例1：
本实例公开一种基于多项式展开技术的高阶容错卫星轨道确定方法，具体实现步骤如下：
步骤一：建立无量纲的卫星轨道状态方程为
Figure GDA0003339982680000081
其中
Figure GDA0003339982680000082
分别表示无量纲的时间，卫星的无量纲位置矢量、速度矢量和角速度值，t，r，v，ω表示相应的有量纲值，
Figure GDA0003339982680000083
表示卫星无量纲位置矢量的范数；另外，卫星的标称角速度
Figure GDA0003339982680000084
表示卫星在半长轴为a＝42164km的轨道上运行的平均的轨道角速度，μe表示地球的引力常数。
步骤二：建立非线性测量方程，包括径向距离，赤经和赤纬，
Figure GDA0003339982680000085
其中u＝[u1,u2,u3]T表示测量噪声矢量，径向距离的测量噪声满足零均值标准差为1米的高斯分布，赤经和赤纬的测量噪声满足零均值标准差为1.745×10-6rad的高斯分布。同时，我们假设每个轨道可以获得12个等间距的测量值。
步骤三：非线性状态预测：
假设tk时刻航天器的标称状态为
Figure GDA0003339982680000086
初始偏差为δxk，因此初值
Figure GDA0003339982680000087
的一个邻域可表示为
Figure GDA0003339982680000088
被称为多项式状态变量，在区间[tk,tk+1]上进行积分后得到微分方程在tk+1时刻的解为
Figure GDA0003339982680000089
对该解在标称状态
Figure GDA00033399826800000810
处做泰勒展开可得到高阶展开式
Figure GDA00033399826800000811
其中
Figure GDA00033399826800000812
表示状态t1状态邻域[xk+1]对初始邻域大小δxk的非线性依赖关系，其γ阶多项式近似解可表示为
Figure GDA0003339982680000091
其中，i表示状态矢量的索引值，
Figure GDA0003339982680000092
表示标称状态；δxk＝[δxk,1,…,δxk,n]T表示初始偏差，
Figure GDA0003339982680000093
表示相应的泰勒展开式系数，γ＝[γ1,…,γn]T表示各状态偏差分量的指数。
进一步的，在tk+1时刻，使用γ阶多项式形式的解，即方程(5)，预测状态均值
Figure GDA0003339982680000094
和协方差矩阵
Figure GDA0003339982680000095
如下
Figure GDA0003339982680000096
Figure GDA0003339982680000097
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure GDA0003339982680000098
表示过程噪声的协方差矩阵；方程(7)中的系数
Figure GDA0003339982680000099
与
Figure GDA00033399826800000910
除了
Figure GDA00033399826800000911
和
Figure GDA00033399826800000912
之外，
Figure GDA00033399826800000913
Figure GDA00033399826800000914
与对应的
Figure GDA00033399826800000915
均相等。
类似的，根据测量方程(2)，tk+1的测量方程可表示为
zk+1＝h(xk+1,tk+1)+uk+1 (3)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，uk+1表示tk+1时刻的观测噪声，在tk+1时刻，将方程(3)在标称状态
Figure GDA00033399826800000916
处做泰勒展开可得到高阶展开式
Figure GDA00033399826800000917
其中p表示m维测量矢量的索引值，
Figure GDA00033399826800000918
表示对应于标称观测值，系数
Figure GDA00033399826800000919
表示泰勒展开系数，可通过将方程(5)代入方程(3)计算得到，m表示观测方程的数量；在tk+1时刻，使用γ阶多项式形式的解，即方程(4)，预测观测值的均值
Figure GDA00033399826800000920
如下：
Figure GDA0003339982680000101
Figure GDA0003339982680000102
Figure GDA0003339982680000103
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure GDA0003339982680000104
表示观测噪声的协方差矩阵，系数
Figure GDA0003339982680000105
与
Figure GDA0003339982680000106
除了
Figure GDA0003339982680000107
和
Figure GDA0003339982680000108
之外，
Figure GDA0003339982680000109
Figure GDA00033399826800001010
与对应的
Figure GDA00033399826800001011
均相等。
步骤四：滤波阶切换策略
定义测量新息矩阵vk+1和相应的协方差矩阵Sk+1为
Figure GDA00033399826800001012
Figure GDA00033399826800001013
Figure GDA00033399826800001014
进一步，qk+1定义为
Figure GDA00033399826800001015
定义时间平均的新息平方归一化值(NIS)
Figure GDA00033399826800001016
为
Figure GDA00033399826800001017
其中，L表示数据采集的窗口宽度，默认为3。
为了设计自适应的滤波阶切换策略，先提出两个假设
1)H0：滤波器在滤波阶为β时正常工作；
2)H1：滤波器需要保持滤波阶β。
在假设H0成立时，新息平方归一化平均值
Figure GDA00033399826800001018
是一个具有m个自由度卡方分布
Figure GDA00033399826800001019
其中重要性系数α＝99％，因此我们可以通过卡方分布的概率表得到对应于参数m和α的临界值
Figure GDA0003339982680000111
当
Figure GDA0003339982680000112
时，假设H0不成立而假设H1成立，即滤波器不正常工作，滤波器自动增加滤波阶。当
Figure GDA0003339982680000113
时，假设H1不成立而假设H0成立，即滤波器正常工作，在这种情况下，我们需要进一步验证是否存在降低滤波阶的可能性。类似的，我们提出以下两个假设：
1)
Figure GDA0003339982680000114
滤波器可以尝试降低滤波阶γ；
2)
Figure GDA0003339982680000115
滤波器在滤波阶为γ时不正常工作。
为了验证
Figure GDA0003339982680000116
和
Figure GDA0003339982680000117
我们需要围绕
Figure GDA0003339982680000118
定义一个严格的区间
Figure GDA0003339982680000119
一般情况，αl＝10％，αu＝75％，因此我们可以通过卡方分布的概率表得到对应于参数αl和αu的临界值
Figure GDA00033399826800001110
和
Figure GDA00033399826800001111
当
Figure GDA00033399826800001112
表明当前滤波器的一致性非常号，因此滤波器将尝试降低滤波阶γ，否则将保持滤波阶不变。
步骤五：非线性状态更新
通过观测设备得到tk+1时刻新的观测值
Figure GDA00033399826800001113
时，将其融合到预测值中，可以得到卡尔曼滤波器的增益，以及状态的估计值
Figure GDA00033399826800001114
和其协方差矩阵
Figure GDA00033399826800001115
Figure GDA00033399826800001116
Figure GDA00033399826800001117
Figure GDA00033399826800001118
步骤五：反复执行上述算法，得到任意时刻航天器的轨道状态。
对本发明方法进行验证：
地球同步轨道卫星的初始状态为
Figure GDA00033399826800001119
初始位置误差为100km，初始速度误差为0.01km/s。同时我们假设地球同步轨道卫星只有在晚上的12个小时可以进行观察，且每天可以均匀地获取十二个测量量。
图1和图2展示了分别由标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的位置和速度误差。仿真结果表明二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器具有相同的估计精度，最终的位置误差为4.9m，速度误差为0.0004m/s，相反，一阶非线性扩展卡尔曼滤波器由于在初始滤波阶段对状态误差的高估导致了精度的损失，如图3和图4所示。
另外，表1展示了三种滤波器消耗的计算成本(CPU时间)以及滤波器在平稳阶段的速度和位置误差。
表1：仿真得到的平均误差和CPU时长
滤波器 位置误差(千米) 速度误差(米/秒) CPU时长(秒)
HEKF-1 0.1994 0.0146 24.3
HEKF-2 0.0049 0.0004 38.7
AHEKF 0.0049 0.0004 26.1
表1的结果表明，本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器，在计算效率方面只比一阶非线性扩展卡尔曼滤波器慢一点，却获得了和二阶非线性扩展卡尔曼滤波器几乎相同的估计精度。
综上所述，为了改善现有的非线性扩展卡尔曼滤波器在低阶滤波和高阶滤波时存在的缺点，本发明通过动态地监测滤波过程中测量新息，验证滤波一致性特征，并以此为依据提出了一种自适应自主切换滤波阶的高效策略，该方法不仅保持原有非线性扩展卡尔曼滤波采用高阶时的高精度特点，同时继承了滤波器在低阶运行时的高效率特点。
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

Claims (10)

1.一种基于自主滤波阶切换的非线性卫星轨道确定方法，其特征在于，包括如下步骤：
步骤一，建立方程：分别建立航天器轨道动力学方程和观测方程；
步骤二，非线性状态和测量预测：分别对步骤一的轨道动力学方程和观测方程进行泰勒展开逼近，得到轨道状态的高阶多项式映射和观测方程的高阶多项式映射，并根据上一时刻状态统计信息计算得到状态预测值；
步骤三，计算测量新息，并进行卡方检测，根据结果对步骤二的两个高阶多项式映射过程进行自主阶切换；当系统非线性状态较强时，自主阶切换算法采用高阶算法进行滤波；当系统非线性状态较弱或滤波过程进入稳定状态时，自主阶切换算法采用低阶算法进行滤波；
步骤四，状态更新：设置时间参数、初始状态和当前时刻测量值，采用步骤三选择的自主阶切换算法，对步骤二中得到的状态预测值进行更新，得到目标时刻航天器的轨道状态。
2.根据权利要求1所述的非线性卫星轨道确定方法，其特征在于，步骤一建立的轨道动力学方程为：
Figure FDA0003339982670000011
式(1)中，x表示n维状态变量，t表示时间；假设t0时刻航天器的状态为x0，该常微分方程的解表示为x(t)＝Φ(t；t0,x0)；
航天器的观测方程为：
z＝h(x,t)+u (2)
式(2)中，z表示t时刻的观测量，x表示t时刻的状态值，u表示t时刻的观测噪声。
3.根据权利要求1所述的非线性卫星轨道确定方法，其特征在于，步骤二所述的非线性状态预测包括两个，分别为
1)对轨道运动进行非线性状态预测，具体为：
假设tk时刻航天器的标称状态为
Figure FDA0003339982670000021
初始偏差为δxk，此时多项式形式的状态变量表示为
Figure FDA0003339982670000022
将该多项式变量在区间[tk,tk+1]上进行积分，得到微分方程在tk+1时刻的解为
Figure FDA0003339982670000023
对该解在标称状态
Figure FDA0003339982670000024
处做泰勒展开，得到高阶展开式
Figure FDA0003339982670000025
其中，
Figure FDA0003339982670000026
表示状态邻域[xk+1]对初始邻域大小δxk的非线性依赖关系；
2)对航天器的观测进行非线性预测，具体为：
时间为tk+1时的测量方程为：
zk+1＝h(xk+1,tk+1)+uk+1 (3)
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，uk+1表示tk+1时刻的观测噪声，在tk+1时刻，将方程(3)在标称状态
Figure FDA0003339982670000027
处做泰勒展开，得到高阶展开式，
Figure FDA0003339982670000028
式(4)中，p表示m维测量矢量的索引值，
Figure FDA0003339982670000029
表示tk+1时刻的标称观测矢量，系数
Figure FDA00033399826700000213
表示泰勒展开系数，
Figure FDA00033399826700000210
表示tk时刻第i个状态偏差分量，其阶数为γi；γ＝[γ1,…,γn]T，表示各状态偏差分量的阶数。
4.根据权利要求3所述的非线性卫星轨道确定方法，其特征在于，高阶展开式
Figure FDA00033399826700000211
的γ阶多项式的近似解表示为：
Figure FDA00033399826700000212
式(5)中，i表示n维状态矢量的索引值；
Figure FDA0003339982670000031
表示标称状态；δxk＝[δxk,1，…,δxk,n]T，表示初始偏差；
Figure FDA0003339982670000032
表示相应的泰勒展开式系数，γ＝[γ1,…,γn]T表示各状态偏差分量的指数。
5.根据权利要求4所述的非线性卫星轨道确定方法，其特征在于，步骤二所述的非线性状态预测过程中，轨道运动的预测状态均值
Figure FDA0003339982670000033
和协方差矩阵
Figure FDA0003339982670000034
为：
Figure FDA0003339982670000035
Figure FDA0003339982670000036
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure FDA0003339982670000037
表示过程噪声的协方差矩阵；式(7)中，
Figure FDA0003339982670000038
和
Figure FDA0003339982670000039
其余的
Figure FDA00033399826700000310
和
Figure FDA00033399826700000311
与对应的
Figure FDA00033399826700000312
均相等。
6.根据权利要求5所述的非线性卫星轨道确定方法，其特征在于，预测观测值的均值
Figure FDA00033399826700000313
和协方差矩阵Pk+1如下：
Figure FDA00033399826700000314
Figure FDA00033399826700000315
Figure FDA00033399826700000316
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure FDA00033399826700000317
表示观测噪声的协方差矩阵；
Figure FDA00033399826700000318
表示状态矢量和测量矢量的协态矩阵，
Figure FDA00033399826700000319
表示测量值协方差矩阵；E[]表示期望值；式(9)中，
Figure FDA00033399826700000320
其余的
Figure FDA00033399826700000321
均与对应的
Figure FDA00033399826700000322
相等。
7.根据权利要求6所述的非线性卫星轨道确定方法，其特征在于，步骤三所述的滤波阶切换具体为：
定义测量新息矢量vk+1和相应的协方差矩阵Sk+1为
Figure FDA0003339982670000041
Figure FDA0003339982670000042
Figure FDA0003339982670000043
定义tk+1时刻的卡方指标归一化值为
Figure FDA0003339982670000044
定义tk+1时刻的卡方指标的时间平均值
Figure FDA0003339982670000045
为
Figure FDA0003339982670000046
其中，uk+1表示tk+1时刻的观测噪声；Rk+1表示tk+1时刻观测噪声的协方差矩阵，E[]表示期望值，L表示数据采集的窗口宽度，L等于3。
8.根据权利要求7所述的非线性卫星轨道确定方法，其特征在于，步骤三所述的滤波阶切换的具体过程包括两种状态，分别为：
1)滤波器需要保持滤波阶γ，此时
Figure FDA0003339982670000047
滤波器不能正常工作，滤波器自动增加滤波阶，其中
Figure FDA0003339982670000048
可根据m和α值查卡方值表得到；
2)滤波器在滤波阶为γ时正常工作，此时
Figure FDA0003339982670000049
滤波器正常工作，之后进行降低滤波阶存在的可能验证过程。
9.根据权利要求8所述的非线性卫星轨道确定方法，其特征在于，降低滤波阶存在的可能性验证过程，包括：当
Figure FDA00033399826700000410
滤波器将尝试降低滤波阶γ；否则将保持滤波阶不变；其中，
Figure FDA00033399826700000411
和
Figure FDA00033399826700000412
通过卡方分布的概率表得到。
10.根据权利要求6所述的非线性卫星轨道确定方法，其特征在于，步骤四所述的状态更新的具体过程为：
通过观测设备得到tk+1时刻新的观测值
Figure FDA00033399826700000413
时，将其融合到预测值中，得到卡尔曼滤波器的增益Kk+1、状态的估计值
Figure FDA0003339982670000051
和协方差矩阵
Figure FDA0003339982670000052
Figure FDA0003339982670000053
Figure FDA0003339982670000054
Figure FDA0003339982670000055
CN202011166646.9A 2020-10-27 2020-10-27 一种基于自主滤波阶切换的非线性卫星轨道确定方法 Active CN112319859B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011166646.9A CN112319859B (zh) 2020-10-27 2020-10-27 一种基于自主滤波阶切换的非线性卫星轨道确定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011166646.9A CN112319859B (zh) 2020-10-27 2020-10-27 一种基于自主滤波阶切换的非线性卫星轨道确定方法

Publications (2)

Publication Number Publication Date
CN112319859A CN112319859A (zh) 2021-02-05
CN112319859B true CN112319859B (zh) 2021-12-31

Family

ID=74296993

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011166646.9A Active CN112319859B (zh) 2020-10-27 2020-10-27 一种基于自主滤波阶切换的非线性卫星轨道确定方法

Country Status (1)

Country Link
CN (1) CN112319859B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116331523B (zh) * 2023-05-29 2023-08-25 哈尔滨工业大学 带大惯量旋转载荷卫星的未知参数辨识方法、装置及介质

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102305630B (zh) * 2011-05-17 2013-04-03 哈尔滨工业大学 基于扩展卡尔曼滤波的sar卫星自主定轨方法
US8938413B2 (en) * 2012-09-12 2015-01-20 Numerica Corp. Method and system for predicting a location of an object in a multi-dimensional space
CN103268407B (zh) * 2013-05-10 2016-12-28 航天东方红卫星有限公司 一种基于拉格朗日插值及卡尔曼滤波的轨道数据插值方法
CN103500455B (zh) * 2013-10-15 2016-05-11 北京航空航天大学 一种基于无偏有限冲击响应滤波器（ufir）的改进机动目标跟踪方法
US20190251215A1 (en) * 2018-02-15 2019-08-15 Regents Of The University Of Minnesota Accurate estimation of upper atmospheric density using satellite observations

Also Published As

Publication number Publication date
CN112319859A (zh) 2021-02-05

Similar Documents

Publication Publication Date Title
Fraser et al. Adaptive extended Kalman filtering strategies for spacecraft formation relative navigation
CN109974714B (zh) 一种Sage-Husa自适应无迹卡尔曼滤波姿态数据融合方法
CN109032176B (zh) 一种基于微分代数的地球同步轨道确定和参数确定方法
CN112319859B (zh) 一种基于自主滤波阶切换的非线性卫星轨道确定方法
CN111366156A (zh) 基于神经网络辅助的变电站巡检机器人导航方法及系统
Sabzevari et al. INS/GPS sensor fusion based on adaptive fuzzy EKF with sensitivity to disturbances
CN111428369A (zh) 一种空间目标碰撞预警结果置信度计算方法
CN112632756A (zh) 基于太阳敏感器的卫星地影自主预报方法及系统
CN110595434B (zh) 基于mems传感器的四元数融合姿态估计方法
Reihs et al. A method for perturbed initial orbit determination and correlation of radar measurements
Lin et al. Adaptive kernel auxiliary particle filter method for degradation state estimation
Garcia et al. Unscented Kalman filter and smoothing applied to attitude estimation of artificial satellites
US10889302B2 (en) Methods for monitoring the output performance of state estimators in navigation systems
Qian et al. An INS/DVL integrated navigation filtering method against complex underwater environment
CN111623764B (zh) 微纳卫星姿态估计方法
Valizadeh et al. Improvement of navigation accuracy using tightly coupled kalman filter
Qiu et al. Modified multiplicative quaternion cubature Kalman filter for attitude estimation
CN108981689B (zh) 基于dsp tms320c6748的uwb/ins组合导航系统
Li et al. Fault detection approach applied to inertial navigation system/air data system integrated navigation system with time‐offset
CN112415559A (zh) 一种基于多项式展开技术的高阶容错卫星轨道确定方法
CN115327587A (zh) 基于gnss定位信息的低轨卫星轨道误差修正方法及系统
Su et al. Quasi-consistent fusion navigation algorithm for DSS
CN113587926A (zh) 一种航天器空间自主交会对接相对导航方法
Zanette et al. RealSysId: A software tool for real-time aircraft model structure selection and parameter estimation
Nanda et al. Performance comparison of EKF and UKF for estimation of COE using Kozai mechanism

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