CN101435863A - 一种导航卫星实时精密定轨的方法 - Google Patents

Abstract

本发明公开了一种基于滑动数据窗口的导航卫星实时精密定轨的方法，该方法的主要思想是：在对状态X_k，包括卫星初始状态、力模型参数、地球自转参数、大气参数、测站坐标等进行估计时，只利用离k时刻最近的前N个量测，而把之前其余的量测完全甩掉，这里N是预先规定的滑动窗口长度，再通过建立，合并，求解法方程等步骤获得精确的第N弧段的卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户。采用本方法进行卫星定轨具有数据处理速度快，轨道精度高，数值稳定等优点。

Description

一种导航卫星实时精密定轨的方法

技术领域

本发明涉及一种卫星导航技术，具体的说是一种导航卫星实时精密定轨的方法。

背景技术

导航卫星精密定轨作为卫星导航系统的关键技术，一直受到国内外众多研究机构的重视，而成为相关研究领域的前沿课题。自从1994年国际GNSS地球动力学服务组织成立以来，精密GPS轨道与钟差产品在精度、时间延迟上不断的改进。目前，IGS提供的事后精密星历轨道精度优于5cm，钟差精度优于0.1ns，可满足高精度用户的事后定位需求。随着精密单点定位技术应用的广泛深入，已有越来越多的GPS用户对导航星座轨道的实时性和精度提出了更高要求。为促进GPS高精度实时应用领域的发展，IGS组织、美国的JPL、欧洲的BKG、GFZ、ESA等知名研究机构都在开展相关研究。

导航卫星实时轨道改进最早在广域差分系统中提出，早期通过几何法对广播星历计算的卫星位置进行实时定轨，虽然此方法简单、易于实现，也能够在一定程度上改善轨道的精度，但由于没有考虑卫星轨道的运动特性，其噪声大而且易受观测的几何形状和数据丢失的影响，以致轨道预报的精度难以保证，因此，我们对于高精度的导航卫星实时定轨应充分考虑卫星的动力学信息。导航卫星定轨是一个大型复杂的数据处理系统，涉及到复杂的动力学模型和观测模型。由于导航卫星实时精密轨道生成过程中数据传输、数据处理以及实时轨道产品发播都存在时间延迟，实时用户需求的卫星位置一般通过实时更新的精密卫星轨道初值和力模型参数向外积分卫星运动方程得到。因而，实时轨道的精度与卫星轨道初值和力模型参数精度、轨道预报模型精度、外推时间长短三者相关。又由于GNSS卫星轨道预报模型及方法的改进空间很小，因此导航卫星实时定轨技术的关键在于快速精确地更新轨道初值与力模型参数。

目前，实时精密定轨参数一般采用线性滤波的方法，并且逐个历元更新卫星轨道。对于高度在20000多公里的导航卫星，在非机动及无故障的正常情况下，卫星轨道变化是平滑的。因此，短时间的轨道预报精度很高。实践表明，基于1天对IGS事后精密轨道进行拟合结果并外推1小时轨道的精度达到5cm左右，与估计轨道精度基本相当。在采用逐历元滤波过程中，为提高解算速度以满足实时需求，系统中估计参数的个数与观测值的个数要求尽量少，这就会影响到系统模型以及参数解的精度。因而，对于实时定轨而言，是否有必要采用逐历元估计的方法是值得进一步研究的。

一般而言，滤波方法不可避免地存在初始化、数值稳定性、数据饱和以及处理速度等棘手问题。初始化是指观测数据累积到一定的时间跨度后参数估计才能提供稳定可靠的结果。在实时定轨中这个问题非常突出，因为初始化观测时间要大于12小时，初始化所需要的时间至少等于处理12小时数据的计算时间。同时滤波中稳健的实时数据质量控制是非常困难的问题，由此产生的滤波发散将导致估计系统崩溃，系统需要重新初始化。

美国JPL的RTG软件采用平方根信息滤波的方法实现实时精密定轨，也同样存在以上这些方面的问题，其参与实时解算的测站数量在70个左右，定轨精度在20cm左右，甚至差于IGS的快速轨道IGU预报部分精度；并且实时定轨系统需要多个备份系统支持，确保在处理系统崩溃的状态下迅速恢复提供用户实时轨道。

