CN101608921B - 一种脉冲星/cns组合导航方法 - Google Patents

Abstract

一种脉冲星/CNS组合导航方法，属于航天器自主轨道确定领域。该方法设计了UKF/H_∞滤波器，融合天文导航和脉冲星导航***(CNS)提供的导航定位数据，并在脉冲星导航***中，提出了历元间差分法以消除时钟误差和脉冲星位置误差。本发明能够获得高精度的航天器位置信息，为航天器提供长时间、高精度的自主导航，使其定位精度不受星载时钟和测量设备的限制。

Description

一种脉冲星/CNS组合导航方法

技术领域

本发明涉及一种航天飞行器自主导航方法。

背景技术

目前，以卫星为代表的航天飞行器的轨道确定和预报需要依靠地面测控网的支持，随着进入太空的卫星数目的急剧增加，地面测控站的负担越来越重。为了降低地面支持成本，或不依赖地面支持，就要求卫星实现自主导航。自主导航技术是实现卫星自主的一个重要方面，是当今卫星控制技术的发展趋势，它对减轻地面测控负担、降低卫星运行费用、提高卫星生存能力尤其是深空探测器的生存能力等方面都有重要意义。

在现有的的航天器导航方式中，只有天文导航***(CNS，CelestialNavigation System)和脉冲星导航是完全自主导航方式。CNS依靠天体敏感器观测天体的方位信息进行定位导航，其定位误差与时间无关，滤波周期短，但是由于地平仪精度低，该方法定位性能较差。脉冲星导航是基于脉冲星发射的脉冲信号和已知的脉冲星位置信息实现航天器轨道确定的技术，定位精度高，但是，受量测手段的限制，脉冲星位置信息存在误差，该误差严重影响脉冲星导航***性能，另外该方法定位精度与脉冲到达时间(Time of Arrival，TOA)精度相关，而TOA精度又受时钟误差影响，时钟误差随着时间流逝而增大。

CNS和脉冲星导航各有优缺点，如果把这两种以及其它多种导航***组合起来，应用最优估计理论，可以使组合导航***在精度方面有所提高。在申请号为200710191527.7的中国专利文件公布了一种基于多信息融合的组合导航方法，该方法融合了四种不同导航信息，其中就包含有CNS和脉冲星导航信息。该方法在一定程度上提高了定位精度，但是没有考虑脉冲星的脉冲星位置误差和时钟误差对定位的影响。

发明内容

本发明的目的在于提供一种脉冲星/CNS组合导航方法，有效消除脉冲星位置误差和时钟误差，提高了定位精度。

一种脉冲星/CNS组合导航方法，具体为：在历元t，将历元t-1得到的卫星位置估计值和历元t时的星光仰角观测值作为无迹卡尔曼滤波器的输入，无迹卡尔曼滤波器的输出为中间卫星位置估计值；判断历元t是否等于N_pk，N_p为脉冲星导航与天文导航的采样周期比值，k为自然数，若不等于，则中间卫星位置估计值即为最终卫星位置估计值，否则，将中间卫星位置估计值和在历元N_pk时脉冲星导航观测的脉冲到达时间差作为H_∞滤波器的输入，H_∞滤波器的输出即为最终卫星位置估计值；

在历元N_pk时脉冲星导航观测的脉冲到达时间差 ${\stackrel{\‾}{Z}}_2\left(N_{p}k\right)={\tilde{Z}}_2(N_{p}k-N_p)+Z_2\left(N_{p}k\right)-Z_2(N_{p}k-N_p),$ 历元N_pk-N_p时脉冲星导航的量测参考值 ${\tilde{Z}}_2(N_{p}k-N_p)=h_2(\tilde{X}(N_{p}k-N_p),N_{p}k-N_p),$ h₂( )为脉冲星导航量测方程，为对历元(N_pk-N_p+1)，(N_pk-N_p+2)，……，(N_pk-1)时获取的星光仰角作平滑处理后得到的值，Z₂(N_pk)和Z₂(N_pk-N_p)分别为脉冲星导航在历元N_pk和N_pk-N_p实际量测的脉冲到达时间差。

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

