一种环火探测器精密同步定位守时方法及系统

摘要

一种环火探测器精密同步定位守时方法及系统，首先输入初始数据，包括模拟的标称轨道和用于滤波初始化的相关参数；进行观测数据模拟，根据脉冲星测量方程模拟脉冲星观测量，该观测量涉及的测量噪声包括钟差，其中钟差由钟差模型模拟；根据环火探测器自主定位和守时的滤波状态方程和观测方程，进行自适应卡尔曼滤波，得到探测器的轨道和钟差。本发明在建立X射线脉冲星观测模型基础上，将观测模型、探测器动力学模型和星载原子钟的钟差模型进行有机结合，同时估算了探测器的位置、速度和钟差参数，在实现钟差校正的同时也削弱了钟差对定位精度的影响，提高了自主定位精度。

说明书

技术领域

本发明涉及火星探测技术领域，尤其涉及一种环火探测器精密同步定位守时方法及系统。

背景技术

对于执行火星探索任务的火星轨道探测器而言，实时、高精度地获取探测器的位置、速 度和姿态等相关信息，不仅是对探测器进行精密导航控制的基础，也是确保各种预定科学实 验和探测任务顺利实施的决定性因素之一。

火星和地球之间的遥远距离和探测器所面临的复杂太空环境对基于地面测控通讯站的传 统导航方法而言是一个巨大挑战，航天器的自主导航将是突破该瓶颈的关键技术之一。X射 线脉冲星导航是一种新型自主导航方法，在火星探测器自主导航方面极具应用潜力。

在利用X射线脉冲星观测量对环火探测器进行定位时涉及到环火探测器的动力学模型和 脉冲星观测模型，由环火探测器的动力学模型通过数值积分可以外推但仅能获得探测器的近 似轨道信息，此外，由于技术条件限制，探测器的初始轨道状态也存在误差。

X射线脉冲星导航系统中光子到达探测器的时间是由星载原子钟测量得到的，因此星载 原子钟的精度及其稳定性直接关系到导航精度的高低。而由于火星探测器飞行时间长，探测 器上搭载的原子钟，受到各种因素的限制，若不定期校正将难以长期保持较高精度，这样就 回对导航系统的精度造成损害。

目前已有的X射线脉冲星导航是一种新型自主导航方法，具备可靠性高、自主性强、精 确性良好和使用范围广等特点，在深空探测器自主导航领域应用潜力巨大。目前国际上许多 学者和专家都在这种新型自主导航方法进行研究和验证，国内外有大量的相关文献对其进行 了介绍，Suneel I S.2005.The Use of Variable Celestial X-ray Sources for Spacecraft Navigation [D]:[Ph.D.].Maryland:Department of Aerospace Engineering University of Maryland.Wei Erhu, Jin Shuanggen,Zhang Qi,et al.2013Autonomous navigation of Mars probe using X-ray pulsars: Modeling and results[J].Advances in Space Research,51(2013):849–857.Emadzadeh A A,Robert  A,Speyer J L,et al.2011.Relative Navigation between Two Spacecraft Using X-ray Pulsars[J]. IEEE Transactions on Control Systems Technology,19(5):1021-1035.Yang Tinggao.2008. Determination of X-ray pulsar pulse time of arrival at spacecraft[J].Chin.J.Space Sci., 28(4):330-334.Chester T J,Butman S A.1981.Navigation Using X-ray Pulsars[R].NASA Technical Reports N81-27129,21-25.

利用X射线脉冲星观测量进行探测器自主定位的原理如图1所示。在基准点(太阳系质 心SSB)坐标系中，空间飞行器Probe上的X射线探测设备接收脉冲信号，进而获得信号到 达时间观测量，然后将其与基准点处的时间相位模型的计算值较差，差值τ反映了同一脉冲 信号到达航天器和基准点处的时间延迟，τ与光速c的乘积为c·τ。该差分观测量可以表示 为航天器位置矢量r的函数，在已知导航脉冲星方向矢量n的前提下，结合航天器的轨道动 力学方程即可滤波估算出航天器在基准点坐标系中的位置向量。

通常情况下脉冲星时间相位模型通常建立在太阳系质心SSB处，在J2000.0太阳系质心 坐标系中：

某一观测历元，脉冲星i的辐射信号至航天器SC和火星质心M的直接距离观测量ρ_SCi、 ρ_Mi为：

ρ_SCi＝||r_SC-D_i||+d_Rel,SCi      (1-1)

ρ_Mi＝||r_M-D_i||+d_Rel,Mi      (1-2)

