CN103913173A - 单x射线脉冲星导航观星序列选择方法 - Google Patents

Abstract

本发明提供了一种基于可观性分析的单X射线脉冲星导航观星序列选择方法，去除在航天器导航过程中会发生周期跃变的脉冲星、导航过程中会被天体遮挡的脉冲星、在太阳遮蔽角内的脉冲星后，选择使特定的指标函数达到最大值的脉冲星序列。本发明提出的方法原理简单，可操作性强，易推广和使用，大大提高了单X射线脉冲星导航的可靠性和稳定性，为其推广提供了条件。

Description

单X射线脉冲星导航观星序列选择方法

技术领域

本发明涉及导航方法领域，特别地，涉及一种用于X射线脉冲星导航的观星序列选择方法，适用于不同飞行任务的单X射线脉冲星导航任务。

背景技术

X射线脉冲星导航是一种新兴的天文导航方法。同传统的天文导航方法相比，X射线脉冲星导航可提供良好的外部参考时间频率基准。同以GPS（Global Positioning System）为例的卫星导航***相比，X射线脉冲星导航不局限于近地空间，可服务于深空探测器。同其他波段相比，X射线脉冲星在X射线频段的特征明显，可避免空间各种信号的干扰。在满足流量和时空分辨率的前提下可降低探测器的体积和功耗。此外，X射线脉冲星的数目较少，信号结构简单，降低了信号识别的难度。

单X射线脉冲星导航方法仅要求航天器搭载一颗探测器，可大幅降低X射线脉冲星导航***的体积和功耗。为了保证导航性能，单X射线脉冲星导航***要求航天器分时段观测不同的脉冲星。因此，观星序列的选择会直接影响单X射线脉冲星导航***的性能。为了选择最佳的导航脉冲星，就必须提出相应的指标函数。论文“毛悦,宋小勇,柴飞,脉冲星TOA测量误差及几何精度分析[J],测绘科学技术学报,2009,26(2):140-143.”和“褚永辉,王大轶,黄翔宇,脉冲星导航中最优脉冲星组合选取方法[J],中国空间科学技术,2011,5:64-69.”均提出了各自的选星指标函数。这两种指标函数借鉴了GPS导航中的GDOP概念，但仅考虑了脉冲星自身的几何分布和相应的脉冲到达时间测量精度的影响。申请号为201110095075.9的专利提出了一种综合考虑脉冲星星表位置误差、几何分布、脉冲星到达时间测量精度的选星策略。申请号为201310117103.1的专利提出了一种基于Fisher信息阵的脉冲星导航选星方法。这两份专利提出的方法也仅考虑了脉冲星自身特征对导航精度的影响。然而，单X射线脉冲星导航***必须通过滤波技术融合航天器轨道动力学信息和单X射线脉冲星测量信息才能实现导航定位。在该种情况下，若仅考虑脉冲星自身特征的影响而忽略航天器轨道的影响，难以准确选择最佳的导航脉冲星。

发明内容

本发明要解决的技术问题为：为了提高单X射线脉冲星导航***的性能，本专利综合考虑航天器轨道信息和脉冲星特征信息，提出一套可行的单X射线脉冲星导航观星序列选择方法。

为解决上述技术问题，本发明采用以下技术方案：

一种基于可观性分析的单X射线脉冲星导航观星序列选择方法，步骤为：

（1）根据天文观测数据，去除在航天器导航过程中会发生周期跃变的脉冲星；

（2）去除在导航过程中会被天体遮挡的脉冲星；

（3）分析航天器和太阳连线与航天器和脉冲星连线的夹角（即太阳遮蔽角），根据航天器的热控***要求，去除在太阳遮蔽角内的脉冲星；

（4）在通过了（1）-（3）步筛选的脉冲星中，利用下式计算每颗脉冲星的脉冲到达时间的估计精度：

${\mathrm{\σ}}_{\mathrm{TOA}}=\frac{\sqrt{\frac{1}{4}{T}_{50}^2+T_b^2}}{\sqrt{A\·\mathrm{\Δt}}}\·\frac{\sqrt{{\mathrm{\λ}}_n+{\mathrm{\λ}}_p}}{{\mathrm{\λ}}_p}---\left(1\right)$

其中，T₅₀是X射线脉冲信号半流量密度持续时间；T_b是探测器时间分辨率；λ_p和λ_n分别是脉冲信号和背景噪声的平均流量密度，A是探测器有效面积，Δt是观测持续时间；

（5）在通过了（1）-（3）步筛选的脉冲星中，选择使如下的指标函数达到最大值的脉冲星序列。

