CN104048664A - 一种导航卫星星座自主定轨的方法 - Google Patents

Abstract

本发明属于卫星自主定轨领域，涉及一种导航卫星星座自主定轨的方法。其包括如下步骤：(1)：自主定轨开始，***初始化；(2)：获取卫星之间的星间测距；(3)：计算星间测距的观测矩阵；(4)：利用轨道动力学模型进行轨道预报；(5)：利用卡尔曼滤波算法分别进行量测更新和状态估计；(6)：判断自主定轨是否结束，若未结束，则再次运行步骤(2)～步骤(6)；反之，则结束退出自主定轨程序。本发明相对于传统的导航卫星星座自主定轨的方法，具有以下优点：1、可以解决导航星座仅利用星间测距进行自主定轨时的秩亏问题；2、长期自主定轨精度较高。

Description

一种导航卫星星座自主定轨的方法

技术领域

本发明属于卫星自主定轨领域，具体涉及一种导航卫星星座自主定轨的方法。

背景技术

自主定轨能力对卫星导航***的发展将具有越来越重要的意义。拥有自主定轨能力的卫星导航***不需要通过全球布站来提供支持，同时可以减轻有限的地面站的负担，提高卫星导航***的安全性。

卫星导航***的自主定轨既可以利用单星自主定轨方法，通过***中每颗卫星独立的自主定轨实现整个导航***的自主定轨。也可以利用卫星之间的相对测量实现卫星导航***的整网定轨。相对测量包括导航***中卫星间的相对角度、相对距离和相对速度等信息。相对于卫星与一个表面不规则的天体之间的测量，两颗卫星之间的相对测量精度要高的多。因此利用导航卫星间的相对测量信息进行自主定轨是一个可行的方式。早期研究指出只需要利用星间测角信息就可以进行自主定轨，仅用测角的自主定轨精度为几百米。但联合使用星间测距和测角的自主定轨精度为米级。这说明在利用星间测量信息进行自主定轨时，星间测距比测角的作用更大。然而仅利用星间测距对近地卫星导航***进行自主定轨存在秩亏问题。

发明内容

发明目的：为解决上述技术问题，本发明提供了一种导航卫星星座自主定轨的方法。该方法中的需要至少有一颗运行在地月系拉格朗日轨道上的导航卫星。

技术方案：一种导航卫星星座自主定轨的方法，包括如下步骤：

(1)：自主定轨开始，***初始化；

(2)：获取卫星之间的星间测距；

(3)：计算星间测距的观测矩阵；

(4)：利用轨道动力学模型进行轨道预报；

(5)：利用卡尔曼滤波算法分别进行量测更新和状态估计；

(6)：判断自主定轨是否结束，若未结束，则再次运行步骤(2)～步骤(6)；反之，则结束退出自主定轨程序。

作为本发明中导航卫星星座自主定轨的方法的优选方案：步骤(2)中所述测距包括导航卫星之间的测距和导航卫星与拉格朗日卫星之间的测距。

作为本发明中导航卫星星座自主定轨的方法的优选方案：若步骤(2)中的测距为导航卫星之间的测距矩阵，则按(1)式计算导航卫星之间的观测矩阵：

$H_{k+1}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂[{\left({\mathrm{\σ}}_i\right)}_{k+1}^T,{\left({\mathrm{\σ}}_j\right)}_{k+1}^T]}{|}_{x_{k+1}={\hat{x}}_{k+1/k}}---\left(1\right)$

ρ_ij为星间测距，σ_i,σ_j为两颗导航卫星的轨道根数；

(1)式中

其中，

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_i}^T}=\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\\ \frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_j}^T}=-\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\end{array}\right.---\left(3\right)$

[x_Ei  y_Ei  z_Ei]^T和[x_Ej  y_Ej  z_Ej]^T分别为近地卫星i和近地卫星j在惯性系下的位置坐标。

作为本发明中导航卫星星座自主定轨的方法的优选方案：若步骤(2)中的测距矩阵为导航卫星与拉格朗日卫星之间的测距，则按(4)式计算导航卫星与拉格朗日卫星之间的观测矩阵：

其中，

且

${\left[\begin{array}{ccc}{x}_{{\mathrm{EL}}_k}& y_{{\mathrm{EL}}_k}& z_{{\mathrm{EL}}_k}\end{array}\right]}^T$ 为拉格朗日轨道导航卫星在地心惯性系下的位置向量。