式中，r_M和r_SC分别表示火星质心和探测器的位置矢量；D_i表示观测脉冲星的位置矢量；d_Rel,Mi和d_Rel,SCi分别表示相应的脉冲信号延迟观测误差。

式(1-2)减去(1-1)得到单差观测量δρ_i＝ρ_Mi-ρ_SCi和“单差测量方程”为：

δρ_i＝(||r_M-D_i||-||r_SC-D_i||)+(d_Rel,Mi-d_Rel,SCi)      (1-3)

因为脉冲星距离参考点的距离D_i的量级远大于航天器参考点的距离，因此对式(1-3)进 行展开，并忽略无穷小量以上的高阶项可得

δρ_i＝n_i·(r_SC-r_M)+(d_Rel,Mi-d_Rel,SCi)      (1-4)

火星轨道探测器相对于火星质心的位置矢量r_SC/M＝r_SC-r_M,脉冲星的方向矢量为n_i，则包 含误差项的单差分观测方程为：

$\delta {\tilde{\rho }}_i=n_i·\delta r_{\mathrm{SC}/M}+{\mathrm{c\delta t}}_{\mathrm{SC}}+{\eta }_i---(1-5)$

式中，单差改正量$\delta {\tilde{\rho }}_i=c(t_{\mathrm{Mi}}-{\tilde{t}}_{\mathrm{SCi}})-n_i·{\tilde{r}}_{\mathrm{SC}/M}-d_{\mathrm{Rel},\mathrm{Mi}}+d_{\mathrm{Rel},\mathrm{SCi}},$和δr_SC/M分别表示探 测器相对于火星质心的位置向量及其改正数；t_Mi表示火星脉冲信号到达火星质心的时间；理 想时钟读数为t_SC；和δt_SC分别为星载原子钟读出的时间和钟差改正量；c为光速；η_i代表 随机测量噪声。

而在环火探测器自主定位中为了简化数据处理，通常用火星质心代替太阳系质心。

测量模型(1-5)中包含了原子钟的钟差参数，而研究分析表明0.5us的原子钟钟差会带 来约300m的位置误差和0.03m/s的速度误差，而在目前的技术条件下应用于卫星和航天器的 原子钟的理论稳定指标为应用于卫星和航天器的原子钟的理论稳定指标为5×10^-12/s，对于任 务周期较长的深空探测器来说，若不定期进行钟差校正，原子钟的钟差将累计至数十微秒甚 至更多，这将给探测器的导航精度造成严重损害。

脉冲星与火星的距离遥远，脉冲信号微弱，X射线光子的累积一般需要8到10钟的时间， 因此在新的观测数据未被采集期间，需要对探测器轨道进行精确预报。常用数值积分的方法 进行探测器的轨道外推，该方法须要考虑探测器受到的中心引力和作用于探测器上的各项摄 动力。

在J2000火星质心地球赤道坐标系S_ECI中，探测器的动力学模型可表示为：

$\stackrel{··}{r}=-\frac{{\mu }_{M}r}{r^3}+a_r+W---(1-6)$

式中，r＝[x y z]^T为探测器在S_ECI中的位置矢量，为其二阶导数，r＝||r||；μ_M是 火星引力常数；W为满足高斯分布的随机噪声；a_r则代表探测器受到的各种摄动力所产生 的摄动加速度总和。研究表明，对于中高轨道火星探测器而言火星非球形形状摄动中的J₂项 为最主要的影响。若仅考虑J₂项摄动，则有：

$a_r=a_{J_2}---(1-7)$

式中

$a_{J_2}=\left(\begin{array}{c}-\frac{{\mu }_{M}x}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(1-5\frac{z^2}{r^2})\\ -\frac{{\mu }_{M}y}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(1-5\frac{z^2}{r^2})\\ -\frac{{\mu }_{M}z}{r^3}{J}_2{\left(\frac{R_M}{r}\right)}^2\frac{3}{2}(3-5\frac{z^2}{r^2})\end{array}\right)---(1-8)$

上式中，R_M为火星平均半径。

若获得某一时刻的探测器轨道状态参数，便可利用建立的轨道动力学模型，通过数值积 分的方法外推下一个时刻的轨道值。

综上所述，利用环火探测器轨道动力学模型进行轨道外推时，外推轨道的精度直接受轨 道动力学模型的精确性和轨道积分方法的影响，此外初始轨道参数若存在误差，外推轨道将 很快发散。从式(1-8)中可以看到，若不定期进行钟差校正，原子钟的钟差将对定位精度造 成损害。