$J\left(\mathrm{pulsar}\right)=\mathrm{tr}\left(\underset{i=1}{\overset{6}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_{\mathrm{TOA},i}}^2}{\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)}^T{H_i}^T{H}_i\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)\right)---\left(2\right)$

在式(2)中，σ_TOA,i是第i颗观测的脉冲星的脉冲到达时间估计精度，Φ_k+1,k为航天器从k时刻到k+1时刻的状态转移矩阵，H_i为第i颗脉冲星的观测矩阵，tr(·)表示对矩阵的迹函数。

与现有技术相比，本发明的优点在于：

（1）本发明提出的单X射线脉冲星导航观星序列选择方法综合考虑了航天器轨道信息和脉冲星自身特征信息对导航性能的影响。由式（2）可知，航天器轨道的影响，可通过Φ_k+1,k的变化来反映。脉冲星的自身特征信息的影响可通过σ_TOA,i和H_i的变化来反映。

现有的选星方法均仅考虑了脉冲星自身特征的影响。然而，在论文“王奕迪.深空探测中的X射线脉冲星导航方法研究[D].国防科学技术大学,2011.”中，发明人充分论证了如下结论：在观测同一组脉冲星的情况下，航天器处于不同高度的轨道，单X射线脉冲星导航***的导航性能不相同。因此，对单X射线脉冲星导航***而言，若仅考虑脉冲星自身特征影响难以选择到合适的导航脉冲星。而本发明提出的观星序列选择方法，可综合考虑航天器轨道信息和脉冲星自身特征信息的影响，更为贴切地衡量脉冲星导航的性能，可为事先选择单X射线脉冲星导航的观星序列提供理论指导。在相同的条件下，提高滤波收敛速度，提高导航精度。

（2）本发明提出的方法可操作性强，易推广和使用，大大提高了单X射线脉冲星导航的可靠性和稳定性，为其推广提供了条件。

附图说明

图1是本发明具体应用实例的流程示意图。

图2是本发明具体应用实例中天体对X射线脉冲星的可用性影响示意图。

图3是本发明具体应用实例中太阳规避角对X射线脉冲星的可用性影响示意图。

图4是本发明具体应用实例中的选星与不选星的导航结果对比图。

其中，1、地球，2、阴影区，3、脉冲星，4、太阳，5、地球卫星，6、轨道。

具体实施方式

以下将结合附图和具体实施例对本发明作进一步详细说明。

在具体应用实例中，所采用的航天器为一种地球卫星。

如图1所示，本发明选星方法的具体步骤如下：

（1）根据天文观测数据，去除在航天器导航过程中会发生周期跃变的脉冲星。

虽然X射线脉冲星自转周期的长期稳定性可媲美当前的原子钟。但一些“年轻”脉冲星的自转周期存在不规则的变化。当脉冲星发生周期跃变，该脉冲性提供的测量信息不可靠。因此，首先需要基于天文观测数据预报脉冲星信号的稳定性，去除航天器运行期间导航星数据库中可能发生周期跃变的X射线脉冲星。

例如，天文台存有大量Ｘ射线脉冲星观测数据且定期发布处理结果，可以通过天文台官方网站查询导航Ｘ射线脉冲星数据库中正在或者在短期内可能发生周期跃变的X射线脉冲星，在导航过程中对这些脉冲星不予使用。

（2）去除在导航过程中会被天体遮挡的脉冲星。

尽管X射线脉冲星距离太阳系非常遥远，但是当航天器进入某天体的阴影区内时，即任何经过航天器与X射线脉冲星之间的天体都会阻挡航天器探测器对X射线脉冲星的可见性。要根据航天器的飞行任务轨道预报脉冲星的可用性，即考虑星体遮挡的影响，根据航天器的飞行任务轨道预报脉冲星的可用性，去除不能使用的脉冲星。

如图2所示，当卫星进入地球阴影区内时，地球会阻挡Ｘ射线探测器接收X射线脉冲星信号。图中角度α为航天器相对于地球的位置矢量r和X射线脉冲星的方向n之间的夹角，地球半径为R，当α₁＜α＜α₂时，航天器位于地球的阴影中，此时该脉冲星不可用，航天器位于遮挡天体阴影处时要满足：

$\mathrm{\π}-\arccos \left(\frac{\sqrt{{\left|r\right|}^2-R^2}}{\left|r\right|}\right)\≤\arccos (n\·r)\≤\mathrm{\π}+\arccos \left(\frac{\sqrt{{\left|r\right|}^2-R^2}}{\left|r\right|}\right)---\left(3\right)$

