CN103674032B - 融合脉冲星辐射矢量和计时观测的卫星自主导航***及方法 - Google Patents

Abstract

本发明公开了一种融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位***及方法。***利用时辨探测器配合准直器获取脉冲到达时间和精确获取脉冲星辐射方向矢量，并利用红外或紫外地平敏感器获取地心方向矢量；利用脉冲星脉冲到达时间、辐射方向矢量和地心方向矢量，导航计算机分别完成地心坐标系下脉冲星计时观测导航算法、基于角距测量的脉冲星定位算法和基于角距测量的几何定向算法，并通过几种算法的融合实时获取航天器位置、速度、姿态以及在轨运动方向，完成轨道卫星自主导航，输出高精度的位置、姿态和速度信息。本发明适用于地球或其它行星的轨道卫星或星座的自主运控。

Description

融合脉冲星辐射矢量和计时观测的卫星自主导航***及方法

技术领域

本发明属于导航技术领域，特别涉及融合X射线脉冲星辐射矢量和计时观测的卫星自主导航方法，用于为地球或其它近天体的轨道卫星，星际航行航天器，或者无稠密大气星体表面巡游车提供高精度自主导航定位服务。

背景技术

基于X射线脉冲星的导航定位是一种利用X射线脉冲星辐射脉冲到达时间测量作为信息输入的天文导航方法，能为近地、深空及星际空间航天器提供位置、速度、时间、姿态等丰富的导航信息。它通过测量同一脉冲到达航天器和到达SSB的时间差来计算航天器在该脉冲星方向相对于SSB的位移，再结合多颗脉冲星方向上的位移即可确定航天器在SSB坐标系下航天器位置。其中，时间差测量通常使用星载探测器接收X射线光子，在星载钟支持下经过周期累积得到累积脉冲轮廓，再与标准脉冲轮廓比对相位得到。

基于脉冲星的导航技术最早由美国喷气推进实验室的德恩斯博士于1974年提出；1981年，美国通信***研究所的切斯特和巴特曼提出利用脉冲星X射线源为航天器导航的构想；1999年，搭载USA试验设备的美国空军”先进研究与全球观测卫星”被发射进入预定轨道，USA开展空间科学试验研究中包括了一项有关X射线导航的主题研究；2004年8月，美国航空航天局和海军天文台等多家单位着手拟定和启动脉冲星导航的研究计划，同时X射线脉冲星导航已纳入国防部长期发展战略规划纲要，并逐年增加项目研究经费，持续开展脉冲星导航的理论方法研究、关键技术攻关和原理样机研制等方面的研究工作。我国在空间科学十一五规划中也包括了发射用于实现1-250keV宽频段X射线成像巡天的天文卫星HXMT，其中的分***-低能X射线，其频段为1.0-15keV，该X射线望远镜用于探索脉冲星导航。

脉冲星导航有能力具有多种优越的导航特性包括：(1)自主导航能力增强(2)从近地轨道、深空至星际空间飞行的持续导航能力(3)提供良好的外部时间频率基准(4)航天器编队飞行星座导航或星间跟踪(5)有利于减少维护成本和提高抗干扰能力。但是，脉冲星信号极其微弱，在目前的技术和认知水平下，很难达到理想的导航定位精度，特别是对于轨道卫星定轨应用，与传统的基于地基测控的定轨方法相比，除了在自主导航和抗干扰能力上具有明显优势外，其定轨精度还有较大差距。如果在实现脉冲星自主导航的同时，能进一步提高其定位精度，不仅可以拓展脉冲星应用空间，还十分有利于挖掘更多潜在的应用。而从目前脉冲星导航的原理来看，其仅仅利用脉冲星辐射脉冲稳定的计时特性测量信号光经过航天器到达太阳系质心的时间，再通过换算得到距离。基于这一原理，导航精度受诸多因素的限制，如相对论修正精度、脉冲星计时模型精度、位相测量精度等等。传统增强脉冲星导航精度可以通过优化信息获取，如精化相位预测模型、改进信号的时间测量等，但这类方法只能使定位精度接近其无偏估计的下界，并不增加可用信息量，因而也不能从本质上提升脉冲星导航的性能。因此如何在有限的观测技术水平下有效增强导航性能对X射线脉冲星导航理论研究和工程应用都具有十分重要的意义。

发明内容

针对上述问题，本发明的主要目的在于：针对传统X射线脉冲星导航方法中，仅仅利用脉冲星辐射脉冲的计时观测作为导航信息输入，造成脉冲星导航能力利用不足、轨道卫星定位精度需要进一步提高的问题，提供一种融合脉冲星辐射矢量和计时观测的卫星自主导航方法，以达到充分利用脉冲星导航能力，提高脉冲星导航精度，实现轨道卫星自主高精度轨道确定、时间同步和姿态测量。

本发明的另一个目的还在于能够为编队卫星或者卫星星座提供定向手段，以辅助解决星座整体旋转和地球自传误差积累。