(1)本发明在获取脉冲星导航观测的脉冲到达时间差时，利用了当前历元和前一历元的实际量测脉冲作差分运算，差分结果使得时钟误差和脉冲星位置误差各自抵消，既消除了时钟误差影响，使脉冲星导航不受星载时钟精度的限制，降低了星载时钟的成本；又消除了脉冲星位置误差的影响，使脉冲星导航不受量测技术手段的限制，使其在现有技术水平下就能达到较高定位精度。

(2)本发明采用UKF/H_∞滤波器融合天文导航和脉冲星导航提供的定位信息，在保证***鲁棒性的前提下达到精确定位的效果，具有计算量小，实时性强，设计简单，易于工程实现等优点。

附图说明

图1为本发明脉冲星/CNS导航流程图。

具体实施方式

下面以卫星为例对本发明作进一步说明。

卫星轨道动力学模型表示为

$\stackrel{\·}{X}\left(t\right)=f(X\left(t\right),t)+w\left(t\right)---\left(1\right)$

其中，t为历元，状态矢量X(t)＝[x y z v_x v_y v_z]^T，v_x在X轴方向的位置x处的速度，v_y在Y轴方向的位置y处的速度，v_z在Z轴方向的位置z处的速度。 $w\left(t\right)={\left[\begin{array}{cccccc}{w}_x& w_y& w_z& w_{v_x}& w_{v_y}& w_{v_z}\end{array}\right]}^T$ 为X(t)＝[x y z v_x v_y v_z]^T对应的状态处理噪声，其协方差矩阵为Q，T表示矢量转置。该模型具体表达式如下：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x+w_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y+w_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z+w_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{{\left|r\right|}^3}[1-J_2{\left(\frac{R_e}{\left|r\right|}\right)}^2(7.5\frac{z^2}{{\left|r\right|}^2}-1.5)]+{\mathrm{\ΔF}}_x+w_{v_x}\\ \frac{{\mathrm{dv}}_t}{\mathrm{dt}}=-\mathrm{\μ}\frac{y}{{\left|r\right|}^3}[1-J_2{\left(\frac{R_e}{\left|r\right|}\right)}^2(7.5\frac{z^2}{{\left|r\right|}^2}-1.5)]+{\mathrm{\ΔF}}_y+w_{v_y}\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mathrm{\μ}\frac{z}{{\left|r\right|}^3}[1-J_2{\left(\frac{R_e}{\left|r\right|}\right)}^2(7.5\frac{z^2}{{\left|r\right|}^2}-4.5)]+{\mathrm{\ΔF}}_z+w_{v_z}\end{array}\right.---\left(2\right)$

其中，r＝[x，y，z]^T为卫星相对于地球的位置矢量，则航天器与地球质心之间的距离 $\left|r\right|=\sqrt{x^2+y^2+z^2},$ μ是地球引力常数，J₂为二阶带谐项系数，R_e为地球平均半径，ΔF_x，ΔF_y，ΔF_z分别为摄动力在X、Y、Z轴三个方向上的影响，摄动力涉及有地球非球形摄动的高阶摄动项、日月摄动、太阳光压摄动以及大气摄动。

对卫星的组合导航步骤具体为：

1)构建天文导航量测模型

天文导航量测模型如下：

$Z_1=h_1(X\left(t\right),t)+v_1\left(t\right)=\arccos (-\frac{r\·s}{\left|r\right|})-\arcsin \left(\frac{R_e}{\left|r\right|}\right)+v_1\left(t\right)---\left(3\right)$

其中，z₁为星光仰角观测量，s是恒星的方向矢量，v₁(t)为天文导航量测噪声，天文导航量测方程h₁(X(t)，t)表达式如下：

$h_1(X\left(t\right),t)=\arccos (-\frac{r\·s}{\left|r\right|})-\arcsin \left(\frac{R_e}{\left|r\right|}\right)---\left(4\right)$

