基于记忆衰减滤波的X射线脉冲星自主导航

doi:10.6052/1000-0879-19-086

基于记忆衰减滤波的X射线脉冲星自主导航

胡腾戈, 武迪^,¹⁾

清华大学航天航空学院,北京 100084

AUTONOMOUS NAVIGATION OF PROBE BY X-RAY PULSARS BASED ON FADING MEMORY FILTERING

HU Tengge, WU Di^,¹⁾

School of Aerospace Engineering, Tsinghua University, Beijing 100084, China

通讯作者: 1) E-mail:wud17@mails.tsinghua.edu.cn

收稿日期: 2019-03-5 网络出版日期: 2019-08-27

Received: 2019-03-5 Online: 2019-08-27

作者简介 About authors

摘要

针对扩展Kalman滤波(Extend Kalman Filter, EKF)在处理X射线脉冲星的自主导航定轨问题上易发散的问题,本文提出了基于记忆衰减滤波的抑制发散自主导航算法。在研究了自主导航滤波中误差变化的特征的基础上,从原理上分析了预报误差随迭代次数增加而增大的原因。针对EKF存在的这一缺陷,本文应用记忆衰减滤波控制误差发散,推迟了误差发散的时间并减小了发散的幅值,使得滤波结果稳定于更精确的结果。本文的结论均通过数值计算的方法验证

关键词： 自主导航 ; X射线脉冲星 ; 记忆衰减滤波

Abstract

In view of the fact that the Extend Kalman Filter (EKF) is prone to be divergent in the autonomous navigation by X-ray pulsars, this paper proposes an algorithm utilizing the fading memory filter. Based on the study of the error pattern of the EKF with the increase of iterations, some possible reasons for the divergency are analyzed. The fading memory filter can be used to reduce the peak of error, while delaying the divergence

Keywords： autonomous navigation ; X-ray pulsar ; fading memory filtering