发明内容

针对以上问题，本发明在建立X射线脉冲星观测模型基础上，将观测模型、探测器动力 学模型和星载原子钟的钟差模型进行有机结合，同时估算了探测器的位置、速度和钟差参数， 在实现钟差校正的同时也削弱了钟差对定位精度的影响，提高了自主定位精度。

本发明的技术方案提供一种环火探测器精密同步定位守时方法，包括以下步骤，

步骤1，输入初始数据，包括模拟的标称轨道和用于滤波初始化的相关参数；

步骤2，进行观测数据模拟，根据脉冲星测量方程模拟脉冲星观测量，该观测量涉及的测量 噪声包括钟差，其中钟差由钟差模型模拟，

设导航脉冲星组合中包括p颗脉冲星，环火探测器的位置矢量r＝[x y z]^T、速度矢量 v＝[v_x v_y v_z]^T，加速度矢量a＝[a_x a_y a_z]^T，

所述脉冲星测量方程如下式，

$Z=\left(\begin{array}{c}\Delta t_1\\ ·\\ ·\\ ·\\ \Delta t_i\\ ·\\ ·\\ ·\\ \Delta t_p\end{array}\right)=H·X+\eta$

其中，Δt_i表示第i颗脉冲星时延观测量，i的取值为1,2,…,p，p颗脉冲星观测量组成观测向 量Z，η为测量噪声，H为观测方程系数矩阵，状态参数X＝[x y z v_x v_y v_z x₁ x₂ x₃]^T；

所述钟差模型如下式，

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+ϵ·x_2\left(t_k\right)+\frac{{ϵ}^2}{2}·x_3\left(t_k\right)+W\left(t_k\right)$

其中，ε为采样间隔，x₁(t_k)、x₂(t_k)和x₃(t_k)分别表示时刻t_k原子钟的钟差、频率漂移和频 率漂移变化率，组成了钟差参数[x₁ x₂ x₃]；W(t_k)为时刻t_k相应的系统噪声，x₁(t_k+1)表示 时刻t_k+1原子钟的钟差；

步骤3，根据环火探测器自主定位和守时的滤波状态方程和观测方程，进行自适应卡尔曼滤 波，

所述观测方程采用脉冲星测量方程，

设时刻t的状态参数X记为X(t)，X(t)的一阶导数为所述滤波状态方程如下，

$\stackrel{·}{X}\left(t\right)=F\left(t\right)·X\left(t\right)+W\left(t\right)$

其中，W(t)为系统状态噪声矩阵，

状态方程系数矩阵F(t)为，$F\left(t\right)=\left(\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y}& \frac{\partial a_x}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_z}{\partial x}& \frac{\partial a_z}{\partial y}& \frac{\partial a_z}{\partial z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$

步骤4，输出步骤3所得滤波结果，包括滤波所得探测器的轨道和钟差。

而且，步骤3进行自适应卡尔曼滤波的方式如下，

1)滤波的相关参数初始化，包括计算初始时刻的状态参数x_k及其方差矩阵P_k，状态误差 δx_k；

2)进行一步预测值x_k+1,k和修正值δx_k+1,k计算如下，

x_k+1,k＝Φ_k+1,k·x_k

δx_k+1,k＝0

其中，Φ_k+1,k为状态转移矩阵；

3)进行自适应因子λ_k+1的计算如下，

$M_{k+1}=H_{k+1}{\Phi }_{k+1,k}{P}_k{\Phi }_{k+1,k}^T{H}_{k+1}^T+R_{k+1}$

δz_k+1＝z_k+1-H_k+1x_k+1,k

$N_{k+1}=\delta z_{k+1}\delta z_{k+1}^T$

${\lambda }_{k+1}=\mathrm{Max}\{1,\sqrt{\mathrm{sp}\left(N_{k+1}\right)/\mathrm{sp}\left(M_{k+1}\right)}\}$

式中，H_k+1表示观测方程系数矩阵；符号M_k+1和N_k+1用于计算自适应因子；z_k+1为观测值向 量；预测残差为δz_k+1；R_k+1为观测值误差方差矩阵；sp(.)为求矩阵迹符号；

4)进行一步预测估计误差方差阵P_k+1,k的计算如下，

$P_{k+1,k}={\lambda }_k{\Phi }_{k+1,k}{P}_k{\Phi }_{k+1,k}^T+Q_k$

式中，Q_k为系统状态噪声方差矩阵；

5)进行滤波增益K_k+1的计算如下，