2)构建脉冲星导航量测模型

(21)构建传统的脉冲星导航量测模型

选择Crab(蟹状星云)脉冲星作为导航星。卫星量测的脉冲TOA为t_SC，其在SSB(Solar System Barycentre，太阳系质心)处的相应时间t_b。脉冲到达时间差t_b-t_SC与卫星位置的关系为：

$t_b-t_{\mathrm{SC}}=\frac{1}{c}n\·r_{\mathrm{SC}}+\frac{1}{2c{D}_0}[-r_{\mathrm{SC}}^2+{(n\·r_{\mathrm{SC}})}^2-2b\·r_{\mathrm{SC}}+2(n\·b)(n\·r_{\mathrm{SC}})]+\frac{2{\mathrm{\μ}}_{\mathrm{Sum}}}{c^3}\ln |\frac{n\·r_{\mathrm{SC}}+r_{\mathrm{SC}}}{n\·b+b}+1|---\left(5\right)$

其中，c为光速，D₀为脉冲星到太阳系质心的距离，b为SSB相对于太阳的位置矢量，μ_Sum为太阳引力常数，n是脉冲星的实际方向矢量，r_SC是卫星相对于SSB的位置矢量，r_E为地球相对于SSB的位置矢量。

卫星相对于SSB的位置矢量r_SC与卫星相对于地球的位置矢量r满足方程式(6)：

r＝r_SC-r_E                     (6)

结合方程式(5)和(6)得到传统的脉冲星导航量测模型，表示为：

z₂＝h₂(X(t)，t)+V₂(t)＝t_b-t_SC+v₂(t)          (7)

其中，Z₂为脉冲到达时间差t-t_SC，v₂(t)为脉冲星导航量测噪声，脉冲星导航量测方程h₂(X(t)，t)表达式如下：

$h_2(X\left(t\right),t)=\frac{1}{c}\hat{n}\·r_{\mathrm{SC}}+\frac{1}{2c{D}_0}[-r_{\mathrm{SC}}^2+{(\hat{n}\·r_{\mathrm{SC}})}^2-2b\·r_{\mathrm{SC}}+2(\hat{n}\·b)(\hat{n}\·r_{\mathrm{SC}})]+\frac{2{\mathrm{\μ}}_{\mathrm{Sum}}}{c^3}\ln |\frac{\hat{n}\·r_{\mathrm{SC}}+r_{\mathrm{SC}}}{\hat{n}\·b+b}+1|---\left(8\right)$

为量测的脉冲星方向矢量。

(22)消除脉冲星导航量测模型的脉冲星位置误差和时钟误差。

脉冲星导航量测噪声v₂(t)可表示为：

v₂(t)＝v_c(t)+v_p(t)+v_m(t)                     (9)

其中，v_c(t)为钟差，v_m(t)为TOA量测噪声，其方差为σ_R ²。方差可按Taylor(泰勒)方法进行计算。v_p(t)为脉冲星位置误差引起的***偏差，可近似表示为

$v_p\left(t\right)=(\hat{n}-n)\·r_{\mathrm{SC}}/c---\left(10\right)$

从式(9)可以看出，脉冲星导航量测噪声由三部分组成。其中，v_m(t)为白噪声，可用滤波方法进行处理，v_c(t)和v_p(t)在短时间内可看作常值偏差，不能用滤波的方法消除其影响，本发明的核心处就是提出历元间差分法以消除v_c(t)和v_p(t)。

为便于说明，先介绍采样周期比的含义。采样周期比N_p按下式计算：

N_p＝P_p/P_CNS                                  (11)

其中，P_p和P_CNS分别为脉冲星导航和CNS的采样周期。选择合适的P_p和P_CNS，使N_p为整数。也就是说，每个历元都可观测星光仰角，只在历元N_pk(k为自然数)才能观测到脉冲到达时间差。

在历元间差分法中，首先计算前一次量测的估计值，再加上本次与前一次实际量测之差，即可得到一个不含v_c(t)和v_p(t)的新量测，具体如下：

