CN112319859B - 一种基于自主滤波阶切换的非线性卫星轨道确定方法 - Google Patents
一种基于自主滤波阶切换的非线性卫星轨道确定方法 Download PDFInfo
- 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
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时,非线性扩展卡尔曼滤波器将退化为常规的扩展卡尔曼滤波器,显然,低阶算法可以获得良好的计算效率却有可能遭遇精度损失,相反,当选择滤波阶较高时,非线性扩展卡尔曼滤波器可以获得较高的精度却面临着巨大的计算压力。因此,如何克服非线性扩展卡尔曼滤波器在不同滤波阶时遭遇的问题而同时获得低级算法和高阶算法的优势,具有非常重要的实际意义。
发明内容
本发明的目的在于克服上述现有技术中,非线性扩展卡尔曼滤波器不能同时获得低级算法和高阶算法的优势的缺点,提供一种基于自主滤波阶切换的非线性卫星轨道确定方法
为了达到上述目的,本发明采用以下技术方案予以实现:
一种基于自主滤波阶切换的非线性卫星轨道确定方法,包括如下步骤:
步骤一,建立方程:分别建立航天器轨道动力学方程和观测方程;
步骤二,非线性状态和测量预测:分别对步骤一的轨道动力学方程和观测方程进行泰勒展开逼近,得到轨道状态的高阶多项式映射和观测方程的高阶多项式映射,并根据上一时刻状态统计信息计算得到状态预测值;
步骤三,计算测量新息,并进行卡方检测,根据结果对步骤二的两个高阶多项式映射过程进行自主阶切换;当系统非线性状态较强时,自主阶切换算法采用高阶算法进行滤波;当系统非线性状态较弱或滤波过程进入稳定状态时,自主阶切换算法采用低阶算法进行滤波;
步骤四,状态更新:设置时间参数、初始状态和当前时刻测量值,采用步骤三选择的自主阶切换算法,对步骤二中得到的状态预测值进行更新,得到目标时刻航天器的轨道状态。
优选地,步骤一建立的轨道动力学方程为:
式(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)对轨道运动进行非线性状态预测,具体为:
2)对航天器的观测进行非线性预测,具体为:
时间为tk+1时的测量方程为:
zk+1=h(xk+1,tk+1)+uk+1 (3)
优选地,步骤三所述的滤波阶切换具体为:
定义测量新息矢量vk+1和相应的协方差矩阵Sk+1为
定义tk+1时刻的卡方指标归一化值为
其中,uk+1表示tk+1时刻的观测噪声;Rk+1表示tk+1时刻观测噪声的协方差矩阵,E[]表示期望值,L表示数据采集的窗口宽度,L等于3。
优选地,步骤三所述的滤波阶切换的具体过程包括两种状态,分别为:
优选地,步骤四所述的状态更新的具体过程为:
与现有技术相比,本发明具有以下有益效果:
本发明公开了一种基于自主滤波阶切换的非线性卫星轨道确定方法,利用高阶方法适合于非线性较强时的精确状态估计,及低阶方法适合于滤波器进入平稳阶段的高效率滤波,通过设计自适应滤波阶切换策略,在保证精度的前提下,在一个估计过程中根据需要动态改变滤波阶。在滤波过程中,通过计算测量新息并检验滤波一致性特性,动态地设计了一个自适应自主滤波阶切换策略,当系统非线性较强时,采用高阶算法保证精度;当滤波过程进入稳定状态时,采用低阶算法保证效率。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法是一种混合算法,同时获取了低阶算法和高阶算法的优点。确保非线性扩展卡尔曼滤波器在保证精度的前提下,加快滤波过程的执行速度。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法,在保证精度的前提下,可以降低滤波器在滤波平稳阶段的滤波阶,显著减少计算成本,提高计算效率。克服了传统的非线性扩展卡尔曼滤波器在一个估计过程中无法改变滤波阶的困境,不仅保持原有非线性扩展卡尔曼滤波采用高阶时的高精度特点,同时继承了低阶运行时的高效率。本发明提出的基于自主滤波阶切换的非线性卫星轨道确定方法普遍适用于卫星、空间碎片等的轨道确定,可加速高精度卫星轨道确定的执行速度。
进一步地,通过对滤波过程中测量新息的监测在保证滤波一致性的要求下,动态地调整滤波平稳阶段的滤波阶,大大降低了滤波器的计算成本,进一步优化了卫星轨道确定的执行速度。
附图说明
图1为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的位置误差示意图;
图2为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的速度误差示意图;
图3为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的位置分量标准差估计值示意图;
图4为标准一阶、二阶非线性扩展卡尔曼滤波器和本发明提出的基于滤波阶切换的非线性扩展卡尔曼滤波器得到的速度分量标准差估计值示意图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分的实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都应当属于本发明保护的范围。
需要说明的是,本发明的说明书和权利要求书及上述附图中的术语是用于区别类似的对象,而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换,以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外,术语“包括”和“具有”以及他们的任何变形,意图在于覆盖不排他的包含,例如,包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元,而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
下面结合附图对本发明做进一步详细描述:
实施例1:
本实例公开一种基于多项式展开技术的高阶容错卫星轨道确定方法,具体实现步骤如下:
步骤一:建立无量纲的卫星轨道状态方程为
其中分别表示无量纲的时间,卫星的无量纲位置矢量、速度矢量和角速度值,t,r,v,ω表示相应的有量纲值,表示卫星无量纲位置矢量的范数;另外,卫星的标称角速度表示卫星在半长轴为a=42164km的轨道上运行的平均的轨道角速度,μe表示地球的引力常数。
步骤二:建立非线性测量方程,包括径向距离,赤经和赤纬,
其中u=[u1,u2,u3]T表示测量噪声矢量,径向距离的测量噪声满足零均值标准差为1米的高斯分布,赤经和赤纬的测量噪声满足零均值标准差为1.745×10-6rad的高斯分布。同时,我们假设每个轨道可以获得12个等间距的测量值。
步骤三:非线性状态预测:
假设tk时刻航天器的标称状态为初始偏差为δxk,因此初值的一个邻域可表示为被称为多项式状态变量,在区间[tk,tk+1]上进行积分后得到微分方程在tk+1时刻的解为对该解在标称状态处做泰勒展开可得到高阶展开式其中表示状态t1状态邻域[xk+1]对初始邻域大小δxk的非线性依赖关系,其γ阶多项式近似解可表示为
类似的,根据测量方程(2),tk+1的测量方程可表示为
zk+1=h(xk+1,tk+1)+uk+1 (3)
其中p表示m维测量矢量的索引值,表示对应于标称观测值,系数表示泰勒展开系数,可通过将方程(5)代入方程(3)计算得到,m表示观测方程的数量;在tk+1时刻,使用γ阶多项式形式的解,即方程(4),预测观测值的均值如下:
步骤四:滤波阶切换策略
定义测量新息矩阵vk+1和相应的协方差矩阵Sk+1为
进一步,qk+1定义为
其中,L表示数据采集的窗口宽度,默认为3。
为了设计自适应的滤波阶切换策略,先提出两个假设
1)H0:滤波器在滤波阶为β时正常工作;
2)H1:滤波器需要保持滤波阶β。
在假设H0成立时,新息平方归一化平均值是一个具有m个自由度卡方分布其中重要性系数α=99%,因此我们可以通过卡方分布的概率表得到对应于参数m和α的临界值当时,假设H0不成立而假设H1成立,即滤波器不正常工作,滤波器自动增加滤波阶。当时,假设H1不成立而假设H0成立,即滤波器正常工作,在这种情况下,我们需要进一步验证是否存在降低滤波阶的可能性。类似的,我们提出以下两个假设:
为了验证和我们需要围绕定义一个严格的区间一般情况,αl=10%,αu=75%,因此我们可以通过卡方分布的概率表得到对应于参数αl和αu的临界值和当表明当前滤波器的一致性非常号,因此滤波器将尝试降低滤波阶γ,否则将保持滤波阶不变。
步骤五:非线性状态更新
步骤五:反复执行上述算法,得到任意时刻航天器的轨道状态。
对本发明方法进行验证:
地球同步轨道卫星的初始状态为
初始位置误差为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.一种基于自主滤波阶切换的非线性卫星轨道确定方法,其特征在于,包括如下步骤:
步骤一,建立方程:分别建立航天器轨道动力学方程和观测方程;
步骤二,非线性状态和测量预测:分别对步骤一的轨道动力学方程和观测方程进行泰勒展开逼近,得到轨道状态的高阶多项式映射和观测方程的高阶多项式映射,并根据上一时刻状态统计信息计算得到状态预测值;
步骤三,计算测量新息,并进行卡方检测,根据结果对步骤二的两个高阶多项式映射过程进行自主阶切换;当系统非线性状态较强时,自主阶切换算法采用高阶算法进行滤波;当系统非线性状态较弱或滤波过程进入稳定状态时,自主阶切换算法采用低阶算法进行滤波;
步骤四,状态更新:设置时间参数、初始状态和当前时刻测量值,采用步骤三选择的自主阶切换算法,对步骤二中得到的状态预测值进行更新,得到目标时刻航天器的轨道状态。
3.根据权利要求1所述的非线性卫星轨道确定方法,其特征在于,步骤二所述的非线性状态预测包括两个,分别为
1)对轨道运动进行非线性状态预测,具体为:
2)对航天器的观测进行非线性预测,具体为:
时间为tk+1时的测量方程为:
zk+1=h(xk+1,tk+1)+uk+1 (3)
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)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116331523B (zh) * | 2023-05-29 | 2023-08-25 | 哈尔滨工业大学 | 带大惯量旋转载荷卫星的未知参数辨识方法、装置及介质 |
Family Cites Families (5)
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 |
-
2020
- 2020-10-27 CN CN202011166646.9A patent/CN112319859B/zh active Active
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 | |
CN109032176B (zh) | 一种基于微分代数的地球同步轨道确定和参数确定方法 | |
CN112319859B (zh) | 一种基于自主滤波阶切换的非线性卫星轨道确定方法 | |
Leonard et al. | Absolute orbit determination and gravity field recovery for 433 eros using satellite-to-satellite tracking | |
CN110595434B (zh) | 基于mems传感器的四元数融合姿态估计方法 | |
Wang et al. | An EMD-MRLS de-noising method for fiber optic gyro Signal | |
Reihs et al. | A method for perturbed initial orbit determination and correlation of radar measurements | |
Garcia et al. | Unscented Kalman filter and smoothing applied to attitude estimation of artificial satellites | |
Lin et al. | Adaptive kernel auxiliary particle filter method for degradation state estimation | |
KR20200074660A (ko) | 위성 탑재용 소프트웨어에 탑재 가능한 위성 이벤트 예측 기법 | |
Qian et al. | An INS/DVL integrated navigation filtering method against complex underwater environment | |
CN111623764B (zh) | 微纳卫星姿态估计方法 | |
Qiu et al. | Modified multiplicative quaternion cubature Kalman filter for attitude estimation | |
CN110912535B (zh) | 一种新型无先导卡尔曼滤波方法 | |
Su et al. | Variational Bayesian adaptive high‐degree cubature Huber‐based filter for vision‐aided inertial navigation on asteroid missions | |
Wang et al. | Implementation of solution separation-based Kalman filter integrity monitoring against all-source faults for multi-sensor integrated navigation | |
CN108981689B (zh) | 基于dsp tms320c6748的uwb/ins组合导航系统 | |
Valizadeh et al. | Improvement of navigation accuracy using tightly coupled kalman filter | |
Li et al. | Fault detection approach applied to inertial navigation system/air data system integrated navigation system with time‐offset | |
CN112415559A (zh) | 一种基于多项式展开技术的高阶容错卫星轨道确定方法 | |
Sabzevari et al. | Observability analysis and design of two nested filters for the satellite attitude estimation with magnetometer‐only | |
Su et al. | Quasi-consistent fusion navigation algorithm for DSS | |
Reifler et al. | Multi-target ensemble Gaussian mixture tracking with sparse observations | |
Nanda et al. | Performance comparison of EKF and UKF for estimation of COE using Kozai mechanism | |
Martin et al. | The Physics-Informed Neural Network Gravity Model: Generation III |
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 |