CN104864875A - 一种基于非线性h∞滤波的航天器自主定位方法 - Google Patents

Abstract

一种基于非线性H∞滤波的航天器自主定位方法，选择航天器的位置矢量和速度矢量在地心惯性系的投影作为状态变量，将航天器轨道动力学方程作为系统模型，将脉冲星时间转换方程作为测量模型，将脉冲星星表误差描述为模型中的有界误差，基于模型设计用于状态估计的非线性H∞滤波算法；进而，采用所设计的非线性H∞滤波算法，处理脉冲到达时间观测量序列，通过递推计算估计出航天器的位置和速度。本发明所述方法能够增强X射线脉冲星导航系统克服星表误差影响的能力，改善航天器自主定位精度。

Description

一种基于非线性H∞滤波的航天器自主定位方法

技术领域

本发明涉及一种基于非线性H∞滤波的航天器自主定位方法，属于航天器自主导航技术领域。

背景技术

X射线脉冲星导航作为一项新兴的航天器自主导航技术，是当前国内外航天技术领域的研究热点之一。由NASA(美国航空航天局)和DARPA(国防高级研究计划局)资助的研究团队所取得的成果展示了X射线脉冲星在导航领域的巨大潜力，国内中科院高能所、西光所和中国空间技术研究院等正在进行用于空间导航的X射线探测器的研制。

典型的脉冲星导航方式是基于航天器轨道动力学模型建立状态方程，基于脉冲星角位置信息和脉冲星时间模型建立观测方程，采用适当的滤波算法处理由X射线探测器获得的脉冲到达时间(TOA)观测量，估计得到航天器的位置矢量和速度矢量。受目前X射线脉冲星天文观测水平和模型精度的限制，用于建立观测方程的脉冲星角位置数据不可避免会存在误差。脉冲星角位置误差通常称为星表误差，星表误差是影响脉冲星导航系统性能的主要因素之一。

事实上，星表误差可看作模型不确定性，存在模型不确定性的情况下，基于最优滤波理论设计的卡尔曼滤波(KF)及其改进算法，实际上都是次优的。为了解决这一问题，提出一种基于非线性H∞滤波的脉冲星导航方法，该方法能够消弱星表误差的影响，提高X射线脉冲星导航精度。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提出一种基于非线性H∞滤波的航天器自主定位方法，将有关星表误差的先验知识用于优化滤波增益阵的设计，使包含在系统模型和观测量中的信息得到充分合理的运用，从而抑制星表误差对估计精度的不利影响，提高地球轨道航天器自主导航精度。

本发明的技术解决方案是：

一种基于非线性H∞滤波的航天器自主定位方法，步骤如下：

(1)选择参与导航的地球轨道航天器的位置矢量和速度矢量在地心惯性坐标系的分量作为状态变量；

(2)在所述地球轨道航天器上配置X射线探测器，通过X射线探测器测量脉冲星辐射的脉冲信号，获取脉冲到达时间观测量；

(3)利用非线性H∞滤波算法处理步骤(2)中得到的脉冲到达时间观测量，获得状态变量的估计值，即地球轨道航天器的位置矢量和速度矢量的估计值，从而实现了地球轨道航天器的自主定位。

所述步骤(1)中状态变量为：

$x_t={\left[\begin{array}{cc}{r}^T& {\stackrel{.}{r}}^T\end{array}\right]}^T$

其中，

r＝[r_x r_y r_z]^T

$\stackrel{.}{r}={\left[\begin{array}{ccc}{v}_x& v_y& v_z\end{array}\right]}^T$

r和表示地球轨道航天器的位置矢量和速度矢量，r_x,r_y,r_z表示位置矢量的三个分量，v_x,v_y,v_z表示速度矢量的三个分量。

所述步骤(2)中脉冲到达时间观测量为：

y_t＝h(x_t)+v_t

其中，

$h\left(x_t\right)=t_b^{\left(p\right)}-\frac{1}{c}{n}^{\left(p\right)}\·r_{\mathrm{obs}}+\frac{1}{2c{D}_0^{\left(p\right)}}[r_{\mathrm{obs}}^2-{(n^{\left(p\right)}\·r_{\mathrm{obs}})}^2]-\frac{2{\mathrm{\μ}}_s}{c^3}\ln |\frac{n^{\left(p\right)}\·r_{\mathrm{obs}}+r_{\mathrm{obs}}}{n^{\left(p\right)}\·b+b}+1|$