（3）分析航天器和太阳连线与航天器和脉冲星连线的夹角，根据航天器的热控***要求，去除在太阳遮蔽角内的脉冲星。

根据航天器的热控需求，会规定不同给的太阳规避角。航天器不会接收到在太阳规避角中的脉冲星信号。图2给出了太阳规避角和脉冲星、太阳、航天器之间的关系。其中，β为太阳和航天器的连线与脉冲星和航天器连线的夹角，可由式(4)计算。若β小于规定的太阳规避角，则该颗脉冲星不可用。

β＝arccos(n·R_Sun/sc)          (4)

在式(4)中，R_Sun/sc为太阳相对于航天器的方向矢量。

（4）在通过了（1）-（3）步筛选的脉冲星中，计算每颗脉冲星的脉冲到达时间估计精度。

探测器测量到的光子到达事件是服从泊松分布的，依据泊松分布均值等于方差的性质，可以使用下列模型来计算脉冲到达时间的估计精度：

${\mathrm{\σ}}_{\mathrm{TOA}}=\frac{\sqrt{{\mathrm{HWHM}}^2+T_b^2}}{S/\sqrt{S+B}}---\left(5\right)$

其中，HWHM是脉冲信号半流量密度持续时间T₅₀的一半；S是观测期间探测到的来自脉冲星的光子数；B是探测到的背景噪声的光子数；T_b是探测器时间分辨率。

S＝Aλ_pΔt      (6)

B＝Aλ_nΔt      (7)

其中，λ_p和λ_n分别是脉冲信号和背景噪声的平均流量密度，A是探测器有效面积，Δt是观测持续时间。

将式(6)和(7)代入式(5)，可得到X射线脉冲星的脉冲到达时间估计精度为：

${\mathrm{\σ}}_{\mathrm{TOA}}=\frac{\sqrt{\frac{1}{4}{T}_{50}^2+T_b^2}}{\sqrt{A\·\mathrm{\Δt}}}\·\frac{\sqrt{{\mathrm{\λ}}_n+{\mathrm{\λ}}_p}}{{\mathrm{\λ}}_p}---\left(8\right)$

其中，T₅₀是X射线脉冲信号半流量密度持续时间；T_b是探测器时间分辨率；λ_p和λ_n分别是脉冲信号和背景噪声的平均流量密度，A是探测器有效面积，Δt是观测持续时间。

（5）选择使指标函数达到最大值的脉冲星序列。

在地心惯性系J2000.0中，建立近地卫星的轨道动力学方程为

$\left[\begin{array}{c}\stackrel{\·}{r}\\ \stackrel{\·}{v}\end{array}\right]=\left[\begin{array}{c}v\\ a\end{array}\right]+\begin{array}{c}\left[\begin{array}{c}{w}_r\\ w_v\end{array}\right]\end{array}---\left(9\right)$

其中，x=[r^T,v^T]^T是航天器的状态矢量，w=[w_r ^Tw_v ^T]^T为动力学模型噪声，可建模为零均值高斯白噪声。该噪声的方差阵为Q。a＝a_TB+a_NS+a_T+a_H.O.T是航天器受到的加速度，包含以下几项。

1）a_TB＝-μ_Er/|r|³是航天器受到的地球二体引力加速度，其中μ_E是地球的引力常数。

2）是地球的非球形引力摄动。U_NSE可以表示为

$U_{\mathrm{NSE}}=-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_n{P}_n\sin \mathrm{\φ}-\frac{{\mathrm{\μ}}_E}{\left|r\right|}\underset{n=2}{\overset{\∞}{\mathrm{\Σ}}}\underset{m=1}{\overset{n}{\mathrm{\Σ}}}{\left(\frac{R_e}{\left|r\right|}\right)}^n{J}_{n,m}{P}_{\mathrm{nm}}\sin \mathrm{\φ}\cos m(\mathrm{\λ}-{\mathrm{\λ}}_{n,m})---\left(10\right)$

其中，R_e是地球的半径，φ和λ分别是经度和纬度，P_n和P_nm是勒让德多项式，J_n是带谐项系数，λ_n,m是田谐项，J_n,m为田谐项系数。

3）是三体摄动加速度。μ_i是第i个天体的引力常数，r_i是第i个天体相对于地球的位置矢量。

4）a_H.O.T是影响航天器加速度的高阶项。相对于已建模的摄动加速度，这些高阶项的影响可以忽略。

式(9)可简记为