发明内容

本发明的目的就是为了克服上述背景技术的不足之处，提供一种导航卫星实时精密定轨的方法，它能有效解决实时定轨中初始化、数值稳定性、数据饱和以及处理速度等方面的问题。

为了实现上述目的，本发明提供的一种导航卫星精密定轨的方法，该方法包括以下步骤：

(1)将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为N*T小时；

(2)处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数据，生成该弧段上的短弧法方程，其短弧法方程的公式为N_k，kδx_k＝W_k，k，其中δx_k为该弧段的轨道初始位置与力模型初值相应的改正参数，N_k，k为法方程矩阵，k为短弧段标记；

(3)将滑动数据窗口内前N-1个短弧段生成的短弧法方程进行合并；其短弧合并利用的公式如下

两个相连弧段的轨道合并公式

$\left[\begin{array}{cc}{N}_{k,k}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{(k+1,k)}& -{\mathrm{\Φ}}^T{P}_{(k+1,k)}\\ -P{\mathrm{\Φ}}_{(k+1,k)}& N_{k+1,k+1}+P_w^{k+1,k}\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}{x}_k\\ \mathrm{\δ}{x}_{k+1}\end{array}\right]=\left[\begin{array}{c}{W}_{k,k}+{\mathrm{\Φ}}_{k+1,k}^T{W}_w^{(k+1,k)}\\ W_{k+1,k+1}-W_w^{(k+1,k)}\end{array}\right]$                     式1

n个相连弧段的轨道合并公式

