2018, Vol. 38
GNSS卫星实时精密轨道有2种确定方式：事后预报和实时滤波。近年来，众多国内外学者对BDS卫星事后轨道进行了确定，目前精度为：GEO卫星为dm~m量级，IGSO和MEO卫星为10~30 cm量级[1-2]。与此同时，BDS卫星实时精密轨道估计也取得了一些研究成果[3-4]，其中事后预报6 h轨道，GEO卫星优于1 m，IGSO与MEO卫星优于20 cm[4]。BDS实时滤波轨道也达到了较高的精度水平[5]。
相比于事后预报方式，基于滤波的实时轨道估计具有较为灵活的操作性及系统异常实时检测能力(例如机动卫星的探测)[5]。但当动力学模型(尤其是光压模型)不准确时，卫星预报轨道精度会显著降低，从而严重影响实时轨道的可靠性[6]。由于BDS卫星采用动偏与零偏相结合的姿态控制策略[7]，传统光压模型(如ECOM模型)对于零偏姿控或者动偏与零偏姿态切换期间卫星定轨会带来较大的误差影响[8]，因此对于BDS卫星实时精密定轨，事后预报方式也存在一定的局限性。然而，滤波轨道往往在滤波启动后需要经过较长的收敛时间才能达到比较稳定的状态[5]，这极大降低了实时导航定位的用户体验。鉴于此，本文提出利用预报轨道约束的实时滤波精密定轨方法，并分析该方法在滤波轨道收敛期间的改善作用，评估BDS卫星实时滤波轨道精度。
1 数据处理模型与方法 1.1 BDS实时滤波精密定轨数学模型及轨道参数约束方程假定BDS卫星运动方程可表示为一阶微分方程：
| $ \begin{align} &{{{\mathit{\boldsymbol{\dot{X}}}}}^{\text{s}}}\text{=}\mathit{f}\left( {{\mathit{\boldsymbol{X}}}^{\text{s}}}\text{, }\mathit{k} \right) \\ &{{\left. {{\mathit{\boldsymbol{X}}}^{\text{s}}} \right|}_{\mathit{k}\text{-1}}}\text{=}\mathit{\boldsymbol{X}}_{_{\mathit{k}\text{-1}}}^{\text{s}} \\ \end{align} $ | (1) |
式中，
| $ \mathit{\boldsymbol{x}}{\rm{}}_\mathit{k}^{\rm{s}}{\rm{ = }}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{s}}}\left( {\mathit{k}{\rm{, }}\mathit{k}{\rm{ - 1}}} \right){\rm{\cdot}}\mathit{\boldsymbol{x}}_{_{\mathit{k}{\rm{ - 1}}}}^{\rm{s}} $ | (2) |
式中，xks表示k时刻相对于参考(积分)轨道的改正数，Φs(k, k-1)表示从(k-1)至k时刻的状态转移矩阵(或称时间更新矩阵)。
利用BDS卫星连续弧段的实测观测值即可实时确定卫星在每个历元的轨道状态参数。将原始相位与伪距观测值进行消电离层组合，并在参考轨道处线性化后得到误差方程，其在历元k时可表示为：
| $ \begin{array}{l} {\mathit{V}_{\mathit{P}_{\mathit{k}{\rm{, }}\mathit{r}{\rm{, if}}}^{\rm{s}}}}{\rm{ = - }}\mathit{\boldsymbol{l}}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}{\mathit{\boldsymbol{x}}_{\mathit{k}{\rm{, }}\mathit{r}}}{\rm{ + }}\mathit{\boldsymbol{l}}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}\mathit{\boldsymbol{x}}{\rm{}}_\mathit{k}^{\rm{s}}{\rm{ + }}\\ \mathit{c}{\rm{(}}{\mathit{t}_{\mathit{k}{\rm{, }}\mathit{r}}}{\rm{ - }}\mathit{t}_\mathit{k}^{\rm{s}}{\rm{) + }}\mathit{c}{\rm{(}}{\mathit{b}_{\mathit{k}{\rm{, }}\mathit{r}{\rm{, if}}}}{\rm{ - }}\mathit{b}_{\mathit{k}{\rm{, if}}}^{\rm{s}}{\rm{) + }}\mathit{M}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}{\mathit{T}_{\mathit{k}{\rm{, }}\mathit{r}}}\\ {\mathit{V}_{\mathit{L}_{\mathit{k}{\rm{, }}\mathit{r}{\rm{, if}}}^{\rm{s}}}}{\rm{ = - }}\mathit{\boldsymbol{l}}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}{\mathit{\boldsymbol{x}}_{\mathit{k}{\rm{, }}\mathit{r}}}{\rm{ + }}\mathit{\boldsymbol{l}}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}\mathit{\boldsymbol{x}}{\rm{}}_\mathit{k}^{\rm{s}}{\rm{ + }}\mathit{c}{\rm{(}}{\mathit{t}_{\mathit{k}{\rm{, }}\mathit{r}}}{\rm{ - }}\mathit{t}_\mathit{k}^{\rm{s}}{\rm{) + }}\\ \mathit{N}_{\mathit{k}{\rm{, }}\mathit{r}{\rm{, if}}}^{\rm{s}}{\rm{ + }}{\mathit{\varphi }_{\mathit{k}{\rm{, }}\mathit{r}{\rm{, if}}}}{\rm{ - }}\mathit{\varphi }_{\mathit{k}{\rm{, }}{\rm{ if}}}^{\rm{s}}{\rm{ + }}\mathit{M}_{\mathit{k}{\rm{, }}\mathit{r}}^{\rm{s}}{\mathit{T}_{\mathit{k}{\rm{, }}\mathit{r}}}\end{array} $ | (3) |
式中，VPk, r, ifs与VLk, r, ifs分别表示伪距与相位消电离层组合OMC，其中r表示测站编号，lk, rs表示站星方向余弦向量，xk, r表示测站坐标改正数，xks含义与式(2)相同，tk, r与tks分别表示接收机与卫星钟差，c为光速，bk, r, if与bk, ifs分别表示测站与卫星伪距偏差，Nk, r, ifs表示相位模糊度，φk, r, if与φk, ifs分别表示相位偏差，Tk, r表示天顶对流层延迟，Mk, rs表示对流层投影函数。
在参数估计中，由于接收机钟差与卫星钟差线性相关导致方程秩亏，因此需要提供一个基准钟，通常固定一个测站钟作为基准钟。伪距偏差与相位偏差将分别被钟差与模糊度吸收，因此无需处理[9]。
在滤波启动后的一段时间内，由于参数先验精度、几何结构等因素的限制，往往需要经过较长时间收敛其轨道参数才能达到较高的精度。考虑到轨道参数受动力学模型约束，其具有较强的规律性和较高的预报精度，因此可以利用事后解的预报轨道对轨道参数进行约束，以增强方程解的强度，其约束方程可以表示为：
| $ {\mathit{\boldsymbol{V}}_{\mathit{x}_\mathit{k}^{\rm{s}}}}{\rm{ = }}\mathit{\boldsymbol{x}}_\mathit{k}^{\rm{s}}{\rm{ - }}\mathit{\boldsymbol{\bar x}}_\mathit{k}^{\rm{s}} $ | (4) |
式中，xks表示先验预报轨道在k历元处的轨道状态参数，其方差为DVxks。由此可见，通过构建先验轨道约束方程可以使滤波轨道参数估值接近于先验轨道参数值，其接近程度由DVxks决定，从而达到改善滤波启动后轨道参数精度的目的。
1.2 先验轨道约束方法与步骤不同于事后批处理模式，基于滤波模式的定轨不估计参考时刻的轨道状态，而是逐历元更新卫星的轨道状态参数及动力学参数。在滤波启动后，利用预报轨道进行约束，主要分为3个步骤：1)利用一定时长(如6 h)的预报轨道进行动力学拟合，获取较为精确的初始轨道参数，并以此通过动力学积分得到先验参考轨道；2)在滤波启动之后的一定时长(如6 h)内，从先验参考轨道实时获取当前历元先验轨道状态参数，并利用式(4)在式(3)上添加约束方程，经过参数估计得到施加了约束的当前历元轨道状态参数估值；3)在滤波启动一定时长(如6 h)之后，取消先验轨道约束，由滤波器无约束运行。
2 BDS卫星实时精密定轨参数估计策略为了验证本文方法的有效性，利用多模GNSS试验跟踪网(multi-GNSS experiment，MGEX)全球分布的51个测站2015年年积日36~42的实测数据对BDS卫星的精密轨道进行实时滤波估计。测站分布及BDS卫星星下点轨迹见图 1。由于联合GPS观测数据有利于提高BDS卫星定轨精度[10]，本文同时估计GPS与BDS卫星轨道，因此需要在式(3)基础上增加估计系统间偏差参数。
|
图 1 51个测站分布及BDS卫星星下点轨迹 Fig. 1 The distribution of 51 selected MGEX stations and sub-satellite points of BDS satellites |
实时轨道估计时，采用平方根信息滤波(SRIF)逐历元更新轨道状态参数和光压模型参数，观测值采样率为30 s，卫星截止高度角为7°，待估参数估计策略见表 1。定轨中主要的动力学模型及误差修正模型选用IGS推荐的策略(http://acc.igs.org/reprocess2.html)，其中光压模型采用ECOM5参数模型[11]，BDS卫星姿态采用动偏和零偏切换模型[6]，电离层延迟则通过无电离层组合进行消除。
|
表 1 SRIF实时精密定轨参数估计策略 Tab. 1 The parameter estimation strategy of real-time precise orbit determination based on SRIF |
利用3 d(72 h)连续弧长进行事后轨道确定，并通过积分外推得到预报轨道。事后轨道精度评估时，将不同天定轨弧段最后24 h与中间24 h的重叠弧段部分求差，其差异的RMS可以用来评估事后定轨精度；预报轨道精度评估时，将预报弧段与事后3 d解中间24 h的重叠弧段部分进行求差，其差异的RMS可以用来评估预报轨道的精度。考虑到目前超快速星历更新频率为6 h，因此本文也主要统计6 h的轨道预报精度。C13卫星在2015年已经停止服务，因此未参与统计。
图 2为7 d时间内BDS 3 d解事后轨道实测部分24 h重叠弧段差异RMS。可以看出，BDS卫星事后解轨道具有较高的精度，其中3个方向上GEO卫星均优于30 cm，IGSO和MEO卫星均优于15 cm。图 3为7 d时间内BDS卫星3 d解事后预报部分6 h重叠弧段差异RMS。可以看出，3个方向上GEO卫星均优于80 cm，IGSO和MEO卫星均优于30 cm。BDS不同轨道类型卫星3 d解重叠弧段差异RMS的平均值见表 2。可以看出，BDS卫星轨道重叠弧段差异主要体现在切向上，这与其他GNSS系统具有相同的规律；IGSO和MEO卫星径向上的精度最高，法向次之，而GEO卫星切向和法向较其他方向明显偏大，这主要是由于GEO的静地特性导致其在切向和法向上几何变化较小所致。对于预报6 h轨道，IGSO与MEO卫星径向与法向重叠弧段与实测部分相当，说明这两类卫星预报精度较高; GEO卫星预报精度较实测部分有所衰减，但径向上依然小于20 cm。综上所述，BDS卫星6 h预报轨道具有较高的精度，因此在实际应用中可以将其作为比较精确的先验轨道对滤波启动后的轨道参数加以约束。
|
图 2 BDS 3 d解事后轨道实测部分24 h重叠弧段差异RMS Fig. 2 The orbit overlap RMS of 24 h estimation part of a 3-day solution obtained in a batch estimation |
|
图 3 BDS卫星3 d解事后预报部分6 h重叠弧段差异RMS Fig. 3 The orbit overlap RMS of 6 h prediction part of a 3-day solution obtained in a batch estimation |
|
|
表 2 BDS不同轨道类型卫星3 d解实测与预报轨道重叠弧段差异RMS的平均值 Tab. 2 The mean RMS of estimation part and prediction part orbit overlap for different orbit types of BDS satellites |
设置2套方案进行实验：方案1，利用广播星历为初始轨道；方案2，利用预报星历为初始轨道并逐历元添加轨道约束。方案1中，卫星位置、速度和光压参数先验精度分别设为10 m、10 mm/s以及0.1；方案2中，为使滤波轨道不受预报轨道过强影响，同时顾及参数解的强度，卫星位置、速度先验精度以及约束方程精度均设较为宽松的值，而光压参数则采用较强约束，其中GEO卫星位置和速度参数分别设为1.5 m、5 mm/s，IGSO和MEO卫星分别设为0.5 m、0.5 mm/s，光压参数均设为0.000 1。
图 4、5分别为方案1、2情况下GEO、IGSO与MEO每颗卫星解算的实时滤波轨道与3 d解事后轨道中间24 h三维位置差异时序。为了更直观地说明本文方法对提高滤波收敛期间轨道精度的效果，图中只显示滤波启动后72 h内的轨道结果。对比2幅图可以看出，方案1情况下，滤波轨道需要一定的时间收敛才能达到较稳定的精度水平。分别以6 m、1 m以及0.6 m作为GEO、IGSO与MEO卫星轨道收敛指标，GEO卫星需要约13 h收敛，IGSO卫星需要约18 h收敛，MEO卫星需要约14 h收敛，说明利用广播星历作为先验轨道需要较长的收敛时间，这极大增加了实时精密服务初始化时间。而方案2中，由于在滤波启动后6 h时间内添加了预报轨道约束，IGSO与MEO卫星未出现方案1中的长期收敛过程，并且在整个运行期间轨道变化较为稳定。这不仅说明本文提出的方法能够很好地改善滤波启动后十几h内的轨道精度，也反映出仅需6 h时间的约束就可以让轨道参数达到稳定状态，并在6 h之后无约束运行中都处于稳定的精度水平。然而，GEO卫星虽然有效改善了轨道收敛期间的估值精度，但轨道三维位置差异随时间推移有变大的趋势，这一方面是因为GEO卫星误差本身较大，添加约束方程使轨道估值结果偏向于预报轨道，而随着约束的取消以及滤波的推进会逐渐恢复到无约束情况下的误差大小；另一方面是由于GEO卫星在切向上与事后轨道存在系统性偏差，而文中约束方程(式(4))是建立在服从零均值正态分布的误差理论基础上，其无偏特性并不完全适用于GEO卫星切向系统性偏差的情况，具体表现为约束期间轨道过于偏向预报轨道，因此在下阶段研究中可以考虑增加估计切向上的常量偏差来补偿系统性偏差的影响。此外，从图中还可以发现，在跨天的时候(24 h及48 h)，滤波轨道时序出现了不连续的跳跃现象，这是由于事后轨道在跨天处不连续所致，这也体现了滤波轨道在提供实时连续轨道方面的优势。综上所述，本文方法充分利用预报轨道改善滤波收敛的同时，也保持了滤波轨道连续性、操作灵活性以及系统异常实时监测等优势。
|
图 4 利用广播星历为初始轨道的北斗卫星实时滤波轨道与3 d解事后轨道三维位置差异时序图 Fig. 4 The BDS satellites' 3D position difference between the 3-day solution obtained in a batch estimation and the real-time filtering orbit using the initial orbit parameter obtained from broadcast ephemeris |
|
图 5 利用本文方法的北斗卫星实时滤波轨道与3 d解事后轨道三维位置差异时序图 Fig. 5 The BDS satellites' 3D position difference between the 3-day solution obtained in a batch estimation and the real-time filtering orbit based on the method proposed in this paper |
图 6为7 d时间内本文方法的实时滤波轨道在切向、法向和径向上与事后轨道中间24 h差异的RMS。为更好地显示所有方向上的精度，未全刻度显示GEO卫星在切向上的值。可以看出，除GEO卫星切向外，其余卫星各个方向的RMS均优于30 cm(切向、法向和径向上的平均值见表 3)，说明本文方法可以精确确定BDS卫星实时滤波轨道。进一步观察发现，GEO卫星法向精度明显优于预报轨道，主要是因为轨道参数过程噪声的引入在一定程度上减小了动力学模型误差的影响；而IGSO和MEO卫星各个方向RMS较预报轨道略大，这可能与轨道参数过程噪声选取的合理性有关，并且与2种数据处理模式之间的差异也有一定的关系，其原因有待进一步深入研究。
|
图 6 利用本文方法的北斗卫星实时滤波轨道与3 d解事后轨道差异RMS Fig. 6 The RMS of orbit difference between the 3-day solution obtained in a batch estimation and the real-time filtering orbit based on the method proposed in this paper |
|
|
表 3 BDS不同轨道类型卫星实时滤波轨道与3 d解事后轨道差异RMS的平均值 Tab. 3 The mean RMS of the difference between the 3-day solution and the real-time filtering orbit for different orbit types of BDS satellites |
高精度的实时轨道是北斗系统提供实时精密服务的基础。利用滤波方式计算北斗卫星实时精密轨道是北斗系统实时精密服务领域的一个重要发展方向。由于参数先验精度、几何结构等因素限制，利用广播星历为初始轨道时，卫星轨道在滤波启动后需要较长时间收敛才能达到较稳定的精度水平。鉴于此，本文提出基于北斗精密预报轨道约束的实时滤波定轨方法，并利用MGEX跟踪网51个测站7 d的实测数据进行验证。实验结果表明：1)与事后轨道相比，BDS预报6 h轨道三维精度GEO卫星优于90 cm，IGSO卫星优于25 cm，MEO卫星优于20 cm；2)利用以广播星历为初始轨道的传统方式进行实时轨道计算时，BDS卫星在滤波启动后平均需要约15 h收敛才能达到较为稳定的精度水平，而本文方法在这段时间内轨道变化较为平稳，未出现明显的收敛现象，说明本文方法能够有效解决传统方法轨道收敛时间长的问题；3)与事后轨道相比，实时滤波轨道GEO卫星径向上的RMS为21 cm，IGSO卫星为10 cm，MEO卫星为5 cm，说明本文方法可以精确确定BDS卫星的实时轨道，而通过不同类型卫星选取合理的轨道参数过程噪声有待进一步研究。
| [1] |
Guo J, Xu X L, Zhao Q L, et al. Precise Orbit Determination for Quad-Constellation Satellites at Wuhan University: Strategy, Result Validation, and Comparison[J]. Journal of Geodesy, 2016, 90(2): 143-159 DOI:10.1007/s00190-015-0862-9
(  0) |
| [2] |
Guo F, Li X X, Zhang X H, et al. Assessment of Precise Orbit and Clock Products for Galileo, Beidou, and QZSS from IGS Multi-GNSS Experiment (MGEX)[J]. GPS Solutions, 2017, 21(1): 279-290 DOI:10.1007/s10291-016-0523-3
( 0) |
| [3] |
Lou Y D, Liu Y, Shi C, et al. Precise Orbit Determination of Beidou Constellation: Method Comparison[J]. GPS Solutions, 2016, 20(2): 259-268 DOI:10.1007/s10291-014-0436-y
( 0) |
| [4] |
Li X X, Ge M R, Dai X L, et al. Accuracy and Reliability of Multi-GNSS Real-Time Precise Positioning: GPS, GLONASS, Beidou, and Galileo[J]. Journal of Geodesy, 2015, 89(6): 607-635 DOI:10.1007/s00190-015-0802-8
( 0) |
| [5] |
戴小蕾.基于平方根信息滤波的GNSS导航卫星实时精密定轨理论与方法[D].武汉: 武汉大学, 2016 (Dai Xiaolei. Real-Time Precise GNSS Satellite Orbit Determination Using the SRIF Method: Theory and Implemencation[D]. Wuhan : Wuhan University, 2016) http://cdmd.cnki.com.cn/Article/CDMD-10486-1016124266.htm
( 0) |
| [6] |
Rodriguez-Solano C J, Hugentobler U, Steigenberger P, et al. Improving the Orbits of GPS Block ⅡA Satellites during Eclipse Seasons[J]. Advances in Space Research, 2013, 52(8): 1 511-1 529 DOI:10.1016/j.asr.2013.07.013
( 0) |
| [7] |
Dai X L, Ge M R, Lou Y D, et al. Estimating the Yaw-Attitude of BDS IGSO and MEO Satellites[J]. Journal of Geodesy, 2015, 89(10): 1 005-1 018 DOI:10.1007/s00190-015-0829-x
( 0) |
| [8] |
郭靖.姿态、光压和函数模型对导航卫星精密定轨影响的研究[D].武汉: 武汉大学, 2014 (Guo Jing. The Impact of Attitude, Solar Radiation and Function Model on Precise Orbit Determination for GNSS Satellites[D]. Wuhan: Wuhan University, 2014) http://cdmd.cnki.com.cn/Article/CDMD-10486-1015528772.htm
( 0) |
| [9] |
范磊, 钟世明, 李子申, 等. 跟踪站分布对非组合精密单点定位提取GPS卫星差分码偏差的影响[J]. 武汉大学学报:信息科学版, 2016, 41(3): 316-321 (Fan Lei, Zhong Shiming, Li Zishen, et al. Effect of Tracking Stations Distribution on the Estimation of Differential Code Biases by GPS Satellites Based on Uncombined Precise Point Positioning[J]. Geomatics and Information Science of Wuhan University, 2016, 41(3): 316-321)
( 0) |
| [10] |
张睿, 杨元喜, 张勤, 等. BDS/GPS联合定轨的贡献分析[J]. 武汉大学学报:信息科学版, 2017, 42(5): 600-608 (Zhang Rui, Yang Yuanxi, Zhang Qin, et al. Contribution Analysis of BDS/GPS Combined Orbit Determination[J]. Geomatics and Information Science of Wuhan University, 2017, 42(5): 600-608)
( 0) |
| [11] |
Springer T A, Beutler G, Rothacher M. A New Solar Radiation Pressure Model for GPS[J]. Advances in Space Research, 1999, 23(4): 673-676 DOI:10.1016/S0273-1177(99)00158-1
( 0) |