$\stackrel{\·}{x}=f\left(x\right)+w---\left(11\right)$

脉冲星导航的观测方程为

Z＝h_P ^j(x)+υ        (12)

其中，V_p为观测噪声，h_P ^j(x)为观测第j个脉冲星的观测方程，其表达式为

$\begin{array}{c}{h_P}^j\left(x\right)={\tilde{n}}^j\·(r+{\tilde{r}}_E)+\frac{1}{{{2D}_0}^j}\{{[{\tilde{n}}^j\·(r+{\tilde{r}}_E)]}^2-{\left|\right|r+{\tilde{r}}_E\left|\right|}^2+2({\tilde{n}}^j\·b)[{\tilde{n}}^j\·(r+{\tilde{r}}_E)]-2[b\·(r+{\tilde{r}}_E)\left]\right\}\\ +\frac{{2\mathrm{\μ}}_S}{c^2}\ln |\frac{{\tilde{n}}^j\·(r+{\tilde{r}}_E)+\left|\right|r+{\tilde{r}}_E\left|\right|}{{\tilde{n}}^j\·b+\left|\right|b\left|\right|}+1|\end{array}---\left(13\right)$

其中，为第j个脉冲星方向矢量的测量值，为由星历预报的地球位置，D₀ ^j为第j个脉冲星相距太阳质心的距离，b为太阳质心相对于太阳系质心的位置，μ_S为太阳的引力常数，c为光速。

式(11)和(13)即为X射线脉冲星导航的***方程。将该***方程在预估位置x^*处线性化，可得

$\left\{\begin{array}{c}{\mathrm{\Δx}}_{k+1}={\mathrm{\Φ}}_{k+1,k}{\mathrm{\Δx}}_k+W_k\\ {\mathrm{\ΔZ}}_{k+1}=H_{k+1}{\mathrm{\Δx}}_{k+1}+V_k\end{array}\right.---\left(14\right)$

其中，

Φ_k+1,k＝e^FΔt(15)

$F_k=\frac{\∂f\left(x\right)}{\∂x}{|}_{x=x*},H_{k+1}=\frac{{\∂h}_p^j\left(x\right)}{\∂x}{|}_{x=x*}---\left(16\right)$

建立观测性方程

$\left[\begin{array}{c}\mathrm{\Δ}{Z}_k\\ {\mathrm{\ΔZ}}_{k+1}\\ \mathrm{\Δ}{Z}_{k+2}\\ {\mathrm{\ΔZ}}_{k+3}\\ \mathrm{\Δ}{Z}_{k+4}\\ {\mathrm{\ΔZ}}_{k+5}\end{array}\right]=M_k{\mathrm{\Δx}}_k+\left[\begin{array}{c}{V}_k\\ V_{k+1}\\ V_{k+2}\\ V_{k+3}\\ V_{k+4}\\ V_{k+5}\end{array}\right]---\left(17\right)$

其中，

$M_k=\left[\begin{array}{c}{H}_k\\ H_{k+1}{\mathrm{\Φ}}_{k+1,k}\\ H_{k+2}{\mathrm{\Φ}}_{k+1,k}^2\\ H_{k+3}{\mathrm{\Φ}}_{k+1,k}^3\\ H_{k+4}{\mathrm{\Φ}}_{k+1,k}^4\\ H_{k+5}{\mathrm{\Φ}}_{k+1,k}^5\end{array}\right]---\left(18\right)$

式(17)可通过最小二乘求解。该解的方差阵为

$P_k={\left({M_k}^T{R}_k^{-1}{M}_k\right)}^{-1}---\left(19\right)$

其中， 为第i个观测的脉冲星的脉冲到达时间估计精度。

式(19)展开，可写成

$P_k={\left(\underset{i=1}{\overset{6}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_{\mathrm{TOA},i}}^2}{\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)}^T{H_i}^T{H}_i\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)\right)}^{-1}---\left(20\right)$

由方差阵的意义可知，P_k减小，对应着脉冲星导航精度的提高。为此，可取P_k ^-1的迹函数作为选星的指标函数，即

$J\left(\mathrm{pulsar}\right)=\mathrm{tr}\left(\underset{i=1}{\overset{6}{\mathrm{\Σ}}}\frac{1}{{{\mathrm{\σ}}_{\mathrm{TOA},i}}^2}{\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)}^T{H_i}^T{H}_i\left({\mathrm{\Φ}}_{k+1,k}^{i-1}\right)\right)---\left(21\right)$

其中，tr(·)表示对矩阵的迹函数。

