-
X射线脉冲星导航是利用X射线敏感器[1]收集X射线脉冲星的微弱辐射信号,累积出脉冲星轮廓,并估计出脉冲到达时间(Time-Of-Arrival,TOA)来实现导航的[2-5]。但是,航天器高速飞行会改变接收到的脉冲星辐射信号周期,导致脉冲轮廓出现畸变,严重影响脉冲TOA估计的准确性。如何在高速飞行条件下确保X射线脉冲星导航性能已成为一个研究热点。解决该问题的方法可分为两类:脉冲TOA偏差建模法和脉冲轮廓畸变检测法。
脉冲TOA偏差建模法的解决思路如下:研究脉冲TOA偏差的形成机理,建立航天器速度估计误差对其的影响机制模型,将该影响机制模型融入到X射线脉冲星导航量测模型中,构成新的量测模型。基于新量测模型,优化卡尔曼(Kalman)滤波器,如:EKF-CMBSEE(Extended Kalman Filter- Correlated Measurement Bias and State Estimation Error)[6]、闭环卡尔曼滤波器[7]以及跟踪滤波器[8]等,以达到抑制脉冲星信号中的多普勒效应的目的。在此类方法中,航天器速度并未被直接测量,而是由脉冲TOA偏差来体现。本文将脉冲TOA偏差建模法又称为间接测速法。
脉冲轮廓畸变检测法按历元折叠次数可分为两类:多次历元折叠和单次历元折叠。多次历元折叠的脉冲轮廓畸变检测法按多个周期折叠信号,可得多个畸变度的脉冲累积轮廓,最小畸变度对应的周期视为固有周期。时域χ2验证法、H验证法[9]、快速折叠法[10]、傅里叶频谱法[11]、循环平稳信号相干统计量法[12]、最大相关方差搜索法[13]、基于Lomb的χ2法[14]、脉冲轮廓熵法[15]、轮廓特征法[16]、蝶形算法[17-18]均属此类。但星载计算机计算量有限,无法承受多次历元折叠。因而,此类方法难以在轨实现。单次历元折叠的脉冲轮廓畸变检测法则是根据畸变度直接估计出固有周期,避免了多次历元折叠。基于FFT-CS(Fast Fourier Transform- Compressed Sensing)的脉冲星周期估计方法[19],基于EMD-CS(Empirical Mode Decomposition - Compressed Sensing)的脉冲星周期估计方法[20],直接提供了速度信息,本文将其简称为直接测速法。由于脉冲轮廓的畸变是一种微弱特征,当X射线脉冲星的信噪比高时,脉冲轮廓畸变检测法可正常工作;反之不能。在众多脉冲星中,仅Crab脉冲星的光子流量高于X射线背景流量,可用于脉冲轮廓畸变检测法,即直接测速法。
直接测速法可直接提供速度,但仅对Crab脉冲星有作用;间接测速法虽对X射线脉冲星光子流量无要求,但是无法直接提供速度信息。X射线脉冲星导航系统至少要观测3颗脉冲星才能实现三维定位[21],除Crab脉冲星外,还需两颗低流量脉冲星,因此,对Crab脉冲星采用直接测速法,低流量脉冲星则采用间接测速法。
鉴于直接与间接测速法二者的互补性,本文利用直接测速与间接测速方法轮流对卡尔曼滤波器更新,以提供高精度导航信息。从单次历元折叠的脉冲轮廓畸变检测法的基本原理出发,提出了一种面向直接测速测距的畸变轮廓匹配估计方法。该方法能一次性解算出速度和距离等导航信息。
-
根据脉冲星畸变轮廓形成机理[22],利用一个脉冲星畸变轮廓直接解算出距离和速度信息,提出一种面向直接测速测距的畸变轮廓匹配估计方法。该方法包括脉冲星畸变轮廓集的构建和超分辨率估计。
将门函数与脉冲星标准轮廓循环互相关即可得到脉冲星畸变轮廓。设畸变度为m的脉冲星畸变轮廓hm为
$$ {h}_{m}={\mathit{G}}_{m}\star s $$ (1) 其中:
$\star $ 表示循环相关;m={0,1, ···,MD–1};MD为最大畸变度;s(N×1)是脉冲星标准轮廓;Gm(N×1)是门函数。门函数为$$ {{\boldsymbol{G}}_m}(n) = \left\{ \begin{gathered} \frac{1}{m},0 \leqslant n \leqslant m \hfill \\ 0,m{\text{ + }}1 \leqslant n \leqslant N - 1 \hfill \\ \end{gathered} \right. $$ (2) 将第m个脉冲星畸变轮廓循环移位构成第m个脉冲星畸变轮廓子集
$ {{\boldsymbol{\varphi }}_m}(N \times P) $ ,公式为$$ \begin{gathered} {{\boldsymbol{\varphi }}_m} = \left\{ {{{\boldsymbol{\varphi }}_{mi}}(t)|{{\boldsymbol{\varphi }}_{mi}}(t) = {{{{\boldsymbol{h}}_m}(t \pm i)} \mathord{\left/ {\vphantom {{{{\boldsymbol{h}}_m}(t \pm i)} {\left\| {{{\boldsymbol{h}}_m}} \right\|}}} \right. } {\left\| {{{\boldsymbol{h}}_m}} \right\|}}} \right\}\quad \hfill \\ {\text{ }}i = 0,1,2 ,\cdots ,(P - 1)/2 \hfill \\ \end{gathered} $$ (3) 其中:
$ ||\cdot ||$ 为2范数;P是最大相位偏移量;$ {{\boldsymbol{\varphi }}_{mi}}(t) $ 是子集$ {{\boldsymbol{\varphi }}_m} $ 中第i个原子。可将$ {{\boldsymbol{\varphi }}_m} $ 构成一个三维集Ψ(N×P×MD)。将脉冲星累积轮廓
$ {\boldsymbol{\tilde h}} $ (N×1)与三维集Ψ相乘,可得到匹配矩阵C。$$ {\boldsymbol{C}} = {\boldsymbol{\tilde h}} \cdot {\boldsymbol{\varPsi }} $$ (4) 设C的最大值对应脉冲星累积轮廓相位和畸变度的索引值分别
$ \tilde p $ 、$ \tilde m $ 。$$ [\tilde p,\tilde m] = \mathop {\arg \max }\limits_{i,j} [{\boldsymbol{C}}(i,j)] $$ (5) 利用超分辨率估计,可得高精度频率估计为
$$ \hat f = \frac{1}{{N \times {T_{\text{obs}}}}} \times \left\{ {\tilde m - \frac{{{\text{0}}{\text{.5[max(}}{{\boldsymbol{C}}_{\tilde m + 1}}{\text{)}} - {\text{max}}({{\boldsymbol{C}}_{\tilde m - 1}}{\text{)]}}}}{{{\text{max(}}{{\boldsymbol{C}}_{\tilde m + 1}}{\text{)}} + {\text{max}}({{\boldsymbol{C}}_{\tilde m - 1}}{\text{)}} - 2{\boldsymbol{C}}(\tilde p,\tilde m)}}} \right\} $$ (6) 其中:
$ {{\boldsymbol{C}}_{\tilde m + 1}} $ 和$ {{\boldsymbol{C}}_{\tilde m - 1}} $ 分别表示匹配矩阵C中畸变度索引值分别是$ \tilde m + 1 $ 和$ \tilde m - 1 $ 时的最大元素值;Tobs是观测时间。将频率估计转化为高精度速度为
$$ \begin{gathered} \hat v = c \times {T_ {P}} \times \hat f = \hfill \\ \frac{{c \times {T_{P}}}}{{N \times {T_{\text {obs}}}}} \times \left\{ {\tilde m - \frac{{{\text{0}}{\text{.5[max(}}{{\boldsymbol{C}}_{\tilde m + 1}}{\text{)}} - {\text{max}}({{\boldsymbol{C}}_{\tilde m - 1}}{\text{)]}}}}{{{\text{max(}}{{\boldsymbol{C}}_{\tilde m + 1}}{\text{)}} + {\text{max}}({{\boldsymbol{C}}_{\tilde m - 1}}{\text{)}} - 2{\boldsymbol{C}}(\tilde p,\tilde m)}}} \right\} \hfill \\ \end{gathered}$$ (7) 其中:TP为脉冲星周期。
超分辨率估计可得脉冲星到达时间t为
$$ \begin{gathered} t= I \times{T_ {P}} + \hfill \\ \frac{{{T_ {P}}}}{N}\left\{ {\tilde p - \frac{{{\text{0}}{\text{.5[}}{\boldsymbol{C}}{\text{(}}\tilde p + 1,\tilde m{\text{)}} - {\boldsymbol{C}}{\text{(}}\tilde p - 1,\tilde m{\text{)]}}}}{{{\boldsymbol{C}}{\text{(}}\tilde p + 1,\tilde m{\text{) + }}{\boldsymbol{C}}{\text{(}}\tilde p - 1,\tilde m{\text{)}} - 2{\boldsymbol{C}}(\tilde p,\tilde m)}} - \hat m/2} \right\} \hfill \\ \end{gathered} $$ (8) 其中:I为整数个脉冲星周期,由航天器导航系统估算得到。值得一提的是,脉冲星信号的频移会引起时移[7]。若要得到高精度的脉冲星到达时间,需对频移引起的时移量作补偿。
${{\hat m \times {T_P}} / {2N}}$ 为频移引起的时移补偿量。面向直接测速测距的畸变轮廓匹配估计方法仅适用于高流量的Crab脉冲星。低流量脉冲星仅需提供到达时间。鉴于传统的互相关法具有较好的抗噪能力,本文采用互相关法估计低流量脉冲星的到达时间,详见文献[23]。
面向直接测速测距的畸变轮廓匹配估计方法的计算复杂度计算量主要集中在式(4),即脉冲星累积轮廓与三维集Ψ相乘,约为N×P×MD次乘法运算与(N–1)×P×MD次加法运算。可以看出,面向直接测速测距的畸变轮廓匹配估计方法的计算复杂量与脉冲星轮廓间隔数、最大相位偏移量和最大畸变度均成正比关系。
互相关法的计算量约为N×P次乘法运算与(N–1)×P次加法运算。与面向直接测速测距的畸变轮廓匹配估计方法相比,互相关法的计算量大幅减小。以Crab脉冲星为例,分析面向直接测速测距的畸变轮廓匹配估计方法。此时,N= 33 400,P=MD= 81。计算量约为2.2×108次乘法运算与2.2 × 108次加法运算。以B1821-24为例,分析互相关法。此时,N= 3 050,P= 701。计算量约为2.1×106次乘法运算与2.1×106次加法运算。若以B1937+21为例,互相关法计算量约为106次乘法运算与106次加法运算。互相关法的计算量远小于面向直接测速测距的畸变轮廓匹配估计方法。对于Crab脉冲星,导航系统以计算量为代价,获得了测速信息;对于其它脉冲星,既然畸变轮廓匹配估计方法无法提供测速信息,不如采用小计算量的互相关法提供定速信息。从计算量上考虑,直接/间接混合测速是一种合理方式。
-
本节建立间接测速和直接测速的脉冲星量测模型,在此基础上,利用闭环卡尔曼滤波器实现X射线脉冲星导航。为便于表达,本文将高流量脉冲星(Crab脉冲星,B0531+21)和两颗低流量脉冲星(B1821-24、B1937+21)编号。PSR B0531+21、B1821-24、B1937+21的编号分别为1、2、3。
-
直接测速模型是直接测量多普勒效应引起的速度,适用于Crab脉冲星。脉冲星到达时间模型可表示为
$$ \begin{aligned}[b] & c\left({t}_{\text {b}}^{1}-{t}^{1}\right)={\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{r}}+\\ & \frac{1}{2{D}_{0}}\left[-{||{\boldsymbol{r}} ||}^{2}+{\left({\boldsymbol n}_{1}^{\text{T}}\cdot r\right)}^{2}-2{\boldsymbol b}^{\text{T}}\cdot{\boldsymbol{r}}+2\left({\boldsymbol n}_{1}^{\text{T}}\cdot{\boldsymbol{b}}\right)\left({\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{r}}\right)\right]+\\ & \frac{2{\mu }_{\text {Sun}}}{{c}^{2}}\mathrm{ln}\left|\frac{{\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{r}}+{\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{b}}+||{\boldsymbol{r+b}}||}{{\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{b}}+||{\boldsymbol{b}}||}+1\right| \end{aligned} $$ (9) 其中:c为光速;n1为Crab脉冲星方位矢量;r为航天器位置;
$ {D}_{0} $ 为脉冲星到太阳系质心的距离;${\boldsymbol b}$ 为太阳系质心(Solar System Barycentre,SSB)相对于太阳的位置矢量;${\mu _{{\text{Sun}}}}$ 为太阳引力常数;tb和t分别为SSB处和航天器上的脉冲星到达时间。在直接测速与到达时间模型中,量测量
$ {{\boldsymbol{Y}}_{\text{D}}} $ 包括达到时间$ t_b^1 - {t^1} $ 与速度$\boldsymbol{v}$ ,如式(9)所示。$$ {{\boldsymbol{Y}}_{\text{D}}} = \left[ \begin{gathered} c\left( {t_b^1 - {t^1}} \right) \hfill \\ \boldsymbol {v} \hfill \\ \end{gathered} \right] $$ (10) 直接测速与到达时间模型可表示为
$$ {{\boldsymbol{Y}}_{\text{D}}} = {h_{\text{D}}}\left( {{{\boldsymbol{X}}},t} \right) + {{\boldsymbol{\xi }}_{\text{D}}} $$ (11) $$ \begin{gathered} {{\boldsymbol{h}}_{\text{D}}}\left( {{{\boldsymbol{X}}},t} \right) = \hfill \\ \left[ \begin{gathered} {\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{r}} + \frac{1}{{2D_0}}\left[ { - {{|| {\boldsymbol{r}} ||}^2} + {{\left( {{\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{r}}} \right)}^2} - 2{{\boldsymbol{b}}^{\text{T}}}\cdot{\boldsymbol{r}} + 2\left( {{\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{b}}} \right)\left( {{\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{r}}} \right)} \right] + \hfill \\ {{ }}\frac{{2{\mu _{\rm Sun}}}}{{{c^2}}}\ln \Bigg| {\frac{{{\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{r}} +{\boldsymbol n}_{1}^{\text{T}}\cdot {\boldsymbol{b}}+ || {\boldsymbol{r+b}} ||}}{{{\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{b}} + || {\boldsymbol{b}} ||}} } \Bigg| \hfill \\ {\boldsymbol{n}}_1^{\text{T}} \cdot {\boldsymbol{v}} \hfill \\ \end{gathered} \right] \hfill \\ \end{gathered} $$ (12) 其中:X为状态矢量,包括航天器位置r和速度v这两个分量。
$$ {\boldsymbol{X}}{\text{ = }}\left[ \begin{gathered} {\boldsymbol{r}} \hfill \\ {\boldsymbol{v}} \hfill \\ \end{gathered} \right] $$ (13) 测量噪声矢量
$ {{\boldsymbol{\xi }}_{\text{D}}} $ 包括达到时间与速度噪声两个分量$$ {{\boldsymbol{\xi }}_{\text{D}}} = {\left[ {\varsigma _{{\text{TOA}}}^1,{\varsigma _{\text{v}}}} \right]^{\text{T}}} $$ (14) 其中:
$ \varsigma _{{\text{TOA}}}^1 $ 和$ {\varsigma _{\text{v}}} $ 分别为达到时间噪声与速度噪声。 -
间接测速模型建立航天器速度误差与脉冲星到达时间之间的关系。考虑脉冲星TOA偏差
${\beta ^j}$ 的脉冲到达时间模型[7]可表示为$$ \begin{aligned}[b] & c\left({t}_{b}^{j}-{t}^{j}\right)={\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol r}+\\ & \frac{1}{2{D}_{0}^{j}}\left[-{||{\boldsymbol r}||}^{2}+{\left({\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol r}\right)}^{2}-2{\boldsymbol b}^{\text{T}}\cdot{\boldsymbol r}+2\left({\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol{b}}\right)\left({\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol r}\right)\right]+\\ & \frac{2{\mu }_{\text{Sun}}}{{c}^{2}}\mathrm{ln}\Bigg|\frac{{\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol r}+{\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol{b}}+||{\boldsymbol r+\boldsymbol b}||}{{\boldsymbol n}_{j}^{\text{T}}\cdot {\boldsymbol{b}}+||{\boldsymbol{b}}||}+1\Bigg|-c{\beta }^{j} \end{aligned}$$ (15) 值得一提的是,式(8)与式(14)之间的差异在于脉冲星TOA偏差
${\beta ^j}$ 。在间接测速模型中,量测量
$ {{\boldsymbol{Y}}_{\text{I}}} $ 为两个弱脉冲星的到达时间$ t_b^2 - {t^2} $ 和$ t_b^3 - {t^3} $ 之一,如式(15)所示。$$ {{\boldsymbol{Y}}_{\text{I}}} = \left[ {c\left( {t_b^j - {t^j}} \right)} \right] $$ (16) 其中:j=2,3。
间接测速模型可表示为
$$ {{\boldsymbol{Y}}_{\text{I}}} = h_{\text{I}}^j\left( {{{\boldsymbol{X}}},t} \right) - c{\beta ^j} + {\boldsymbol{\xi }}_{\text{I}}^j $$ (17) $$\begin{aligned}[b] &{h}_{\text{I}}^{j}\left({{\boldsymbol{X}}}_,t\right)={{{\boldsymbol{n}}}_{j}^{\rm{T}}}\cdot{\boldsymbol{r}}+\\ &\frac{1}{2{D}_{0}^{j}}\left[-{||{\boldsymbol{r}}||}^{2}+{\left({{\boldsymbol n}_{j}^{\rm T}}\cdot{\boldsymbol{r}}\right)}^{2}-2{\boldsymbol b}^{\text{T}}\cdot{\boldsymbol{r}}+2\left({{{\boldsymbol{n}}}_{j}^{\rm{T}}}\cdot{\boldsymbol{b}}\right)\left({{\boldsymbol n}_{j}^{\text{T}}}\cdot{\boldsymbol{r}}\right)\right]+\\ &\frac{2{\mu }_{\text {Sun}}}{{c}^{2}}\mathrm{ln}\Bigg|\frac{{{\boldsymbol n}_{j}^{\text{T}}}\cdot{\boldsymbol{r}}+{{{\boldsymbol{n}}}_{j}^{\rm{T}}}\cdot{\boldsymbol{b}}+||{\boldsymbol{r+b}}||}{{{\boldsymbol n}_{j}^{\text{T}}}\cdot{\boldsymbol{b}}+||{\boldsymbol{b}}||}\Bigg| \end{aligned}$$ (18) 测量噪声
$ {\boldsymbol{\xi }}_{\text{I}}^j $ 是弱脉冲星的到达时间噪声。脉冲星TOA偏差$\; {\beta ^j} $ 可表示为$$ - c{\beta ^j} = {\boldsymbol{n}}_j^{\text{T}}\left[ {\frac{{5T_{{\text{obs}}}^2}}{6}{\boldsymbol{S}},{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{{T_{{\text{obs}}}}}}{2}{{\boldsymbol{I}}_{3 \times 3}} + \frac{{T_{{\text{obs}}}^3}}{3}{\boldsymbol{S}}} \right]\delta {\boldsymbol{X}} $$ (19) 其中:
$ \delta {\boldsymbol{X}} $ 为状态误差;S为轨道转移矩阵。 -
考虑到TOA测量与状态误差相关,采用误差状态EKF。在这种滤波器中,状态预测误差被作为新的状态估计。采用反馈修正的状态预测误差滤波,是将状态预测误差估计反馈到轨道动力学模型中,对状态进行修正。X射线敏感器载重大,因此,X射线脉冲星导航系统只安装一个X射线敏感器,对3颗导航脉冲星轮流观测[24]。
设X射线脉冲星系统的状态转移模型和量测模型分别为
$$ \delta {{\boldsymbol{X}}_k} = {{\boldsymbol{\varPhi }}_k}\delta {{\boldsymbol{X}}_{k - 1}} + {{\boldsymbol{\omega }}_{k - 1}} $$ (20) $$ \delta {{\boldsymbol{Y}}_k} = {{\boldsymbol{H}}_k}\delta {{\boldsymbol{X}}_k} + {{\boldsymbol{\delta }}_k} $$ (21) 其中:
$ {\boldsymbol{\varPhi }} $ 为状态转移矩阵,通过对轨道动力学模型微分得到;$ {\boldsymbol{\omega }} $ 为状态处理噪声,其协方差矩阵为Q;量测矩阵为H;测量噪声为${\boldsymbol{\xi }}$ ;其对应的协方差矩阵为R。对于弱脉冲星,测量误差
$ \delta {{\boldsymbol{Y}}_k} $ 可表示为$$\begin{aligned}[b] & \delta {\boldsymbol Y}_{k}={\boldsymbol Y}_{\text{I}}-{\boldsymbol h}_{\text{I}}^{j}\left({\tilde{\boldsymbol X}}_{k},k\right)=\\ & \left[{\boldsymbol n}_{j}^{\text{T}},{\boldsymbol 0}_{3\times 3}\right]\delta {X}_{k}-c{\beta }^{j}+{\boldsymbol {\xi} }_{\text{I}}^{j}=\\ & {\boldsymbol n}_{j}^{\text{T}}\left[{\boldsymbol I}_{3\times 3}+\frac{5{T}^{2}}{6}{\boldsymbol S}_{k-1},\frac{T}{2}{\boldsymbol I}_{3\times 3}+\frac{{T}^{3}}{3}{\boldsymbol S}_{k-1}\right]\delta {\boldsymbol X}_{k}+{\boldsymbol {\xi} }_{\text{I}}^{j}=\\ & {\boldsymbol H}_{\text{I}}\delta {\boldsymbol X}_{k}+{\boldsymbol {\xi} }_{\text{I}}^{j} \end{aligned}$$ (22) 其中:
${{\boldsymbol{H}}_{\text{I}}}$ 为间接测速量测矩阵。$$ {{\boldsymbol{H}}_{\text{I}}} = {\boldsymbol{n}}_j^{\text{T}}\left[ {{{\boldsymbol{I}}_{3 \times 3}} + \frac{{5{T^2}}}{6}{{\boldsymbol{S}}_{k - 1}},{\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{T}{2}{{\boldsymbol{I}}_{3 \times 3}} + \frac{{{T^3}}}{3}{{\boldsymbol{S}}_{k - 1}}} \right] $$ (23) 对于Crab脉冲星,测量误差
$ \delta {{\boldsymbol{Y}}_k} $ 可表示为$$\begin{aligned}[b] & \delta {\boldsymbol Y}_{k}={\boldsymbol Y}_{\text{D}}-{h}_{\text{D}}\left({\tilde{\boldsymbol X}}_{k},k\right)=\\ &\left[{\boldsymbol n}_{1}^{\text{T}},{\boldsymbol n}_{1}^{\text{T}}\right]\delta {\boldsymbol X}_{k}+{\xi }_{\text{D}}={\boldsymbol H}_{\text{D}}\delta {\boldsymbol X}_{k}+{\xi }_{\text{D}} \end{aligned}$$ (24) 其中:
${{\boldsymbol{H}}_{\text{D}}}$ 为直接测速测距量测矩阵。$$ {{\boldsymbol{H}}_{\text{D}}} = \left[ {{\boldsymbol{n}}_j^{\text{T}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\boldsymbol{n}}_j^{\text{T}}} \right] $$ (25) 误差状态EKF的原理如图1所示。轨道动力学模型需要利用当前状态预测下一个状态。Kalman滤波器提供状态预测误差,用于修正轨道动力学模型中的状态。各量测量轮流进入卡尔曼滤波器对状态进行修正,实现了数据融合。
误差状态EKF的滤波公式为
$$ {\boldsymbol{P}}_k^ - = {{\boldsymbol{\varPhi }}_{k - 1}}{\boldsymbol{P}}_{k - 1}^ + {\boldsymbol{\varPhi }}_{k - 1}^\text {T} + {{\boldsymbol{Q}}_{k - 1}} $$ (26) $$ {{\boldsymbol{K}}_k} = \left( {{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}}} \right){\left( {{{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - {\boldsymbol{H}}_k^{\text{T}} + {{\boldsymbol{R}}_k}} \right)^{ - 1}} $$ (27) $$ \delta {{\boldsymbol{X}}_k} = {{\boldsymbol{K}}_k}\delta {{\boldsymbol{Y}}_k} $$ (28) $$ {\boldsymbol{P}}_k^ + = {\boldsymbol{P}}_k^ - - {{\boldsymbol{K}}_k}{{\boldsymbol{H}}_k}{\boldsymbol{P}}_k^ - $$ (29) $$ {\boldsymbol{X}}_k^ + = {\boldsymbol{X}}_k^ - + \delta {{\boldsymbol{X}}_k} $$ (30) X射线脉冲星导航系统利用一个X射线敏感器轮流观测3颗X射线脉冲星。闭环EKF的滤波中的测量误差
$ \delta {{\boldsymbol{Y}}_k} $ 、量测矩阵H以及测量噪声协方差矩阵R需要与脉冲星对应。 -
为了验证基于直接/间接混合测速的X射线脉冲星导航的可行性和有效性,本节通过仿真将其与基于间接测速的X射线脉冲星导航比较。
-
本文的直接测速法是畸变轮廓匹配法。仅Crab脉冲星可使用该方法,其仿真参数如表1所示。表1同时给出了两颗低流量脉冲星的仿真参数。其中脉冲星周期由参考文献[21]提供,光子流量密度来自文献[25],噪声流量密度、分辨率、探测器有效面积均来自文献[26]。航天器轨道采用环火轨道,如表2所示。基于直接/间接混合测速的X射线脉冲星导航滤波器参数如表3所示。X射线脉冲星定位和定速精度由3.2节估计得到。Crab脉冲星的定位和定速的相关系数为–0.626,由仿真得到。
表 1X射线脉冲星仿真参数
Table 1.Simulation parameters of X-ray pulsars
参数 数值 脉冲星 B0531+21 B1821-24 B1937+21 周期/s 0.033 40 0.003 05 0.001 56 光子流量密度/(ph·cm–2·s–1) 1.54 1.93×10–4 4.99×10–5 噪声流量密度/(ph·cm–2·s–1) 0.005 0.005 0.005 分辨率/μs 1 1 1 探测器有效面积/m2 1 1 1 观测时间/s 1 000 1 000 1 000 畸变轮廓数 81 — — 最大相位偏移量 81 701 701 表 2环火轨道初始六要素
Table 2.Initial orbital elements of the Mars
轨道六要素 数值 半长轴/km 6 794 偏心率 0 轨道倾角/(°) 45 升交点赤经/(°) 0 近地点幅角/(°) 0 真近点角/(°) 0 表 3导航滤波器参数
Table 3.Parameters of navigation filter
参数 数值 X射线敏感器数量 1 观测周期 1 000 s 初始轨道误差 $\begin{aligned} & \delta {\boldsymbol{X} }\left( 0 \right) = [6\;000\;{\text{m} },6\;000\;{\text{m} },6\;000\;{\text{m,} } \hfill \\& 20\;{\text{m/s} },20\;{\text{m/s} },15\;{\text{m/s} }{]^{\text{T} } } \hfill \end{aligned}$ 状态处理噪声协方差矩阵 ${\boldsymbol{Q } } = {\text{diag} }\left( {q_1^2,q_1^2,q_1^2,q_2^2,q_2^2,q_2^2} \right)$,
${q}_{1}=6\; \text{m},$$ {q_2} = 0.04\;{\text{m/s}} $测量噪声协方差矩阵 B0531+21:${\boldsymbol{R} } = \left[ \begin{gathered} \sigma _{\text{D} }^2 \;\;\;\;\;\;- 0.626{\sigma _{\text{D} } }{\sigma _{\text{v} } } \hfill \\ - 0.626{\sigma _{\text{D} } }{\sigma _{\text{v} } }\;\;\;\;\;\;\sigma _\text v^2 \hfill \\ \end{gathered} \right]$
$ {\sigma _{\text{D}}}{\text{ = }}89\;{\text{m}},{\sigma _{\text{v}}} = 0.1\;{\text{m/s}} $
B1821-24:${\boldsymbol{R}}=\left[{\sigma }_{\text{D} }^{2}\right],$$ {\sigma _{\text{D}}}{\text{ = 624}}\;{\text{m}} $
B1937+21:${\boldsymbol{ R}}=\left[{\sigma }_{\text{D} }^{2}\right],$$ {\sigma _{\text{D}}}{\text{ = 1}}\;{\text{463}}\;{\text{m}} $假设导航滤波的速度估计误差为
${\sigma _\text v}$ , 其引起的畸变度为$\left( {{\sigma _{\text v}} \times {T_{{\text{obs}}}}} \right)/\left( {{\varDelta _{{\text{bin}}}} \times c} \right)$ ,其中,${\varDelta _{{\text{bin}}}}$ 为时间分辨率。为确保Crab脉冲星的测速估计值在字典范围内,考虑正负号,最大畸变度应设置为$10\left( {{\sigma _{\text v}} \times {T_{{\text{obs}}}}} \right)/\left( {{\varDelta _{{\text{bin}}}} \times c} \right)$ 。假设定位估计误差为${\sigma _{\text{P}}}$ , 其引起的相位偏移量为${\sigma _{\text{P}}}/\left( {{\varDelta _{{\text{bin}}}} \times c} \right)$ 。为确保Crab脉冲星定位估计值在字典范围内,考虑正负号,最大相位偏移量应为${\sigma _{\text{P}}}/\left( {{\varDelta _{{\text{bin}}}} \times c} \right)$ 的10倍。对于低流量脉冲星,最大相位偏移量为${\sigma _{\text{p}}}/\left( {{\varDelta _{{\text{bin}}}} \times c} \right)$ 的4倍。根据后文的仿真试验得,Crab脉冲星的测速误差$\sigma_{\text v}$ 最大约为1.2 m/s。Crab脉冲星的定位误差$\sigma_{\rm{p}} $ 最大为2.1 km;对于低流量脉冲星,$\sigma_{\rm{p}} $ 最大为50 km。 -
畸变轮廓匹配估计方法适用于强流量密度的X射线脉冲星即Crab脉冲星,可同时估计出位置和速度信息。而互相关算法对高流量和低流量脉冲星均有效。本文将互相关算法用于低流量脉冲星的定位估计。图2给出了3颗脉冲星累积轮廓。从图2中可以看出,Crab脉冲星轮廓信噪比高,能从中检测出轮廓畸变;而两颗弱脉冲星轮廓信噪比低,仅能检测出相移。这是为何对弱脉冲星仅采用间接测速方法的原因。
畸变轮廓匹配法和互相关法的仿真结果如表4所示。运行平台为Intel笔记本,CPUi7-8750H主频2.20 GHz,8 G内存。
表 4畸变轮廓匹配估计
Table 4.Matching estimation of distorted profiles
脉冲星 方法 定位精度/m 定速精度/(m·s–1) 计算时间/s B0531+21 畸变轮廓匹配 89 0.097 6 1.413 量子CS 102 0.098 8 1.621 B1821-24 互相关 624 / 0.001 4 B1937+21 互相关 1463 / 0.000 9 从表4中可以看出,对于Crab脉冲星(B0531+21),畸变轮廓匹配法提供的定位和定速误差分别在100 m和0.1 m/s以下,而B1821-24和B1937+21仅能提供1 km量级的定位精度,无法提供速度信息。为体现畸变轮廓匹配法的优越性,本文将其与量子压缩感知法比较。畸变轮廓匹配法的定位定速精度和计算量均优于量子压缩感知法。
Crab脉冲星的计算时间远大于另两颗低流量脉冲星,因为低流量脉冲星无需计算速度,计算量大幅下降。此外,3颗脉冲星的匹配估计的计算时间小于1.5 s。航天器计算机的主频约为100 MHz,而本文计算机的主频为2.2 GHz。假设计算时间与主频成反比,那么3颗脉冲星的匹配估计在航天器计算机上的计算时间约为33 s,远小于脉冲星观测时间。这表明匹配估计具有较好的实时性。
-
本文将基于直接/间接混合测速的X射线脉冲星导航方法与基于间接测速的X射线脉冲星导航方法[6]对比。二者滤波结果如图3所示。
从图3可看出,两种导航方式均收敛。这表明两种导航方式均能正常工作。此外,基于直接/间接混合测速的X射线脉冲星导航误差略低于基于间接测速的X射线脉冲星导航误差。
表5列出了两种导航方式100次蒙特卡洛试验的结果。从表5可以看出,与基于间接测速的X射线脉冲星导航相比,基于直接/间接混合测速的X射线脉冲星导航的定位和定速精度分别提高了16.72%和14.59%。这表明直接/间接混合测速方式有利于提高导航精度。究其原因,直接测速方法提供了速度信息,为X射线脉冲星导航滤波器提供了新的观测量。
表 5两种导航方式的对比
Table 5.Comparison of two navigation methods
导航方式 定位精度/m 定速精度/(m·s-1) 间接测速导航 6 871 2.823 直接/间接混合测速导航 5 722 2.411 表6给出了100 s和500 s观测时间下的导航性能。可见在不同观测时间下,基于直接/间接混合测速的X射线脉冲星导航方法都优于基于间接测速的X射线脉冲星导航方法。
表 6不同观测时间
Table 6.Different observation times
导航方式 观测时间100 s 观测时间500 s 定位精度/
m定速精度/
(m·s–1)定位精度/
m定速精度/
(m·s–1)间接测速导航 12 718 4.436 7 549 2.991 直接/间接
混合测速导航12 482 4.370 7 164 2.771 与表5所示的1 000 s观测时间相比,100 s和500 s观测时间的导航精度略有下降。究其原因,B1937+21和B1821-21定位精度随观测时间的缩短而大幅下降,进而导致整个导航性能下降。理论上,脉冲星的定位误差与时间的平方根成反比。而仿真中PSR B1937+21和B1821-21的定位误差随时间的下降速度快于理论值。
X射线探测器有效面积与导航之间的关系如图4所示。由图4可知,X射线探测器有效面积的增大能提高导航精度。混合测速导航方法明显优于间接测速导航方法。究其原因,混合测速利用了Crab脉冲星测速信息。当X射线探测器有效面积较小时,混合测速导航方法优势明显。
-
鉴于脉冲星直接测速和间接测速二者各有优势,将这两种测速方法复合使用,提出了基于直接/间接混合测速的X射线脉冲星导航方法。直接测速法精度高,但仅适用于Crab脉冲星;间接测速方法则适用于低流量脉冲星。混合利用直接测速法与间接测速法,可充分发挥二者各自的优势,在抑制多普勒效应的同时有效提高精度。此外,为进一步提高实时性与精度,提出了面向直接测速测距的畸变轮廓匹配估计方法。该方法具有简单、实时、高精度等优点。
基于直接/间接混合测速的X射线脉冲星导航具有以下优点:高精度,与间接测速脉冲星导航相比,基于直接/间接混合测速的X射线脉冲星导航方法的精度提高了15%,因发挥了两种测速法各自的优势;小计算量,匹配估计的计算时间小于1.5 s,可以满足实时性要求。基于直接/间接混合测速的X射线脉冲星导航具有实时高精度的特点,能满足深空探测任务对导航的要求。
X-Ray Pulsar Navigation Based on Hybrid Velocimetry Method
-
摘要:为抑制由航天器速度引起的脉冲星信号多普勒效应,提出了基于直接/间接混合测速的X射线复合脉冲星导航方法。直接测速方法根据脉冲星畸变轮廓直接提供速度信息,提出了一种面向直接测速测距的畸变轮廓匹配估计方法。直接测速精度高,但仅适用于高流量脉冲星,即Crab脉冲星。间接测速方法通过建立速度与脉冲星到达时间之间的关系模型,由卡尔曼滤波器间接提供速度信息。低流量脉冲星畸变轮廓信噪比低,无法体现脉冲星畸变这一微弱特征。因此,间接测速方法适用于低流量脉冲星。在此基础上,利用直接测速与间接测速方法轮流对滤波器更新,卡尔曼滤波器提供高精度导航信息。仿真结果表明,与间接测速脉冲星导航相比,基于直接/间接混合测速的X射线脉冲星导航方法的精度提高了15%。Abstract:In order to suppress the Doppler effects in pulsar signals caused by spacecraft velocity, a novel method of X-ray pulsar navigation based on direct/indirect hybrid velocimetry was proposed. The direct velocimetry method provided velocity information directly according to the distorted pulsar profile. It had high accuracy, but was only suitable for high flux pulsar, namely Crab pulsar. A distorted profile matching estimation method for direct velocity and distance measurement was proposed. The indirect velocimetry method was to establish the relationship model between velocity and pulsar time-of-arrival, and provide velocity information by Kalman filter indirectly. The signal-to-noise ratio of low flux pulsar distortion profile was very low, which could not reflect the weak pulsar distortion. Therefore, the indirect velocimetry method was suitable for low flux pulsars. Based on this, the direct velocimetry method and the indirect one were used to update the Kalman filter in turn. And the Kalman filter was used to provide high-accuracy navigation information. Simulation results show that the accuracy of X-ray pulsar navigation method based on direct/indirect hybrid velocimetry is 15% higher than that of the indirect pulsar navigation.
-
Key words:
- velocimetry/
- navigation/
- Kalman filters/
- TOA/
- X-ray pulsars
Highlights● The direct velocimetry method and the indirect velocimetry method are suitable for the Crab pulsar and low flux pulsars, respectively. ● The direct velocimetry method and the indirect velocimetry method are combined with each other to suppress the Doppler effects. ● A distorted profile matching estimation method for direct velocity and distance measurement is proposed. ● The accuracy of X-ray pulsar navigation method based on direct/indirect velocimetry is 15% higher than that of the indirect pulsar navigation. -
表 1X射线脉冲星仿真参数
Table 1Simulation parameters of X-ray pulsars
参数 数值 脉冲星 B0531+21 B1821-24 B1937+21 周期/s 0.033 40 0.003 05 0.001 56 光子流量密度/(ph·cm–2·s–1) 1.54 1.93×10–4 4.99×10–5 噪声流量密度/(ph·cm–2·s–1) 0.005 0.005 0.005 分辨率/μs 1 1 1 探测器有效面积/m2 1 1 1 观测时间/s 1 000 1 000 1 000 畸变轮廓数 81 — — 最大相位偏移量 81 701 701 表 2环火轨道初始六要素
Table 2Initial orbital elements of the Mars
轨道六要素 数值 半长轴/km 6 794 偏心率 0 轨道倾角/(°) 45 升交点赤经/(°) 0 近地点幅角/(°) 0 真近点角/(°) 0 表 3导航滤波器参数
Table 3Parameters of navigation filter
参数 数值 X射线敏感器数量 1 观测周期 1 000 s 初始轨道误差 $\begin{aligned} & \delta {\boldsymbol{X} }\left( 0 \right) = [6\;000\;{\text{m} },6\;000\;{\text{m} },6\;000\;{\text{m,} } \hfill \\& 20\;{\text{m/s} },20\;{\text{m/s} },15\;{\text{m/s} }{]^{\text{T} } } \hfill \end{aligned}$ 状态处理噪声协方差矩阵 ${\boldsymbol{Q } } = {\text{diag} }\left( {q_1^2,q_1^2,q_1^2,q_2^2,q_2^2,q_2^2} \right)$,
${q}_{1}=6\; \text{m},$$ {q_2} = 0.04\;{\text{m/s}} $测量噪声协方差矩阵 B0531+21:${\boldsymbol{R} } = \left[ \begin{gathered} \sigma _{\text{D} }^2 \;\;\;\;\;\;- 0.626{\sigma _{\text{D} } }{\sigma _{\text{v} } } \hfill \\ - 0.626{\sigma _{\text{D} } }{\sigma _{\text{v} } }\;\;\;\;\;\;\sigma _\text v^2 \hfill \\ \end{gathered} \right]$
$ {\sigma _{\text{D}}}{\text{ = }}89\;{\text{m}},{\sigma _{\text{v}}} = 0.1\;{\text{m/s}} $
B1821-24:${\boldsymbol{R}}=\left[{\sigma }_{\text{D} }^{2}\right],$$ {\sigma _{\text{D}}}{\text{ = 624}}\;{\text{m}} $
B1937+21:${\boldsymbol{ R}}=\left[{\sigma }_{\text{D} }^{2}\right],$$ {\sigma _{\text{D}}}{\text{ = 1}}\;{\text{463}}\;{\text{m}} $表 4畸变轮廓匹配估计
Table 4Matching estimation of distorted profiles
脉冲星 方法 定位精度/m 定速精度/(m·s–1) 计算时间/s B0531+21 畸变轮廓匹配 89 0.097 6 1.413 量子CS 102 0.098 8 1.621 B1821-24 互相关 624 / 0.001 4 B1937+21 互相关 1463 / 0.000 9 表 5两种导航方式的对比
Table 5Comparison of two navigation methods
导航方式 定位精度/m 定速精度/(m·s-1) 间接测速导航 6 871 2.823 直接/间接混合测速导航 5 722 2.411 表 6不同观测时间
Table 6Different observation times
导航方式 观测时间100 s 观测时间500 s 定位精度/
m定速精度/
(m·s–1)定位精度/
m定速精度/
(m·s–1)间接测速导航 12 718 4.436 7 549 2.991 直接/间接
混合测速导航12 482 4.370 7 164 2.771 -
[1] 房建成, 宁晓琳, 刘劲等. 航天器自主天文导航原理与方法(第2版)[M]. 北京: 国防工业出版社, 2017. [2] 焦荣,甘伟,肖志红,等. 利用增量相位和星光仰角增强XNAV[J]. 宇航学报,2019,40(6):666-672.JIAO R,GAN W,XIAO Z H,et al. Enhancement of XNAV utilizing incremental phase and star elevation[J]. Journal of Astronautics,2019,40(6):666-672. [3] RODIN A E,ORESHKO V V,POTAPOV V A,et al. Principles of pulsar space navigation[J]. Astronomy Reports,2020,64(6):499-525.doi:10.1134/S1063772920070057 [4] XU Q,FAN X H,ZHAO A G,et al. Pre-correction X-ray pulsar navigation algorithm based on asynchronous overlapping observation method[J]. Advances in Space Research,2021,67(1):583-596.doi:10.1016/j.asr.2020.10.019 [5] WANG Y S,WANG Y D,ZHENG W. On-orbit pulse phase estimation based on CE-adam algorithm[J]. Aerospace,2021,8(4):95.doi:10.3390/aerospace8040095 [6] LIU J,FANG J C,KANG Z W,et al. Novel algorithm for X-ray pulsar navigation against doppler effects[J]. IEEE Transactions on Aerospace and Electronic Systems,2015,51(1):228-241.doi:10.1109/TAES.2014.130463 [7] LIU J,FANG J C,NING X L,et al. Closed-loop EKF-based pulsar navigation for Mars explorer with Doppler effects[J]. Journal of Navigation,2014,67(5):776-790.doi:10.1017/S0373463314000216 [8] HUANG L W,LIANG B,ZHANG T. Pulse phase and doppler frequency estimation of X-ray pulsars under conditions of spacecraft and binary motion and its application in navigation[J]. Science China,Physics,Mechanics and Astronomy,2013,56(4):848-858. [9] XUE M F,PENG D L,SUN H F,et al. X-ray pulsar navigation based on two-stage estimation of Doppler frequency and phase delay[J]. Aerospace Science and Technology,2021,110:106470.doi:10.1016/j.ast.2020.106470 [10] LYNE A G,GRAHAM-SMITH F. Book review:pulsar astronomy.- 2nd ed. / Cambridge U Press,1998[J]. Journal of the British Astronomical Association,1999,109:145. [11] BURNS W R,CLARK B G. Pulsar search techniques[J]. Astronomy & Astrophysics,1969,2:280-287. [12] 李建勋,柯熙政. 基于循环平稳信号相干统计量的脉冲星周期估计新方法[J]. 物理学报,2010,59(11):8304-8310.doi:10.7498/aps.59.8304LI J X,KE X Z. Period estimation method for weak pulsars based on coherent statistic of cyclostationary signal[J]. Acta Physica Sinica,2010,59(11):8304-8310.doi:10.7498/aps.59.8304 [13] 李建勋,柯熙政,赵宝升. 一种脉冲星周期的时域估计新方法[J]. 物理学报,2012,61(6):537-543.LI J X,KE X Z,ZHAO B S. A new time-domain estimation method for period of pulsars[J]. Acta Physica Sinica,2012,61(6):537-543. [14] 周庆勇,姬剑锋,任红飞. 非等间隔计时数据的X射线脉冲星周期快速搜索算法[J]. 物理学报,2013,62(1):525-532.ZHOU Q Y,JI J F,REN H F. The quick search algorithm of pulsar period based on unevenly spaced timing data[J]. Acta Physica Sinica,2013,62(1):525-532. [15] ZHANG H,XU L P,XIE Q. Modeling and doppler measurement of X-ray pulsar[J]. Science China,2011,54(6):1068-1076. [16] 谢强,许录平,张华. 基于轮廓特征的X射线脉冲星信号多普勒估计[J]. 宇航学报,2012,33(9):1301-1307.XIE Q,XU L P,ZHANG H,et al. Doppler estimation of X-ray pulsar signals based on profile feature[J]. Journal of Astronautics,2012,33(9):1301-1307. [17] 张新源,帅平,黄良伟. 脉冲星导航轮廓折叠失真与周期估计算法[J]. 宇航学报,2015,36(9):1056-1060.ZHANG X Y,SHUAI P,HUANG L W. Profile folding distortion and period estimation for pulsar navigation[J]. Journal of Astronautics,2015,36(9):1056-1060. [18] LIU L,LIU L L,NING X L,et al. Fast butterfly epoch folding-based X-ray pulsar period estimation with a few distorted profiles[J]. IEEE Access,2020,8:4211-4219.doi:10.1109/ACCESS.2019.2962977 [19] LIU J,YANG Z H,KANG Z W,et al. Fast CS-based pulsar period estimation method without tentative epoch folding and its CRLB[J]. Acta Astronautica,2019,160:90-100.doi:10.1016/j.actaastro.2019.04.023 [20] 刘劲,韩雪侠,宁晓琳,等. 基于EMD-CS的脉冲星周期超快速估计[J]. 航空学报,2020,41(8):623486.LIU J,HAN X X,NING X L,et al. Ultra-fast estimation of pulsar period based on EMD-CS[J]. Acta Aeronauticaet Astronautica Sinica,2020,41(8):623486. [21] SHEIKH S I,PINES D J,RAY P S,et al. Spacecraft navigation using X-ray pulsars[J]. Journal of Guidance,Control,and Dynamics,2006,29(1):49-63. [22] WU D L,WU J,LIU J,et al. Quick X-ray pulsar positioning and velocimetry approach based on quantum CS[J]. Optik - International Journal for Light and Electron Optics,2021,241:166649.doi:10.1016/j.ijleo.2021.166649 [23] EMADZADEH A A,SPEYER J L. X-ray pulsar-based relative navigation using epoch folding[J]. IEEE Transactions on Aerospace & Electronic Systems,2011,47(4):2317-2328. [24] LIANG H, ZHAN Y F, YIN H L. New observation strategy for X-ray pulsar navigation using a single detector[J]. IET Radar Sonar and Navigation, 2016, 10(6): 1107-1111. [25] CRAVEN P H, COLLINS J T, et al. Spacecraft navigation X-ray pulsars[C]//7th International ESA Conference on Guidance. [S. l. ]: ESA, 2008. [26] DENNIS W W. The use of X-ray pulsars for aiding GPS satellite orbit determination[D]. Colorado Springs, CO: Air Force Institute of Technology, 2005. -