作为本发明中导航卫星星座自主定轨的方法的优选方案：步骤(4)中所述利用轨道动力学模型进行轨道预报包括以下步骤：

(41)将前一时刻的定轨结果作为轨道预报的初值；

(42)将初值带入卫星的轨道动力学模型，其中轨道动力学模型考虑了大气阻力、太阳光压、地球非球形引力、第三体引力、潮汐力的影响；

(43)利用RKF78数值积分方法得到轨道预报的结果。

作为本发明中导航卫星星座自主定轨的方法的优选方案：步骤(5)中利用卡尔曼滤波算法分别进行量测更新和状态估计，包括如下步骤：

(51)按(7)式和(8)式对被估状态进行一步预测

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1/k}{X}_{k/k}+U_k=f({\hat{X}}_k,k)---\left(7\right)$

$P_{k+1/k}={\mathrm{\Φ}}_{k+1/k}{P}_{k/k}{\mathrm{\Φ}}_{k+1/k}^T+{\mathrm{\Γ}}_k{Q}_k{\mathrm{\Γ}}_k^T---\left(8\right)$

(52)将步骤(1)的结果作为先验信息，按(9)式计算增益矩阵K_k+1

$K_{k+1}=P_{k+1/k}{H}_{k+1}^T{(H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R_{k+1})}^{-1}---\left(9\right)$

(53)利用K_k+1，按(10)式得到新的状态估计值

${\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-Z^{*})---\left(10\right)$

(54)按(11)式计算新的协方差矩阵，为下一步滤波做准备

$\begin{array}{c}{P}_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T\\ +K_{k+1}{R}_{k+1}{K}_{k+1}^T\end{array}---\left(11\right).$

有益效果：本发明相对于传统的导航卫星星座自主定轨的方法，具有以下优点：

1、可以解决导航星座仅利用星间测距进行自主定轨时的秩亏问题；

2、长期自主定轨精度较高。

附图说明

图1为圆型限制性三体问题的示意图

图2为圆型限制性三体问题中五个拉格朗日点位置示意图

图3为本发明的方法流程图

图4对本发明进行仿真验证示意图

图5为仿真实验中GPS导航卫星PRN9的位置误差和用户伪距误差

图6为仿真实验中GPS导航卫星PRN30的位置误差和用户伪距误差

图7为仿真实验中GPS导航卫星PRN31的位置误差和用户伪距误差

具体实施方式：

为了更好地理解本发明相对于现有技术所作出的改进，在对本发明的具体实施方式进行详细说明之前，先对发明内容中的地月系拉格朗日轨道进一步补充。

对于运行在拉格朗日轨道上的导航卫星，可以描述其动力学特征的最简单的模型为如图1所示圆型限制性三体问题模型——设两个大天体m₁和m₂仅在相互引力作用下运动，且它们之间相互绕行的轨道为半径为r₁₂的圆。设一非惯性坐标系xyz(会合坐标系)，坐标原点为两个大天体的质心，x轴由m₁指向m₂，z轴垂直于m₂绕m₁的运行轨道平面，y轴与x、z轴满足右手定则。在此坐标系中，m₁和m₂看起来静止不动，现在引入质量为m的第三个物体，其质量与m₁和m₂相比可以忽略不计，m在m₁和m₂引力场作用下的运动就是所谓的圆型限制性三体问题。

在会合坐标系中有五个拉格朗日点，如图2所示，其中三个共线拉格朗日点是不稳定的，两个三角拉格朗日点是稳定的，不管是稳定的拉格朗日点还是不稳定的拉格朗日点，其附近都存在周期轨道，本申请文件中的导航卫星星座正是分布在拉格朗日点周期轨道上的。

如图3所示，一种导航卫星星座自主定轨的方法，包括如下步骤：

(1)：自主定轨开始，***初始化；

(2)：获取卫星之间的星间测距；

(3)：计算星间测距的观测矩阵；

(4)：利用轨道动力学模型进行轨道预报；

(5)：利用卡尔曼滤波算法分别进行量测更新和状态估计；

(6)：判断自主定轨是否结束，若未结束，则再次运行步骤(2)～步骤(6)；反之，则结束退出自主定轨程序。

步骤(2)中所述测距包括导航卫星之间的测距和导航卫星与拉格朗日卫星之间的测距。

若步骤(2)中的测距为导航卫星之间的测距，则按(1)式计算导航卫星之间的观测矩阵：