PDF (1386KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

胡腾戈, 武迪. 基于记忆衰减滤波的X射线脉冲星自主导航. 力学与实践[J], 2019, 41(4): 382-384 DOI:10.6052/1000-0879-19-086

HU Tengge, WU Di. AUTONOMOUS NAVIGATION OF PROBE BY X-RAY PULSARS BASED ON FADING MEMORY FILTERING. MECHANICS IN ENGINEERING[J], 2019, 41(4): 382-384 DOI:10.6052/1000-0879-19-086

近年来,随着中国对于火星、小行星等的探测项目不断开展,深空探测逐渐受到更多的关注。在深空探测任务中,由于探测器与地球间遥远的距离,依赖于地面测控的导航方法会出现较大的通信时延,无法根据探测器所处的具体情况即时地给出相应的控制^[1];同时深空探测器巡航阶段持续时间长,常规的导航方式将占用大量的基站资源^[2]。因此,深空探测项目的独特性,使得其发展不需要地面设备支持,仅利用星载设备的自主导航系统成为了必然选择。

脉冲星是一种中子星,其磁极能够发射出一定波长的电磁波束。因星体高速自转的缘故,地球处观测到的电磁波信号呈现出脉冲的特性,因而得名。由于其自转周期具有较好的稳定性,脉冲星作为天然时钟参与导航的想法自20世纪70年代起便受到了关注^[1]。Maryland大学的Sheikh等^[3]提出了考虑相对论效应的改进到达时间计算公式,目前得到了广泛的认可与使用。以此公式为基础观测方程,基于X射线脉冲星的研究正受到广泛的关注。毛悦等^[4]提出脉冲星导航与传统导航方法的联合应用,王奕迪等^[5]研究了单脉冲星导航方法的可行性。同时,基于X射线脉冲星的误差分析问题也同样值得关注^[6-7]。

在X射线脉冲星自主导航中,扩展Kalman滤波(Extend Kalman Filter, EKF)是一种较常见的位置估计算法^[2,8]。但是,EKF滤波在处理X射线脉冲星自主导航时容易发散,这使得人们开始研究更适合的滤波算法,其中记忆衰减滤波是一种比较简便的抑制发散方法^[9]。但EKF滤波发散的原理与特征缺乏有效的分析,记忆衰减滤波抑制发散的原理与效果也较少被研究。

因此,本文将主要探讨基于EKF的X射线脉冲星自主导航的滤波发散特征,并研究记忆衰减滤波在抑制EKF发散时发挥的作用。研究主要分为三个方面：首先,建立探测器绕火星的精密动力学模型,应用Sheikh的计算模型建立脉冲星的测量方程,并基于EKF算法实现X射线脉冲星的自主导航。其次,针对EKF的误差发散的特征,从算法原理上分析出现问题的原因,并在此基础上,应用记忆衰减滤波避免发散,从而提高导航效果。最后,选取合适的参数给出数值仿真。

1 导航模型

1.1 探测器轨道方程

探测器在火心赤道惯性系中的运动方程写为

(1) $\begin{equation} \label{eq1} \ddot{r}={a}_{0} +{a}_{\rm{ns}}+{a}_{\rm{nbody}}+{a}_{\rm{drag}} \end{equation}$

式中,${r}$为探测器在火心赤道惯性系中的位置矢量。${a}_{0}$为火星对探测器的球形引力加速度,为平方反比形式

(2)$\label{eq2} {a}_{0}=-\frac{\mu_{\rm{m}}}{r^{3}}{r} $

式中$\mu_{\rm m}=GM$,其中$G$为引力常数,$M$为火星质量。${a}_{\rm{ns}}$为火星非球形引力项。由于火星并非严格的球体,可以假设真实引力为球体质心的万有引力(即${a}_{0}$项)与非球形引力项的叠加。通常情况下近球形天体的引力场能用球谐系数展开表示,根据火星引力场的GMM-1模型^[10],产生非球形引力的引力势可表示为

(3) $ \label{eq3} U=\frac{\mu_{\rm m}}{r} \bigg[1-\sum_{i=1}^n \Big( \frac{R}{r}\Big)^{i} J_{i} P_{i}(\sin\phi)\bigg]$

式中,$P_{i}$为各阶勒让德多项式,$J_{i}$为带谐项各阶系数,$R$为火星半径。本文取的阶数为4阶,应用此模型,即可得到相应的带谐项系数;将引力势求偏导,即可得到火星的非球形引力。${a}_{\rm{nbody}}$ 为其他天体引力项,其计算公式为

(4) $ \label{eq4} {a}_{\rm{nbody}} =- \frac{\mu_{\rm s}}{r_{\rm s}^{3}}{r}_{\rm{s}} + \frac{\mu_{\rm m}}{r_{\rm m}^{3}}{r}_{\rm{m}} $

在火心坐标系中,探测器受到太阳,火星的两颗天然卫星火卫一与火卫二,以及太阳系中其他大行星的作用力,与坐标系自身的加速度抵消,起到了类似潮汐力的作用。根据文献^[11]的结果,其他大行星的引力摄动量级远小于其余项,故本文仅考虑太阳与火星卫星的引力摄动,而其他大行星的引力摄动忽略不计。式(4)中下标s代表太阳(或火卫),下标m 代表火星。${a}_{\rm{drag}}$为大气阻力项。使用文献^[12]中引用的Joel Benito 2008火星大气模型,有

(5)$ \label{eq5} \left.\begin{array}{*{20}l} & \rho =\rho_{0}\exp \Big(\dfrac{h-h_0}{h_0}\Big)\\[3mm] & {a}_{\rm{drag}}=-\dfrac{1}{2}C_{\rm{d}}\dfrac{A}{m}\rho v\cdot {v} \end{array}\right\}$

式(5)中$\rho_{0}$为参考大气密度,$h_{0}$为参考大气高度,${v}$为探测器相对于火星地面的速度,即探测器在火心赤道系中的速度与当地自转速度之差。$C_{\rm{d}}$ 为探测器的阻力系数,本文中取为${1.0}$。$A$为探测器速度方向截面积,本文中取${A}/{m} =0.02$。

1.2 观测方程

自主导航算法的核心思想是建立探测器位置与某一观测量之间的函数关系,即观测方程。在X射线脉冲星导航中,选取的观测量是同一脉冲信号到达探测器与某一参考点之间的时间差(time of arrival, TOA)。参考点一般选取太阳系质心SSB,每一脉冲到达此点的时间可以预先计算出存储于探测器中。在观测时,通过信号处理方法分析X射线探测器的数据,得到脉冲信号到达探测器的时间。这样即可测量到同一信号的TOA。Sheikh等^[3]综合考虑了电磁信号在空间中的传播时间、引力对信号传播的相对论效应等影响因素,给出了量度TOA与探测器坐标的精确方程,其表达式为