为实现上述目的，本发明的实现方法和技术解决方案为：融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，其特点在于包括X射线探测器、准直器、X射线脉冲星定时模型和参数数据库、太阳系行星参数库，星载原子钟、星载计算机，红外(或紫外)地平敏感、导航算法库、星载万向支架等；星载计算机根据当前记录的航天器姿态和X射线脉冲星参数数据库，控制星载万向架支撑星载X射线探测器指向目标脉冲星；X射线探测器配合准直器，通过测量到达探测器的X射线脉冲星的光子流量强度计算脉冲星指向与脉冲星辐射矢量的偏差，该偏差反馈到星载计算机用于增强姿态控制，精确获取脉冲星辐射方向矢量；在探测器准确指向目标脉冲星的同时，X射线探测器记录X射线光子到达时间，在星载脉冲星参数库和太阳系行星参数库支持下将光子到达时间改正转换到太阳系质心；转换后的光子按照周期统计其在采样时间分辨单位内的重复频率直方图，进而形成累积轮廓；该累积轮廓与脉冲星标准平均轮廓比相，获取累计轮廓位相，该位相再与星载原子钟保持的脉冲星相位演化模型做差，所得相差转换成距离差，送至导航计算机中的基于地心坐标系的脉冲星计时观测导航算法；同时地平敏感仪敏感地平得到地平角距，并计算地心方向矢量，该角距与上述脉冲星辐射方向矢量共同送至星载计算机中基于角距测量的脉冲星定位算法，同时该地心方向矢量与上述脉冲星辐射方向矢量共同送至基于角距测量的几何定向算法；星载计算机调用自主导航算法库中的融合导航滤波算法对基于脉冲星计时观测的导航算法、基于角距测量的脉冲星定位算法和基于角距测量的几何定向算法融合，输出高精度的自主导航参数信息，包括位置、姿态和速度信息。

所述的X射线探测器为时间分辨探测器，所述的准直器安装于探测器前方，准直器由一系列平行的空心柱体紧密排列，只允许与柱体平行的光线通过，到达探测器，依次达到背景噪声滤波、辐射矢量识别的目的。

所述的重复频率直方图的生成方法为：假设为一段X射线脉冲星观测信号，它由N个采样构成，采样间隔为δ_t，起始时刻为t₀，周期为T_P，Nδt＝mT_P，m为观测信号中的周期数，z_i(i＝1，...，N)表示每个采样间隔内光子计数，写成向量形式为z＝[z₁，z₂，...，z_N]^T，令为了区别观测信号，用x＝[x₁，x₂，...，x_m]^T表示累积后的信号，令p_i＝|x_i|/||z||，重复频率直方图的生成方法依次调整变量i，使x的熵：

$H\left(x\right)=-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{p}_i\log p_i$

最小。

所述的脉冲星计时观测导航算法的实现方法为：星载计算机利用太阳系行星参数数据库和脉冲星计时模型和参数数据库计算地球坐标系中心相对于太阳系质心的位置和速度，并通过光行时转换和时间转换将建立在太阳系质心和太阳系质心力学时下的脉冲星相位演化模型转换到地球坐标系中心和地球质心力学时下；同时星载计算机在星载原子钟支持下，维持地心坐标系下的脉冲星相位演化模型；星载X射线探测器记录X射线脉冲星光子到达时间，整合出累积脉冲轮廓，该轮廓与所述的地心坐标系下的脉冲星相位演化模型比相获得相差，进而利用相差计算距离误差，该距离误差送入导航算法库。

所述的地平角距为地心坐标系下卫星到地球两侧表面切线方向矢量的夹角。

所述的地心方向矢量为地球坐标系下卫星中心到地球中心连线方向矢量，计算方法为地平角距的中心线方向。

所述的脉冲星辐射矢量的识别方法为：以万向支架支撑点为中心，在一个圆球面上以往返形式等幅度步进探测器指向，当探测器接收到的目标脉冲星辐射流量强度最大时即认为探测器指向与脉冲星辐射方向矢量平行。

所述的脉冲星辐射流量强度表示方法为对探测器记录的光子序列做FFT变换，统计目标脉冲星信号在其特征频率点的平均功率，用于标定流量强度。

所述的基于角距测量的脉冲星定位算法的实现方法为：星载计算机同时接收X射线探测器输出的脉冲星矢量信息和地平敏感仪输出的地心方向矢量以及地平仪输出的地平角距，首先建立地平角距与卫星位置矢量的函数关系式，命名为关系式A，再建立地心矢量方向与脉冲星辐射方向矢量夹角与卫星位置的函数关系式，命名为关系式B；利用关系式A和关系式B建立观测方程送入导航滤波算法。

所述的关系式A为 $\left|{\stackrel{\→}{r}}_{U/S}\right|\sin \frac{\mathrm{\α}}{2}=R_e.$

所述的关系式B为

$\mathrm{\β}=\arccos (-\frac{{\stackrel{\→}{r}}_{U/S}\·{\hat{r}}_s}{\left|{\stackrel{\→}{r}}_{U/S}\right|})$

所述的观测方程为：设卫星对地心矢量方向为它与脉冲星视线方向的矢量的夹角为β，地平仪观测的两侧地平方向夹角为α，地球半径为R_e。

$\mathrm{\α}=\widetilde{\mathrm{\α}}+{\left[\frac{\∂\mathrm{\α}}{\∂{\stackrel{\→}{r}}_{U/S}}\right]}^T\mathrm{\δX}+n_s=\widetilde{\mathrm{\α}}-\frac{R_e}{\sqrt{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2-R_e^2}}\frac{{\stackrel{\→}{r}}_e^T}{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2}\mathrm{\δX}+n_s=\widetilde{\mathrm{\α}}-F_s\mathrm{\δX}+n_s---\left(1\right)$

其中n_s包括泰勒展开的高阶项和噪声，简单起见，假设噪声服从高斯分布近似有E[n_s(t)]＝0， $E\left[n_s\left(t\right)n_s^T\left(t\right)\right]={\mathrm{\δ}}_s^2.$

$\mathrm{\β}=\widetilde{\mathrm{\β}}+{\left[\frac{\∂\mathrm{\β}}{\∂{\stackrel{\→}{r}}_{U/S}}\right]}^T\mathrm{\δX}+n_h=\widetilde{\mathrm{\β}}+\frac{{\stackrel{\→}{r}}_s^T}{\sqrt{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2-{\stackrel{\→}{r}}_{U/S}\·r_s}}\mathrm{\δX}+n_h=\widetilde{\mathrm{\α}}-F_h\mathrm{\δX}+n_h---\left(2\right)$

