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

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

Info

Publication number
CN112319859A
CN112319859A CN202011166646.9A CN202011166646A CN112319859A CN 112319859 A CN112319859 A CN 112319859A CN 202011166646 A CN202011166646 A CN 202011166646A CN 112319859 A CN112319859 A CN 112319859A
Authority
CN
China
Prior art keywords
order
state
filtering
time
nonlinear
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.)
Granted
Application number
CN202011166646.9A
Other languages
English (en)
Other versions
CN112319859B (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 BDA0002746010560000021
式(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 BDA0002746010560000031
初始偏差为δxk，此时多项式形式的状 态变量表示为
Figure BDA0002746010560000032
将该多项式变量在区间[tk,tk+1]上进行积分，得到微分 方程在tk+1时刻的解为
Figure BDA0002746010560000033
对该解在标称状态
Figure BDA0002746010560000039
处做泰勒展开，得到高阶展开式
Figure BDA0002746010560000034
其中，δxk+1＝Txk+1(δxk)，表示状态邻域[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 BDA0002746010560000035
处做泰勒展开，得到 高阶展开式，
Figure BDA0002746010560000036
式(4)中，p表示m维测量矢量的索引值，
Figure BDA0002746010560000037
表示tk+1时刻的 标称观测矢量，系数
Figure BDA00027460105600000310
表示泰勒展开系数，
Figure BDA0002746010560000038
表示tk时刻第i个状态偏差 分量，其阶数为γi；γ＝[γ1,L,γn]T，表示各状态偏差分量的阶数。
优选地，高阶展开式
Figure BDA0002746010560000041
的γ阶多项式的近似解表示为：
Figure BDA0002746010560000042
式(5)中，i表示n维状态矢量的索引值；
Figure BDA0002746010560000043
表示标称状态； δxk＝[δxk,1，L,δxk,n]T，表示初始偏差；
Figure BDA0002746010560000044
表示相应的泰勒展开式系数， γ＝[γ1,L,γn]T表示各状态偏差分量的指数。
优选地，步骤二所述的非线性状态预测过程中，轨道运动的预测状态均值
Figure BDA0002746010560000045
和协方差矩阵
Figure BDA0002746010560000046
为：
Figure BDA0002746010560000047
Figure BDA0002746010560000048
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure RE-GDA0002826276090000048
表示过程噪声的协方差矩阵；式(7)中，
Figure RE-GDA0002826276090000049
和
Figure RE-GDA00028262760900000410
其余的
Figure RE-GDA00028262760900000411
和
Figure RE-GDA00028262760900000412
与对应的
Figure RE-GDA00028262760900000413
均相等。
优选地，预测观测值的均值
Figure BDA00027460105600000415
和协方差矩阵Pk+1如下：
Figure BDA00027460105600000416
Figure BDA00027460105600000417
Figure BDA00027460105600000418
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure BDA00027460105600000419
表示观 测噪声的协方差矩阵；
Figure BDA00027460105600000420
表示状态矢量和测量矢量的协态矩阵，
Figure BDA00027460105600000421
表示 测量值协方差矩阵；E[]表示期望值；式(9)中，
Figure BDA00027460105600000422
其 余的
Figure BDA00027460105600000423
均与对应的
Figure BDA00027460105600000424
相等。
优选地，步骤三所述的滤波阶切换具体为：
定义测量新息矢量vk+1和相应的协方差矩阵Sk+1为
Figure BDA0002746010560000051
Figure BDA0002746010560000052
Figure BDA0002746010560000053
定义tk+1时刻的卡方指标归一化值为
Figure BDA0002746010560000054
定义tk+1时刻的卡方指标的时间平均值
Figure BDA0002746010560000055
为
Figure BDA0002746010560000056
其中，uk+1表示tk+1时刻的观测噪声；Rk+1表示tk+1时刻观测噪声的协方差矩阵， E[]表示期望值，L表示数据采集的窗口宽度，L等于3。
优选地，步骤三所述的自适应滤波切换的具体过程包括两种状态，分别为：
1)滤波器需要保持滤波阶γ，此时
Figure BDA0002746010560000057
滤波器不能正常工作，滤波器自 动增加滤波阶，其中
Figure BDA0002746010560000058
可根据m和α值查卡方值表得到；
2)滤波器在滤波阶为γ时正常工作，此时
Figure BDA0002746010560000059
滤波器正常工作，之后进 行降低滤波阶存在的可能验证过程。
优选地，降低滤波阶存在的可能性验证过程，包括：当
Figure BDA00027460105600000510
滤 波器将尝试降低滤波阶γ；否则将保持滤波阶不变；其中，
Figure BDA00027460105600000511
和
Figure BDA00027460105600000512
通过卡方 分布的概率表得到。
优选地，步骤四所述的非线性状态更新的具体过程为：
通过观测设备得到tk+1时刻新的观测值
Figure BDA00027460105600000513
时，将其融合到预测值中，得到卡 尔曼滤波器的增益Kk+1、状态的估计值
Figure BDA00027460105600000514
和协方差矩阵
Figure BDA00027460105600000515
Figure BDA0002746010560000061
Figure BDA0002746010560000062
Figure BDA0002746010560000063
与现有技术相比，本发明具有以下有益效果：
本发明公开了一种基于自主滤波阶切换的非线性卫星轨道确定方法，利用高 阶方法适合于非线性较强时的精确状态估计，及低阶方法适合于滤波器进入平稳 阶段的高效率滤波，通过设计自适应滤波阶切换策略，在保证精度的前提下，在 一个估计过程中根据需要动态改变滤波阶。在滤波过程中，通过计算测量新息并 检验滤波一致性特性，动态地设计了一个自适应自主滤波阶切换策略，当系统非 线性较强时，采用高阶算法保证精度；当滤波过程进入稳定状态时，采用低阶算 法保证效率。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法是一 种混合算法，同时获取了低阶算法和高阶算法的优点。确保非线性扩展卡尔曼滤 波器在保证精度的前提下，加快滤波过程的执行速度。本发明提出的基于自主滤 波阶切换的非线性卫星轨道确定方法，在保证精度的前提下，可以降低滤波器在 滤波平稳阶段的滤波阶，显著减少计算成本，提高计算效率。克服了传统的非线 性扩展卡尔曼滤波器在一个估计过程中无法改变滤波阶的困境，不仅保持原有非 线性扩展卡尔曼滤波采用高阶时的高精度特点，同时继承了低阶运行时的高效率。 本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法普遍适用于卫星、 空间碎片等的轨道确定，可加速高精度卫星轨道确定的执行速度。
进一步地，通过对滤波过程中测量新息的监测在保证滤波一致性的要求下， 动态地调整滤波平稳阶段的滤波阶，大大降低了滤波器的计算成本，进一步优化 了卫星轨道确定的执行速度。
附图说明
图1为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶 切换的非线性扩展卡尔曼滤波器得到的位置误差示意图；
图2为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶 切换的非线性扩展卡尔曼滤波器得到的速度误差示意图；
图3为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶 切换的非线性扩展卡尔曼滤波器得到的位置分量标准差估计值示意图；
图4为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶 切换的非线性扩展卡尔曼滤波器得到的速度分量标准差估计值示意图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例 中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述 的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的 实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实 施例，都应当属于本发明保护的范围。
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语是用于区 别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数 据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图 示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任 何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、 方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没 有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
下面结合附图对本发明做进一步详细描述：
实施例1：
本实例公开一种基于多项式展开技术的高阶容错卫星轨道确定方法，具体实 现步骤如下：
步骤一：建立无量纲的卫星轨道状态方程为
Figure BDA0002746010560000081
其中
Figure BDA0002746010560000082
分别表示无量纲的时间，卫星的无量纲位置矢量、速度矢量和角速度值，t，r，v，ω表示相应的有量纲值，
Figure BDA00027460105600000812
表示卫星 无量纲位置矢量的范数；另外，卫星的标称角速度
Figure BDA0002746010560000083
表示卫星在半长轴为 a＝42164km的轨道上运行的平均的轨道角速度，μe表示地球的引力常数。
步骤二：建立非线性测量方程，包括径向距离，赤经和赤纬，
Figure BDA0002746010560000084
其中u＝[u1,u2,u3]T表示测量噪声矢量，径向距离的测量噪声满足零均值标准 差为1米的高斯分布，赤经和赤纬的测量噪声满足零均值标准差为1.745×10-6 rad的高斯分布。同时，我们假设每个轨道可以获得12个等间距的测量值。
步骤三：非线性状态预测：
假设tk时刻航天器的标称状态为
Figure BDA0002746010560000085
初始偏差为δxk，因此初值
Figure BDA0002746010560000086
的一个邻 域可表示为
Figure BDA0002746010560000087
被称为多项式状态变量，在区间[tk,tk+1]上进行积分后得 到微分方程在tk+1时刻的解为
Figure BDA0002746010560000088
对该解在标称状态
Figure BDA0002746010560000089
处做泰勒展开 可得到高阶展开式
Figure BDA00027460105600000810
其中
Figure BDA00027460105600000811
表示状态t1状态邻域[xk+1] 对初始邻域大小δxk的非线性依赖关系，其γ阶多项式近似解可表示为
Figure BDA0002746010560000091
其中，i表示状态矢量的索引值，
Figure BDA0002746010560000092
表示标称状态； δxk＝[δxk,1,L,δxk,n]T表示初始偏差，
Figure BDA0002746010560000093
表示相应的泰勒展开式系数， γ＝[γ1,L,γn]T表示各状态偏差分量的指数。
进一步的，在tk+1时刻，使用γ阶多项式形式的解，即方程(5)，预测状态均 值
Figure BDA0002746010560000094
和协方差矩阵
Figure BDA0002746010560000095
如下
Figure BDA0002746010560000096
Figure BDA0002746010560000097
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure RE-GDA0002826276090000098
表示过程噪声的协方差矩阵；方程(7)中的系数
Figure RE-GDA0002826276090000099
与
Figure RE-GDA00028262760900000910
除了
Figure RE-GDA00028262760900000911
和
Figure RE-GDA00028262760900000912
之外，
Figure RE-GDA00028262760900000913
与对应的
Figure RE-GDA00028262760900000914
均相等。
类似的，根据测量方程(2)，tk+1的测量方程可表示为
Figure BDA00027460105600000915
其中，zk+1表示tk+1时刻的观测量，xk+1表示tk+1时刻的状态预测值，uk+1表示 tk+1时刻的观测噪声，在tk+1时刻，将方程(3)在标称状态
Figure BDA00027460105600000916
处做泰勒展开可得 到高阶展开式
Figure BDA00027460105600000917
其中p表示m维测量矢量的索引值，
Figure BDA00027460105600000918
表示对应于标称观 测值，系数
Figure BDA00027460105600000919
表示泰勒展开系数，可通过将方程(5)代入方程(3)计算得到， m表示观测方程的数量；在tk+1时刻，使用γ阶多项式形式的解，即方程(4)，预 测观测值的均值
Figure BDA00027460105600000920
如下：
Figure BDA0002746010560000101
Figure BDA0002746010560000102
Figure BDA0002746010560000103
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure BDA0002746010560000104
表示观 测噪声的协方差矩阵，系数
Figure BDA0002746010560000105
与
Figure BDA0002746010560000106
除了
Figure BDA0002746010560000107
和
Figure BDA0002746010560000108
之外，
Figure BDA0002746010560000109
与对应的
Figure BDA00027460105600001010
均相等。
步骤四：滤波阶切换策略
定义测量新息矩阵vk+1和相应的协方差矩阵Sk+1为
Figure BDA00027460105600001011
Figure BDA00027460105600001012
Figure BDA00027460105600001013
进一步，qk+1定义为
Figure BDA00027460105600001014
定义时间平均的新息平方归一化值(NIS)
Figure BDA00027460105600001015
为
Figure BDA00027460105600001016
其中，L表示数据采集的窗口宽度，默认为3。
为了设计自适应的滤波阶切换策略，先提出两个假设
1)H0：滤波器在滤波阶为γ时正常工作；
2)H1：滤波器需要保持滤波阶γ。
在假设H0成立时，新息平方归一化平均值
Figure BDA00027460105600001017
是一个具有m个自由度卡方分 布
Figure BDA00027460105600001018
其中重要性系数α＝99％，因此我们可以通过卡方分布的概率表得到对应 于参数m和α的临界值
Figure BDA0002746010560000111
当
Figure BDA0002746010560000112
时，假设H0不成立而假设H1成立，即滤波 器不正常工作，滤波器自动增加滤波阶。当
Figure BDA0002746010560000113
时，假设H1不成立而假设H0成立，即滤波器正常工作，在这种情况下，我们需要进一步验证是否存在降低滤 波阶的可能性。类似的，我们提出以下两个假设：
1)
Figure BDA0002746010560000114
滤波器可以尝试降低滤波阶γ；
2)
Figure BDA0002746010560000115
滤波器在滤波阶为γ时不正常工作。
为了验证
Figure BDA0002746010560000116
和
Figure BDA0002746010560000117
我们需要围绕
Figure BDA0002746010560000118
定义一个严格的区间
Figure BDA0002746010560000119
一般情况，αl＝10％，αu＝75％，因此我们可以通过卡方分布的概率表得到对应于 参数αl和αu的临界值
Figure BDA00027460105600001110
和
Figure BDA00027460105600001111
当
Figure BDA00027460105600001112
表明当前滤波器的一致性非常号，因此滤波器将尝试降低滤波阶γ，否则将保持滤波阶不变。
步骤五：非线性状态更新
通过观测设备得到tk+1时刻新的观测值
Figure BDA00027460105600001113
时，将其融合到预测值中，可以得 到卡尔曼滤波器的增益，以及状态的估计值
Figure BDA00027460105600001114
和其协方差矩阵
Figure BDA00027460105600001115
Figure BDA00027460105600001116
Figure BDA00027460105600001117
Figure BDA00027460105600001118
步骤五：反复执行上述算法，得到任意时刻航天器的轨道状态。
对本发明方法进行验证：
地球同步轨道卫星的初始状态为
Figure BDA00027460105600001119
初始位置误差为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 FDA0002746010550000011
式(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 FDA0002746010550000021
初始偏差为δxk，此时多项式形式的状态变量表示为
Figure FDA0002746010550000022
将该多项式变量在区间[tk，tk+1]上进行积分，得到微分方程在tk+1时刻的解为
Figure FDA0002746010550000023
对该解在标称状态
Figure FDA0002746010550000024
处做泰勒展开，得到高阶展开式
Figure FDA0002746010550000025
其中，
Figure FDA0002746010550000026
表示状态邻域[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 FDA0002746010550000027
处做泰勒展开，得到高阶展开式，
Figure FDA0002746010550000028
式(4)中，p表示m维测量矢量的索引值，
Figure FDA0002746010550000029
表示tk+1时刻的标称观测矢量，系数
Figure FDA00027460105500000210
表示泰勒展开系数，
Figure FDA00027460105500000211
表示tk时刻第i个状态偏差分量，其阶数为γi；γ＝[γ1，L，γn]T，表示各状态偏差分量的阶数。
4.根据权利要求3所述的非线性卫星轨道确定方法，其特征在于，高阶展开式
Figure FDA00027460105500000212
的γ阶多项式的近似解表示为：
Figure FDA00027460105500000213
式(5)中，i表示n维状态矢量的索引值；
Figure FDA0002746010550000031
表示标称状态；δxk，＝[δxk，1，L，δxk，n]T，表示初始偏差；
Figure FDA0002746010550000032
表示相应的泰勒展开式系数，γ＝[γ1，L，γn]T表示各状态偏差分量的指数。
5.根据权利要求4所述的非线性卫星轨道确定方法，其特征在于，步骤二所述的非线性状态预测过程中，轨道运动的预测状态均值
Figure RE-FDA0002826276080000033
和协方差矩阵
Figure RE-FDA0002826276080000034
为：
Figure RE-FDA0002826276080000035
Figure RE-FDA0002826276080000036
其中，E[]表示期望值，γ＝[γ1,…,γn]T和l＝[l1,…,ln]T表示变量偏差分量的指数，
Figure RE-FDA0002826276080000037
表示过程噪声的协方差矩阵；式(7)中，
Figure RE-FDA0002826276080000038
和
Figure RE-FDA0002826276080000039
其余的
Figure RE-FDA00028262760800000310
和
Figure RE-FDA00028262760800000311
与对应的
Figure RE-FDA00028262760800000312
均相等。
6.根据权利要求5所述的非线性卫星轨道确定方法，其特征在于，预测观测值的均值
Figure FDA00027460105500000313
和协方差矩阵Pk+1如下：
Figure FDA00027460105500000314
Figure FDA00027460105500000315
Figure FDA00027460105500000316
其中p和q表示测量矢量的索引值，i表示状态的索引值，
Figure FDA00027460105500000317
表示观测噪声的协方差矩阵；
Figure FDA00027460105500000318
表示状态矢量和测量矢量的协态矩阵，
Figure FDA00027460105500000319
表示测量值协方差矩阵；E[]表示期望值；式(9)中，
Figure FDA00027460105500000320
其余的
Figure FDA00027460105500000321
均与对应的
Figure FDA00027460105500000322
相等。
7.根据权利要求6所述的非线性卫星轨道确定方法，其特征在于，步骤三所述的滤波阶切换具体为：
定义测量新息矢量vk+1和相应的协方差矩阵Sk+1为
Figure FDA00027460105500000413
Figure FDA0002746010550000041
Figure FDA0002746010550000042
定义tk+1时刻的卡方指标归一化值为
Figure FDA0002746010550000043
定义tk+1时刻的卡方指标的时间平均值
Figure FDA0002746010550000044
为
Figure FDA0002746010550000045
其中，uk+1表示tk+1时刻的观测噪声；Rk+1表示tk+1时刻观测噪声的协方差矩阵，E[]表示期望值，L表示数据采集的窗口宽度，L等于3。
8.根据权利要求7所述的非线性卫星轨道确定方法，其特征在于，步骤三所述的自适应滤波切换的具体过程包括两种状态，分别为：
1)滤波器需要保持滤波阶γ，此时
Figure FDA0002746010550000046
滤波器不能正常工作，滤波器自动增加滤波阶，其中
Figure FDA0002746010550000047
可根据m和α值查卡方值表得到；
2)滤波器在滤波阶为γ时正常工作，此时
Figure FDA0002746010550000048
滤波器正常工作，之后进行降低滤波阶存在的可能验证过程。
9.根据权利要求8所述的非线性卫星轨道确定方法，其特征在于，降低滤波阶存在的可能性验证过程，包括：当
Figure FDA0002746010550000049
滤波器将尝试降低滤波阶γ；否则将保持滤波阶不变；其中，
Figure FDA00027460105500000410
和
Figure FDA00027460105500000411
通过卡方分布的概率表得到。
10.根据权利要求6所述的非线性卫星轨道确定方法，其特征在于，步骤四所述的非线性状态更新的具体过程为：
通过观测设备得到tk+1时刻新的观测值
Figure FDA00027460105500000412
时，将其融合到预测值中，得到卡尔曼滤波器的增益Kk+1、状态的估计值
Figure FDA0002746010550000051
和协方差矩阵
Figure FDA0002746010550000052
Figure FDA0002746010550000053
Figure FDA0002746010550000054
Figure FDA0002746010550000055
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 true CN112319859A (zh) 2021-02-05
CN112319859B 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)