(6)$ t_{\rm{SSB}}-t_{\rm{SC}} = \frac{1}{c}\left| {D}-{b} \right|-\frac{1}{c}\left| {D}-{p} \right| + \notag\\{} \sum_1^N \frac{2\mu_{k}}{c^{3}} \ln \Big| {n}_{\rm{SC}}\cdot {p}_{k} +\left\| {p}_{k} \right\| \Big|- \notag \\ \frac{2\mu_{\rm{s}}^{2}}{c^{5}D_{y}^{2}}\bigg\{ \left| {D}-{p} \right| \bigg[\Big( \frac{{n}_{\rm{SC}}\cdot {D}}{\left\| {D} \right\|} \Big)^{2} +1 \bigg] + \notag\\ 2\left( {n}_{\rm{SC}}\cdot {D} \right) \left( \frac{\left\| {p} \right\|}{\left\| {D} \right\|} - 1 \right) + \notag\\ D_{y} \Big[\arctan \Big(\frac{{n}_{\rm{SC}}\cdot {D}}{D_{y}} \Big) - \arctan \Big(\frac{{n}_{\rm{SC}}\cdot {p}}{D_{y}} \Big) \Big] \bigg\}$

式中${D}$为脉冲星相对太阳质心的位置矢量,${b}$为参考点相对太阳质心的位置矢量,${n}_{\rm{SC}}$为探测器相对脉冲星的单位位置矢量,${p}$为探测器到太阳质心的位置矢量,${p}_{k}$为探测器到第$k$颗行星的位置矢量。

考虑到式(5)对于计算的要求比较高,而选取太阳系质心SSB作为参考点的传统简化算法无法满足定位精度的要求,Ma等^[13]提出了将参考点选取在火星质心的计算方法,即

(7)$ t_{\rm{mars}}-t_{\rm{SC}}=\frac{1}{c}{n\cdot r}$

式中,${n}$为火星到脉冲星的单位位置矢量,${r}$为火星到探测器的距离矢量。火星处的参考时间可以利用式(5)预先算出内置于探测器中。利用此方法,是由于火星与探测器之间的距离较近,相对论效应项产生的影响几乎一致,因此在两时间相减$t_{\rm{mars}}-t_{\rm{SC}}$时,其影响相较距离项可以忽略不计,故仅需考虑二者的在脉冲星方向的距离差在TOA中的贡献。

这一简化算法的合理性可以通过数值计算的方法给出。仿真计算的结果如图1所示,其中观测值选取为TOA与光速的乘积$\Delta z= c(t_{\rm{mars}}-t_{\rm{SC}})$。

计算的参数与下文一致。由此可见,近似公式带来的位置误差量级远在1 m 以下,可以满足定位精度的需求。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 精确计算公式与近似计算公式之间误差

2 轨道估计方法

2.1 扩展Kalman滤波

Kalman滤波是一种根据时间序列的状态转移矩阵与观测方程,通过已有的估计值与实时更新的观测值,给出时间序列中的最优估计问题的算法。对于线性的、离散的系统,满足式

(8)$ {X}={\varPhi X} + {\varGamma W} \\ {Z}={HX} +{V} $

其中${X}$为系统的状态向量,${\varPhi}$为状态转移矩阵;${W}$为过程噪声,${\varGamma}$为噪声输入矩阵。${Z}$为观测向量,${V}$为观测噪声。各噪声满足的条件

(9)$ E[ {W}_{k} ]=0,~~ E[{W}_{k}{W}_{j}^{\rm{T}}]={Q}_{k}\delta_{kj} \\ E[ {V}_{k} ]=0,~~ E [{V}_{k}{V}_{j}^{\rm{T}}]={R}_{k}\delta_{kj} $

则对于此系统,有状态预测

(10)$ {P}_{k\vert k-1}= {\varPhi } {P}_{k-1}{\varPhi }^{\rm{T}} + {\varGamma Q}{\varGamma }^{\rm T} \\ {K}={P}_{k\vert k-1}{H}^{\rm T} [{H}{P}_{k\vert k-1}{H}^{\rm T} + {R}]^{-1} \\ \hat{X}_{k}=\hat{X}_{k\vert k-1}+{K}[{Z}_{k}-h(\hat{X}_{k\vert k-1})] \\ {P}_{k}=[{I}-{{KH}}]{P}_{k\vert k-1} $

在式中,$\hat{X}_{k}$为$k$时刻对向量${X}_k$给出的预测。${P}_{k}$ 为$k$时刻的误差方差矩阵估计,在滤波过程中不断被更新。

在实际应用中,探测器状态的转移方程与观测方程均不是线性、离散的,但可以通过线性化、离散化处理得到如式(8)的形式;这其中线性化处理即使用EKF进行。利用Kalman滤波器,可以利用脉冲星观测数据,迭代计算探测器状态,给出探测器状态的最优估计。