其中，n_h包括泰勒展开的高阶项和地心矢量观测噪声，E[n_h]＝0， $E\left[n_h\left(t\right){n_h}^T\left(t\right)\right]=$ ${\mathrm{\δ}}_h^2.$

所述的基于角距测量的几何定向算法实现方法为：设轨道运动方向用脉冲星辐射方向矢量为n_p，t₀时刻的地心方向矢量分别为和则可表示为

${\widehat{\mathrm{\θ}}}_{t_0}=\arccos n_p\·n_{{\mathrm{rt}}_1}-{\arccos n}_p\·n_{{\mathrm{rt}}_0}---\left(3\right)$

所述的导航滤波算法观测矩阵由，(1)式、(2)式和(22)式所示的测量方程组成，设组合量测向量表示为Y＝[Δτ，Δα，Δβ]^T，其中 $\mathrm{\Δ\τ}=\mathrm{\Δ}{t}_{\mathrm{TOA},N}-\mathrm{\Δ}{\tilde{t}}_{\mathrm{TOA},N},$ $\mathrm{\Δ\α}=\mathrm{\α}-\widetilde{\mathrm{\α}},$ $\mathrm{\Δ\β}=\mathrm{\β}-\widetilde{\mathrm{\β}},$ 此时观测方程为：

$Y=\mathrm{H\σX}={[F_s^T,F_h^T,G^T]}^T\mathrm{\δX}+\mathrm{\η}\left(t\right)---\left(4\right)$

其中H为观测矩阵，η(t)为零均值高斯观测噪声且η(t)＝[n_s(t)，n_h(t)，w_N(t)]^T，方差阵为：

$R=E\left({\mathrm{\η\η}}^T\right)=\left(\begin{array}{ccc}{n}_s{\left(t\right)}^2& & \\ & n_h{\left(t\right)}^2& \\ & & w_N{\left(t\right)}^2\end{array}\right)$

其w_N(t)²为脉冲星时间测量误差转换得到的位置量测误差方差，i＝(1，...，N)，表示探测器对第i颗脉冲星时间测量误差，W为脉冲半宽度， 其中F_x为X射线脉冲星辐射光子流，A表示探测器面积，p_f为一个脉冲周期内的脉冲辐射流量与平均辐射流量比，Δt_obs表示观测时间，B_x表示X射线背景辐射噪声，d＝W/P，P为脉冲周期。为了使提高TOA测量的精度，通常要使用累积的方法得到高信噪比的脉冲轮廓，此时Δt_obs反映了脉冲叠加段的长度。

所述的导航滤波算法的特征在于：在(2)式和(1)式的基础上，多星观测矩阵是单星观测变量的有序组合，此时矢量观测矩阵表示为

Y₁＝H₁δX+η＝[F_h，1，F_h，2，...，F_h，n]^TδX+η₁(5)

其中，Y₁＝[Δβ₁，Δβ₂，...，Δβ_n]^T，噪声均值，

同理到达时间观测矩阵表示为：

$Y_2=H_2\mathrm{\δX}+{\mathrm{\η}}_2={[G_1^T,G_2^T,...,G_N^T]}^T\mathrm{\δX}+{\mathrm{\η}}_2---\left(6\right)$

其中，Y₂＝[Δτ₂，Δτ₂，...，Δτ_n]^T，噪声均值，E(η₂)＝0，

所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，其特征在于：基于计时观测的导航或者基于角距观测均可以独立完成近地轨道卫星导航定位，二者也可以共同用于导航以增加导航精度，适合于近地轨道卫星或者卫星星座高精度自主导航。

所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法其特征在于，导航***的参考坐标中心均在地心。

本发明与现有基于地基测控和基于星间链路的卫星和星座导航技术相比，优点有：(1)能够为卫星或星座提供长时间的高精度自主导航能力，有利于缓解地面测控站的压力；(2)由于采用的地基时空基准，可以克服星间链路自主导航方法中星座整体旋转问题；(3)***被动接收天文X射线信号，在地面上无法进行干扰。

此外，与传统的单纯利用计时观测的脉冲星导航方法相比，本发明增加了角距观测，因此对脉冲星导航信息的利用更为充分，可以显著提高脉冲星导航精度。

附图说明

图1为本发明融合脉冲星辐射矢量和计时观测的卫星自主导航***组成图；

图2为本发明的在地球质心参考坐标系中的脉冲星、太阳质心、导航卫星、地平角距之间的几何关系以及基于计时观测导航算法和基于角距测量的导航算法原理图；

图3为本发明的在地球质心参考坐标系中地球中心与脉冲星辐射方向矢量之间的几何关系和基于角距测量的几何定向算法原理图；

图4为本发明的融合计时观测和角距测量的脉冲星导航融合算法框图。

具体实施方式