$H_{k+1}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂[{\left({\mathrm{\σ}}_i\right)}_{k+1}^T,{\left({\mathrm{\σ}}_j\right)}_{k+1}^T]}{|}_{x_{k+1}={\hat{x}}_{k+1/k}}---\left(1\right)$

ρ_ij为星间测距，σ_i,σ_j为两颗导航卫星的轨道根数；

(1)式中

其中，

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_i}^T}=\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\\ \frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_j}^T}=-\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\end{array}\right.---\left(3\right)$

[x_Ei  y_Ei  z_Ei]^T和[x_Ej  y_Ej  z_Ej]^T分别为近地卫星i和近地卫星j在惯性系下的位置坐标。

若步骤(2)中的测距为导航卫星与拉格朗日卫星之间的测距，

则按(4)式计算导航卫星与拉格朗日卫星之间的观测矩阵：

其中，

且

${\left[\begin{array}{ccc}{x}_{{\mathrm{EL}}_k}& y_{{\mathrm{EL}}_k}& z_{{\mathrm{EL}}_k}\end{array}\right]}^T$ 为拉格朗日轨道导航卫星在地心惯性系下的位置向量。

步骤(4)中所述利用轨道动力学模型进行轨道预报包括以下步骤：

(41)将前一时刻的定轨结果作为轨道预报的初值；

(42)将初值带入卫星的轨道动力学模型，其中轨道动力学模型考虑了大气阻力、太阳光压、地球非球形引力、第三体引力、潮汐力的影响；

(43)利用RKF78数值积分方法得到轨道预报的结果。

步骤(5)中利用卡尔曼滤波算法分别进行量测更新和状态估计，包括如下步骤：

(51)按(7)式和(8)式对被估状态进行一步预测

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1/k}{X}_{k/k}+U_k=f({\hat{X}}_k,k)---\left(7\right)$

(52)将步骤(1)的结果作为先验信息，按(9)式计算增益矩阵K_k+1

$K_{k+1}=P_{k+1/k}{H}_{k+1}^T{(H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R_{k+1})}^{-1}---\left(9\right)$

(53)利用K_k+1，按(10)式得到新的状态估计值

${\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-Z^{*})---\left(10\right)$

(54)按(11)式计算新的协方差矩阵，为下一步滤波做准备

$\begin{array}{c}{P}_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T\\ +K_{k+1}{R}_{k+1}{K}_{k+1}^T\end{array}---\left(11\right).$

仿真验证方案

如图4所示，为了验证本发明的方法的正确性，采用仿真验证的方法进行验证，步骤如下：

(1)、利用GPS卫星的精密星历产生导航卫星之间的星间测距；

(2)、利用拉格朗日卫星的精密星历结合GPS卫星的精密星历产生拉格朗日卫星与GPS卫星之间的星间测距，将两类测距中加入适当的测距误差；

(3)、利用卡尔曼滤波方法结合两类测距进行自主定轨；

(4)、将定轨结果与精密星历比较，对定轨精度进行评估。

仿真起始时刻选为2000/4/13，23:59:47(UTC)，仿真时长180天，地球引力场模型采用10×10阶WGS-84模型，第三体(日月)引力摄动采用DE200计算日月星历，太阳光压摄动采用新的RPR模型。星间测距的观测间隔设为1小时，假设近地导航卫星间测距的***误差为0.1m，随机差0.1m，近地导航卫星与拉格朗日导航卫星间测距***误差0.5m，随机差0.5m；拉个朗日导航卫星之间测距***差1m，随机差1m。LU为地月距离。

表1：拉格朗日导航卫星的初始状态

图5～图7给出了拉格朗日卫星星间测距误差1m，初始位置误差0.1m条件下，利用扩展kalman滤波算法对拉格朗日点轨道上的卫星与GPS卫星组成的导航星座仅利用星间测距进行自主定轨的仿真结果。三个坐标轴上的最大误差见下表：

表2：近地导航卫星的RNT坐标最大定轨误差

表3:拉格朗日导航卫星的最大定轨误差

Claims (6)

1.一种导航卫星星座自主定轨的方法，包括如下步骤：

(1)、自主定轨开始，***初始化；

(2)、获取卫星之间的星间测距；

(3)、计算星间测距的观测矩阵；

(4)、利用轨道动力学模型进行轨道预报；

(5)、利用卡尔曼滤波算法分别进行量测更新和状态估计；

