CN103542854B

CN103542854B - 基于星载处理器的自主定轨方法

Info

Publication number: CN103542854B
Application number: CN201310531747.5A
Authority: CN
Inventors: 冯旭哲; 杨俊�; 陈建云; 黄文德; 胡梅; 李鑫
Original assignee: National University of Defense Technology
Current assignee: National University of Defense Technology
Priority date: 2013-11-02
Filing date: 2013-11-02
Publication date: 2016-03-23
Anticipated expiration: 2033-11-02
Also published as: CN103542854A

Abstract

本发明提供了一种基于星载处理器的自主定轨方法。该方法利用星间测距，克服导航卫星定轨在脱离地面站支持的条件下定轨精度发散的难题，完成导航卫星之间精密测量数据转换、不同时间测量值归算、几何定轨、动力学定轨、自主星历生成等多项步骤，实现脱离地面站数据上注支持的导航卫星自主定轨算法。本方法在脱离地面站支持的条件下，即没有地面站测量数据，也没有地面站上注数据，导航卫星仅仅凭借星际双向测量的结果，并结合导航卫星上预存的60天长期预报星历，通过导航卫星自主运算，得到导航卫星的实时轨道预报结果，并可生成广播星历下发。

Description

基于星载处理器的自主定轨方法

技术领域

本发明主要涉及一种基于星载处理器的自主定轨方法。

背景技术

全球卫星导航系统（GlobalNavigationSatelliteSystem，GNSS）能够为地球及近地空间的任意地点提供全天候的精密位置和时间信息。全球卫星导航系统可支持各种需要精确位置与时间信息的战术操作，与通信、计算机和情报监视等其他系统构成多兵种协同作战指挥系统，已成为武装力量的支撑系统和战斗力倍增器，是我国基于信息系统体系作战能力形成过程中必不可少的核心关键系统之一。除此之外，卫星导航系统已广泛应用于国民经济各个领域，不仅为经济发展提供了强大的动力，且已成为关系国计民生的重要基础设施。

传统的导航卫星定轨方法为：多个已知精确坐标的地面站对导航卫星进行大量、长时间的测量，测得的数据在地面站计算并最终得到导航卫星轨道数据。导航卫星提供给地面用户的广播星历，需要地面站通过星历上注的方式获取。如果失去地面站支持，则导航卫星无法得到自身轨道以及提供给地面用户的广播星历。即使导航卫星通过预存保留了一段时间的预报星历，在地面站无法定期上注更新参数的情况下，预报星历因为精度发散而无法使用。

基于星间测距的导航卫星自主定轨方法，主要是导航卫星在脱离地面站支持的前提下，仅依靠星间精密测量的数据，对卫星轨道进行定位，并生成广播星历下发，使导航卫星在脱离地面站支持的情况下仍具备一定的功效，保持了导航能力。

发明内容

本发明要解决的技术问题是克服导航卫星定轨在脱离地面站支持的条件下定轨精度发散的难题，提供了一种基于星间测距的导航卫星自主定轨方法，完成导航卫星之间精密测量数据转换、不同时间测量值归算、几何定轨、动力学定轨、自主星历生成等多项步骤，实现了脱离地面站数据上注支持的导航卫星自主定轨算法。

本发明的技术方案是：通过星间（星地）精密测量得到星间（星地）精密测量值，即距离观测量，对星间距离观测量进行时间归化，得到同一时刻的星间测距观测量；经过一段时间的累积，得到卫星的几何定位值；在卫星几何定位的基础上，进行卫星动力学定轨；轨道动力学外推然后输出预报星历。具体步骤如下：

下列步骤中的卫星编号并非人为指定，而仅仅为叙述方便，卫星编号和卫星之间的对应关系并不影响本发明的技术方案。

步骤1：星间双向观测量生成

星间观测量并非连续获得，而是根据星间链路的节拍分时获取，假设节拍周期为，起始时刻为，则在2时间内完成一次双向测量。具体为：时刻A星发B星收，+时刻A星收B星发；+2时刻A星发C星收，+3时刻A星收C星发；以此类推。在一个测量周期内的第二个节拍，B星在完成测量的同时，将第一个节拍得到的测量数据发送回A星，即A星在一个测量周期内得到了AB星的双向测量数据。

以下以AB两颗星为例详细描述得到星间观测量的过程。

（1）起始时刻为，A星发B星收，测得的伪距存储在B星。

