专利名称：一种导航卫星实时精密定轨的方法
技术领域：
本发明涉及一种卫星导航技术，具体的说是一种导航卫星实时精密定轨的方法。
背景技术：
导航卫星精密定轨作为卫星导航系统的关键技术，一直受到国内外众多研究 机构的重视，而成为相关研究领域的前沿课题。自从1994年国际GNSS地球动力 学服务组织成立以来，精密GPS轨道与钟差产品在精度、时间延迟上不断的改进。 目前，IGS提供的事后精密星历轨道精度优于5cm，钟差精度优于0.1ns，可满 足高精度用户的事后定位需求。随着精密单点定位技术应用的广泛深入，已有越 来越多的GPS用户对导航星座轨道的实时性和精度提出了更高要求。为促进GPS 高精度实时应用领域的发展，IGS组织、美国的JPL、欧洲的BKG、 GFZ、 ESA等 知名研究机构都在开展相关研究。
导航卫星实时轨道改进最早在广域差分系统中提出，早期通过几何法对广播 星历计算的卫星位置进行实时定轨，虽然此方法简单、易于实现，也能够在一定 程度上改善轨道的精度，但由于没有考虑卫星轨道的运动特性，其噪声大而且易 受观测的几何形状和数据丢失的影响，以致轨道预报的精度难以保证，因此，我 们对于高精度的导航卫星实时定轨应充分考虑卫星的动力学信息。导航卫星定轨 是一个大型复杂的数据处理系统，涉及到复杂的动力学模型和观测模型。由于导 航卫星实时精密轨道生成过程中数据传输、数据处理以及实时轨道产品发播都存 在时间延迟，实时用户需求的卫星位置一般通过实时更新的精密卫星轨道初值和 力模型参数向外积分卫星运动方程得到。因而，实时轨道的精度与卫星轨道初值 和力模型参数精度、轨道预报模型精度、外推时间长短三者相关。又由于GNSS 卫星轨道预报模型及方法的改进空间很�。虼说己轿佬鞘凳倍ü旒际醯墓丶� 于快速精确地更新轨道初值与力模型参数。
目前，实时精密定轨参数一般采用线性滤波的方法，并且逐个历元更新卫星 轨道。对于高度在20000多公里的导航卫星，在非机动及无故障的正常情况下， 卫星轨道变化是平滑的。因此，短时间的轨道预报精度很高。实践表明，基于l 天对IGS事后精密轨道进行拟合结果并外推1小时轨道的精度达到5cm左右，与 估计轨道精度基本相当。在采用逐历元滤波过程中，为提高解算速度以满足实时 需求，系统中估计参数的个数与观测值的个数要求尽量少，这就会影响到系统模型以及参数解的精度。因而，对于实时定轨而言，是否有必要采用逐历元估计的 方法是值得进一步研究的。
一般而言，滤波方法不可避免地存在初始化、数值稳定性、数据饱和以及处 理速度等棘手问题。初始化是指观测数据累积到一定的时间跨度后参数估计才能 提供稳定可靠的结果。在实时定轨中这个问题非常突出，因为初始化观测时间要 大于12小时，初始化所需要的时间至少等于处理12小时数据的计算时间。同时 滤波中稳健的实时数据质量控制是非常困难的问题，由此产生的滤波发散将导致 估计系统崩溃，系统需要重新初始化。
美国JPL的RTG软件采用平方根信息滤波的方法实现实时精密定轨，也同样 存在以上这些方面的问题，其参与实时解算的测站数量在70个左右，定轨精度 在20cm左右，甚至差于IGS的快速轨道IGU预报部分精度；并且实时定轨系统 需要多个备份系统支持，确保在处理系统崩溃的状态下迅速恢复提供用户实时轨 道。
发明内容
本发明的目的就是为了克服上述背景技术的不足之处，提供一种导航卫星实 时精密定轨的方法，它能有效解决实时定轨中初始化、数值稳定性、数据饱和以 及处理速度等方面的问题。
为了实现上述目的，本发明提供的一种导航卫星精密定轨的方法，该方法包
括以下步骤
(1) 将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段
的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为
(2) 处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数 据，生成该弧段上的短弧法方程，其短弧法方程的公式为a^々&,=『i4，其中^a 为该弧段的轨道初始位置与力模型初值相应的改正参数，iV^为法方程矩阵，k 为短弧段标记；
(3) 将滑动数据窗口内前N-1个短弧段生成的短弧法方程进行合并；其短
弧合并利用的公式如下
两个相连弧段的轨道合并公式
<formula>formula see original document page 5</formula>W个相连弧段的轨道合并公式
<formula>formula see original document page 6</formula> 式2
其中，O(UJ是状态转移矩阵；(D描述了函数关系模型上的连接关系，Pw 描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程
度；
根据短弧轨道间参数类型的不同相连弧段间连接时①与&的设置分为以下 三类，
(a) 对于前后弧段相同或具有确定性的状态转移关系的参数，在轨道连接 处是相同的，设尺为无穷大实现，
(b) 对于前后两个弧段不相关的分段参数或白噪声参数，设0) = 0， Pw=0
实现，
(c) 对于随机过程噪声参数，通过其过程噪声的标准差给定相应的^值实
现；
(4)将当前短弧段的法方程与由前N-1个短弧段法方程合并成的一个法方 程进行合并；
(5) 求解由整个滑动数据窗口内短弧合并而成的法方程，得到卫星轨道初 值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给 用户；
(6) 将滑动窗口向下移动一个短弧段，返回步骤(2)。在上述技术方案中，步骤(1)中所述的滑动数据窗口的长度即N*T的值 不超过24小时。
在上述技术方案中，步骤(2)中所述的短弧段上获得的实时数据包括卫星 初始状态、力模型参数、地球自转参数、大气参数、测站坐标。
在上述技术方案中，步骤(3)中所述的短弧合并公式推导过程如下
(A)卫星运动状态连接，在导航卫星运动的实际过程中，在正常情况下， 卫星运动轨�：退俣仁橇�。如图2所示，在相邻的两个短弧的交界处即同一 时刻，卫星的状态与受到的摄动力是相同的，可用以下的关联条件表示
式3
其中，^(U、 4(4+,)、 ^是第A弧段^,时刻的卫星位置、速度与力模型参数； 4+1ft+1)、 4+1(U、《w是第A: + l弧段L,时刻的卫星位置、速度与力模型参数。
将上述关联条件通过两个相连弧段各自的轨道初值与力模型参数可表示为 如下式4， 5。
<formula>formula see original document page 7</formula>
其中，
^U。、 4+1(U。、 &+1(4+1)。为"1弧段4+1时刻的轨道与力模型参数初值； &U、 &U、 、+1("+1)为"1弧段~+1时刻轨道初始位置与力模型初
值相应的改正参数；
"U。、 4(U。、 A(U。为A弧段U寸刻轨道与力模型参数初值；
5 ("，54ft) ， ^Aft)为6弧段4时刻轨道初始位置与力模型初值相应的改正参数；
O(UJ是状态转移矩阵。
(B) 观测模型参数连接，与卫星连续运动相同，导航卫星的观测在正常情 况下也是连续的，由此相连两个弧段的观测模型参数也是关联的。观测模型参数 包括对流层参数、卫星钟差参数、接收机参数、地球自转参数、地面站坐标、模 糊度参数等，这些待估参数可归类为随机过程噪声参数与状态参数。对于状态参 数同样可以建立与上述卫星运动状态连接类似的关联方程，只是状态转移矩阵具 有不同的定义。对于随机过程噪声参数， 一般可以通过以下公式表示
HA+^ 式6 由于v^可以认为是一独立随机过程噪声，并服从正态分布，上式可以转换为
<formula>formula see original document page 8</formula> 式7
其中Q为相应过程噪声的标准差，考虑到^为零均值，^A可设为零，A、 A+1 为随机过程参数，M^^为映射矩阵，这样随机过程噪声参数与状态参数一样建 立类似的关联方程。
(C) 由此将上述两种数学模型统一起来表述，将不同的参数都用义表示，这 样关联条件则表示为
<formula>formula see original document page 8</formula> 式8
约束方程表述为
<formula>formula see original document page 8</formula>
式9
对于短弧轨道间不同参数的连接或合并，可以通过设置参数前后弧段的映射 因子o与动态噪声&来实现。实际上，o描述了函数关系模型上的连接关系，A
描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程 度，况为观测值残差，以下讨论不同类型参数相连弧段间连接时O与A的设置。
对于前后弧段相同或具有确定性的状态转移关系的参数，如卫星状态参数、测站坐标参数、连续的�：炔问龋诠斓懒哟κ窍嗤模梢陨鑎为无 穷大实现；
对于前后两个弧段不相关的分段参数或白噪声参数，如卫星钟差、接收机钟 差等，可以设0 = 0， &=0实现；
对于随机过程噪声参数，可以通过其过程噪声的标准差给定相应的&值实现。
(D)基于相连两个弧段各自的短弧法方程信息，与边界约束条件方程，推 导短弧合并的统一数学模型。
假设相连弧段+1对应的法方程为
<formula>formula see original document page 9</formula>式10
边界约束条件方程为
<formula>formula see original document page 9</formula> 式11
将上式io假设为先验约束，式ii作为虚拟观测方程，则轨道合并的问题就 可以转换为附有先验信息的平差问题，采用最小方差估计，要求满足性能函数
J(x)为最小
<formula>formula see original document page 9</formula> 式12
则^M-0可推导得到 <formula>formula see original document page 9</formula>简化表达，设<formula>formula see original document page 9</formula>上式13可简化为:<formula>formula see original document page 10</formula>式1
以上推导了两个相连弧段的轨道合并公式，下面直接给出一般性的w个弧段
的轨道合并公式
<formula>formula see original document page 10</formula>
式2
相比较现有技术而言，运用本方法实现导航卫星实时精密定轨具有如下优 点 一、处理速度快，由于法方程合并中大量的参数都可以预先消去，如除了最 后一个弧段外的其它弧段上的轨道参数、�：炔问龋庋喜⒌募扑懔烤� 可以大大减小，从而提高数据处理速度。二、快速初始化， 一旦实时系统需要重 新启动，此方法只需处理当前短弧段上的观测数据，然后与己经存在的法方程合 并就可以实时更新轨道，从而实现快速初始化。三、数值稳定，由于使用滑动数 据窗口的方法，每次轨道更新只用当前窗口中的数据，陈旧的观测信息被去除， 从而可以完全克服"数据饱和"导致的动力学模型参数估计精度降低甚至错误的 问题。四、能有效实现数据质量的控制，大大提高实时定轨的精度。
图1为本发明一种导航卫星实时精密定轨的方法中滑动数据窗口定轨示意
图2为本发明方法中短弧连接示意图。
图3为本发明方法中实施例试算每小时GPS实时轨道与IGS事后轨道比较 结果图。
具体实施方式
下面结合附图对本发明作进一步的描述。
如图1所示， 一种导航卫星实时精密定轨的方法，主要是将实时定轨通过 建立短弧法方程和短弧法方程合并两个并行进程实现，该方法包括以下步骤
(1) 将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段 的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为 N+T小时，其N+T的值不超过24小时；
(2) 处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数
据，建立该弧段上包括卫星初始状态、力模型参数、地球自转参数、大气参数、
测站坐标的短弧法方程，其短弧法方程的公式为7V^&^ =『^，其中5A为该弧 段的轨道初始位置与力模型初值相应的改正参数，7V^为法方程矩阵，k为短弧 段标记；
(3)短弧法方程合并进程将前N-1个弧段生成的法方程进行合并；其短弧 合并利用的公式如下 两个相连弧段的轨道合并公式
—①^
一尸O
「化-
式i
"个弧段的轨道合并公式
-尸①
(2,1)<formula>formula see original document page 11</formula>
式2
其中，O(UJ是状态转移矩阵，①描述了函数关系模型上的连接关系，Pw 描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度；
根据短弧轨道间参数类型的不同相连弧段间违接时0与^的设置分为以下 三类，
(a) 对于前后弧段相同或具有确定性的状态转移关系的参数，如卫星状态 参数、测站坐标参数、连续的�：炔问龋诠斓懒哟κ窍嗤模璩呶� 无穷大实现；
(b) 对于前后两个弧段不相关的分段参数或白噪声参数，如卫星钟差、接 收机钟差等，设0 = 0， &=0实现；
(c) 对于随机过程噪声参数，通过其过程噪声的标准差给定相应的&值实
现；
(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左右。
权利要求
1. 一种导航卫星实时精密定轨的方法，其特征在于该方法包括以下步骤(1)将卫星轨道按观测时间分割成连续等时长的短弧段，设定每个短弧段的时长为T小时，设定处理轨道数据的滑动数据窗口长度为N个短弧段时长为N*T小时；(2)处理当前短弧段，即滑动数据窗口内最后一个短弧段上获得的实时数据，生成该弧段上的短弧法方程，其短弧法方程的公式为Nk，kδxk＝Wkk，其中δxk为该弧段的轨道初始位置与力模型初值相应的改正参数，Nk，k为法方程矩阵，k为短弧段标记；(3)将滑动数据窗口内前N-1个短弧段生成的短弧法方程进行合并；其短弧合并利用的公式如下两个相连弧段的轨道合并公式式1n个相连弧段的轨道合并公式 式2其中，Φ(tk+1，tk)是状态转移矩阵；Φ描述了函数关系模型上的连接关系，Pw描述了随机模型统计特性上的连接，体现了前后弧段对应参数映射关系的约束程度；根据短弧轨道间参数类型的不同相连弧段间连接时Φ与Pw的设置分为以下三类，(a)对于前后弧段相同或具有确定性的状态转移关系的参数，在轨道连接处是相同的，设Pw为无穷大实现，(b)对于前后两个弧段不相关的分段参数或白噪声参数，设Φ＝0，Pw＝0实现，(c)对于随机过程噪声参数，通过其过程噪声的标准差给定相应的Pw值实现；(4)将当前短弧段的法方程与由前N-1个短弧段法方程合并成的一个法方程进行合并；(5)求解由整个滑动数据窗口内短弧合并而成的法方程，得到卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户；(6)将滑动窗口向下移动一个短弧段，返回步骤(2)。
2. 根据权利要求1所述的一种导航卫星实时精密定轨的方法，其特征是 步骤(1)中所述的滑动数据窗口的长度即N^^T的值不超过24小时。
3. 根据权利要求1所述的一种导航卫星实时精密定轨的方法，其特征是 步骤(2)中所述的短弧段上获得的实时数据包括卫星初始状态、力模型参数、 地球自转参数、大气参数、测站坐标。
全文摘要
本发明公开了一种基于滑动数据窗口的导航卫星实时精密定轨的方法，该方法的主要思想是在对状态X<sub>k</sub>，包括卫星初始状态、力模型参数、地球自转参数、大气参数、测站坐标等进行估计时，只利用离k时刻最近的前N个量测，而把之前其余的量测完全甩掉，这里N是预先规定的滑动窗口长度，再通过建立，合并，求解法方程等步骤获得精确的第N弧段的卫星轨道初值和力模型参数，通过轨道积分得到下一时段的预报轨道，作为实时轨道发送给用户。采用本方法进行卫星定轨具有数据处理速度快，轨道精度高，数值稳定等优点。
文档编号G01S1/04GK101435863SQ20081023742
公开日2009年5月20日 申请日期2008年12月25日 优先权日2008年12月25日
发明者刘经南, 唐卫明, 闯 施, 楼益栋, 葛茂荣, 赵齐乐 申请人:武汉大学
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
