1、航空学报Acta Aeronautica et Astronautica SinicaFeb.15 2023 Vol.44 No.3ISSN 1000-6893 CN 11-1929/V526120-1考虑有色噪声影响的脉冲星导航 2级强跟踪差分滤波器许强*，崔洪亮，丁邦平，赵爱罡，赵阳青州高新技术研究所 测试控制系，潍坊 262500摘 要：为提高 X 射线脉冲星导航对有色噪声及太阳系内星历误差的鲁棒性，设计了 2 级强跟踪差分滤波器（TSTDKF）。首先在分析导航原理基础上，导出了中心天体星历误差对导航结果的误差传递关系，并利用扩展卡尔曼滤波器（EKF）进行了仿真验证；在同一运行轨道上，又
2、结合引力摄动模型对第三天体引力摄动数据进行了分析，证明该部分噪声为有色噪声。根据以上结论，将普通 2级卡尔曼滤波器（TKF）的无偏滤波器设计为一种改进的差分卡尔曼滤波器以降低有色噪声对导航系统的影响，同时又在其独立偏差滤波器中根据观测残差构建了多重自适应调节因子以增强其跟踪性能，两者共同构成 TSTDKF 的 2 个并行滤波器。通过仿真实验证明，TSTDKF 的位置误差性能最大可比EKF 和 TKF 改进 56.49%和 35.18%，速度误差性能改进 27.66%和 17.07%；对星历误差的跟踪效果也整体好于TKF。关键词：星历误差；有色噪声；2级卡尔曼滤波器；差分卡尔曼滤波器；强跟踪卡尔
3、曼滤波器中图分类号：V324.2+4 文献标识码：A 文章编号：1000-6893（2023）03-526120-13X 射线脉冲星是一种能够周期性发射稳定 X射线脉冲信号的中子星1。利用这一脉冲信号实现定位、授时功能的导航方法称为 X 射线脉冲星导航。它相比于传统的惯性导航、无线电导航和星光导航等导航方式具有累积误差较小、潜在应用情景广泛、信号不易受地面干扰源和天体运动影响、兼具授时功能等优点。目前已经有大量的国内外学者研究了其在近地卫星导航任务和地月轨道转移、火星探测等行星际探测任务中的应用问题2-4。根据相关原理，X 射线脉冲星导航需要使用到包括中心天体在内的众多太阳系天体的星历数据5。
4、而以近地卫星为例，目前已发表的相关文献证明 1 km 的地球轨道误差就会对导航结果带来将近千米级的定位误差6-7。所以 X 射线脉冲星导航技术的应用对太阳系内天体的星历数据精度提出了较高要求。此外，目前绝大多数的研究文献均将导航滤波器模型中的过程噪声与观测噪声认为是高斯白噪声8-11。过程噪声主要为除中心天体外其他天体对航天器的摄动影响，而观测噪声为探测器对 X 射线光子的探测噪声。实际上，这 2 部分噪声在天体运动等各方面因素的影响下不仅很难保证为白噪声，而且两者间还可能存在相关的情况。如果采用基于白噪声假设的数据模型进行建模解算，那么同样会影响计算精度。http:/ 引用格式：许强，崔洪亮
5、，丁邦平，等.考虑有色噪声影响的脉冲星导航 2 级强跟踪差分滤波器 J.航空学报，2023，44（3）：526120.XU Q，CUI H L，DING B P，et al.Two-stage strong tracking differential Kalman filter for X-ray pulsar navigation with coloured noise J.Acta Aeronautica et Astronautica Sinica，2023，44（3）：526120（in Chinese）.doi：10.7527/S1000-6893.2021.26120收稿日期：20
6、21-07-15；退修日期：2021-07-29；录用日期：2021-08-05；网络出版时间：2021-08-26 10：37网络出版地址：https：/ 长 基 线 干 涉 测 量 技 术（Very Long Baseline Interferometry，VLBI）技术和不同空间任务返回的观测数据，测定精度着实有限。以目前应用较为广泛的由美国喷气推进实验室发布的 DE（Development Ephemerides）星历为例，地球的位置精度只能精确到千米以内，木星更是仅能到几十千米12-13。因此，如何降低太阳系内星历误差对X射线脉冲星导航的影响已经成为一项亟待解决的实际问题。有学者研究
7、了增广估计算法以排除星历误差的干扰7。但是该方法会导致系统的状态量维数增加，带来较重的运算负担和数值病态问题。而且，增广估计算法将太阳系内星历误差认为是不变的。考虑到实际星历误差会存在缓慢的变化，增广估计算法长时间运行时的跟踪性能无法保证。其他学者还研究了不同形式的差分导航算法14-16。但是这些导航算法大多没有考虑差分计算过程中会导致过程噪声与观测噪声相关的问题。虽然文献 14 考虑了该问题，但它与其他研究文献一样，对原始模型中的过程噪声和观测噪声都非常理想地认为是高斯白噪声。而文献 16采用的方法是多航天器间观测数据的差分计算。该方法虽然具有一定效果，但并不适合单个航天器的导航系统。本文针
8、对以上问题提出了适合绕某一中心天体轨道运行的 2 级强跟踪差分卡尔曼滤波器（Two-stage Strong Tracking Differential Kalman Filter，TSTDKF）。首先将变化较为缓慢的中心天体星历误差作为独立偏差滤波器的状态量予以实时估计。考虑到星历的周年变化，构建了多重自适应调节因子以保证算法的跟踪精度。为解决原始模型中可能存在的有色噪声问题，将第 1级滤波器设计为一种改进的差分卡尔曼滤波器，提高其对有色噪声的鲁棒性能，拓宽 X 射线脉冲星的应用范围。1X射线脉冲星导航原理以绕地球运行的航天器为例，X 射线脉冲星导航的原理可描述为如图 1 所示。其中，c为光
9、速；rsat为卫星在国际天球参考坐标系（International Celestial Reference System，ICRS）17中的位置矢量；rE和E分别为地球在 ICRS 坐标系中的位置和速度矢量；r和分别为航天器在地球质心坐标系（Earth Centered Inertial，ECI）中的位置矢量和速度矢量，且满足rsat=rE+r；tSSB和tSC分别为太阳系质心 SSB（Solar System Barycenter，SSB）和航天器处的脉冲信号到达时间TOA（Time-of-Arrival，TOA）；n为 脉 冲 星 在ICRS坐标系中的单位方向矢量。安装在航天器上的 X 射
10、线探测器配合相应的信号恢复算法，可以获得tSC。而tSSB可以根据对应脉冲星的相位时间模型进行预测18。tSC与tSSB之间的差值t便是脉冲信号在航天器与 SSB之间沿脉冲星方向传播的时间，也就等价于两点间距离在该脉冲星方向的投影。如果选用 3颗不同方向的脉冲星，则通过坐标变换即可得到航天器在 ICRS 坐标系中的绝对位置。配合卡尔曼滤波等迭代算法即可获得稳定持续的导航信息。根据几何关系，tSSB和tSC之间的时间差t与rsat满足：tSSB-tSC=t=nrsatc（1）考虑到脉冲信号在传播过程中会受到引力延迟效应和周年视差效应等因素的影响，其高阶模型实际上可近似表示为5 地球XSSBZSS
11、BYSSBrErrsatn脉冲数航天器X射线脉冲星cttSSBtSC123国际天球参考坐标系SSBE图 1X射线脉冲星导航原理Fig.1X-ray pulsar navigation principle航空学报526120-3tSSB-tSC=t=nrsatc+2GMsunc3ln|1+nrsat+rsatnb+b|+12cD()nrsat2-r2sat-2brsat+2()nb()nrsat （2）式 中：D为 脉 冲 星 到 SSB 的 距 离；b是 SSB 在ICRS 坐标系中相对于太阳质心的位置矢量；b和rsat分别为b和rsat的长度模量；G是万有引力常量；Msun是太阳质量。式（2
12、）便是目前 TOA 计算应用最为广泛的一种高阶模型，也是脉冲星导航卡尔曼滤波器中的观测模型。而状态模型则可用航天器绕中心天体飞行所满足的动力学方程。对近地航天器而言，其在 ECI坐标系中满足的动力学方程为f(rv)=r?v?=vxvyvz-rxr3 1-J2()REarthr2()7.5rz2r2-1.5+Fx-ryr3 1-J2()REarthr2()7.5rz2r2-1.5+Fy-rzr3 1-J2()REarthr2()7.5rz2r2-4.5+Fz（3）式 中：r=rx，ry，rzT；v=vx，vy，vzT；=GMEarth是地球引力常数；MEarth为地球质量；J2为地球的二阶带谐项
13、；REarth为地球半径。Fx、Fy、Fz为其他天体对航天器的引力摄动，是过程噪声的主要组成部分，f()为式（3）中的函数关系。若将系统的状态量x选为航天器在 ECI坐标系中的位置r与速度v，观测量确定为tSSB和tSC之间的时间差t，则系统的状态空间方程为x?=f(x)+w（4）z=h(x)+V（5）式中：w为系统的过程噪声，除其他天体的引力扰动外，还包括太阳光压，地球潮汐等摄动影响；x=rkvkT为状态量；z为观测量，即不同脉冲星观测的时间差t；h()为式（2）等号右侧的函数关系式；V为观测噪声。为满足线性滤波器的应用，常将式（4）利用雅克比矩阵进行离散并线性化，将式（5）采取近似省略和泰
14、勒展开方式进行线性化。离散线性化后的状态空间方程可表示为19-20Xk=k|k-1Xk-1+GkWk（6）Zk=t1-w1()kt2-w2()kt3-w3()k=HkXk+HkEk+Vk（7）wi()k=2GMsunc3ln|1+nirsatk-1+rsatk-1nib+b|-2GMsunc3nirsatk-1nib+b+rsatk-1+nirsatk-1 i=1，2，3（8）Hk=H1()k01 3H2()k01 3H3()k01 3（9）HTi()k=nic+2GMsunc3ninib+b+rsatk-1+nirsatk-1+12cDi(nirsatk-1)ni-rsatk-1-2b+2(
15、nib)ni i=1，2，3（10）式中：k|k-1、Gk、Xk、Wk、Zk、Vk分别为离散化后的状态转移矩阵、噪声驱动矩阵、状态量、过程噪声、观 测 量 及 观 测 噪 声，且Xk=rkvkT=rxk，ryk，rzk，vxk，vyk，vzkT；ti为脉冲星i在航天器与 SSB 之间的时间差；rsatk-1为离散化后ICRS 坐标系中航天器在 k-1 时刻的位置矢量，rsatk-1为其标量形式；Ek=rEkvEkT为离散化后ICRS坐标系中地球的位置及速度矢量；rk、vk为离散化后 ECI坐标系中航天器的位置及速度矢量，且满足rsatk=rEk+rk；ni与Di分别为不同脉冲星的单位方向矢量和
16、到SSB的距离；Hk为观测矩阵。2误差分析2.1星历误差在 国 际 地 球 自 转 服 务 组 织（International Earth Rotation and Reference Systems Service，航空学报526120-4IERS）2003、2010 规 范 中，分 别 将 DE405 和DE421历表作为 ICRS坐标系中的动力实现21-22。受不同时期观测技术的限制，DE系列星历的精度也是在逐渐提升的。2008年发布的DE421相对于1995年的 DE405而言，其定位精度就得到了极大地改善23。虽然无法获得天体真实的位置数据，但是通过不同星历数据间的比较作差也能较为直观地展示出精度上的提升。以地球为例，分别利用DE405和 DE421计算其从 2015-10-17 04：00：00到2016-10-16 04：00：00在ICRS坐标系中的位置变化情况，采样点间隔设置为5 h。2种星历计算结果间的差值如图 2所示。通 过 图 2，基 本 上 可 认 为 DE421 相 对 于DE405 轨道误差降低了数十千米。但是 DE421仍然是不完美的。其中金星、地球和火