（2）时刻+，A星收B星发，测得的伪距存储在A星。

（3）仍然在时刻+，B星将步骤1测得的伪距从通信通道发至A星。

（4）A星将伪距和配对，做好观测量归化计算准备。

步骤2：星间双向观测量归算

由于星间链路直接观测量为不同时刻的测量值，因此在用于自主定轨之前需要进行预处理，即将观测量归化到同一时刻。输入为本地卫星A与建链卫星B的多个周期接收时刻伪距观测值，输出值为将各周期接收时刻伪距归化到指定接收时刻对应的伪距；

具体做法为：

（1）首先确定本星需要建立双向测量的卫星数量，假设为，则一个测量周期为。

（2）确定不同测量值需要规化的时刻点，一般为某条测量链路起始时刻到测量周期的中点，即时刻；这时计算其它链路对应到此时刻点的值，由于测量链路为时分体制，因此对于不同的测量链路此时刻点不相同。

（3）每条测量链路分别进行归化计算。一般数据积累点为测量周期的10倍以上，确保测量周期以及归化时刻点的前后数据大体相当，此时总的数据点数量假设为。

（4）计算插值系数矩阵，并根据系数矩阵和归化时刻点，计算归化时刻点测量值。

步骤3：几何定位

（1）给定卫星初始位置

此时，对于不同的卫星来说，可以建立不同的观测量方程组。

（1）

对式(1)进行求解，即可得到卫星的位置坐标。但式(1)是一组含观测噪声的非线性方程，直接求解不易，一般需要进行线性化。假设通过某些途径（如测控系统提供的初始轨道，即初定轨结果），获得卫星的概略位置坐标，则只需对概略位置进行相应的改正即可得到卫星位置坐标，改正量用表示。此时得到的卫星位置坐标应为估计值（不是真值），记为，于是有

（2）

（2）对式（1）在卫星概略位置坐标处进行线性化；

对式（1）进行线性化，并忽略2阶以上高阶项，可得

（3）

其中，为地面站到卫星距离的近似值

（4）

为地面站到卫星的观测矢量的方向余弦，即

（5）

于是，令

，，，

则式（3）可写成矩阵形式

（6）

令为式（6）的解算原则（最小二乘原则），即

（7）

则式（7）的最优解（无偏估计）为

（8）

（1）当观测的卫星只有3个时，式（8）的解为

（9）

（2）当观测的卫星数量大于3个时，式（8）的解为

（10）

需要注意的是，由于概略位置可能具有较大的误差，从而导致式（3）存在较大的线性化误差，因此，一般需要迭代数次才能获得精度较高的定位解。另外，上述定位解算是在等精度观测的假设下得到的，对于不是等精度观测的情况，一般需要进行加权求解。相关推导与上述过程类似，此处不再赘述。直接对（10）式给出加权最小二乘解为

（11）

其中，为观测权矩阵。

（3）用上一步计算结果更新卫星位置；

（4）比较两次计算结果，判断是否满足收敛标准。如果满足则直接输出卫星位置估计值；否则用这次计算的卫星位置作为初始值，返回步骤（2）继续计算；

步骤4：卫星动力学定轨

（1）给定卫星初始位置

动力学定轨的输入值为几何定位输出的结果，这些几何定位的结果是一个个独立的点，彼此间没有有机的联系，动力学定轨即利用动力学模型，得到卫星的轨迹曲线。本发明采用的动力学离散模型为：

（12）

其中，、、、、、、分别为：当前历元轨道修正值、前一历元轨道修正值、模型噪声矩阵、观测值、观测矩阵、状态传递矩阵、观测误差矩阵。

应用扩展Kalman滤波（EKF）方法对递归问题（12）进行求解，其步骤如下：

（2）计算得到时刻的积分参考轨道和状态转移矩阵。

（3）从而得到预报状态

（13）

（4）预测状态协方差矩阵

（14）

其中，为模型噪声的协方差矩阵,为上一状态的轨道估计值。

（5）计算新息向量及其协方差矩阵

（15）

（16）

其中，为测量噪声的协方差矩阵，为新息向量，为新息向量的协方差矩阵。

（6）计算增益矩阵

（17）

（7）求解新的状态估计值

（18）

（8）更新状态协方差矩阵

（19）

从而可以对时刻的积分参考轨道进行修正

（20）

（9）如果满足设定门限则直接输出结果，否则返回步骤（1）继续计算。