如图1所示，本发明的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位***，包括X射线探测器3、准直器2、X射线脉冲星定时模型和特征参数数据库7、太阳系行星参数库6，星载原子钟17、星载计算机，红外(或紫外)地平敏感11、导航算法库8、和卫星自主控制平台等；来自X射线脉冲星的辐射光子流1经过准直器2做选通处理，当准直器2的指向与来自X射线脉冲星的辐射光子流方向矢量一致或在选通角范围内时，X射线光子流1到达X射线时变探测器3，否则X射线脉冲星辐射光子流无法通过准直器到达X射线时变探测器3；X射线脉冲星时变探测器3接收X射线脉冲星辐射的光子，送入记录光子到达时间模块9，记录光子到达时间，统计光子到达流量强度，当光子流量强度达到最大时认为准直器与X射线脉冲星辐射矢量方向平行，即准直器此时所指方向为脉冲星在天球坐标系中的位置矢量；该矢量经过坐标转换送至辐射矢量识别模块10；同时X射线探测器3探测的光子信号送入记录光子到达时间模块9记录的X射线脉冲星光子到达时间，该光子到达时间首先通过光行时修正和时间尺度修正转换到太阳系质心，得到太阳系质心处的光子到达时间，利用太阳系行星参数数据库6中的行星参数将太阳系质心处的光子到达时间转换到地心坐标系中心的光子到达时间，然后再利用脉冲星计时模型和特征参数数据库7中的脉冲星特征参数，进行光子轮廓整合，得到累积脉冲星脉冲轮廓；该累积脉冲星脉冲轮廓与星载计算机维持的地心坐标系下的X射线脉冲星辐射相位演化模型进行比对，得到相差；该相差送入基于计时观测的导航模块12，该模块记录多颗脉冲星相差并将相差转换成位置误差，从导航算法库模块8中调用基于计时观测的脉冲星导航算法完成基于计时观测的脉冲星导航；在进行计时观测导航的同时，辐射矢量识别模块10对转换后的脉冲星辐射方向矢量进行辨识，辨识过程不仅要识别该辐射矢量所属的脉冲星，还要根据准直器模型和流量强度的变化规律修正准直器指向与辐射矢量方向的偏差，得到高精度的辐射矢量；该矢量送至基于角距观测导航模块13和基于角距观测的几何定向模块14；此外，地平仪同时敏感地球表面切线方向矢量，送至地平角距测量模块11，地平角距测量模块根据地球表面切向方向矢量计算地心角距，得到地心方向矢量；该地心方向矢量和地球表面切线方向矢量同时送入基于角距观测的导航模块13和基于角距观测的几何定向模块14；基于角距观测的导航模块13利用地球表面切线方向矢量和地心方向矢量与多颗脉冲星辐射方向矢量的夹角建立观测方程，并从导航算法库调用基于角距观测的导航算法完成基于角距观测的导航；基于角距观测的几何定向模块14利用卫星在轨道中的不同位置时，地心方向矢量与脉冲星辐射方向矢量的夹角的变换情况进行卫星运动方向的定向；基于计时观测的导航模块12、基于角距观测的模块13和基于角距观测的几何定向模块14将各自导航状态和参数送至融合导航算法15，输出高精度的导航信息送至卫星自主控制平台。

在图1中，基于计时观测导航模块12中的轮廓累积采用统计重复频率直方图的方法生成，生成方法为：假设x为一段X射线脉冲星观测信号，它由N个采样构成，采样间隔为δ_t，起始时刻为t₀，周期为T_P，Nδt＝mT_P，m为观测信号中的周期数，z_i(i＝1，...，N)表示每个采样间隔内光子计数，写成向量形式为z＝[z₁，z₂，...，z_N]^T，令为了区别观测信号，用x＝[x₁，x₂，...，x_m]^T表示累积后的信号，令p_i＝|x_i|/||z||，重复频率直方图的生成方法依次调整变量i，使x的熵：

$H\left(x\right)=-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{p}_i\log p_i---\left(7\right)$

最小。

在图1中基于计时观测的导航模块12，其实现方法为：基于星载计时原子钟准确预测脉冲星信号到达选定惯性参考系原点的到达时间；探测器实测的脉冲星信号到达航天器时间并转化到参考系原点，与预测的到达时间比对得到时间差，用来修正星载计算***记录的航天器位置。这里所述的选定的惯性参考系原点原则上可以是惯性坐标系中的任意惯性点，为了简化计算过程，通常该参考原点选定为常用的参考系原点，如太阳系质心SSB或地球质心等。SSB坐标系下，脉冲星脉冲相位到达SSB模型可表示为：

其中，代表参考历元t₀时刻的脉冲相位，v代表脉冲星自转频率，分别为v的一阶和二阶导数。式(8)也是脉冲星计时模型和特征参数数据库7中的脉冲星计时模型。该模型中的参数可通过对脉冲星的长期观测得到。脉冲星信号到达航天器时间转化到参考系原点的过程，实际上是实测TOA的固有时通过各项修正转化为质心力学时(TDB)时过程。在假设无引力无介质的平直空间中，该过程可用一阶方程表示为其中为SSB到脉冲星视线方向的矢量。而工程上实现精确的脉冲星定时，必须考虑相对论等效应对TOA测量的影响：

Δt_TOA＝Δ_RΘ+Δ_EΘ+Δ_SΘ+Δ_f，N+δ_t(9)

式中：N表示第N颗脉冲星；Δ_RΘ为Roemer延迟，包括一阶Doppler延迟和Parallax效应；Δ_EΘ为Einstein修正，太阳系中可忽略；Δ_SΘ为太阳系Shapiro延迟；为第N颗脉冲星色散效应误差，由于X射线频率非常高色散量很小，该项忽略；δ_t卫星时钟偏差。本发明中式(8)表示的脉冲星相位演化模型和式(9)表示的太阳系质心光行时转换模型建立在地球质心，星载计算机将式(8)表示的相位演化模型通过式(9)表示的光行时修正方法转换到地球中心，并将导航过程中由于地球运动造成的时间尺度变化，在相位演化模型转换中进行修正。

图1中，脉冲星计时模型和特征参数数据库还包括X射线脉冲星标准轮廓，X射线脉冲星高精度计时模型，X射线脉冲星星历星表、X射线脉冲星辐射流量和背景噪声强度、X射线脉冲星自行参数、X射线脉冲星周期等特征参数。

图1中，太阳系星行星参数库包括太阳系行星星历、太阳及行星引力常数、太阳系质心、地球质心以及其它在维持时空基准和脉冲星光行时转换中必需的参数。

图1中，导航算法库包括常用的最小二乘方法，扩展kalman方法，UnscentedKalman方法，粒子滤波方法和中心差分kalman滤波方法等。最小二乘和扩展kalman滤波方法稳定性较好，但在处理非高斯非线性状态方程时性能较弱，星载计算机优先考虑使用UnscentedKalman等非线性滤波方法，但同时将最小二乘和扩展kalman滤波作为备份滤波器使用，以在得到较好滤波性能的同时保证***的稳定可靠。但这要求星载计算机具有优秀的数据处理能力。