$K_{k+1}=P_{k+1,k}{H}_{k+1}^T{[H_{k+1}{P}_{k+1,k}{H}_{k+1}^T+R_{k+1}]}^{-1}$

6)进行测量更新如下，

δx_k+1＝δx_k+1,k+K_k+1[z_k+1-H_k+1x_k+1,k]

$P_{k+1}=[I-K_{k+1}{H}_{k+1}]P_{k+1,k}{[I-K_{k+1}{H}_{k+1}]}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T$

式中，P_k+1为滤波状态值x_k+1的方差矩阵，I表示相应阶数的单位矩阵；

7)进行滤波结果x_k+1输出如下，

x_k+1＝x_k+1,k+δx_k+1

而且，所述一步预测值x_k+1,k由轨道积分器按滤波状态方程进行数值积分计算。

而且，预先根据脉冲星的可见性和优选指标值选择导航脉冲星组合。

而且，模拟的标称轨道根据STK仿真得到。

而且，将步骤3滤波所得探测器的轨道与模拟的标称轨道较差，将步骤3滤波所得钟差 与步骤2模拟所得钟差较差，得到精度评定结果。

本发明还相应提供一种环火探测器精密同步定位守时系统，包括以下模块，

输入模块，用于输入初始数据，包括模拟的标称轨道和用于滤波初始化的相关参数；

观测数据模拟模块，用于进行观测数据模拟，根据脉冲星测量方程模拟脉冲星观测量，该观 测量涉及的测量噪声包括钟差，其中钟差由钟差模型模拟，

设导航脉冲星组合中包括p颗脉冲星，环火探测器的位置矢量r＝[x y z]^T、速度矢量 v＝[v_x v_y v_z]^T，加速度矢量a＝[a_x a_y a_z]^T，

所述脉冲星测量方程如下式，

$Z=\left(\begin{array}{c}\Delta t_1\\ ·\\ ·\\ ·\\ \Delta t_i\\ ·\\ ·\\ ·\\ \Delta t_p\end{array}\right)=H·X+\eta$

其中，Δt_i表示第i颗脉冲星时延观测量，i的取值为1,2,…,p，p颗脉冲星观测量组成观测向 量Z，η为测量噪声，H为观测方程系数矩阵，状态参数X＝[x y z v_x v_y v_z x₁ x₂ x₃]^T；

所述钟差模型如下式，

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+ϵ·x_2\left(t_k\right)+\frac{{ϵ}^2}{2}·x_3\left(t_k\right)+W\left(t_k\right)$

其中，ε为采样间隔，x₁(t_k)、x₂(t_k)和x₃(t_k)分别表示时刻t_k原子钟的钟差、频率漂移和频 率漂移变化率，组成了钟差参数[x₁ x₂ x₃]；W(t_k)为时刻t_k相应的系统噪声，x₁(t_k+1)表示 时刻t_k+1原子钟的钟差；

自适应滤波模块，用于根据环火探测器自主定位和守时的滤波状态方程和观测方程，进行自 适应卡尔曼滤波，

所述观测方程采用脉冲星测量方程，

设时刻t的状态参数X记为X(t)，X(t)的一阶导数为所述滤波状态方程如下，

$\stackrel{·}{X}\left(t\right)=F\left(t\right)·X\left(t\right)+W\left(t\right)$

其中，W(t)为系统状态噪声矩阵，