Cited By (1)

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

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102305630A (zh) * 2011-05-17 2012-01-04 哈尔滨工业大学 基于扩展卡尔曼滤波的sar卫星自主定轨方法
CN103268407A (zh) * 2013-05-10 2013-08-28 航天东方红卫星有限公司 一种基于拉格朗日插值及卡尔曼滤波的轨道数据插值方法
CN103500455A (zh) * 2013-10-15 2014-01-08 北京航空航天大学 一种基于无偏有限冲击响应滤波器（ufir）的改进机动目标跟踪方法
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
US20190251215A1 (en) * 2018-02-15 2019-08-15 Regents Of The University Of Minnesota Accurate estimation of upper atmospheric density using satellite observations

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102305630A (zh) * 2011-05-17 2012-01-04 哈尔滨工业大学 基于扩展卡尔曼滤波的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
CN103268407A (zh) * 2013-05-10 2013-08-28 航天东方红卫星有限公司 一种基于拉格朗日插值及卡尔曼滤波的轨道数据插值方法
CN103500455A (zh) * 2013-10-15 2014-01-08 北京航空航天大学 一种基于无偏有限冲击响应滤波器（ufir）的改进机动目标跟踪方法
US20190251215A1 (en) * 2018-02-15 2019-08-15 Regents Of The University Of Minnesota Accurate estimation of upper atmospheric density using satellite observations