如图2所示，为本发明的在地球质心参考坐标系中的脉冲星、太阳质心、导航卫星、地平角距之间的几何关系以及基于计时观测导航算法和基于角距测量的导航算法原理图。从图2中可以看到地平角距夹角是卫星轨道位置的函数，地平角距可以使用红外地平敏感或紫外地平敏感设备获取。地平观测方向与地心和脉冲星矢量处于同一平面内，如图2所示，设卫星对地心矢量方向为它与脉冲星视线方向的矢量的夹角为β，地平仪观测的两侧地平方向夹角为α，地球半径为R_e。利用地平仪观测得到式(10)

$\left|{\stackrel{\→}{r}}_{U/S}\right|\sin \frac{\mathrm{\α}}{2}=R_e---\left(10\right)$

地心矢量方向可以利用地平观测经过α/2角度的旋转得到，它与脉冲星视线方向矢量如式(11)所示的关系

$\mathrm{\β}=\arccos (-\frac{{\stackrel{\→}{r}}_{U/S}\·{\hat{r}}_s}{\left|{\stackrel{\→}{r}}_{U/S}\right|})---\left(11\right)$

其中由探测器对脉冲星辐射矢量观测得到。

图2所示的基于计时观测的导航方法，TOA观测方程如式(9)所示，简化表示为：

Δt_TOA，N＝F_N(X，t)+W_N(12)

式中，Δt_TOA，NW为星载探测器测得的第N颗脉冲星到达用户星和SSB的时间差，X为用户星状态矢量，F_N为转换方程，W_N为测量噪声。令Δτ_N为第N颗脉冲星TOA观测与估计之间残差，对(12)式用泰勒级数展开，取一阶项，得到：

$\mathrm{\Δ}{t}_{\mathrm{TOA},N}=\mathrm{\Δ}{\tilde{t}}_{\mathrm{TOA},N}+\frac{\∂{\mathrm{\Δt}}_{\mathrm{TOA},N}}{\∂{\stackrel{\→}{r}}_{U/S}}\mathrm{\σX}+w_N\left(t\right)=G_N\mathrm{\σX}+w_N\left(t\right)---\left(13\right)$

其中w_N(t)为第N颗脉冲星的TOA量测噪声，E[w_N(t)]＝0， G_N为测量矩阵，如(14)式：

$G_N={\left[\begin{array}{c}{\stackrel{\→}{n}}_N+\frac{1}{D_0}\left[\begin{array}{c}({\stackrel{\→}{n}}_N\·r_{U/S})\left({\stackrel{\→}{n}}_N\right)-r_{U/S}\\ +\left(\mathrm{V\Δ}{t}_M\right)-({\stackrel{\→}{n}}_N\·\mathrm{V\Δ}{t}_M)\left({\stackrel{\→}{n}}_N\right)\\ -b+({\stackrel{\→}{n}}_N\·b)\left({\stackrel{\→}{n}}_N\right)\end{array}\right]\\ +\frac{2{\mathrm{\μ}}_{\mathrm{sun}}}{c^3}\ln |\frac{{\stackrel{\→}{n}}_N+(r_{U/S}/r_{U/S})}{({\stackrel{\→}{n}}_N\·r_{U/S}+r_{U/S})+{\stackrel{\→}{n}}_N\·b+b}+1|\end{array}\right]}^T---\left(14\right)$

其中D₀表示脉冲星相对于SSB的距离，表示SSB到第N脉冲星视线方向的单位矢量，b表示太阳系质心相对于太阳质心的位置矢量，V表示脉冲星的固有运动，Δt_M＝t_M-t₀表示从脉冲星初始位置开始到辐射第M个脉冲的时间间隔，由于V很小，D₀＞＞VΔt，仿真计算时VΔt可以忽略，r_U/S与r_U/S分别表示用户星相对于太阳系质心的位置矢量和它的模，μ_sun表示为太阳引力常数，c表示光速。

图2所示的基于角距观测的导航方法，其观测方程由式(12)得到观测角

$\mathrm{\α}=2\arcsin \frac{R_e}{{\stackrel{\→}{r}}_E}---\left(15\right)$

再进行一阶泰勒展开

$\mathrm{\α}=\widetilde{\mathrm{\α}}+{\left[\frac{\∂\mathrm{\α}}{\∂{\stackrel{\→}{r}}_{U/S}}\right]}^T\mathrm{\δX}+n_s$

$=\widetilde{\mathrm{\α}}-\frac{R_e}{\sqrt{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2-R_e^2}}\frac{{\stackrel{\→}{r}}_e^T}{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2}\mathrm{\δX}+n_s---\left(16\right)$

$=\widetilde{\mathrm{\α}}-F_s\mathrm{\δX}+n_s$

其中n_s包括泰勒展开的高阶项和噪声，简单起见，假设噪声高斯分布近似有E[n_s(t)]＝0，同样对(13)式做一阶泰勒近似，

$\mathrm{\β}=\widetilde{\mathrm{\β}}+{\left[\frac{\∂\mathrm{\β}}{\∂{\stackrel{\→}{r}}_{U/S}}\right]}^T\mathrm{\δX}+n_h$

$=\widetilde{\mathrm{\β}}+\frac{{\stackrel{\→}{r}}_s^T}{\sqrt{{\left|{\stackrel{\→}{r}}_{U/S}\right|}^2-{\stackrel{\→}{r}}_{U/S}\·r_s}}\mathrm{\δX}+n_h---\left(17\right)$

$=\widetilde{\mathrm{\α}}-F_h\mathrm{\δX}+n_h$