步骤221：依据天文导航方式提供的星光仰角量测信息，采用Unscented(无迹)Rauch-Tung-Striebel平滑器(见参考文献S.Sarkka，“Unscented Rauch-Tung-Striebel Smoother”，IEEE Transaction on AutomaticControl，Vol.53，No.3，pp845-849，2008.)处理历元(N_pk-N_p+1)～(N_pk-1)时的状态估计，可得到卫星位置的参考状态

和参考状态方差

从而获得历元N_pk-N_p时脉冲星导航的量测参考值

${\tilde{Z}}_2(N_{p}k-N_p)=h_2(\tilde{X}(N_{p}k-N_p),N_{p}k-N_p)$

$=h_2(X(N_{p}k-N_p),N_{p}k-N_p)+v_{\mathrm{process}}(N_{p}k-N_p)---\left(12\right)$

其中，v_process是均值为0，方差为R_est的噪声。

$R_{\mathrm{est}}=E\left[v_p(N_{p}k-N_p)v_p(N_{p}k-N_p)\right]=H^T(N_{p}k-N_p)\tilde{P}(N_{p}k-N_p)H(N_{p}k-N_p)---\left(13\right)$

其中，H( )为量测矩阵，E[ ]表示数学期望。

$H(N_{p}k-N_p)=\frac{{\∂h}_2(X\left(t\right),t)}{\∂X\left(t\right)}{|}_{X\left(t\right)=\tilde{X}(N_{p}k-N_p)}---\left(14\right)$

步骤222：构造脉冲星导航新的量测模型 ${\stackrel{\‾}{Z}}_2\left(N_{p}k\right)={\tilde{Z}}_2(N_{p}k-N_p)+Z_2\left(N_{p}k\right)-Z_2(N_{p}k-N_p)$

其中，z₂(N_pk-N_p)和Z₂(N_pk)分别表示前一次与本次实际量测值。进一步推导，得到

${\stackrel{\‾}{Z}}_2\left(N_{p}k\right)={\tilde{Z}}_2(N_{p}k-N_p)+Z_2\left(N_{p}k\right)-Z_2(N_{p}k-N_p)$

$=h_2(X(N_{p}k-N_p),N_{p}k-N_p)+(t_b\left(N_{p}k\right)-t_{\mathrm{SC}}\left(N_{p}k\right))-(t_b(N_{p}k-N_p)-t_{\mathrm{SC}}(N_{p}k-N_p))+{\stackrel{\‾}{v}}_2\left(N_{p}k\right)$

$=t_b(N_{p}k-N_p)-t_{\mathrm{SC}}(N_{p}k-N_p)+t_b\left(N_{p}k\right)-t_{\mathrm{SC}}\left(N_{p}k\right)-(t_b(N_{p}k-N_p)-t_{\mathrm{SC}}(N_{p}k-N_p))+{\stackrel{\‾}{v}}_2\left(N_{p}k\right)$

$=h_2(X\left(N_{p}k\right),N_{p}k)+{\stackrel{\‾}{v}}_2\left(N_{p}k\right)---\left(15\right)$

新的量测噪声为有色噪声，其表达式如下：

${\stackrel{\‾}{v}}_2\left(N_{p}k\right)=v_{\mathrm{process}}(N_{p}k-N_p)+v_2\left(N_{p}k\right)-v_2(N_{p}k-N_p)$

${=v}_{\mathrm{process}}(N_{p}k-N_p)+v_m\left(N_{p}k\right)-v_m(N_{p}k-N_p)$

可以看出，新的量测噪声

与钟差v_c(t)和脉冲星位置误差v_p(t)无关，克服了钟差和脉冲星位置误差对整个导航的影响。

本发明采用的脉冲星导航量测为

脉冲星导航噪声方差 $R_2={2\mathrm{\σ}}_R^2+R_{\mathrm{est}}.$

(4)采用UKF/H_∞滤波器来融合导航信息