Cited By (2)

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

Also Published As

Publication number Publication date
CN112319859B (zh) 2021-12-31

Similar Documents

Publication Publication Date Title
US20090132164A1 (en) System and method for intelligent turning of kalman filters for ins/gps navigation applications
CN110553653B (zh) 基于多源数据驱动的航天器轨道确定方法
CN113670337A (zh) 一种用于gnss/ins组合导航卫星缓变故障检测方法
CN112319859B (zh) 一种基于自主滤波阶切换的非线性卫星轨道确定方法
CN110395297B (zh) 列车定位方法
CN110006427A (zh) 一种低动态高振动环境下的bds/ins紧组合导航方法
CN110275193A (zh) 一种基于因子图的集群卫星协同导航方法
CN110595434B (zh) 基于mems传感器的四元数融合姿态估计方法
CN111561930A (zh) 一种车载mems陀螺仪随机漂移误差的抑制方法
EP3598070A1 (en) Methods for monitoring the output performance of state estimators in navigation systems
CN111156986A (zh) 一种基于抗差自适应ukf的光谱红移自主组合导航方法
Wang et al. An EMD-MRLS de-noising method for fiber optic gyro Signal
CN106249590A (zh) 一体化自适应纳卫星姿态确定的方法
CN115356754A (zh) 一种基于gnss和低轨卫星的组合导航定位方法
CN113587926A (zh) 一种航天器空间自主交会对接相对导航方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
Reihs et al. A method for perturbed initial orbit determination and correlation of radar measurements
Qian et al. An INS/DVL integrated navigation filtering method against complex underwater environment
KR20200074660A (ko) 위성 탑재용 소프트웨어에 탑재 가능한 위성 이벤트 예측 기법
CN115543637B (zh) 空间目标的关联方法、装置和存储介质
Noureldin et al. Adaptive neuro-fuzzy module for inertial navigation system/global positioning system integration utilising position and velocity updates with real-time cross-validation
CN112415559B (zh) 一种基于多项式展开技术的高阶容错卫星轨道确定方法
Valizadeh et al. Improvement of navigation accuracy using tightly coupled kalman filter
CN112987054B (zh) 标定sins/dvl组合导航系统误差的方法和装置
CN108050997A (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