其中，n_h包括泰勒展开的高阶项和地心矢量观测噪声，E[n_h]＝0， $E\left[n_h\left(t\right){n_h}^T\left(t\right)\right]=$ ${\mathrm{\δ}}_h^2.$

如图2所示，融合计时观测导航算法和基于角距测量的导航算法中，仅利用计时观测或者仅利用矢量观测均存在秩亏问题，避免该问题通常需要多个独立的观测。在基于计时观测的脉冲星导航中，在原子钟足够稳定的前提下，至少需要观测3颗脉冲星信号。为了实现单脉冲星导航，这里我们用单脉冲星计时观测和矢量观测以及地平仪的视场角观测构成3个独立的观测，用于单脉冲星迭代定轨。综合考虑，(13)式、(16)式和(17)式所示的测量方程，设组合量测向量表示为Y＝[Δτ，Δα，Δβ]^T，其中 $\mathrm{\Δ\τ}=\mathrm{\Δ}{t}_{\mathrm{TOA},N}-\mathrm{\Δ}{\tilde{t}}_{\mathrm{TOA},N},$ $\mathrm{\Δ\α}=\mathrm{\α}-\widetilde{\mathrm{\α}},\mathrm{\Δ\β}=\mathrm{\β}-\widetilde{\mathrm{\β}},$ 此时观测方程为：

$Y=\mathrm{H\σX}={[F_s^T,F_h^T,G^T]}^T\mathrm{\δX}+\mathrm{\η}\left(t\right)---\left(18\right)$

其中H为观测矩阵，η(t)为零均值高斯观测噪声且η(t)＝[n_s(t)，n_h(t)，w_N(t)]^T，方差阵为： $R=E\left({\mathrm{\η\η}}^T\right)=\left(\begin{array}{ccc}{n}_s{\left(t\right)}^2& & \\ & n_h{\left(t\right)}^2& \\ & & w_N{\left(t\right)}^2\end{array}\right),$ 其w_N(t)²为脉冲星时间测量误差转换得到的位置量测误差方差，表示探测器对第i颗脉冲星时间测量误差，W为脉冲半宽度，

$\mathrm{SNR}=\frac{F_x{\mathrm{Ap}}_f\mathrm{\Δ}{t}_{\mathrm{obs}}}{\sqrt{[B_x+F_x(1-p_f)]\left(\mathrm{A\Δ}{t}_{\mathrm{obs}}d\right)+F_x{\mathrm{Ap}}_f\mathrm{\Δ}{t}_{\mathrm{obs}}}}---\left(19\right)$

其中F_x为X射线脉冲星辐射光子流，A表示探测器面积，p_f为一个脉冲周期内的脉冲辐射流量与平均辐射流量比，Δt_obs表示观测时间，B_x表示X射线背景辐射噪声，d＝W/P，P为脉冲周期。为了使提高TOA测量的精度，通常要使用累积的方法得到高信噪比的脉冲轮廓，此时Δt_obs反映了脉冲叠加段的长度。对于多脉冲星导航，由于计时观测方法与矢量观测方法均能实现自主迭代定轨，在(14)式和(17)式的基础上，多星观测矩阵是单星观测变量的有序组合，此时矢量观测矩阵表示为

Y₁＝H₁δX+η＝[F_h，1，F_h，2，...，F_h，n]^TδX+η₁(20)

其中，Y₁＝[Δβ₁，Δβ₂，...，Δβ_n]^T，噪声均值，E(η₁)＝0，同理到达时间观测矩阵表示为：

$Y_2=H_2\mathrm{\δX}+{\mathrm{\η}}_2={[G_1^T,G_2^T,...,G_N^T]}^T\mathrm{\δX}+{\mathrm{\η}}_2---\left(21\right)$

其中，Y₂＝[Δτ₁，Δτ₂，...，Δτ_n]^T，噪声均值，E(η₂)＝0，

如图3所示，基于角距测量的几何定向算法实现方法为：设轨道运动方向为脉冲星辐射方向矢量为n_p，t₀时刻的地心方向矢量分别为和则可表示为

${\widehat{\mathrm{\θ}}}_{t_0}=\arccos n_p\·n_{{\mathrm{rt}}_1}-\arccos n_p\·n_{{\mathrm{rt}}_0}---\left(22\right)$

如图4所示，融合计时观测和角距测量的脉冲星导航融合算法具体实现分以下几步完成：(1)初始化P₀，Q₀矩阵，给定航天器初始位置。(2)根据航天器位置和备择脉冲星在天球坐标系中的位置计算航天器可观测到的脉冲星，并调整观测方程。(3)根据(18)式和第2步的结果计算状态转移矩阵H，分两步完成EKF滤波。预测：

R(t_k+1，t_k)＝F(t)R(t_k，t_k)

$\mathrm{\δ}{\hat{X}}_{k+1}=R(t_{k+1},t_k)\mathrm{\δ}{X}_k+{\hat{q}}_k---\left(23\right)$

${\hat{P}}_{k+1}=R{\hat{P}}_k{R}^T+Q_k$

更新：

${\hat{X}}_{k+1}=X_{k+1}+\mathrm{\δ}{X}_{k+1}$

基于Kalman滤波的融合算法如图4所示，设X₁(k)，P₁(k)和X₂(k)，P₂(k)分别为子滤波器1和子滤波器2的局部估计值和估计方差，其全局最优融合方法如(24)式至(26)式：

${\hat{X}}_g\left(k\right)=P_g\left(k\right)\×[P_1^{-1}\left(k\right)X_1\left(k\right)+P_2^{-1}\left(k\right)X_2\left(k\right)]---\left(24\right)$

$P_g\left(k\right)={[P_1^{-1}\left(k\right)+P_2^{-1}\left(k\right)]}^{-1}---\left(25\right)$

${\hat{X}}_i\left(k\right)={\hat{X}}_g\left(k\right)---\left(26\right)$