状态方程系数矩阵F(t)为，$F\left(t\right)=\left(\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y}& \frac{\partial a_x}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_z}{\partial x}& \frac{\partial a_z}{\partial y}& \frac{\partial a_z}{\partial z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$

输出模块，用于输出自适应滤波模块所得滤波结果，包括滤波所得探测器的轨道和钟差。

本发明考虑利用X射线脉冲星观测量，结合探测器的轨道动力学模型和星载原子钟的钟 差模型，通过自适应滤波快速、高精度地估算环火探测器的位置、速度和钟差参数，从而校 正星载原子钟，削弱或消除钟差对定位精度的影响，实现环火探测器的精密同步定位和守时。 本发明还改进了自适应滤波实现方式，进一步提高效率。本发明是在国家自然科学基金支持 下完成的，具有重要的实际推广价值和应用前景，对国民经济的发展和人民生活水平的提高 有不可忽视的作用。

附图说明

图1是现有技术利用X射线脉冲星观测量进行探测器自主定位的原理图；

图2是本发明实施例的流程图。

具体实施方式

以下结合附图和实施例对本发明技术方案进行具体描述。

本发明考虑，X射线脉冲星观测量的精度相对较高，若能通过滤波算法将动力学模型和 测量模型结合起来，利用精度较高的X射线脉冲星观测量对轨道外推状态值进行实时修正， 便可以高精度地估算环火探测器的轨道状态参数。

在环火探测器自主定位中，火星轨道探测器的动力学模型还不能精确建立，在数值积分 进行轨道外推时一些摄动因素的影响通常被忽略，直接导致外推轨道精度较差，而且由于技 术条件限制，初始轨道参数难以精确获得，这也会引起外推轨道的发散。如何利用精度较高 的X射线脉冲星观测量，通过滤波算法实现探测器状态参数的快速、高精度估算是实现精密 定位的首要问题之一。此外，由于火星探测器飞行时间长，其上所搭载的原子钟若不进行定 期校正将难以保持较高的精度，由X射线脉冲星的测量模型可知，探测器的原子钟钟差对定 位精度影响很大，因此如何有效地实现对钟差的估算，是保持原子钟精度和提高探测器自主 定位精度所面临的有一个关键问题。本发明的解决方案就是建立了一种有效融合X射线脉冲 星观测模型、探测器轨道动力学模型和钟差模型的综合模型，设计了一种自适应扩展卡尔曼 滤波算法，实现探测器位置、速度和钟差参数的同步、快速、高精度解算。具体方案如下： (1)建立基于X射线脉冲星的环火探测器精密自主定位和守时模型，包括测量模型、探测 器轨道动力学模型和星载原子钟的钟差模型。

结合环火探测器的具体情况，建立相应的测量模型；要分析环火探测器受力情况，在综 合精度和计算效率的基础上对一些摄动因素进行取舍，然后确定探测器的轨道动力学模型； 了解星载原子钟的工作原理和性能评估方法，在此基础上对钟差进行有效地建模。

(2)设计自适应扩展卡尔曼滤波算法，将测量模型、轨道动力学模型和原子钟的钟差模型有 效地融合利用，将钟差模型参数引入到滤波算法中，同探测器的位置和速度参数同步解算， 实现快速、高精度定位和守时。

(3)具体实现

包括观测数据模拟和自适应扩展卡尔曼滤波解算两大环节，具体设计如下：

具体实施时，可利用STK软件仿真环火探测器轨道作为“标称轨道”，并从导航脉冲星 数据库中选择合适的导航星，根据轨道测量模型(1-5)式模拟脉冲星观测量，模拟过程中要 顾及时钟误差、观测误差、相对论效应延迟等影响。

结合环火探测器的具体情况，建立探测器的X射线脉冲星观测方程。

在J2000火星质心地球赤道坐标系S_ECI中，若同一时段观测了p颗脉冲星，则环火探测 器自主定位和守时系统的测量方程为：

$Z=\left(\begin{array}{c}\Delta t_1\\ ·\\ ·\\ ·\\ \Delta t_i\\ ·\\ ·\\ ·\\ \Delta t_p\end{array}\right)=H·X+\eta ---(2-1)$

式中，Δt_i表示第i颗脉冲星时延观测量，i的取值为1,2,…,p，p颗脉冲星观测量组成观测向 量Z；η为测量噪声；状态参数X＝[x y z v_x v_y v_z x₁ x₂ x₃]^T，观测方程系数矩阵H为：

$H=\left(\begin{array}{ccccccc}\frac{n_1}{c}& 0& 0& 0& 1& 0& 0\\ & & & ·& & & \\ & & & ·& & & \\ & & & ·& & & \\ \frac{n_i}{c}& 0& 0& 0& 1& 0& 0\\ & & & ·& & & \\ & & & ·& & & \\ & & & ·& & & \\ \frac{n_p}{c}& 0& 0& 0& 1& 0& 0\end{array}\right)---(2-2)$

式中，n_i为第i颗脉冲星的方向矢量；c为光速值。

分析环火探测器的受力情况，确定定轨中要考虑的摄动因素，建立环火探测器的轨道动 力学方程，并利用泰勒级数进行线性化展开。

在J2000火星质心地球赤道坐标系S_ECI中，设环火探测器的位置矢量r＝[x y z]^T、 速度矢量v＝[v_x v_y v_z]^T、加速度矢量a＝[a_x a_y a_z]^T；探测器的状态参数向量为 x(t)＝[x y z v_x v_y v_z]^T，x(t)的一阶导数为ω(t)为t时刻状态噪声，则t时刻它的动力 学方程可表达为：

$\stackrel{·}{x}\left(t\right)=f(x\left(t\right),t)+\omega \left(t\right)---(2-3)$

在仅考虑火星中心引力及J₂项摄动的情况下该方程可作如下表示：

$\stackrel{·}{x}=f(x\left(t\right),t)=\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mu }_M\frac{x}{r_3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{2}\right)}^2[\frac{5z^2}{r^2}-1\left]\right]\\ -{\mu }_M\frac{y}{r^3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{r}\right)}^2[\frac{5z^2}{r^2}-1\left]\right]\\ -{\mu }_M\frac{z}{r^3}[1-\frac{3}{2}{J}_2{\left(\frac{R_M}{r}\right)}^2[\frac{5z^2}{r^2}-3\left]\right]\end{array}\right)---(2-4)$