y_t表示脉冲到达时间观测量，h(x_t)为中间变量，v_t表示测量噪声，表示脉冲信号到达太阳系质心SSB的时间，根据脉冲星时间模型预测得到；c表示光速，n^(p)表示脉冲星视线矢量，表示SSB到脉冲星的距离，b是SSB相对于太阳质心的位置矢量，b为SSB相对于太阳质心的距离，μ_s是太阳引力常数，r_obs是地球轨道航天器相对于SSB的位置矢量，其计算公式为

r_obs＝r_e+r

其中，r_e表示地球中心相对于SSB的位置矢量。

所述步骤(3)中通过非线性H∞滤波算法处理脉冲到达时间观测量，获得状态变量的估计值具体为：

${\hat{x}}_{t|t-1}={\hat{x}}_{t-1}+f\left({\hat{x}}_{t-1}\right)T$

${\hat{x}}_t={\hat{x}}_{t|t-1}+K_t[y_t-h\left({\hat{x}}_{t|t-1}\right)]$

其中，和分别为t时刻状态变量的估计值和预测值，为t-1时刻状态变量的估计值，T为滤波周期，为已知的状态转移函数，K_t为非线性H∞滤波算法的增益阵，且有：

$K_t={(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T{[H_t{(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T+D_t^T{\hat{x}}_{t|t-1}{\hat{x}}_{t|t-1}^T{D}_t+R_t]}^{-1};$

其中，I为单位矩阵，P_t|t-1为估计误差方差阵，通过所述非线性H∞滤波算法递推计算得到；γ为可调滤波参数且γ＝(1×10⁴)²；H_t为观测方程的雅克比矩阵，其计算公式为

$H_t=\frac{\∂h\left(x_t\right)}{\∂x_t}{|}_{x_t={\hat{x}}_{t|t-1}}$

D_t是用于描述星表误差影响的不确定性矩阵，其计算公式为

$D_t=\left[\begin{array}{c}-\cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ 0\end{array}\right]b_{\mathrm{\α}}^{\left(p\right)}+\left[\begin{array}{c}-\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ -\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\end{array}\right]b_{\mathrm{\δ}}^{\left(p\right)}$

其中，和分别表示脉冲星星表中带有误差的赤经和赤纬，和分别表示星表中赤经和赤纬的误差界；R_t表示测量噪声方差阵，为正定阵。

本发明与现有技术相比的有益效果是：

传统扩展卡尔曼滤波(EKF)算法的估计精度会受到脉冲星星表误差的不利影响，造成导航精度下降。本发明将星表误差视为模型不确定性的一部分，采用非线性H∞滤波(NRF)技术进行处理，通过优化滤波增益阵的设计，所提算法能够增强X射线脉冲星导航系统克服星表误差影响的能力，获得优于传统扩展卡尔曼滤波的导航精度。此外，本发明所设计的基于非线性H∞滤波的脉冲星导航算法与传统扩展卡尔曼滤波算法在实现形式上的区别仅在于滤波增益阵的计算方式不同，不会显著增大算法的计算量，具有较强的实用性。

附图说明

图1为本发明流程图；

图2为扩展卡尔曼滤波的估计误差曲线；

图3为非线性H∞滤波的估计误差曲线。

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

在脉冲星导航系统中，脉冲星星表误差体现为测量模型的不确定性。在存在模型不确定的系统中，如何设计滤波器，使其估计结果满足一定的设计指标要求，这属于H∞滤波的研究范畴。与传统卡尔曼滤波等最优滤波算法不同，H∞滤波在设计阶段就考虑了模型不确定性的影响，并将有关模型不确定性的先验知识用于优化滤波增益阵的设计，使包含在系统模型和观测量中的信息得到合理的运用，从而抑制模型不确定性对估计精度的不利影响。基于非线性H∞滤波的脉冲星导航方法流程可归纳如下：选择航天器的位置矢量和速度矢量在地心惯性系的投影作为状态变量，将航天器轨道动力学方程作为系统模型，将脉冲星时间转换方程作为测量模型，将脉冲星星表误差描述为模型中的有界误差，基于模型设计用于状态估计的非线性H∞滤波算法；进而，采用所设计的非线性H∞滤波算法，处理脉冲到达时间观测量序列，通过递推计算估计出航天器的位置和速度。

本发明提出一种基于非线性H∞滤波的航天器自主定位方法，如图1所示，步骤如下：

(1)选择参与导航的地球轨道航天器的位置矢量和速度矢量在地心惯性坐标系的分量作为状态变量。状态变量的数学表达式为：

$x_t={\left[\begin{array}{cc}{r}^T& {\stackrel{.}{r}}^T\end{array}\right]}^T$

其中，

r＝[r_x r_y r_z]^T

$\stackrel{.}{r}={\left[\begin{array}{ccc}{v}_x& v_y& v_z\end{array}\right]}^T$

r和表示地球轨道航天器的位置矢量和速度矢量，r_x,r_y,r_z表示位置矢量的三个分量，v_x,v_y,v_z表示速度矢量的三个分量。

(2)在所述地球轨道航天器上配置X射线探测器，通过X射线探测器获取脉冲到达时间观测量；目前在国内，北京控制工程研究所等单位均研制了用于脉冲到达时间测量的X射线探测器。脉冲到达时间观测量的数学表达式为：

y_t＝h(x_t)+v_t

其中，

$h\left(x_t\right)=t_b^{\left(p\right)}-\frac{1}{c}{n}^{\left(p\right)}\·r_{\mathrm{obs}}+\frac{1}{2c{D}_0^{\left(p\right)}}[r_{\mathrm{obs}}^2-{(n^{\left(p\right)}\·r_{\mathrm{obs}})}^2]-\frac{2{\mathrm{\μ}}_s}{c^3}\ln |\frac{n^{\left(p\right)}\·r_{\mathrm{obs}}+r_{\mathrm{obs}}}{n^{\left(p\right)}\·b+b}+1|$

y_t表示脉冲到达时间观测量，v_t表示测量噪声，表示脉冲信号到达太阳系质心(SSB)的时间，可根据脉冲星时间模型精确的预测得到，c表示光速，n^(p)表示脉冲星视线矢量，表示SSB到脉冲星的距离，b是SSB相对于太阳质心的位置矢量，μ_s是太阳引力常数，r_obs是航天器相对于SSB的位置矢量，其计算公式为

r_obs＝r_e+r

其中，r_e表示地球中心相对于SSB的位置矢量。上述两式中的c、n^(p)、b、μ_s和r_e均为已知量。观测方程的推导过程可参见中国宇航出版社2009年出版的由帅平、李明、陈绍龙、黄震编著的《X射线脉冲星导航系统原理与方法》一书。

(3)利用非线性H∞滤波算法处理步骤(2)中得到的脉冲到达时间观测量，获得状态变量的估计值，即地球轨道航天器的位置矢量和速度矢量的估计值，从而实现了地球轨道航天器的自主导航。

扩展卡尔曼滤波算法是通过如下公式进行：

${\hat{x}}_{t|t-1}={\hat{x}}_{t-1}+f\left({\hat{x}}_{t-1}\right)T$

${\hat{x}}_t={\hat{x}}_{t|t-1}+K_t[y_t-h\left({\hat{x}}_{t|t-1}\right)]$

其中，和分别为t时刻状态变量的估计值和预测值，为t-1时刻状态变量的估计值，T为滤波周期，可取为1s，为已知的状态转移函数，K_t为非线性H∞滤波算法的增益阵。

考虑地球的非球形引力摄动，状态转移函数如下所示：

$f\left(x_t\right)=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -\frac{\mathrm{\μ}{r}_x}{r^3}[1+\frac{3}{2}{J}_2{\left(\frac{R_e}{r}\right)}^2(1-5\frac{r_z^2}{r^2})]\\ -\frac{\mathrm{\μ}{r}_y}{r^3}[1+\frac{3}{2}{J}_2{\left(\frac{R_e}{r}\right)}^2(1-5\frac{r_z^2}{r^2})]\\ -\frac{\mathrm{\μ}{r}_z}{r^3}[1+\frac{3}{2}{J}_2{\left(\frac{R_e}{r}\right)}^2(3-5\frac{r_z^2}{r^2})]\end{array}\right]$

其中，μ为地球引力常数，J₂为二阶带谐项系数，R_e为地球平均赤道半径，表示地心距。参数μ、J₂和R_e均为已知量。状态转移函数的推导过程可参见北京航空航天大学出版社1998年出版的由章仁为编著的《卫星轨道姿态动力学与控制》一书。

非线性H∞滤波算法的增益阵计算公式如下所示：

$K_t={(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T{[H_t{(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T+D_t^T{\hat{x}}_{t|t-1}{\hat{x}}_{t|t-1}^T{D}_t+R_t]}^{-1}$

其中，P_t|t-1为估计误差方差阵，可通过滤波算法的递推计算得到，估计误差方差阵的递推计算过程可参考西北工业大学出版社1998出版的由秦永元、张洪钺、汪叔华编写的《卡尔曼滤波与组合导航原理》一书；γ为可调滤波参数，通常取为正常数，对于本实施例，可选择γ＝(1×10⁴)²；H_t为观测方程的雅克比矩阵，其计算公式为

$H_t=\frac{\∂h\left(x_t\right)}{\∂x_t}{|}_{x_t={\hat{x}}_{t|t-1}}$

D_t是用于描述星表误差影响的不确定性矩阵，其计算公式为

$D_t=\left[\begin{array}{c}-\cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ 0\end{array}\right]b_{\mathrm{\α}}^{\left(p\right)}+\left[\begin{array}{c}-\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ -\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\end{array}\right]b_{\mathrm{\δ}}^{\left(p\right)}$

其中，和分别表示脉冲星星表中带有误差的赤经和赤纬，和分别表示星表中赤经和赤纬的误差界，均为已知量；R_t表示测量噪声方差阵，为已知正定阵。

下面，以地球轨道航天器自主导航为例，通过仿真实例验证本发明所述方法的有效性。设航天器在半长轴12275km，偏心率0.0038，轨道倾角109.8°的轨道上环绕地球运动。假定X射线探测器分时段观测3颗脉冲星B0531+21、B1821-24和B1937+21，仿真时间为地球卫星的3个轨道周期，轨道外推数据每1s更新一次。地球轨道航天器的初始位置误差设为1km，初始速度误差设为1m/s。

仿真过程中，分别采用扩展卡尔曼滤波(EKF)算法和非线性H∞滤波(NRF)算法处理脉冲到达时间观测量，估计地球轨道航天器的位置矢量和速度矢量，EKF的位置估计误差曲线如图2所示，NRF的位置估计误差曲线如图3所示。图2和图3中实线均表示位置估计误差曲线，虚线均表示通过估计误差方差阵递推得到的期望定位精度，每幅图均包括从上到下三幅图，分别对应卫星位置矢量的3个分量。纵坐标表示位置估计误差的大小，单位为m，横坐标表示时间，单位为轨道周期。根据图2和图3不难看出，存在星表误差的情况下，采用非线性H∞滤波算法进行航天器自主定位导航，能够获得优于传统扩展卡尔曼滤波的导航精度。通过统计计算可知，EKF的定位精度约为1211.1m，NRF的定位精度约为559.8m。通过对比可知，本发明所述方法精度更高。

仿真结果表明，本发明提出的基于非线性H∞滤波的脉冲星定位方法能够实现地球轨道航天器高精度自主导航。本发明的主要技术内容可用于实现具有较高精度需求的卫星自主导航任务，有助于降低卫星对地面测控的依赖程度，增强卫星系统在紧急情况下的自主生存能力。

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

Claims (4)

1.一种基于非线性H∞滤波的航天器自主定位方法，其特征在于步骤如下：

(1)选择参与导航的地球轨道航天器的位置矢量和速度矢量在地心惯性坐标系的分量作为状态变量；

(2)在所述地球轨道航天器上配置X射线探测器，通过X射线探测器测量脉冲星辐射的脉冲信号，获取脉冲到达时间观测量；

(3)利用非线性H∞滤波算法处理步骤(2)中得到的脉冲到达时间观测量，获得状态变量的估计值，即地球轨道航天器的位置矢量和速度矢量的估计值，从而实现了地球轨道航天器的自主定位。

2.根据权利要求1所述的一种基于非线性H∞滤波的航天器自主定位方法，其特征在于：所述步骤(1)中状态变量为：

$x_t={\left[\begin{array}{cc}{r}^T& {\stackrel{.}{r}}^T\end{array}\right]}^T$

其中，

r＝[r_x r_y r_z]^T

$\stackrel{.}{r}={\left[\begin{array}{ccc}{v}_x& v_y& v_z\end{array}\right]}^T$

r和表示地球轨道航天器的位置矢量和速度矢量，r_x,r_y,r_z表示位置矢量的三个分量，v_x,v_y,v_z表示速度矢量的三个分量。

3.根据权利要求1所述的一种基于非线性H∞滤波的航天器自主定位方法，其特征在于：所述步骤(2)中脉冲到达时间观测量为：

y_t＝h(x_t)+v_t

其中，

$h\left(x_t\right)=t_b^{\left(p\right)}-\frac{1}{c}{n}^{\left(p\right)}\·r_{\mathrm{obs}}+\frac{1}{2c{D}_0^{\left(p\right)}}[r_{\mathrm{obs}}^2-{(n^{\left(p\right)}\·r_{\mathrm{obs}})}^2]-\frac{2{\mathrm{\μ}}_s}{c^3}\ln |\frac{n^{\left(p\right)}\·r_{\mathrm{obs}}+r_{\mathrm{obs}}}{n^{\left(p\right)}\·b+b}+1|$

y_t表示脉冲到达时间观测量，h(x_t)为中间变量，v_t表示测量噪声，表示脉冲信号到达太阳系质心SSB的时间，根据脉冲星时间模型预测得到；c表示光速，n^(p)表示脉冲星视线矢量，表示SSB到脉冲星的距离，b是SSB相对于太阳质心的位置矢量，b为SSB相对于太阳质心的距离，μ_s是太阳引力常数，r_obs是地球轨道航天器相对于SSB的位置矢量，其计算公式为

r_obs＝r_e+r

其中，r_e表示地球中心相对于SSB的位置矢量。

4.根据权利要求1所述的一种基于非线性H∞滤波的航天器自主定位方法，其特征在于：所述步骤(3)中通过非线性H∞滤波算法处理脉冲到达时间观测量，获得状态变量的估计值具体为：

${\hat{x}}_{t|t-1}={\hat{x}}_{t-1}+f\left({\hat{x}}_{t-1}\right)T$

${\hat{x}}_t={\hat{x}}_{t|t-1}+K_t[y_t-h\left({\hat{x}}_{t|t-1}\right)]$

其中，和分别为t时刻状态变量的估计值和预测值，为t-1时刻状态变量的估计值，T为滤波周期，为已知的状态转移函数，K_t为非线性H∞滤波算法的增益阵，且有：

$K_t={(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T{[H_t{(P_{t|t-1}^{-1}-{\mathrm{\γ}}^{-2}I)}^{-1}{H}_t^T+D_t^T{\hat{x}}_{t|t-1}{\hat{x}}_{t|t-1}^T{D}_t+R_t]}^{-1};$

其中，I为单位矩阵，P_t|t-1为估计误差方差阵，通过所述非线性H∞滤波算法递推计算得到；γ为可调滤波参数且γ＝(1×10⁴)²；H_t为观测方程的雅克比矩阵，其计算公式为

$H_t=\frac{\∂h\left(x_t\right)}{\∂x_t}{|}_{x_t={\hat{x}}_{t|t-1}}$

D_t是用于描述星表误差影响的不确定性矩阵，其计算公式为

$D_t=\left[\begin{array}{c}-\cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ 0\end{array}\right]b_{\mathrm{\α}}^{\left(p\right)}+\left[\begin{array}{c}-\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\cos \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ -\sin \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\sin \left({\widehat{\mathrm{\α}}}^{\left(p\right)}\right)\\ \cos \left({\widehat{\mathrm{\δ}}}^{\left(p\right)}\right)\end{array}\right]b_{\mathrm{\δ}}^{\left(p\right)}$

其中，和分别表示脉冲星星表中带有误差的赤经和赤纬，和分别表示星表中赤经和赤纬的误差界；R_t表示测量噪声方差阵，为正定阵。