子滤波器中只反馈状态变量。

本发明以地球卫星为例对所发明的***和方法进行了阐述，而从理论方法上讲，本发明的融合脉冲星辐射矢量和计时观测的卫星自主导航***及方法完全适用于除地球外的其它近天体轨道航天器，仅需要将时间和坐标***转换到目标惯性系下，该转换方法有完备的理论和方法支撑，除此之外其它组成可以不加修改的使用适用。本领域专业人员在不背离本发明权利要求范围和主旨的前提下可以实现多种显而易见的改进，因此上述内容只是本发明借以阐述的实施例，本发明的权利要求范围并不限制于以上论述。

本说明未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (10)

1.融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，其特征在于包括X射线探测器、准直器、X射线脉冲星定时模型和参数数据库、太阳系行星参数库，星载原子钟、星载计算机，红外地平敏感或紫外地平敏感、导航算法库、星载万向支架；星载计算机根据当前记录的航天器姿态和X射线脉冲星参数数据库，控制星载万向架支撑星载X射线探测器指向目标脉冲星；X射线探测器配合准直器，通过测量到达探测器的X射线脉冲星的光子流量强度计算脉冲星指向与脉冲星辐射矢量的偏差，该偏差反馈到星载计算机用于增强姿态控制，精确获取脉冲星辐射方向矢量；在探测器准确指向目标脉冲星的同时，X射线探测器记录X射线光子到达时间，在星载脉冲星参数库和太阳系行星参数库支持下将光子到达时间改正转换到太阳系质心；转换后的光子按照周期统计其在采样时间分辨单位内的重复频率直方图，进而形成累积轮廓；该累积轮廓与脉冲星标准平均轮廓比相，获取累计轮廓位相，该位相再与星载原子钟保持的脉冲星相位演化模型做差，所得相差转换成距离差，送至导航计算机中的基于地心坐标系的脉冲星计时观测导航算法；同时地平敏感仪敏感地平得到地平角距，并计算地心方向矢量，该角距与上述脉冲星辐射方向矢量共同送至星载计算机中基于角距测量的脉冲星定位算法，同时该地心方向矢量与上述脉冲星辐射方向矢量共同送至基于角距测量的几何定向算法；星载计算机调用自主导航算法库中的融合导航滤波算法对基于脉冲星计时观测的导航算法、基于角距测量的脉冲星定位算法和基于角距测量的几何定向算法融合，输出高精度的自主导航参数信息，包括位置、姿态和速度信息。

2.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的X射线探测器为时间分辨探测器，所述的准直器安装于探测器前方，准直器由一系列平行的空心柱体紧密排列，只允许与柱体平行的光线通过，到达探测器，依次达到背景噪声滤波、辐射矢量识别的目的。

3.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的重复频率直方图的生成方法为：假设为一段X射线脉冲星观测信号，它由N个采样构成，采样间隔为δ_t，起始时刻为t₀，周期为T_P，Nδt＝mT_P，m为观测信号中的周期数，z_i(i＝1，...，N)表示每个采样间隔内光子计数，写成向量形式为z＝[z₁，z₂，...，z_N]^T，令为了区别观测信号，用x＝[x₁，x₂，...，x_m]^T表示累积后的信号，令p_i＝|x_i|/||z||，重复频率直方图的生成方法依次调整变量i，使x的熵： $H\left(x\right)=-\underset{i=1}{\overset{m}{\Σ}}{p}_i{\mathrm{logp}}_i$ 最小。

4.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的脉冲星计时观测导航算法的实现方法为：星载计算机利用太阳系行星参数数据库和脉冲星计时模型和参数数据库计算地球坐标系中心相对于太阳系质心的位置和速度，并通过光行时转换和时间转换将建立在太阳系质心和太阳系质心力学时下的脉冲星相位演化模型转换到地球坐标系中心和地球质心力学时下；同时星载计算机在星载原子钟支持下，维持地心坐标系下的脉冲星相位演化模型；星载X射线探测器记录X射线脉冲星光子到达时间，整合出累积脉冲轮廓，该轮廓与所述的地心坐标系下的脉冲星相位演化模型比相获得相差，进而利用相差计算距离误差，该距离误差送入导航算法库；所述的地平角距为地心坐标系下卫星到地球两侧表面切线方向矢量的夹角；所述的地心方向矢量为地球坐标系下卫星中心到地球中心连线方向矢量，计算方法为地平角距的中心线方向；所述的脉冲星辐射矢量的识别方法为：以万向支架支撑点为中心，在一个圆球面上以往返形式等幅度步进探测器指向，当探测器接收到的目标脉冲星辐射流量强度最大时即认为探测器指向与脉冲星辐射方向矢量平行。

5.根据权利要求4所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的脉冲星辐射流量强度表示方法为对探测器记录的光子序列做FFT变换，统计目标脉冲星信号在其特征频率点的平均功率，用于标定流量强度。

6.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的基于角距测量的脉冲星定位算法的实现方法为：星载计算机同时接收X射线探测器输出的脉冲星矢量信息和地平敏感仪输出的地心方向矢量以及地平仪输出的地平角距，首先建立地平角距与卫星位置矢量的函数关系式，命名为关系式A，再建立地心矢量方向与脉冲星辐射方向矢量夹角与卫星位置的函数关系式，命名为关系式B；利用关系式A和关系式B建立观测方程送入导航滤波算法。

7.根据权利要求6所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的关系式A为 $\left|{\stackrel{\→}{r}}_{U/S}\right|sin\frac{\mathrm{\α}}{2}=R_e$ 所述的关系式B为 $\mathrm{\β}=\arccos (-\frac{{\stackrel{\→}{r}}_{U/S}\·{\stackrel{\→}{r}}_s}{\left|{\stackrel{\→}{r}}_{U/S}\right|});$