2.2 EKF性能分析

虽然EKF是自主导航中常用的轨道估计方法,但众多研究者的仿真结果表明^[14-16],EKF算法得出的轨道估计误差容易发散。EKF算法得出估计误差发散的原因主要有：

(1) 建立的轨道模型与真实不符。当自身的状态转移方程与观测方程不精确时,EKF可能出现发散。而若相关方程在滤波过程中由于摄动出现变化,EKF也可能出现发散^[15]。

(2) 算法自身的原因。在式(10)中,误差方差矩阵${P}_{k}$在迭代时,其值与新的观测量${Z}_{k}$没有直接关系,仅由前一状态下的误差方差矩阵${P}_{k-1}$ 修正得到^[16]。因此,当方程的相关初始值设定不恰当时,${P}$ 可能无法对于误差给出较好的估计。由于误差方差阵将影响式(10)中状态估计量$\hat{X}_{k}$ —$\hat{X}_{k}$与观测值${Z}_{k}$的权重,即误差方差较大时前一状态估计量的权重将降低。如果${P}$长时间低于真实的误差方差,估计将会过于依赖前一状态估计量,从而造成误差的积累,而出现滤波的发散。

2.3 记忆衰减滤波

本文引用的记忆衰减滤波是基于EKF做出的改进,其处理的系统仍为满足式(8)、式(9)的线性离散系统。记忆衰减算法主要针对的问题为上面的第2条,即避免误差方差估计${P}$长时间低于真实的误差方差。基于此,算法给出的状态预测为

(11)\begin{align}\left.\begin{array}{P}_{k\vert k-1}=c\cdot {\varPhi }{P}_{k-1}{\varPhi }^{\rm{T}} + {\varGamma Q}{\varGamma }^{\rm{T}} \\{K} = {P}_{k\vert k-1}{H}^{\rm{T}}[{H}{P}_{k\vert k-1}{H}^{\rm{T}}+{R}]^{-1}\\\hat{X}_{k}=\hat{X}_{k\vert k-1}+{K}[{Z}_{k}-h(\hat{X}_{k\vert k-1})]\\{P}_{k}=[{I}-{{KH}}]{P}_{k\vert k-1}\end{array} \right\}\end{align}

相较于式(10),记忆衰减滤波在误差方差矩阵${P}_{k\vert k-1}$估计时乘上了常数$c$。在实际应用中$c$一般为大于1的常数。从直观上,增大误差方差阵将使得已有系统信息的可信度降低,从而系统更少利用已有的系统信息,而更多基于新的观测值给出预测,从而避免误差积累而导致的发散。

3 仿真处理结果

3.1 参数设定

本文中选用的三颗脉冲星参数由表1给出。

表1 脉冲星参数表

新窗口打开| 下载CSV

表中kpc为千秒差距。探测器在火心赤道系中的初始坐标为

(12)$ \begin{equation} \label{eq12} \left. {\begin{array}{*{20}l}{ x} =(3.7\times 10^6 ,0,0)\rm~ (m)\\ { v}=(0,4000,0)\rm ~(m/s) \end{array} } \right\}\end{equation}$

在滤波器中，在动力学模型中三方向加速度的误差方差阵${ Q}$为

(13)$\begin{equation} \label{eq8} { Q} =diag(10^{-16}, 10^{-16}, 10^{-16}) \end{equation}$

观测方程三个方向距离的误差方差阵${ R}$为

(14)$\begin{equation} \label{eq9}{ R} = diag[\rm (300~m)^2, (300~m)^2,(300~m)^2] \end{equation}$

Kalman滤波初始预报的误差方差阵${ P}_{0}$设定为

(15)$\begin{equation} \label{eq10} { P}_{0} =diag\rm \big[10\,km,10\,km,10\,km,1\,m/s,1\,m/s,1\,m/s\big] \end{equation}$

3.2 EKF处理结果

在数值模拟中,计算了探测器一段轨道。一方面,通过数值求解动力学方程式(2),并加入随机作用项,得到真实情况下探测器的坐标。基于此真实坐标,计算出探测器得到的观测值,即真实情况下探测器处理X射线后得到的TOA。将TOA与初始位置、初始预报误差方差阵代入EKF中,计算得到估计的位置坐标。将得到的位置坐标与用受力方程迭代得到的真实位置坐标相比较,即可得到估计位置误差,用此衡量算法的优劣。在动力学方程的求解中,与EKF位置更新的时间间隔均取为15s。