步骤5：星历参数拟合

导航电文星历参数是提供给用户的关于卫星轨道的基准。设星历参考时刻，则相对于定轨时刻是未来时。因此需要将定轨时刻得到的轨道进行外推预报。

设动力学定轨得到的轨道为，为轨道历元，为轨道根数。这里的轨道外推涉及两个过程：一是将时刻的轨道预报到星历参考时刻，得到；二是为了计算星历参数所需的轨道外推，一般在时刻轨道的基础上外推1-2小时。

在时刻轨道的基础上外推数小时的目的是获得星历参数拟合所需的轨道观测值。本文采用卫星位置参数作为星历参数拟合的观测值，可以得到一系列的采样点观测值。因为星历参数一般在地固坐标系下给出，因此需要将观测值转换到地固坐标系下，即转换为。假设外推时长为4h，采样间隔为200s，则共有72个采样点观测值。最终输出导航卫星的16个基本广播星历参数：、、e、、、、、、、、、、、、、。由这16个参数可以计算卫星在地固坐标系中的位置坐标。

采用本发明可以达到以下技术效果：

在脱离地面站支持的条件下，即没有地面站测量数据，也没有地面站上注数据，导航卫星仅仅凭借星际双向测量的结果，并结合导航卫星上预存的60天长期预报星历，通过导航卫星自主运算，得到导航卫星的实时轨道预报结果，并可生成广播星历下发。

附图说明

图1是传统的基于地面站测量的导航卫星定轨方法示意图；

图2是本发明的基于星间测距的导航卫星自主定轨算法结构图；

图3是本发明的星间双向观测量生成结构图；

图4是本发明的星间双向观测量归化算法结构图；

图5是本发明的几何定位算法结构图；

图6是本发明的动力学定轨算法流程图；

图7是本发明的星历参数拟合生成算法流程图；

图8是归化计算长度和测量周期关系示意图。

具体实施方式

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

图1是传统的基于地面站测量的导航卫星定轨方法。

传统的导航卫星定轨方法基于地面站对导航卫星进行长期测量，根据测量结果，在地面进行计算，将数据拟合、外推，最终得到导航卫星预报轨道星历，并通过地面站上注给导航卫星，供其分发广播。

图2是本发明的基于星间测距的导航卫星自主定轨算法结构图，本算法由以下实施步骤组成。

步骤1：星间双向观测量生成。

步骤2：星间双向观测量归化计算。

步骤3：几何定位。

步骤4：动力学定轨

步骤5：星历参数拟合生成。

图3是本发明的星间双向观测量生成算法结构图

星间双向观测量通过星间链路的时空分体制得到，双单向建链测量的伪距值不仅含有距离量、钟差量，还包括电离层效应、对流层效应、收发设备延迟等，必须对其进行伪距误差改正；最后，再对经过时间同步改正和伪距误差改正后的量进行距离与钟差的解耦处理，得到距离和钟差观测量，分别用于定轨和时间同步处理。

图4是本发明的星间双向观测量归化算法结构图

由于星间双向观测量为不同时刻的观测量，要进行几何定位，必须进行时间归化，即将不同时刻的测量值归算到同一个时刻。

步骤1：根据建链数量确定计算周期，如建链数为N，单向测量节拍周期为，则测量周期为。

步骤2：确定不同星间观测量需要归化的目标时刻点，此时刻点一般在周期内，一般尽量选择靠近测量周期的中点，同时，确定目标时刻点在不同星间观测量中的位置。

步骤3：根据积累的星间观测量计算归化算法系数矩阵。

步骤4：根据系数矩阵和目标时刻点得到目标时刻点的伪距值（观测量）。

步骤5：插值节点移位，计算下一周期系数矩阵。

图5是本发明的几何定位算法结构图

（1）给定卫星初始位置。

（2）对定位方程进行线性化，并忽略2阶以上高阶项，可得

（21）

其中，为地面站到卫星距离的近似值

（22）

为地面站到卫星的观测矢量的方向余弦，即

（23）

于是，令

，，，

最小二乘解为：

（24）

其中，为观测权矩阵。

（3）用上一步计算结果更新卫星位置。

（4）比较两次计算结果，判断是否满足收敛标准。如果满足则直接输出卫星位置估计值；否则用这次计算的卫星位置作为初始值，返回（2）继续计算。

图6是本发明的动力学定轨算法流程图

（1）给定卫星初始状态