本发明提出了UKF/H_∞滤波器来融合导航数据。该滤波器由UKF和H_∞滤波器串联而成，如图1所示。在UKF/H_∞滤波器中，H_∞滤波器的状态处理噪声和状态转移矩阵分别为0和单位阵。

CNS和脉冲星导航子***分别采用UKF和H_∞滤波器，即星光仰角和基于历元间差分法的脉冲到达时间差分别作为UKF和H_∞滤波器的量测。具体的滤波流程分两种情况：(1)无脉冲到达时间差。在脉冲观测期间，仅CNS运行并输出状态估计值。(2)有脉冲到达时间差。一旦产生了一个脉冲到达时间差，利用脉冲到达时间差对CNS输出的状态估计值进行更新，从而获得高精度的状态估计值。从以上流程可以看出，天文导航和脉冲星导航提供的导航信息都对状态值进行了修正，实现了信息融合。

UKF/H_∞滤波处理的具体过程如下：

给定***模型：

X(t+1)＝f(X(t)，t)+w(t)                     (16)

Z₁(t)＝h₁(X(t)，t)+v₁(t)                    (17)

${\stackrel{\‾}{Z}}_2\left(t\right)=h_2(X\left(t\right),t)+{\stackrel{\‾}{v}}_t\left(t\right)---\left(18\right)$

t为第t个历元。

步骤41：UKF滤波过程(天文导航)

利用天文导航提供的导航信息，按照式(19)～(34)所示的UKF滤波过程运行，可得到卫星的状态估计。在脉冲TOA观测期间，脉冲星导航无法提供导航信息，无需执行步骤42。

(a)初始化

$\hat{X}\left(0\right)=E\left[X\left(0\right)\right]---\left(19\right)$

$P\left(0\right)=E\left[{(X\left(0\right)-\hat{X}\left(0\right))\left(X\left(0\right)\hat{X}\left(0\right)\right))}^T\right]---\left(20\right)$

其中，E[ ]表示数学期望。

(b)计算希格玛(sigma)点

$\mathrm{\χ}(t-1)=\left[\begin{array}{ccc}\hat{X}(t-1)& \hat{X}(t-1)+\sqrt{m+\mathrm{\τ}}{\left(\sqrt{P(t-1)}\right)}_i& \hat{X}(t-1)-\sqrt{m+\mathrm{\τ}}(\sqrt{P(t-1)}{)}_i\end{array}\right]---\left(21\right)$

W₀＝τ/(m+τ)                           (22)

W_i＝1/[2(m+τ)]                         (23)

W_i+m＝1/[2(m+τ)]                       (24)

其中，τ∈R，是矩阵平方根

的第i列，i＝1，2，…m，m为状态矢量的维数。

(c)时间更新

χ(t/t-1)＝f(χ(t-1)，t-1)              (25)

其中，χ(t/t-1)表示t-1时刻的一步预测希格玛点。

$\hat{X}\left(t\right)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i\mathrm{\χ}(t/t-1)---\left(26\right)$

$P\left(t\right)=\mathrm{\Σ}{W}_i{[{\mathrm{\χ}}_i(t/t-1)-\hat{X}\left(t\right)][{\mathrm{\χ}}_i(t/t-1)-\hat{X}\left(t\right)]}^T+Q\left(t\right)---\left(27\right)$

其中，χ_i(t/t-1)表示χ(t/t-1)的第i列。

Z₁(t/t-1)＝h(χ(t/t-1)，t)

其中，Z₁(t/t-1)表示t-1时刻的一步预测量测值。

${\hat{Z}}_1\left(t\right)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i{Z}_1^i(t/t-1)---\left(29\right)$

其中，Z₁ ⁱ(t/t-1)表示Z₁(t/t-1)的第i列。

(d)量测更新

$P_{\hat{z}\hat{z}}\left(t\right)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i[Z_1^i(t/t-1)-{\hat{Z}}_1\left(t\right)]{[Z_1^i(t/t-1)-{\hat{Z}}_1\left(t\right)]}^T+R_1\left(t\right)---\left(30\right)$