f(x(t),t)表示探测器动力学方程式，即探测器轨道模型，下文公式中简写为f。

式(2-4)为非线性系统，需要在“平衡点”处通过泰勒级数展开进行线性化得：

$\stackrel{·}{x}\left(t\right)=G\left(t\right)·x\left(t\right)+\omega \left(t\right)---(2-5)$

公式(2-5)中t时刻雅克比矩阵G(t)简写为G，其计算方法：

$G=\frac{\partial f}{\partial x}=\frac{\partial }{\partial x}\left(\begin{array}{c}v\\ a\end{array}\right)=\left(\begin{array}{cc}\frac{\partial v}{\partial r}& \frac{\partial v}{\partial v}\\ \frac{\partial a}{\partial r}& \frac{\partial a}{\partial v}\end{array}\right)---(2-6)$

$\frac{\partial v}{\partial r}=\frac{\partial a}{\partial v}=0_{3\times 3}---(2-7)$

$\frac{\partial v}{\partial v}=I_{3\times 3}---(2-8)$

$\frac{\partial a}{\partial }=\left(\begin{array}{ccc}\frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y}& \frac{\partial a_x}{\partial z}\\ \frac{\partial a_y}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z}\\ \frac{\partial a_z}{\partial x}& \frac{\partial a_z}{\partial y}& \frac{\partial a_z}{\partial z}\end{array}\right)---(2-9)$

式中，0_3×3表示零矩阵；I_3×3表示单位矩阵。

建立原子钟的钟差模型，并将钟差状态方程与动力学方程结合组成滤波状态方程，将钟 差模型参数与探测器的位置、速度参数作为滤波参数。

一般地，原子钟的时间偏差可以用一个二阶多项式表示：

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+ϵ·x_2\left(t_k\right)+\frac{{ϵ}^2}{2}·x_3\left(t_k\right)+W\left(t_k\right)---(2-10)$

式中，ε为采样间隔，x₁(t_k)、x₂(t_k)和x₃(t_k)分别表示t_k时刻原子钟的钟差、频率漂移和 频率漂移变化率，它们组成了钟差参数[x₁ x₂ x₃]；W(t_k)为t_k时刻相应的系统噪声，属于 有色噪声；x₁(t_k+1)表示t_k+1时刻原子钟的钟差。

由式(2-10)可得，以时钟钟差x₁(t_k)、时钟频率漂移x₂(t_k)和频率漂移变化率x₃(t_k)组 成状态量${\left(\begin{array}{ccc}{x}_1& x_2& x_3\end{array}\right)}_k^T,$则钟差模型为：

${\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)}_{k+1}={\Phi }_{k+1,k}{\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)}_k+{\omega }_k---(2-11)$

式中，Φ_k+1,k为状态转移矩阵，该矩阵与采样间隔τ有关，表达式为：

${\Phi }_{k+1,k}=\left(\begin{array}{ccc}1& ϵ& {ϵ}^2/2\\ 0.& 1& ϵ\\ 0& 0& 1\end{array}\right)---(2-12)$

式中，ω_k为t_k时刻钟差噪声向量。

设时刻t的状态参数X＝[x y z v_x v_y v_z x₁ x₂ x₃]^T记为X(t)，X(t)的一阶导数为结合环火探测器轨道动力学方程(2-4)和钟差模型(2-11)可得环火探测器脉冲星自主定位 和守时系统的状态方程为：

$\stackrel{·}{X}\left(t\right)=F\left(t\right)·X\left(t\right)+W\left(t\right)---(2-13)$

式中，W(t)为系统状态噪声矩阵，状态方程系数矩阵F(t)为：

$F\left(t\right)=\left(\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y}& \frac{\partial a_x}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_z}{\partial x}& \frac{\partial a_z}{\partial y}& \frac{\partial a_z}{\partial z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)---(2-14)$