(16)$\Delta X=\sqrt{\Delta x^2 + \Delta y^2 + \Delta z^2} $

算例的运算结果如图2与图3所示。在算例计算前10000步时(对应时间为1.7d)尽管滤波的误差较小,但呈现出发散的趋势。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 EKF算例误差: 前10000步

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 EKF算例误差: 最终形态

增加迭代步数至65000步(对应时间为11.4d),进一步计算结果表明尽管最终的误差不再增大,总体误差--时间曲线呈现出近似“S”型。但预报的误差量级在10km左右,仍然较大。Ma等^[13]利用EKF得到的轨道预测误差显示出了类似的“S”型曲线,而其误差稳定值同样为10km量级。

下面分析EKF得到的位置估计值发散的原因。如图4所示,EKF迭代中的方差估计值随着迭代次数增加而逐渐减小,无法给出一个有效的预测,而长时间低于真实的滤波误差平方值。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 误差平方值与估计方差对比, $c=$1

因此在式(10)中,计算最优估计$\hat{X}_{k}$时,计算结果过多依赖于已有的估计值$\hat{X}_{k-1}$而非观测值,最终造成误差的累加。滤波发散的其他原因,可能为计算机舍入误差等。

3.3 记忆衰减滤波处理结果

由式(11),记忆衰减滤波通过常数$c$能够控制方差阵${P}$的大小,减缓其逐渐趋于0的速度,如图5所示。下面讨论取不同常数$c$时的滤波效果。随着参数$c$逐渐增大,当取值较接近于1时,"S"型误差增长曲线到来上升段的时间被推后,如图6 所示。同时,最终稳定段的误差随着$c$的上升而逐渐减小,数据见表2.

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 误差平方值与估计方差对比, $c=1.001$

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 不同参数下滤波误差

表2 不同参数下稳定段误差

新窗口打开| 下载CSV

当$c$较大时,在程序所能计算的迭代次数内(即时间≥11.3d),误差不会经历如图6的“S”型增长过程,而停留在较低的水平,如图7所示。这时的误差数量级在100m左右。而在足够长的时间中误差是否仍会经历“S”型增长,仍有待研究。100m量级的导航定位精度,相比起上文图3中10km量级的误差,以及Ma等^[13]的研究结果,均有了较明显的改善。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 $c=1.0015$时滤波误差

当$c$进一步增大时,误差仍然停留在较低水平,但相较图7,滤波的误差又出现了增长,如图8所示。这应当是由于误差方差阵过大,导致估计时对已有的位置估计过于不信任,从而产生次优的估计结果。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8 $c=1.1$时滤波误差

4 结论

本文用仿真方法实现了X射线脉冲星自主导航算法,验证了其导航原理的合理性,并分析了EKF算法的误差发散特征,即为一条“S”型曲线。在此基础上,针对EKF滤波发散的原理,应用记忆衰减滤波抑制了导航误差的发散,并分析了不同参数下记忆衰减滤波的结果特征,即“S”型曲线的增长更为缓慢且稳定误差有所减小。在选取合适的参数$c$时,X 射线脉冲星自主导航算法得出的误差在100m量级,优于EKF的10km量级结果。此导航精度有较好的应用前景。同时,参数$c$的最优性选取和更长时间的滤波结果评估仍有待进一步的研究。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

郑伟, 王奕迪, 汤国建 . X射线脉冲星导航理论与应用.北京: 科学出版社, 2015

[本文引用: 2]

Zheng

Wei

, Wang

Yidi

, Tang

Guojian

X-ray Pulsar-Based Navigation:Theory and Applications.Beijing China Science Publishing & Media Ltd, 2015

[本文引用: 2]

[2]

李俊峰, 崔文, 宝音贺西 .

深空探测自主导航技术综述

力学与实践, 2012,34(2):1-9

[本文引用: 2]

Junfeng

, Cui

Wen

, Baoyin

Hexi

A survey on autonomous navigation for deep space exploration

Mechanics in Engineering, 2012,34(2):1-9

[本文引用: 2]

[3]

Sheikh

, Pines

, Ray

, et al.

Spacecraft navigation using X-ray pulsars

Journal of Guidance Control & Dynamics, 2006,29(1):49-63

[本文引用: 2]

[4]

毛悦, 宋小勇 .

X射线脉冲星导航几何法确定航天器位置

武汉大学学报:信息科学版, 2009,34(7):790-793