$\left[\begin{array}{ccccc}{N}_{\mathrm{1,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{2,1}\right)}& 0& \·\·\·& 0\\ -P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& N_{\mathrm{2,2}}+P_w^{\mathrm{2,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ 0& -P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& N_{\mathrm{3,3}}+P_w^{\mathrm{3,2}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& N_{n,n}+P_w^{n,n-1}\end{array}\right]$

$\×\left[\begin{array}{c}{\mathrm{\δ}}_1\\ {\mathrm{\δ}}_2\\ {\mathrm{\δ}}_3\\ \·\·\·\\ {\mathrm{\δ}}_n\end{array}\right]=\left[\begin{array}{c}{W}_{\mathrm{1,1}}+{\mathrm{\Φ}}_{\mathrm{2,1}}^T{W}_w^{\mathrm{2,1}}\\ W_{\mathrm{2,2}}-W_w^{\mathrm{2,1}}+{\mathrm{\Φ}}_{\mathrm{3,2}}^T{W}_w^{\mathrm{3,2}}\\ W_{\mathrm{3,3}}-W_w^{\mathrm{3,2}}+{\mathrm{\Φ}}_{\mathrm{4,3}}^T{W}_w^{\mathrm{4,3}}\\ \·\·\·\\ W_{n,n}-W_w^{n,n-1}\end{array}\right]$                            式2

其中，Φ(t_k+1，t_k)是状态转移矩阵；Φ描述了函数关系模型上的连接关系，P_w描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度；

根据短弧轨道间参数类型的不同相连弧段间连接时Φ与P_w的设置分为以下三类，

(a)对于前后弧段相同或具有确定性的状态转移关系的参数，在轨道连接处是相同的，设P_w为无穷大实现，

(b)对于前后两个弧段不相关的分段参数或白噪声参数，设Φ＝0，P_w＝0实现，

(c)对于随机过程噪声参数，通过其过程噪声的标准差给定相应的P_w值实现；

(4)将当前短弧段的法方程与由前N-1个短弧段法方程合并成的一个法方程进行合并；

(5)求解由整个滑动数据窗口内短弧合并而成的法方程，得到卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户；

(6)将滑动窗口向下移动一个短弧段，返回步骤(2)。

在上述技术方案中，步骤(1)中所述的滑动数据窗口的长度即N*T的值不超过24小时。

在上述技术方案中，步骤(2)中所述的短弧段上获得的实时数据包括卫星初始状态、力模型参数、地球自转参数、大气参数、测站坐标。

在上述技术方案中，步骤(3)中所述的短弧合并公式推导过程如下：

(A)卫星运动状态连接，在导航卫星运动的实际过程中，在正常情况下，卫星运动轨迹和速度是连续的。如图2所示，在相邻的两个短弧的交界处即同一时刻，卫星的状态与受到的摄动力是相同的，可用以下的关联条件表示

r_k(t_k+1)＝r_k+1(t_k+1)

${\stackrel{.}{r}}_k\left(t_{k+1}\right)={\stackrel{.}{r}}_{k+1}\left(t_{k+1}\right)$

q_k(t_k+1)＝q_k+1(t_k+1)

                                          式3

其中，r_k(t_k+1)、q_k是第k弧段t_k+1时刻的卫星位置、速度与力模型参数；

r_k+1(t_k+1)、

q_k+1是第k+1弧段t_k+1时刻的卫星位置、速度与力模型参数。

将上述关联条件通过两个相连弧段各自的轨道初值与力模型参数可表示为如下式4，5。

$\left[\begin{array}{c}{r}_{k+1}\left(t_{k+1}\right)\\ {\stackrel{.}{r}}_{k+1}\left(t_{k+1}\right)\\ q_{k+1}\left(t_{k+1}\right)\end{array}\right]=\left[\begin{array}{c}{r}_{k+1}{\left(t_{k+1}\right)}_0\\ {\stackrel{.}{r}}_{k+1}{\left(t_{k+1}\right)}_0\\ q_{k+1}{\left(t_{k+1}\right)}_0\end{array}\right]+\left[\begin{array}{c}\mathrm{\δ}{r}_{k+1}\left(t_{k+1}\right)\\ \mathrm{\δ}{\stackrel{.}{r}}_{k+1}\left(t_{k+1}\right)\\ \mathrm{\δ}{q}_{k+1}\left(t_{k+1}\right)\end{array}\right]$                     式4

$\left[\begin{array}{c}{r}_k\left(t_{k+1}\right)\\ {\stackrel{.}{r}}_k\left(t_{k+1}\right)\\ q_k\left(t_{k+1}\right)\end{array}\right]=\left[\begin{array}{c}{r}_k{\left(t_{k+1}\right)}_0\\ {\stackrel{.}{r}}_k{\left(t_{k+1}\right)}_0\\ q_k{\left(t_{k+1}\right)}_0\end{array}\right]+\mathrm{\Φ}(t_{k+1},t_k)\left[\begin{array}{c}\mathrm{\δ}{r}_k\left(t_k\right)\\ \mathrm{\δ}{\stackrel{.}{r}}_k\left(t_k\right)\\ \mathrm{\δ}{q}_k\left(t_k\right)\end{array}\right]$                    式5

其中，

r_k+1(t_k+1)₀、

q_k+1(t_k+1)₀为k+1弧段t_k+1时刻的轨道与力模型参数初值；

δr_k+1(t_k+1)、

δq_k+1(t_k+1)为k+1弧段t_k+1时刻轨道初始位置与力模型初值相应的改正参数；

r_k(t_k+1)₀、

q_k(t_k+1)₀为k弧段t_k+1时刻轨道与力模型参数初值；

δr_k(t_k)，

δq_k(t_k)为k弧段t_k时刻轨道初始位置与力模型初值相应的改正参数；

Φ(t_k+1，t_k)是状态转移矩阵。

(B)观测模型参数连接，与卫星连续运动相同，导航卫星的观测在正常情况下也是连续的，由此相连两个弧段的观测模型参数也是关联的。观测模型参数包括对流层参数、卫星钟差参数、接收机参数、地球自转参数、地面站坐标、模糊度参数等，这些待估参数可归类为随机过程噪声参数与状态参数。对于状态参数同样可以建立与上述卫星运动状态连接类似的关联方程，只是状态转移矩阵具有不同的定义。对于随机过程噪声参数，一般可以通过以下公式表示：

P_k+1＝M_k+1，kP_k+w_k                       式6

由于w_j可以认为是一独立随机过程噪声，并服从正态分布，上式可以转换为

r_wP_k+1＝M_k+1，kr_wP_k+r_ww_k                 式7

其中r_w为相应过程噪声的标准差，考虑到w_k为零均值，r_ww_k可设为零，P_k、P_k+1为随机过程参数，M_k+1，k为映射矩阵，这样随机过程噪声参数与状态参数一样建立类似的关联方程。

(C)由此将上述两种数学模型统一起来表述，将不同的参数都用X表示，这样关联条件则表示为

X_k(t_k+1)₀+Φ(t_k+1，t_k)δx_k＝X_k+1(t_k+1)₀+δx_k+1+Δ             式8

约束方程表述为

v＝Φ(t_k+1，t_k)δx_k-δx_k+1+X_k(t_k+1)₀-X_k+1(t_k+1)₀             式9

 ＝Φ(t_k+1，t_k)δx_k-δx_k+1-δl                         P_w

对于短弧轨道间不同参数的连接或合并，可以通过设置参数前后弧段的映射因子Φ与动态噪声P_w来实现。实际上，Φ描述了函数关系模型上的连接关系，P_w描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度，δl为观测值残差，以下讨论不同类型参数相连弧段间连接时Φ与P_w的设置。

对于前后弧段相同或具有确定性的状态转移关系的参数，如卫星状态参数、测站坐标参数、连续的模糊度参数等，在轨道连接处是相同的，可以设P_w为无穷大实现；

对于前后两个弧段不相关的分段参数或白噪声参数，如卫星钟差、接收机钟差等，可以设Φ＝0，P_w＝0实现；

对于随机过程噪声参数，可以通过其过程噪声的标准差给定相应的P_w值实现。

(D)基于相连两个弧段各自的短弧法方程信息，与边界约束条件方程，推导短弧合并的统一数学模型。

假设相连弧段k，k+1对应的法方程为

N_k，kδx_k＝W_k，k

                                       式10

N_k+1，k+1δx_k+1＝W_k+1，k+1

边界约束条件方程为

$\begin{array}{cc}{v}_w^{k+1,k}=\mathrm{\Φ}(t_{k+1},t_k)\mathrm{\δ}{x}_k-\mathrm{\δ}{x}_{k+1}-\mathrm{\δ}{l}_w^{k+1,k}& P_w^{k+1,k}\end{array}$               式11

将上式10假设为先验约束，式11作为虚拟观测方程，则轨道合并的问题就可以转换为附有先验信息的平差问题，采用最小方差估计，要求满足性能函数J(x)为最小

$J\left(x\right)=v_k^T{p}_k{v}_k+v_{k+1}^T{p}_{k+1}{v}_{k+1}+v_w^T{p}_w{v}_w=\min$                   式12

则 $\frac{\∂J\left(x\right)}{\∂x}=0$ 可推导得到：

$\left[\begin{array}{cc}{N}_{k,k}+{\mathrm{\Φ}}_{k+1,k}^T{P}_w^{k+1,k}{\mathrm{\Φ}}_{k+1,k}& -{\mathrm{\Φ}}_{k+1,k}^T{P}_w^{k+1,k}\\ -P_w^{k+1,k}{\mathrm{\Φ}}_{k+1,k}& N_{k+1,k+1}+P_w^{k+1,k}\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}{x}_k\\ \mathrm{\δ}{x}_{k+1}\end{array}\right]=\left[\begin{array}{c}{W}_{k,k}+{\mathrm{\Φ}}_{k+1}^T{P}_w^{k+1,k}\mathrm{\δ}{l}_w^{k+1,k}\\ W_{k+1,k+1}-P_w^{k+1,k}\mathrm{\δ}{l}_w^{k+1,k}\end{array}\right]$   式13

简化表达，设

${\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{(k+1,k)}={\mathrm{\Φ}}_{k+1,k}^T{P}_w^{k+1,k}{\mathrm{\Φ}}_{k+1,k}$                          式14

$P{\mathrm{\Φ}}_{(k+1,k)}=P_w^{k+1,k}{\mathrm{\Φ}}_{k+1,k}$                              式15

$W_w^{(k+1,k)}=P_w^{k+1,k}\mathrm{\δ}{l}_w^{k+1,k}$                               式16

上式13可简化为：

$\left[\begin{array}{cc}{N}_{k,k}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{(k+1,k)}& -{\mathrm{\Φ}}^T{P}_{(k+1,k)}\\ -P{\mathrm{\Φ}}_{(k+1,k)}& N_{k+1,k+1}+P_w^{k+1,k}\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}{x}_k\\ \mathrm{\δ}{x}_{k+1}\end{array}\right]=\left[\begin{array}{c}{W}_{k,k}+{\mathrm{\Φ}}_{k+1,k}^T{W}_w^{(k+1,k)}\\ W_{k+1,k+1}-W_w^{(k+1,k)}\end{array}\right]$          式1

以上推导了两个相连弧段的轨道合并公式，下面直接给出一般性的n个弧段的轨道合并公式

$\left[\begin{array}{ccccc}{N}_{\mathrm{1,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{2,1}\right)}& 0& \·\·\·& 0\\ -P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& N_{\mathrm{2,2}}+P_w^{\mathrm{2,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ 0& -P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& N_{\mathrm{3,3}}+P_w^{\mathrm{3,2}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& N_{n,n}+P_w^{n,n-1}\end{array}\right]$

$\×\left[\begin{array}{c}{\mathrm{\δ}}_1\\ {\mathrm{\δ}}_2\\ {\mathrm{\δ}}_3\\ \·\·\·\\ {\mathrm{\δ}}_n\end{array}\right]=\left[\begin{array}{c}{W}_{\mathrm{1,1}}+{\mathrm{\Φ}}_{\mathrm{2,1}}^T{W}_w^{\mathrm{2,1}}\\ W_{\mathrm{2,2}}-W_w^{\mathrm{2,1}}+{\mathrm{\Φ}}_{\mathrm{3,2}}^T{W}_w^{\mathrm{3,2}}\\ W_{\mathrm{3,3}}-W_w^{\mathrm{3,2}}+{\mathrm{\Φ}}_{\mathrm{4,3}}^T{W}_w^{\mathrm{4,3}}\\ \·\·\·\\ W_{n,n}-W_w^{n,n-1}\end{array}\right]$                            式2

相比较现有技术而言，运用本方法实现导航卫星实时精密定轨具有如下优点：一、处理速度快，由于法方程合并中大量的参数都可以预先消去，如除了最后一个弧段外的其它弧段上的轨道参数、模糊度参数等，这样，合并的计算量就可以大大减小，从而提高数据处理速度。二、快速初始化，一旦实时系统需要重新启动，此方法只需处理当前短弧段上的观测数据，然后与已经存在的法方程合并就可以实时更新轨道，从而实现快速初始化。三、数值稳定，由于使用滑动数据窗口的方法，每次轨道更新只用当前窗口中的数据，陈旧的观测信息被去除，从而可以完全克服“数据饱和”导致的动力学模型参数估计精度降低甚至错误的问题。四、能有效实现数据质量的控制，大大提高实时定轨的精度。

附图说明

图1为本发明一种导航卫星实时精密定轨的方法中滑动数据窗口定轨示意图。

图2为本发明方法中短弧连接示意图。

图3为本发明方法中实施例试算每小时GPS实时轨道与IGS事后轨道比较结果图。

具体实施方式

下面结合附图对本发明作进一步的描述。

如图1所示，一种导航卫星实时精密定轨的方法，主要是将实时定轨通过建立短弧法方程和短弧法方程合并两个并行进程实现，该方法包括以下步骤：

(1)将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为N*T小时，其N*T的值不超过24小时；

(2)处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数据，建立该弧段上包括卫星初始状态、力模型参数、地球自转参数、大气参数、测站坐标的短弧法方程，其短弧法方程的公式为N_k，kδx_k＝W_k，k，其中δx_k为该弧段的轨道初始位置与力模型初值相应的改正参数，N_k，k为法方程矩阵，k为短弧段标记；

(3)短弧法方程合并进程将前N-1个弧段生成的法方程进行合并；其短弧合并利用的公式如下

两个相连弧段的轨道合并公式

$\left[\begin{array}{cc}{N}_{k,k}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{(k+1,k)}& -{\mathrm{\Φ}}^T{P}_{(k+1,k)}\\ -P{\mathrm{\Φ}}_{(k+1,k)}& N_{k+1,k+1}+P_w^{k+1,k}\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}{x}_k\\ \mathrm{\δ}{x}_{k+1}\end{array}\right]=\left[\begin{array}{c}{W}_{k,k}+{\mathrm{\Φ}}_{k+1,k}^T{W}_w^{(k+1,k)}\\ W_{k+1,k+1}-W_w^{(k+1,k)}\end{array}\right]$                        式1

n个弧段的轨道合并公式

$\left[\begin{array}{ccccc}{N}_{\mathrm{1,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{2,1}\right)}& 0& \·\·\·& 0\\ -P{\mathrm{\Φ}}_{\left(\mathrm{2,1}\right)}& N_{\mathrm{2,2}}+P_w^{\mathrm{2,1}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ 0& -P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& N_{\mathrm{3,3}}+P_w^{\mathrm{3,2}}+{\mathrm{\Φ}}^{T}P{\mathrm{\Φ}}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& N_{n,n}+P_w^{n,n-1}\end{array}\right]$

$\×\left[\begin{array}{c}{\mathrm{\δ}}_1\\ {\mathrm{\δ}}_2\\ {\mathrm{\δ}}_3\\ \·\·\·\\ {\mathrm{\δ}}_n\end{array}\right]=\left[\begin{array}{c}{W}_{\mathrm{1,1}}+{\mathrm{\Φ}}_{\mathrm{2,1}}^T{W}_w^{\mathrm{2,1}}\\ W_{\mathrm{2,2}}-W_w^{\mathrm{2,1}}+{\mathrm{\Φ}}_{\mathrm{3,2}}^T{W}_w^{\mathrm{3,2}}\\ W_{\mathrm{3,3}}-W_w^{\mathrm{3,2}}+{\mathrm{\Φ}}_{\mathrm{4,3}}^T{W}_w^{\mathrm{4,3}}\\ \·\·\·\\ W_{n,n}-W_w^{n,n-1}\end{array}\right]$                                          式2

其中，Φ(t_k+1，t_k)是状态转移矩阵，Φ描述了函数关系模型上的连接关系，P_w描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度；

根据短弧轨道间参数类型的不同相连弧段间连接时Φ与P_w的设置分为以下三类，

(a)对于前后弧段相同或具有确定性的状态转移关系的参数，如卫星状态参数、测站坐标参数、连续的模糊度参数等，在轨道连接处是相同的，设P_w为无穷大实现；

(b)对于前后两个弧段不相关的分段参数或白噪声参数，如卫星钟差、接收机钟差等，设Φ＝0，P_w＝0实现；

(c)对于随机过程噪声参数，通过其过程噪声的标准差给定相应的P_w值实现；

(4)将当前弧段的法方程与由前N-1个弧段法方程合并成的一个法方程合并；

(5)求解由整个滑动窗口中短弧合并而成的法方程，得到卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户；

(6)将滑动窗口向下移动一个弧段，返回步骤(2)。

本实施例采用全球70个IGS站2006年197天至203天一周的观测数据，所选参考站均提供1小时的准实时观测数据文件，进行每小时更新一次的实时精密轨道试算。具体实施方式如下：通过最小二乘方法对每个小时的观测文件形成法方程，实时定轨的滑动窗口长度设为24小时，前23小时为先验信息，第24小时实时估计轨道参数，基于估计的实时轨道参数外推1小时的轨道作为提供用户的实时轨道，以下分析中的实时轨道均为此外推轨道，从197天开始往前递推获得198至203天的实时轨道，197天轨道作为实时系统启动的初始信息。

图3给出了以每小时为单位6天实时轨道与IGS事后精密轨道比较统计的平均RMS值的时间序列，从图中可以看出90％以上的RMS值小于6cm，98％以上的RMS值小于8cm，最大RMS值在10cm以内，绝大部分卫星平均RMS值在4-8cm之间，由此试算结果表明使用本方法每小时更新的实时定轨与IGS事后精密轨道比较，精度达到5cm左右。

Claims (3)

1.一种导航卫星实时精密定轨的方法，其特征在于该方法包括以下步骤：

(1)将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为N*T小时；

(2)处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数据，生成该弧段上的短弧法方程，其短弧法方程的公式为N_k，kδx_k＝W_kk，其中δx_k为该弧段的轨道初始位置与力模型初值相应的改正参数，N_k，k为法方程矩阵，k为短弧段标记；

(3)将滑动数据窗口内前N-1个短弧段生成的短弧法方程进行合并；其短弧合并利用的公式如下

两个相连弧段的轨道合并公式

$\left[\begin{array}{cc}{N}_{k,k}+{\mathrm{\Φ}}^T{\mathrm{P\Φ}}_{(k+1,k)}& -{\mathrm{\Φ}}^T{P}_{(k+1,k)}\\ -{\mathrm{P\Φ}}_{(k+1,k)}& N_{k+1,k+1}+P_w^{k+1,k}\end{array}\right]\left[\begin{array}{c}{\mathrm{\δx}}_k\\ {\mathrm{\δx}}_{k+1}\end{array}\right]=\left[\begin{array}{c}{W}_{k,k}+{\mathrm{\Φ}}_{k+1,k}^T{W}_w^{(k+1,k)}\\ W_{k+1,k+1}-W_w^{(k+1,k)}\end{array}\right]$     式1

n个相连弧段的轨道合并公式

$\left[\begin{array}{ccccc}{N}_{\mathrm{1,1}}+{\mathrm{\Φ}}^T{\mathrm{P\Φ}}_{\left(\mathrm{2,1}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{2,1}\right)}& 0& \·\·\·& 0\\ -{\mathrm{P\Φ}}_{\left(\mathrm{2,1}\right)}& N_{\mathrm{2,2}}+P_w^{\mathrm{2,1}}{\mathrm{\Φ}}^T{\mathrm{P\Φ}}_{\left(\mathrm{3,2}\right)}& -{\mathrm{\Φ}}^T{P}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ 0& -{\mathrm{P\Φ}}_{\left(\mathrm{3,2}\right)}& N_{\mathrm{3,3}}+P_w^{\mathrm{3,2}}{\mathrm{\Φ}}^T{\mathrm{P\Φ}}_{\left(\mathrm{3,2}\right)}& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& \·\·\·& N_{n,n}+P_w^{n,n-1}\end{array}\right]$

$\×\left[\begin{array}{c}{\mathrm{\δ}}_1\\ {\mathrm{\δ}}_2\\ {\mathrm{\δ}}_3\\ \·\·\·\\ {\mathrm{\δ}}_n\end{array}\right]=\left[\begin{array}{c}{W}_{\mathrm{1,1}}+{\mathrm{\Φ}}_{\mathrm{2,1}}^T{W}_w^{\mathrm{2,1}}\\ W_{\mathrm{2,2}}-W_w^{\mathrm{2,1}}+{\mathrm{\Φ}}_{\mathrm{3,2}}^T{W}_w^{\mathrm{3,2}}\\ W_{\mathrm{3,3}}-W_w^{\mathrm{3,2}}+{\mathrm{\Φ}}_{4,3}^T{W}_w^{\mathrm{4,3}}\\ \·\·\·\\ W_{n,n}-W_w^{n,n-1}\end{array}\right]$

                                              式2

其中，Φ(t_k+1，t_k)是状态转移矩阵；Φ描述了函数关系模型上的连接关系，P_w描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度；

根据短弧轨道间参数类型的不同相连弧段间连接时Φ与P_w的设置分为以下三类，

(a)对于前后弧段相同或具有确定性的状态转移关系的参数，在轨道连接处是相同的，设P_w为无穷大实现，

(b)对于前后两个弧段不相关的分段参数或白噪声参数，设Φ＝0，P_w＝0实现，

(c)对于随机过程噪声参数，通过其过程噪声的标准差给定相应的P_w值实现；

(4)将当前短弧段的法方程与由前N-1个短弧段法方程合并成的一个法方程进行合并；

(5)求解由整个滑动数据窗口内短弧合并而成的法方程，得到卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户；

(6)将滑动窗口向下移动一个短弧段，返回步骤(2)。

2.根据权利要求1所述的一种导航卫星实时精密定轨的方法，其特征是：步骤(1)中所述的滑动数据窗口的长度即N*T的值不超过24小时。

3.根据权利要求1所述的一种导航卫星实时精密定轨的方法，其特征是：步骤(2)中所述的短弧段上获得的实时数据包括卫星初始状态、力模型参数、地球自转参数、大气参数、测站坐标。

Images