$P_{\mathrm{xz}}\left(t\right)=\underset{i=0}{\overset{2m}{\mathrm{\Σ}}}{W}_i{[{\mathrm{\χ}}_i(t/t-1)-\hat{X}\left(t\right)][Z_1^i(t/t-1)-{\hat{Z}}_1\left(t\right)]}^T---\left(31\right)$

$G_1\left(t\right)=P_{\mathrm{xz}}\left(t\right)P_{\hat{z}\hat{z}}^{-1}\left(t\right)---\left(32\right)$

$\hat{X}\left(t\right)=\hat{X}\left(t\right)+G_1\left(t\right)(Z_1\left(t\right)-{\hat{Z}}_1\left(t\right))---\left(33\right)$

$P\left(t\right)=P\left(t\right)-G_1\left(t\right)P_{\hat{z}\hat{z}}\left(t\right)G^T\left(t\right)---\left(34\right)$

其中，Q(t)和R₁(t)分别为处理噪声和天文导航量测噪声协方差阵。当X(t)为高斯分布时，通常取m+τ＝3，T表示矩阵转置。

步骤42：H_∞滤波过程(脉冲星导航)

一旦获得了一个脉冲到达时间差，则执行此步骤。利用脉冲到达时间差对天文导航修正过的状态估计进行更新，即按照式(35)～(38)进行处理。这样，天文导航和脉冲星导航提供的信息都对卫星的状态估计进行了更新，也就实现了信息融合。

估计过程如下：

$G_2\left(t\right)=P\left(t\right){[I-\mathrm{\θP}\left(t\right)+H^T\left(t\right)R_2^{-1}\left(t\right)H\left(t\right)P\left(t\right)]}^{-1}{H}^T\left(t\right)R_2^{-1}\left(t\right)---\left(35\right)$

$P\left(t\right)=P\left(t\right){[I-\mathrm{\θP}\left(t\right)+H^T\left(t\right)R_2^{-1}\left(t\right)H\left(t\right)P\left(t\right)]}^{-1}---\left(36\right)$

$\hat{X}\left(t\right)=\hat{X}\left(t\right)+G_2\left(t\right)({\stackrel{\‾}{Z}}_2\left(t\right)-h_2(\hat{X}\left(t\right),t))---\left(37\right)$

其中，θ是用户设定的阈值，取值范围为10^-1～10^-6。R₂(t)为脉冲星导航量测噪声协方差阵。

为使该方案有解，必须满足以下条件。

$P^{-1}\left(t\right)-\mathrm{\θI}+H^Y\left(t\right)R_2^{-1}\left(t\right)H\left(t\right)>0---\left(38\right)$

其中，I为单位矩阵。

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

Claims (1)

1.一种脉冲星/天文导航***组合导航方法，具体为：在历元t，将历元t-1得到的卫星位置估计值和历元t时的星光仰角观测值作为无迹卡尔曼滤波器的输入，无迹卡尔曼滤波器的输出为中间卫星位置估计值；判断历元t是否等于N_pk，N_p为脉冲星导航与天文导航的采样周期比值，k为自然数，若不等于，则中间卫星位置估计值即为最终卫星位置估计值，否则，将中间卫星位置估计值和在历元N_pk时脉冲星导航观测的脉冲到达时间差作为H_∞滤波器的输入，H_∞滤波器的输出即为最终卫星位置估计值；

在历元N_pk时脉冲星导航观测的脉冲到达时间差

历元N_pk-N_p时脉冲星导航的量测参考值

h₂()为脉冲星导航量测表达式，

为对历元(N_pk-N_p+1)，(N_pk-N_p+2)，......，(N_pk-1)时获取的星光仰角作平滑处理后得到的值，Z₂(N_pk)和Z₂(N_pk-N_p)分别为脉冲星导航在历元N_pk和N_pk-N_p实际量测的脉冲到达时间差。

Images