(6)、判断自主定轨是否结束，若未结束，则再次运行步骤(2)～步骤(6)；反之，则结束退出自主定轨程序。

2.如权利要求1中所述的导航卫星星座自主定规的方法，其特征在于，步骤(2)中所述测距包括导航卫星之间的测距和导航卫星与拉格朗日卫星之间的测距。

3.如权利要求1中所述的导航卫星星座自主定轨的方法，其特征在于，若步骤(2)中的测距为导航卫星之间的测距，则按(1)式计算导航卫星之间的观测矩阵：

$H_{k+1}=\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂[{\left({\mathrm{\σ}}_i\right)}_{k+1}^T,{\left({\mathrm{\σ}}_j\right)}_{k+1}^T]}{|}_{x_{k+1}={\hat{x}}_{k+1/k}}---\left(1\right)$

ρ_ij为星间测距，σ_i,σ_j为两颗导航卫星的轨道根数；

(1)式中

其中，

$\left\{\begin{array}{c}\frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_i}^T}=\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\\ \frac{\∂{\mathrm{\ρ}}_{\mathrm{ij}}}{\∂{r_j}^T}=-\frac{1}{{\mathrm{\ρ}}_{\mathrm{ij}}}\left[\begin{array}{ccc}{x}_{\mathrm{Ei}}-x_{\mathrm{Ej}}& y_{\mathrm{Ei}}-y_{\mathrm{Ej}}& z_{\mathrm{Ei}}-z_{\mathrm{Ej}}\end{array}\right]\end{array}\right.---\left(3\right)$ [x_Ei  y_Ei  z_Ei]^T和[x_Ej  y_Ej  z_Ej]^T分别为近地卫星i和近地卫星j在惯性系下的位置坐标。

4.如权利要求1中所述的导航卫星星座自主定轨的方法，其特征在于，若步骤(2)中的测距为导航卫星与拉格朗日卫星之间的测距，则按(4)式计算导航卫星与拉格朗日卫星之间的观测矩阵：

其中，

且

${\left[\begin{array}{ccc}{x}_{{\mathrm{EL}}_k}& y_{{\mathrm{EL}}_k}& z_{{\mathrm{EL}}_k}\end{array}\right]}^T$ 为拉格朗日轨道导航卫星在地心惯性系下的位置向量。

5.如权利要求1中所述的导航卫星星座自主定轨的方法，其特征在于，步骤(4)中所述利用轨道动力学模型进行轨道预报包括以下步骤：

(41)将前一时刻的定轨结果作为轨道预报的初值；

(42)将初值带入卫星的动力学模型，其中动力学模型的考虑了了大气阻力、太阳光压、地球非球形引力、第三体引力、潮汐力的影响；

(43)利用RKF78数值积分方法得到轨道预报的结果。

6.如权利要求1中所述的导航卫星星座自主定轨的方法，其特征在于，步骤(5)中利用卡尔曼滤波算法分别进行量测更新和状态估计，包括如下步骤：

(51)按(7)式和(8)式对被估状态进行一步预测

${\hat{X}}_{k+1/k}={\mathrm{\Φ}}_{k+1/k}{X}_{k/k}+U_k=f({\hat{X}}_k,k)---\left(7\right)$

$P_{k+1/k}={\mathrm{\Φ}}_{k+1/k}{P}_{k/k}{\mathrm{\Φ}}_{k+1/k}^T+{\mathrm{\Γ}}_k{Q}_k{\mathrm{\Γ}}_k^T---\left(8\right)$

(52)将步骤(1)的结果作为先验信息，按(9)式计算增益矩阵K_k+1

$K_{k+1}=P_{k+1/k}{H}_{k+1}^T{(H_{k+1}{P}_{k+1/k}{H}_{k+1}^T+R_{k+1})}^{-1}---\left(9\right)$

(53)利用K_k+1，按(10)式得到新的状态估计值

${\hat{X}}_{k+1/k+1}={\hat{X}}_{k+1/k}+K_{k+1}(Z_{k+1}-Z^{*})---\left(10\right)$

(54)按(11)式计算新的协方差矩阵，为下一步滤波做准备

$\begin{array}{c}{P}_{k+1/k+1}=(I-K_{k+1}{H}_{k+1})P_{k+1/k}{(I-K_{k+1}{H}_{k+1})}^T\\ +K_{k+1}{R}_{k+1}{K}_{k+1}^T\end{array}---\left(11\right).$