（2）计算时刻的积分参考轨道和状态转移矩阵，

（3）计算预报状态

（25）

（4）计算预测状态协方差矩阵

（26）

其中，为模型噪声的协方差矩阵。

（5）计算新息向量及其协方差矩阵

（27）

（28）

其中，为测量噪声的协方差矩阵。

（6）计算增益矩阵

（29）

（7）求解新的状态估计值

（30）

（8）更新状态协方差矩阵

（31）

（9）如果满足设定门限则直接输出结果，否则返回（1）继续计算。

对时刻的积分参考轨道进行修正

（32）

图7是本发明的星历参数拟合生成算法流程图。

计算卫星位置所用到的基本广播星历参数有16个：、、e、、、、、、、、、、、、、，其中是作为已知量给出的，不需要通过拟合得到。我们将前面16个参数中除外的15个作为待求参数，由此可得相应的待估状态参数向量和观测方程为

（33）

（34）

式中，为参考历元时刻的广播星历参数，为一个含个观测量的观测列向量，一个观测量对应卫星的一个位置分量。

设为估值在第次迭代的初值，将观测方程在所给初值处展开，并舍去二阶和二阶以上的小量后可得

（35）

式中，为用参考历元时刻广播星历参数初值计算的卫星位置，分别为相应广播星历参数的改正值，、…为观测量对广播星历参数的偏导数。令

,,

可得误差方程

（36）

由最小二乘原理有

（37）

则第次迭代后的广播星历参数估值为

（38）

在实际计算时，所选用的迭代结束条件为

（39）

其中，是预先给定的任意小的正数，是第次迭代的单位权方差。

导航卫星的基本广播星历参数有16个：、、e、、、、、、、、、、、、、。由这16个参数可以计算卫星在地固坐标系中的位置坐标，计算步骤如下：

计算半长轴：

计算卫星平均角速度：

计算瞬时历元到参考历元的时间差：

改正平角速度：

计算平近点角：

迭代计算偏近点角：

计算真近点角：

计算纬度参数：

计算周期改正项、、：

计算改正后的纬度参数：

计算改正后的向径：

计算改正后的倾角：

计算卫星在轨道平面内的坐标：

改正升交点的经度：

计算卫星在WGS84坐标系中的坐标：

其中，和分别为地球引力常数和地球自转角速度。

将位置坐标改写为

本发明未具体说明部分属于本领域公知技术。

Claims

1.基于星载处理器的自主定轨方法，其特征在于，通过星间精密测量得到星间精密测量值，即距离观测量，对星间距离观测量进行时间归化，得到同一时刻的星间测距观测量；经过一段时间的累积，得到卫星的几何定位值；在卫星几何定位的基础上，进行卫星动力学定轨；轨道动力学外推然后输出预报星历，具体步骤如下：

步骤1：星间双向观测量生成

星间观测量并非连续获得，而是根据星间链路的节拍分时获取，假设节拍周期为，起始时刻为，则在2时间内完成一次双向测量，具体为：时刻A星发B星收，+时刻A星收B星发；+2时刻A星发C星收，+3时刻A星收C星发；以此类推，在一个测量周期内的第二个节拍，B星在完成测量的同时，将第一个节拍得到的测量数据发送回A星，即A星在一个测量周期内得到了AB星的双向测量数据；

步骤2：星间双向观测量归算

由于星间链路直接观测量为不同时刻的测量值，因此在用于自主定轨之前需要进行预处理，即将观测量归化到同一时刻，输入为本地卫星A与建链卫星B的多个周期接收时刻伪距观测值，输出值为将各周期接收时刻伪距归化到指定接收时刻对应的伪距；

步骤3：几何定位

具体为：

（3.1）给定卫星初始位置

对于不同的卫星，建立不同的观测量方程组,

（1）

对式(1)进行求解，即得到卫星的位置坐标，对式(1)进行线性化，假设已知卫星的概略位置坐标，则只需对概略位置进行相应的改正即可得到卫星位置坐标，改正量用表示，此时得到的卫星位置坐标应为估计值，记为，于是有

（2）

（3.2）对式（1）在卫星概略位置坐标处进行线性化，并忽略2阶以上高阶项，得

（3）

其中，为地面站到卫星距离的近似值，

（4）

为地面站到卫星的观测矢量的方向余弦，即

（5）

于是，令

，，，

则式（3）可写成矩阵形式

（6）