根据如式(21)所示的指标函数，若选择的脉冲星观测序列可使得式(21)达到最大值，则该组观星序列是最优的。

为了表明前述选星策略的可行性，现以高轨卫星进行仿真验证。仿真时间为2001年3月11日至2001年3月12日，脉冲星观测周期为300s，脉冲星的星表误差为0.1毫角秒。初始导航误差为（1km，1km，1km）和（1m/s，1m/s，1m/s）。利用UKF为滤波器。导航所用的脉冲星从microcosm公司提出的八颗脉冲星中选取。这八颗脉冲星分别为：PSRB1937+21、PSRJ0218+4232、PSRB0540-69、PSRB1509-58、PSRB1821-24、PSRJ1814-338、PSRB0531+21、SAXJ1808-3658。

第一步，基于英国JodrellBank天文台的Lonvell望远镜对700多颗脉冲星进行的监测数据，在2001年3月12日，脉冲星PSRB1821-24会发生周期跃变。因此，剔除脉冲星PSRB1821-24。

第二步，结合卫星轨道信息和地球、月球轨道信息，按照式（3）计算遮挡天体阴影区范围。计算表明，SAXJ1808-3658、PSRJ0218+4232被地球遮挡，无法保证全导航区间可见。因此，剔除脉冲星SAXJ1808-3658和PSRJ0218+4232。

第三步，根据一般卫星的热控***性能，设定太阳规避角为15°。按照式（4）计算卫星到脉冲星与卫星到太阳的夹角。计算表明，脉冲星PSRB1937+21对应的夹角全程小于15°。因此，剔除脉冲星PSRB1937+21。

第四步，根据脉冲星观测周期，计算脉冲星PSRB0540-69、PSRJ1814-338、PSRB0531+21、PSRB1509-58的脉冲到达时间估计精度。

第五步，根据式（21）计算指标函数，结果表明，顺序观测PSRB0540-69、PSRJ1814-338、PSRB0531+21得到的指标函数最大。因此，在这种情况下，PSRB0540-69，PSRJ1814-338，PSRB0531+21为最佳导航脉冲星组合。

图4给出了选星与不选星所得到的导航结果差异。从图4可以看出，和利用随意选取的脉冲星进行导航得到的导航结果相比，利用根据选星策略选出的脉冲星进行导航，在相同的仿真条件下，可以提高滤波收敛速度，提高导航精度。两种情况的比较也充分显示了本发明所提出的选星策略的可行性以及在脉冲星导航中采用选星策略的必要性。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应该提出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种单X射线脉冲星导航观星序列选择方法，其特征在于，包括以下步骤： 

A、根据天文观测数据，去除在航天器导航过程中会发生周期跃变的脉冲星； 

B、去除在导航过程中会被天体遮挡的脉冲星； 

C、分析航天器和太阳连线与航天器和脉冲星连线的夹角，根据航天器的热控***要求，去除在太阳遮蔽角内的脉冲星； 

D、在通过了A-C步筛选的脉冲星中，利用下式计算每颗脉冲星的脉冲到达时间的估计精度： 

其中，T₅₀是X射线脉冲信号半流量密度持续时间；T_b是探测器时间分辨率；λ_p和λ_n分别是脉冲信号和背景噪声的平均流量密度，A是探测器有效面积，Δt是观测持续时间； 

E、在通过了A-C步筛选的脉冲星中，选择使如下的指标函数达到最大值的脉冲星序列； 

在式Ⅱ中，σ_TOA,i是第i颗观测的脉冲星的脉冲到达时间估计精度，Φ_k+1,k为航天器从k时刻到k+1时刻的状态转移矩阵，H_i为第i颗脉冲星的观测矩阵，tr(·)表示对矩阵的迹函数。 

2.根据权利要求1所述的一种基于可观性分析的单X射线脉冲星导航观星序列设计方法，其特征在于，步骤B中，航天器位于遮挡天体阴影处时满足： 

。

3.根据权利要求1所述的一种基于可观性分析的单X射线脉冲星导航观星序列设计方法，其特征在于，步骤C中， 

太阳和航天器的连线与脉冲星和航天器连线的夹角β小于规定的太阳规避角时，该颗脉冲星不可用； 

β＝arccos(n·R_Sun/sc)           Ⅳ 。

4.根据权利要求1所述的一种基于可观性分析的单X射线脉冲星导航观星序列设计方法，其特征在于，步骤A中， 

去除导航Ｘ射线脉冲星数据库中正在或者在短期内可能发生周期跃变的X射线脉冲星。