.所述的观测方程为：设卫星对地心矢量方向为它与脉冲星视线方向的矢量的夹角为β，地平仪观测的两侧地平方向夹角为α，地球半径为R_e；

$\mathrm{\α}=\widetilde{\mathrm{\α}}+{\[\frac{\∂\mathrm{\α}}{\∂{\stackrel{\→}{r}}_{U/S}}\]}^T\mathrm{\δ}X+n_s=\widetilde{\mathrm{\α}}-\frac{R_e}{\sqrt{|{\stackrel{\→}{r}}_{U/S}{|}^2-R_e^2}}\frac{{\stackrel{\→}{r}}_e^T}{|{\stackrel{\→}{r}}_{U/S}{|}^2}\mathrm{\δ}X+n_s=\widetilde{\mathrm{\α}}-F_s\mathrm{\δ}X+n_s---\left(1\right)$ 其中n_s包括泰勒展开的高阶项和噪声，简单起见，假设噪声服从高斯分布近似有E[n_s(t)]＝0， $E\[n_s\left(t\right)R_s^T\left(t\right)\]={\mathrm{\δ}}_s^2;$

$\mathrm{\β}=\widetilde{\mathrm{\β}}+{\[\frac{\∂\mathrm{\β}}{\∂{\stackrel{\→}{r}}_{U/S}}\]}^T\mathrm{\δ}X+n_h=\widetilde{\mathrm{\β}}+\frac{{\stackrel{\→}{r}}_s^T}{\sqrt{|{\stackrel{\→}{r}}_{U/S}{|}^2-|{\stackrel{\→}{r}}_{U/S}\·r_s}}\mathrm{\δ}X+n_h=\widetilde{\mathrm{\α}}-F_h\mathrm{\δ}X+n_h---\left(2\right)$ 其中，n_h包括泰勒展开的高阶项和地心矢量观测噪声，E[n_h]＝0，

8.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的基于角距测量的几何定向算法为利用卫星在轨运动过程中通过地平敏感获取的地心矢量与脉冲星辐射方向矢量角距变化获取航天器运动方向；实现方法为：设轨道运动方向用脉冲星辐射方向矢量为n_p，t₀时刻的地心方向矢量分别为和则可表示为

${\widehat{\mathrm{\θ}}}_{t_0}=arc\cos n_p\·n_{{\mathrm{rt}}_1}-arc\cos n_p\·n_{{\mathrm{rt}}_0}---\left(3\right)$

所述的导航滤波算法观测矩阵由，(1)式、(2)式和(3)式所示的测量方程组成，设组合量测向量表示为Y＝[Δτ，Δα，Δβ]^T，其中 $\mathrm{\Δ}\mathrm{\τ}={\mathrm{\Δt}}_{TOA,N}-\mathrm{\Δ}{\tilde{t}}_{TOA,N},\mathrm{\Δ}\mathrm{\α}=\mathrm{\α}-\widetilde{\mathrm{\α}},$ $\mathrm{\Δ}\mathrm{\β}=\mathrm{\β}-\widetilde{\mathrm{\β}}$ 此时观测方程为：

$Y=H\mathrm{\σ}X={\[F_s^T,F_h^T,G^T\]}^T\mathrm{\δ}X+\mathrm{\η}\left(t\right)---\left(4\right)$

其中H为观测矩阵，η(t)为零均值高斯观测噪声且η(t)＝[n_s(t)，n_h(t)，w_N(t)]^T，方差阵为：

$R=E\left({\mathrm{\η\η}}^r\right)=\left(\begin{array}{ccc}{n}_s{\left(t\right)}^2& & \\ & n_h{\left(t\right)}^2& \\ & & w_N{\left(t\right)}^2\end{array}\right)$

其w_N(t)²为脉冲星时间测量误差转换得到的位置量测误差方差，表示探测器对第i颗脉冲星时间测量误差，W为脉冲半宽度， $SNR=\frac{F_x{\mathrm{Ap}}_f{\mathrm{\Δt}}_{obs}}{\sqrt{\[B_x+F_x(1-p_f)\]\left({\mathrm{A\Δt}}_{obs}d\right)+F_x{\mathrm{Ap}}_f{\mathrm{\Δt}}_{obs}}}$ 其中F_x为X射线脉冲星辐射光子流，A表示探测器面积，p_f为一个脉冲周期内的脉冲辐射流量与平均辐射流量比，Δt_obs表示观测时间，B_x表示X射线背景辐射噪声，d＝W/P，P为脉冲周期；为了提高TOA测量的精度，通常要使用累积的方法得到高信噪比的脉冲轮廓，此时Δt_obs反映了脉冲叠加段的长度。

9.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，所述的导航滤波算法的特征在于：在(2)式和(1)式的基础上，多星观测矩阵是单星观测变量的有序组合，此时矢量观测矩阵表示为

Y₁＝H₁δX+η＝[F_h，1，F_h，2，...，F_h，n]^TδX+η₁(5)

其中，Y₁＝[Δβ₁，Δβ₂，...，Δβ_n]^T，噪声均值，E(η₁)＝0，

同理到达时间观测矩阵表示为：

$Y_2=H_2\mathrm{\δ}X+{\mathrm{\η}}_2={\[G_1^T,G_2^T,...,G_N^T\]}^T\mathrm{\δ}X+{\mathrm{\η}}_2---\left(6\right)$

其中，Y₂＝[Δτ₁，Δτ₂，...，Δτ_n]^T，噪声均值，E(η₂)＝0，

10.根据权利要求1所述的融合脉冲星辐射方向矢量和计时观测的航天器自主导航定位方法，其特征在于：基于计时观测的导航或者基于角距观测均可以独立完成近地轨道卫星导航定位，二者也可以共同用于导航以增加导航精度，适合于近地轨道卫星或者卫星星座高精度自主导航。