令为式（6）的解算原则，即

（7）

则式（7）的最优解为

（8）

当观测的卫星只有3个时，式（8）的解为

（9）

当观测的卫星数量大于3个时，式（8）的解为

（10）

对（10）式给出加权最小二乘解为

（11）

其中，为观测权矩阵；

（3.3）用上一步计算结果更新卫星位置；

（3.4）比较两次计算结果，判断是否满足收敛标准，如果满足则直接输出卫星位置估计值；否则用这次计算的卫星位置作为初始值，返回步骤（3.2）继续计算；

步骤4：卫星动力学定轨

（4.1）给定卫星初始位置

动力学定轨即利用动力学模型，得到卫星的轨迹曲线，采用动力学离散模型为：

(12)

其中，、、、、、、分别为：当前历元轨道修正值、前一历元轨道修正值、模型噪声矩阵、观测值、观测矩阵、状态传递矩阵、观测误差矩阵；

应用扩展Kalman滤波方法对（12）式进行求解，其步骤如下：

（4.2）计算得到时刻的积分参考轨道和状态转移矩阵；

（4.3）从而得到预报状态

（13）

（4.4）预测状态协方差矩阵

（14）

其中，为模型噪声的协方差矩阵，为为上一状态的轨道估计值；

（4.5）计算新息向量及其协方差矩阵

（15）

（16）

其中，为测量噪声的协方差矩阵，为新息向量，为新息向量的协方差矩阵，

（4.6）计算增益矩阵

（17）

（4.7）求解新的状态估计值

（18）

（4.8）更新状态协方差矩阵

（19）

从而可以对时刻的积分参考轨道进行修正

（20）

（4.9）如果满足设定门限则直接输出结果，否则返回步骤（4.1）继续计算；

步骤5：星历参数拟合

设星历参考时刻，则相对于定轨时刻是未来时，需要将定轨时刻得到的轨道进行外推预报，

设动力学定轨得到的轨道为，为轨道历元，为轨道根数，这里的轨道外推涉及两个过程：一是将时刻的轨道预报到星历参考时刻，得到；二是为了计算星历参数所需的轨道外推，在时刻轨道的基础上外推1-2个小时。

2.根据权利要求1所述的基于星载处理器的自主定轨方法，其特征在于，所述星间双向观测量归算，具体做法为：

（1）首先确定本星需要建立双向测量的卫星数量，假设为，则一个测量周期为，

（2）确定不同测量值需要规化的时刻点，为某条测量链路起始时刻到测量周期的中点，即时刻；这时计算其它链路对应到此时刻点的值，由于测量链路为时分体制，因此对于不同的测量链路此时刻点不相同，

（3）每条测量链路分别进行归化计算，确保测量周期以及归化时刻点的前后数据大体相当，此时总的数据点数量假设为，

3.根据权利要求1所述的基于星载处理器的自主定轨方法，其特征在于，所述星历参数拟合，计算卫星位置所用到的基本广播星历参数有16个：、、e、、、、、、、、、、、、、，其中是作为已知量给出的，将16个参数中除外的15个作为待求参数，得相应的待估状态参数向量和观测方程为

（33）

（34）

式中，为参考历元时刻的广播星历参数，为一个含个观测量的观测列向量，一个观测量对应卫星的一个位置分量，

（35）

式中，为用参考历元时刻广播星历参数初值计算的卫星位置，分别为相应广播星历参数的改正值，、…为观测量对广播星历参数的偏导数，令

,,

可得误差方程

（36）

由最小二乘原理有

（37）

则第次迭代后的广播星历参数估值为

（38）

在实际计算时，所选用的迭代结束条件为

（39）

其中，是预先给定的任意小的正数，是第次迭代的单位权方差，

导航卫星的基本广播星历参数有16个：、、e、、、、、、、、、、、、、，由这16个参数可以计算卫星在地固坐标系中的位置坐标，计算步骤如下：

计算：

计算卫星平均角速度：

计算瞬时历元到参考历元的时间差：

改正平角速度：

计算平近点角：

迭代计算偏近点角：

计算真近点角：

计算纬度参数：

计算周期改正项、、：

计算改正后的纬度参数：

计算改正后的向径：

计算改正后的迭代次数：

计算卫星在轨道平面内的坐标：

改正升交点的经度：

计算卫星在WGS84坐标系中的坐标：

其中，和分别为地球引力常数和地球自转角速度。