设计自适应扩展卡尔曼滤波器，利用X射线脉冲星观测量修正滤波状态参数预估值，快 速、高精度地解算探测器的位置、速度参数和钟差模型参数。

通过公式(2-1)到公式(2-14)建立了环火探测器自主定位和守时的滤波状态方程和观 测方程，接着需要采用自适应扩展卡尔曼滤波器，进行相关的滤波计算。本发明进一步提出， 在现有自适应扩展卡尔曼滤波器(AEKF滤波器)基础上，还可以改进自适应因子和误差方 差阵的计算方式，以提高滤波收敛速度，抑制滤波发散，参见以下公式(2-18)到公式(2-22)。

线性化后系统状态转移矩阵Φ_k+1,k，观测方程矩阵H_k+1，状态方程和观测方程分别为：

$\left(\begin{array}{c}{x}_{k+1}={\Phi }_{k+1,k}·x_k\\ z_{k+1}=H_{k+1}·x_{k+1}\end{array}\right)---(2-15)$

其中k、k+1对应离散化之前的t_k和t_k+1时刻。

实施例自适应扩展卡尔曼滤波器计算过程设计如下：

●滤波参数初始化

根据现有的技术条件，给出计算初始时刻的状态参数x_k及其方差矩阵P_k，状态误差为 δx_k。

●一步预测值x_k+1,k和修正值δx_k+1,k计算

x_k+1,k＝Φ_k+1,k·x_k      (2-16)

δx_k+1,k＝0      (2-17)

●自适应因子λ_k+1的计算

$M_{k+1}=H_{k+1}{\Phi }_{k+1,k}{P}_k{\Phi }_{k+1,k}^T{H}_{k+1}^T+R_{k+1}---(2-18)$

δz_k+1＝z_k+1-H_k+1x_k+1,k      (2-19)

$N_{k+1}=\delta z_{k+1}\delta z_{k+1}^T---(2-20)$

${\lambda }_{k+1}=\mathrm{Max}\{1,\sqrt{\mathrm{sp}\left(N_{k+1}\right)/\mathrm{sp}\left(M_{k+1}\right)}\}---(2-21)$

式中，H_k+1表示观测方程系数矩阵；符号M_k+1和N_k+1用于计算自适应因子；z_k+1为观测值向 量；预测残差为δz_k+1；R_k+1为观测值误差方差矩阵；sp(.)为求矩阵迹符号。

●一步预测估计误差方差阵P_k+1,k的计算

$P_{k+1,k}={\lambda }_k{\Phi }_{k+1,k}{P}_k{\Phi }_{k+1,k}^T+Q_k---(2-22)$

式中，Q_k为系统状态噪声方差矩阵。

●滤波增益K_k+1的计算

$K_{k+1}=P_{k+1,k}{H}_{k+1}^T{[H_{k+1}{P}_{k+1,k}{H}_{k+1}^T+R_{k+1}]}^{-1}---(2-23)$

●测量更新

δx_k+1＝δx_k+1,k+K_k+1[z_k+1-H_k+1x_k+1,k]      (2-24)

$P_{k+1}=[I-K_{k+1}{H}_{k+1}]P_{k+1,k}{[I-K_{k+1}{H}_{k+1}]}^T+K_{k+1}{R}_{k+1}{K}_{k+1}^T---(2-25)$

式中，P_k+1为滤波状态值x_k+1的方差矩阵，I表示相应阶数的单位矩阵。

●滤波结果x_k+1输出

x_k+1＝x_k+1,k+δx_k+1      (2-26)

在给定探测器初始状态和初始钟差参数的情况下，便可利用上述自适应扩展卡尔曼滤波 器进行定位和守时参数的计算。

本发明所提供方法可采用软件方式实现自动运行。参见图2，实施例的实现流程如下： 步骤1，输入初始数据，包括模拟的标称轨道和用于滤波初始化的相关参数。

模拟的标称轨道可根据STK仿真得到，具体实施时，可以预先准备好STK仿真的环火探测 器轨道数据文件，以作为初始轨道参数供STK模拟得到标称轨道。还可以预先准备候选脉冲 星参数文件，根据导航脉冲星个数，利用现有脉冲星优选程序根据脉冲星的可见性和优选指 标值选择最适合的导航脉冲星组合。

步骤2，进行观测数据模拟，根据脉冲星测量方程模拟脉冲星观测量，该观测量涉及的测量 噪声包括钟差，其中钟差由钟差模型模拟。

根据脉冲星测量方程(2-1)式模拟脉冲星观测量，该观测量一般要加入随机观测误差、 钟差、Roemer延迟和Shapiro延迟，其中钟差由钟差模型(2-10)式模拟，基于预设的时刻t_k的钟差初始参数，即可模拟时刻t_k+1原子钟的钟差。其他误差可根据现有技术预设。

步骤3，根据环火探测器自主定位和守时的滤波状态方程和观测方程，进行自适应卡尔曼滤 波。

然后进入自适应卡尔曼滤波和精度评定模块，滤波状态方程和观测方程分别为式(2-13)和 (2-1)，具体滤波流程按公式(2-16)～(2-26)进行，其中一步预测值的计算(2-16)式可 由轨道积分器按式(2-13)进行数值积分计算，获得测数据预报。轨道积分器可采用现有技 术，本发明不予赘述。

滤波完成后可将计算得到的探测器的轨道与①中STK模拟的标称轨道较差，滤波计算得到的 钟差与③中模拟的钟差较差，完成精度评定。

步骤4，输出滤波结果，包括滤波所得探测器的轨道和钟差。

具体实施时，还可将步骤3滤波所得探测器的轨道与模拟的标称轨道较差，将步骤3滤波所 得钟差与步骤2模拟所得钟差较差，得到精度评定结果。

具体实施时，也可以采用模块化方式提供一种环火探测器精密同步定位守时系统，实施 例所提供系统包括以下模块：

输入模块，用于输入初始数据，包括模拟的标称轨道和用于滤波初始化的相关参数；

观测数据模拟模块，用于进行观测数据模拟，根据脉冲星测量方程模拟脉冲星观测量，该观 测量涉及的测量噪声包括钟差，其中钟差由钟差模型模拟，

设导航脉冲星组合中包括p颗脉冲星，环火探测器的位置矢量r＝[x y z]^T、速度矢量 v＝[v_x v_y v_z]^T，加速度矢量a＝[a_x a_y a_z]^T，

所述脉冲星测量方程如下式，

$Z=\left(\begin{array}{c}\Delta t_1\\ ·\\ ·\\ ·\\ \Delta t_i\\ ·\\ ·\\ ·\\ \Delta t_p\end{array}\right)=H·X+\eta$

其中，Δt_i表示第i颗脉冲星时延观测量，i的取值为1,2,…,p，p颗脉冲星观测量组成观测向 量Z，η为测量噪声，H为观测方程系数矩阵，状态参数X＝[x y z v_x v_y v_z x₁ x₂ x₃]^T；

所述钟差模型如下式，

$x_1\left(t_{k+1}\right)=x_1\left(t_k\right)+ϵ·x_2\left(t_k\right)+\frac{{ϵ}^2}{2}·x_3\left(t_k\right)+W\left(t_k\right)$

其中，ε为采样间隔，x₁(t_k)、x₂(t_k)和x₃(t_k)分别表示时刻t_k原子钟的钟差、频率漂移和频 率漂移变化率，组成了钟差参数[x₁ x₂ x₃]；W(t_k)为时刻t_k相应的系统噪声，x₁(t_k+1)表示 时刻t_k+1原子钟的钟差；

自适应滤波模块，用于根据环火探测器自主定位和守时的滤波状态方程和观测方程，进行自 适应卡尔曼滤波，

所述观测方程采用脉冲星测量方程，

设时刻t的状态参数X记为X(t)，X(t)的一阶导数为所述滤波状态方程如下，

$\stackrel{·}{X}\left(t\right)=F\left(t\right)·X\left(t\right)+W\left(t\right)$

其中，W(t)为系统状态噪声矩阵，

状态方程系数矩阵F(t)为，$F\left(t\right)=\left(\begin{array}{ccccccccc}0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_x}{\partial y}& \frac{\partial a_x}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_x}{\partial x}& \frac{\partial a_y}{\partial y}& \frac{\partial a_y}{\partial z}& 0& 0& 0& 0& 0& 0\\ \frac{\partial a_z}{\partial x}& \frac{\partial a_z}{\partial y}& \frac{\partial a_z}{\partial z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right)$

输出模块，用于输出自适应滤波模块所得滤波结果，包括滤波所得探测器的轨道和钟差。

各模块具体实现与方法步骤相应，本发明不予赘述。

以上实施例仅供说明本发明之用，而非对本发明的限制，有关技术领域的技术人员，在 不脱离本发明的精神和范围的情况下，还可以作出各种变换或变型，因此所有等同的技术方 案也应该属于本发明的范畴之内，应由各权利要求限定。