Real-time orbit determination of Beidou satellite based on square root information filtering
-
摘要: 连续、稳定、高精度的实时卫星轨道产品是北斗国际化、规模化、智能化应用的重要前提. 当前，北斗卫星导航系统(BDS)的实时精密轨道产品多基于“批处理解算+轨道预报”的超快速模式获得，存在连续性较差、稳定性较低、精度不高等问题. 为此，本文采用平方根信息滤波(SRIF)方法对北斗卫星精密轨道进行实时逐历元解算. 实验结果表明：相比于超快速定轨模式，基于实时滤波方法的轨道产品能够有效避免边界跳变，具有更好的连续性和稳定性；同时，实时滤波定轨方法能够显著提高BDS的轨道精度，其中中轨道地球卫星(MEO)和倾斜地球同步轨道卫星(IGSO)的三维轨道误差分别减小了46%和68%，卫星激光测距(SLR)检核精度也普遍优于预报轨道.
-
关键词:
- 实时精密定轨 /
- 北斗卫星导航系统(BDS) /
- 平方根信息滤波 /
- 轨道精度 /
- 卫星激光测距(SLR)检核
Abstract: Continuous, stable and high-precision real-time satellite orbit products are essential for the international, scale and intelligent application of the BeiDou system. Current BDS real-time precise orbit products, which rely on “batch processing solution + orbit prediction”, suffer from poor continuity, low stability and low accuracy. Therefore, we introduce the square root information filtering method to perform real-time epoch-by-epoch solutions for the precise orbit of BeiDou satellites. Experimental results demonstrate that the BDS orbits generated by real-time filtering method outperform the ultra-rapid orbit products by effectively avoiding boundary jumps, improving continuity and stability. Moreover, the real-time filtering precise orbit determination method can significantly improve the orbit accuracy of BeiDou satellites, reducing the 3D orbit errors of MEO and IGSO satellites by 46% and 68%, respectively. Also, the validation by satellite laser ranging shows that the orbit accuracy of real-time filtering method is generally better than that of the predicted orbit. -
0. 引 言
北斗卫星导航系统(BDS)是由我国自主建设、独立运行的卫星导航系统，是为全球用户提供全天候、全天时、高精度的定位、导航和授时(PNT)服务的国家重要基础设施. 按照“三步走”建设发展战略，BDS先后经历了实验验证系统(北斗一号)和区域导航系统(北斗二号)阶段[1-2]，并最终于2020年7月完成了北斗三号全球卫星导航系统(BDS-3)的建设. 目前，北斗系统共有45颗卫星在轨，包括15颗北斗二号(BDS-2)和30颗BDS-3.
北斗卫星精密轨道产品提供了导航定位的空间基准，是北斗系统应用的前提条件之一，其精度与可靠性直接影响北斗服务性能. 随着北斗正式迈入全球服务新时代，以“北斗+”和“+北斗”深度融合为特征的导航定位产业进入到全新发展阶段[3]. 物联网、智慧城市、自动驾驶等新兴技术产业对北斗导航定位服务的实时性和精确性提出了更高的要求，对于高精度实时轨道产品的需求愈加迫切. 目前，实时精密轨道主要采用基于事后批处理解算的轨道预报模式获得[4]. 欧洲定轨中心(CODE)、德国地学中心(GFZ)、武汉大学IGS分析中心和多家国际GNSS监测评估系统(iGMAS)分析中心均基于这种方法提供了GNSS卫星的超快速轨道产品，其更新间隔由1 h至6 h不等. 然而，这种“批处理解算+轨道预报”的超快速模式仍然存在一些问题难以解决. 首先，由于卫星非保守力模型不完善的影响，预报方式获得的轨道产品存在精度损失；其次，超快速轨道采用了分弧段计算的方式，导致弧段边界处存在轨道不连续. 此外，BDS采用了地球静止轨道卫星(GEO)+倾斜地球同步轨道卫星(IGSO)+中轨道地球卫星(MEO)的混合星座，轨道机动更加频繁，而超快速定轨模式无法有效处理轨道机动与轨道异常，实时轨道存在可靠性与可用性风险.
与超快速定轨模式不同，基于滤波模式的卫星定轨方法具有连续解算、实时更新、机动异常处理能力强等显著优势，逐渐得到了国内外学者的关注. 英国纽卡斯尔大学采用扩展卡尔曼滤波(EKF)对GPS卫星进行了实时定轨，取得了收敛后三维方向14 cm的轨道精度[5]；法国国家空间研究中心(CNES)则基于均方根信息滤波(SRIF)得到了一维方向3.7 cm的GPS卫星实时定轨精度[6]. 在国内，Dai等[7]开发了基于SRIF的导航卫星实时精密定轨系统，对GPS和GLONASS卫星分别可以取得6.7 cm和9.3 cm的三维轨道精度；匡开发[8]则使用并行计算技术对实时滤波定轨的解算效率进行了优化，在2 s内即可完成GPS单系统轨道解算. 但目前使用实时滤波模式对北斗卫星，尤其是BSD-3的精密轨道确定的研究较少.
因此，本文采用SRIF方法，对北斗卫星精密轨道进行了实时逐历元解算. 结构安排如下：首先给出北斗实时滤波定轨模型与方法；接着介绍实验数据和处理策略；再在此基础上，分析北斗滤波轨道的时序特性和精度情况；最后对论文进行总结.
1. 北斗实时滤波定轨模型与方法
1.1 北斗实时定轨函数模型
BDS提供了丰富的伪距和载波相位观测值，本文使用的观测值组合为双频非差消电离层组合(IF). 为了与轨道动力学所用参考系保持一致，这里给出地心惯性坐标系(ECI)下线性化后的观测方程[9]：
$$ \begin{split} l_{{\rm{r}},{\rm{lc}}}^{\rm{s}} =& {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}^{\rm{s}}} - {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{R}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}_{\rm{r}}} - {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{d}}{\boldsymbol{r}}{{\boldsymbol{d}}_{{\text{eop}}}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}_{{\text{eop}}}} -\\& c{t^{\rm{s}}} + c{t_{\rm{r}}} + {\lambda _{{\rm{lc}}}} \cdot {B_{{\rm{lc}}}} + m_{\rm{r}}^{\rm{s}} \cdot {T_{\rm{r}}} + \varepsilon _{{\rm{r}},{\rm{lc}}}^{\rm{s}}, \\ l_{{\rm{r}},{\rm{pc}}}^{\rm{s}} =& {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}^{\rm{s}}} - {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{R}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}_{\rm{r}}} - {\boldsymbol{u}}_{\rm{r}}^{\rm{s}} \cdot {\boldsymbol{d}}{\boldsymbol{r}}{{\boldsymbol{d}}_{{\text{eop}}}} \cdot {\boldsymbol{\Delta }}{{\boldsymbol{x}}_{{\text{eop}}}} -\\& c{t^{\rm{s}}} + c{t_{\rm{r}}} + m_{\rm{r}}^{\rm{s}} \cdot {T_{\rm{r}}} + \varepsilon _{{\rm{r}},{\rm{pc}}}^{\rm{s}} . \end{split} $$ (1) 式中：
${l}_{{\rm{r}},{\rm{lc}}}^{{\rm{s}}}$ 和${l}_{{\rm{r}},{\rm{pc}}}^{{\rm{s}}}$ 分别为相位和伪距观测值减去计算初值的残差；上标${\rm{s}}$ 和下标${\rm{r}}$ 分别为观测值对应的卫星和测站；${\boldsymbol{u}}_{{\rm{r}}}^{{\rm{s}}}$ 为测站指向卫星的单位方向向量；${\boldsymbol{\Delta }{{\boldsymbol{x}}}}_{{\rm{r}}}$ 和${{\boldsymbol{\Delta}} {{\boldsymbol{x}}}}^{\mathrm{s}}$ 分别为测站位置和卫星轨道状态初始向量的改正量；矩阵$\boldsymbol{R}$ 为地心地固坐标系(ECEF)到ECI坐标系的转换矩阵；${\boldsymbol{d}}{\boldsymbol{r}}{\boldsymbol{d}}_{\mathrm{e}\mathrm{o}\mathrm{p}}$ 为测站坐标对地球方向参数(EOP)的导数；${{\boldsymbol{\Delta}} \boldsymbol{x}}_{\mathrm{e}\mathrm{o}\mathrm{p}}$ 为EOP初始向量的改正量；${t}^{{\rm{s}}}$ 和${t}_{{\rm{r}}}$ 分别为卫星和测站的钟偏差；$ c $ 为光速；${T}_{{\rm{r}}}$ 为修正后的天顶对流层延迟；${m}_{{\rm{r}}}^{{\rm{s}}}$ 为投影函数[10]；${\lambda }_{{\rm{lc}}}$ 和${B}_{{\rm{lc}}}$ 分别为相位观测值IF组合的波长和模糊度；${\varepsilon }_{{\rm{r}},{\rm{lc}}}^{{\rm{s}}}$ 和${\varepsilon }_{{\rm{r}},{\rm{pc}}}^{{\rm{s}}}$ 分别为相位和伪距的观测噪声.BDS卫星轨道状态参数由卫星位置、速度和太阳辐射压参数(SRP)构成. 对于逐历元更新的北斗滤波定轨，需要实时构建下一历元的卫星轨道状态初值及前后历元间的状态转移矩阵. 其求解可通过数值积分中的单步法和多步法进行：首先基于前一历元的初始轨道参数，利用龙格-库塔(Runge-Kutta)单步积分法获取足够的起算点；在此基础上，采用Adams多步积分法完成后续弧段上积分点的计算[11-12]. 相应的轨道状态转移方程可表示为：
$$ \begin{split} &{\boldsymbol{X}}_i^{{\rm{orb}}} \to {\boldsymbol{X}}{_{i + 1}^{{\rm{orb}}*}}, \\& {\boldsymbol{x}}_{i + 1}^{{\rm{orb}}} = {{\boldsymbol{\varPhi }}_{i + 1,i}} \cdot {\boldsymbol{x}}_i^{{\rm{orb}}} + {{\boldsymbol{w}}_i} , \\& {\boldsymbol{X}}_{i + 1}^{{\rm{orb}}} = {\boldsymbol{X}}_i^{{\rm{orb}}*} + {\boldsymbol{x}}_i^{{\rm{orb}}}. \end{split} $$ (2) 式中：
${\boldsymbol{X}}_{i}^{{\rm{orb}}}$ 为$ {t}_{i} $ 时刻卫星轨道状态最优估计值；${\boldsymbol{X}}_{i+1}^{{\rm{orb}}*}$ 为积分得到的$ {t}_{i+1} $ 时刻的卫星轨道状态初值；${\boldsymbol{x}}_{i}^{{\rm{orb}}}$ 和${\boldsymbol{x}}_{i+1}^{{\rm{orb}}}$ 为$ {t}_{i} $ 和$ {t}_{i+1} $ 时刻待估的卫星轨道状态初值的改正量；矩阵${\boldsymbol{\varPhi }}_{i+1,i}$ 为$ {t}_{i} $ 时刻到$ {t}_{i+1} $ 时刻的状态转移矩阵；${\boldsymbol{w}}_{i}$ 为轨道状态转移方程的过程噪声.1.2 SRIF方法
计算机进行数值计算过程中存在截断误差，导致传统卡尔曼滤波容易因数值计算误差而发散，因此扩展出了一些其他形式的平方根滤波算法. 其中，SRIF因其具有数值稳定性和计算高效性，在导航及控制领域被广泛采用[13]. 相较于传统卡尔曼滤波需不断更新状态量及其协方差，SRIF则是维护了基于平方根信息矩阵构成的信息方程.
SRIF算法中包含时间更新和量测更新两个部分. 假设ti时刻的先验信息方程和观测方程分别为
${{\boldsymbol{\tilde z}}_i} = {{\boldsymbol{\tilde R}}_i}{{\boldsymbol{x}}_i} + {{\boldsymbol{\tilde v}}_i}$ 和${{\boldsymbol{z}}_i} = {{\boldsymbol{A}}_i}{\boldsymbol{x}} + {{\boldsymbol{v}}_i}$ ，式中：$[\tilde{\boldsymbol{R}}_{i}\; \tilde{\boldsymbol{z}}_{i}]$ 为先验的平方根信息矩阵；${\tilde{{\boldsymbol{v}}}}_{i}$ 为方程噪声；${\boldsymbol{z}}_{i}$ 为观测向量；${\boldsymbol{v}}_{i}$ 为对应的观测噪声；${\boldsymbol{A}}_{i}$ 为系数矩阵. 则相应量测更新过程可表示为：$$ \begin{split} &{\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\tilde R}}}_i}} \\ {{{\boldsymbol{A}}_i}} \end{array}} \right]{\boldsymbol{x}} = {\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\tilde z}}}_i}} \\ {{{\boldsymbol{z}}_i}} \end{array}} \right] - {\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\tilde v}}}_i}} \\ {{{\boldsymbol{v}}_i}} \end{array}} \right] \\&\Rightarrow \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat R}}}_i}} \\ {\boldsymbol{0}} \end{array}} \right] {\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat z}}}_i}} \\ {{{\boldsymbol{e}}_i}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat v}}}_i}} \\ {{{\boldsymbol{v}}_e}} \end{array}} \right] . \end{split}$$ (3) 式中：矩阵
$\boldsymbol{T}$ 为QR分解得到的正交变换矩阵；${\hat{\boldsymbol{R}}}_{i}$ 和${\hat{\boldsymbol{z}}}_{i}$ 构成了量测更新后的信息方程；${\boldsymbol{e}}_{i}$ 为正交变换后的残余项. 在北斗实时定轨中，根据历元间的时变特性，状态参数$\boldsymbol{x}$ 可划分为确定性参数$\boldsymbol{y}$ 和不确定性参数${\boldsymbol{p}}_{i}$ . 则式(2)中的状态转移方程可进一步转换为$$ \begin{gathered} {{\boldsymbol{R}}_{\boldsymbol{w}}}{{\boldsymbol{p}}_{i + 1}} = {{\boldsymbol{R}}_{\boldsymbol{w}}}{{\boldsymbol{\varPhi }}_{i + 1,i}}{{\boldsymbol{p}}_i} + {{\boldsymbol{R}}_{\boldsymbol{w}}}{\boldsymbol{w}} {\boldsymbol{R}}_{\boldsymbol{w}}^{ - 1}{\boldsymbol{R}}_{\boldsymbol{w}}^{ -{\rm{ T}}} = {\boldsymbol{D}}({\boldsymbol{w}}) . \end{gathered}$$ (4) 式中，
$\boldsymbol{w}$ 和$\boldsymbol{D}\left(\boldsymbol{w}\right)$ 为过程噪声及其协方差. 结合量测更新后的信息方程和式(4)，可以得到时间更新的表达式为：$$ \begin{array}{c} {\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} { - {{\boldsymbol{R}}_{\boldsymbol{w}}}{{\boldsymbol{\varPhi }}_{i + 1,i}}}&{{{\boldsymbol{R}}_{\boldsymbol{w}}}}&{\bf{0}}\;\;\;\;\;\; \\ {{{{\boldsymbol{\hat R}}}_{{{\boldsymbol{p}}_i}}}}&{\bf{0}}&{{{{\boldsymbol{\hat R}}}_{{{\boldsymbol{p}}_i},{\boldsymbol{y}}}}} \\ {\bf{0}}&{\bf{0}}&{{{{\boldsymbol{\hat R}}}_{\boldsymbol{y}}}} \;\;\;\; \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_i}} \\ {{{\boldsymbol{p}}_{i + 1}}} \\ {\boldsymbol{y}} \end{array}} \right] = {\boldsymbol{T}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{R}}_{\boldsymbol{w}}}{\boldsymbol{w}}} \\ {{{{\boldsymbol{\hat z}}}_{{{\boldsymbol{p}}_i}}}} \\ {{{{\boldsymbol{\hat z}}}_{\boldsymbol{y}}}} \end{array}} \right] \\\Rightarrow \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\tilde R}}}_{{{\boldsymbol{p}}_i}}}}&{{{{\boldsymbol{\tilde R}}}_{{{\boldsymbol{p}}_i},{{\boldsymbol{p}}_{i + 1}}}}}&{\bf{0}}\;\;\;\;\;\;\; \\ {\bf{0}}&{{{{\boldsymbol{\tilde R}}}_{{{\boldsymbol{p}}_{i + 1}}}}}\;\;\;\;&{{{{\boldsymbol{\tilde R}}}_{{{\boldsymbol{p}}_{i + 1}},{\boldsymbol{y}}}}} \\ {\bf{0}}&{\bf{0}}\;\;\;\;\;\;\;\;&{{{{\boldsymbol{\tilde R}}}_{\boldsymbol{y}}}} \;\;\;\;\;\; \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_i}} \\ {{{\boldsymbol{p}}_{i + 1}}} \\ {\boldsymbol{y}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\tilde z}}}_{\boldsymbol{w}}}} \\ {{{{\boldsymbol{\tilde z}}}_{{{\boldsymbol{p}}_i}}}} \\ {{{{\boldsymbol{\tilde z}}}_{\boldsymbol{y}}}} \end{array}} \right]. \end{array}$$ (5) 式中，参数
${\boldsymbol{p}}_{i}$ 可以直接消除以减少后续计算负担，而参数${\boldsymbol{p}}_{i}$ 和$\boldsymbol{y}$ 所对应的正交变换后的平方根信息矩阵则构成了时间更新后的信息方程. 至此，SRIF完成了当前历元的解算，通过循环上述过程不断更新相应的信息方程.2. 实验数据及处理策略
为充分测试基于SRIF方法解算得到的北斗实时轨道精度，本文基于GREAT软件[14]解算了2021年年积日第106—113天的地面测站的北斗观测数据. 其中，选取了全球分布的125个MGEX(Multi-GNSS Experiment)测站，具体分布如图1所示.
表1给出了逐历元解算北斗滤波轨道时的相关模型及处理策略. 由于处理卫星中同时包含BDS-2和BDS-3卫星，相应的IF组合观测值所采用的频率为B1I/B3I. 对于北斗轨道状态参数，这里采用随机游走模型进行估计并设置相应的经验性模型噪声. 而对于非保守力模型，在轨道状态参数估计时主要考虑了太阳光压模型，并采用了ECOM五参数模型，其余非保守力模型则暂不考虑；对于其他保守力模型，则采用了常用的模型进行计算[15]. 另外，为最大程度降低实时模糊度错误固定对轨道精度的影响，采用了对所有模糊度双差基线进行松约束固定的策略[16]，并将固定约束仅作用在当前历元而非传递到后续历元.
表 1 实时滤波轨道解算时的处理策略相关内容 处理策略 观测值 MGEX观测数据；伪距和相位；非差IF组合 观测频率选用 B1I+B3I 观测处理间隔 300 s 截止高度角 7° 观测先验 伪距0.6 m；相位0.02周 对流层干延迟 采用Saastamonien模型和GMF的投影函数 电离层延迟 采用IF组合消除，不考虑高阶延迟 天线相位中心改正 采用igs14.atx天线文件改正 卫星天线相位缠绕 模型改正 地球潮汐影响 对其中固体潮汐，海洋潮汐和极移潮汐进行改正 相对论效应 模型改正 地球重力场 采用了12×12阶次的EGM08模型 N体引力 考虑月亮及太阳系内行星的影响；采用JPL DE405计算星体位置 地球潮汐 采用IERS2010模型和FES2004模型计算 相对论效应 采用IERS2010模型计算 太阳光压模型 采用ECOM2模型 卫星位置/速度参数 采用随机游走模型估计；模型噪声为10−12 卫星太阳光压参数 采用随机游走模型估计；模型噪声为10−12 测站坐标参数 采用IGS的SNX文件进行强约束 测站钟差/卫星钟差参数 采用白噪声模型估计 模糊度固定参数 采用常数模型进行估计 对流层参数 采用分段常数模型估计；时段间隔为2 h 系统间偏差参数 采用分段常数模型估计；时间间隔为24 h 地球自转定向参数 采用分段线性模型估计；时段间隔为24 h 模糊度固定策略 松约束固定所有双差基线 最短卫星共视时长 900 s 宽巷模糊度平滑处理时间间隔 30 s 双差基线长度阈值 3500 km 3. 结果分析
为进一步验证北斗实时滤波轨道的精度水平，本文选用了武汉大学IGS分析中心生成的同时段北斗超快速轨道(Ultra-Rapid)产品作为对比，并将其事后精密产品(Final)作为参考轨道. 两者均与参考轨道进行了轨道比较以及采用SLR数据进行检核.
3.1 轨道时序特性分析
为分析北斗实时滤波轨道收敛期间及收敛后的时序特性，图2和图3给出了北斗实时滤波轨道浮点解(Float)及固定解(Fix)在切向(A)、法向(C)和径向(R)上与参考轨道互差结果的时序图. 可以看到，所有卫星基本在一天内完成了收敛，同时由于实时模糊度固定约束设定为滤波开始后24 h进行，因此浮点解和固定解在收敛期间具有相同的收敛情况. 这里进一步统计和分析不同类型卫星在不同方向上的收敛时间，其中收敛定义具体如下：对于MEO卫星，其A、C和R上的轨道互差绝对值在连续12 h内分别小于25 cm、20 cm和15 cm；对于IGSO卫星，其A、C和R上的轨道互差绝对值小于25 cm[17]. 表2给出了北斗实时滤波轨道浮点解的收敛时间统计结果. MEO卫星在A和C上的收敛时间接近，分别为13.3 h和12.9 h，而R上的收敛时间则最长，为17.9 h. 对于IGSO卫星，三个方向上的收敛时间则分别为18.8 h、22.0 h和26.6 h，其轨道类型决定了观测几何构型相较MEO卫星变化更慢，因此整体收敛时间也更长. 同时，由于R上的观测几何构型相较于其他两个方向上更差，无论是IGSO还是MEO卫星，R方向上的收敛都需要更长的时间.
表 2 北斗卫星实时滤波轨道浮点解的平均收敛时间统计结果h 卫星类型 A C R BDS-MEO 13.3 12.9 17.9 BDS-IGSO 18.8 22.0 26.6 对于收敛后的北斗实时滤波轨道浮点解的时序结果，MEO卫星在A、C和R上的变化幅度基本为±20 cm、±15 cm和±10 cm，而IGSO卫星在A、C和R上的变化幅度则基本为±25 cm. 对于收敛后的固定解时序结果，其变化幅度略小于浮点解结果. 同时由于前述的模糊度固定策略，固定解时序结果中包含了更多的锯齿状变化. 而无论是浮点解和固定解，所有卫星在天边界处都包含了不同程度的跳变，其主要由参考轨道所有的天边界不连续性导致. 进一步地，图4给出了超快轨道产品的预报弧段与WUM产品互差的时序图. 由于超快轨道分时段更新的特性，因此互差结果存在不同程度的周期性跳变，同时由于IGSO卫星轨道特性，其展现了比MEO卫星更为显著的跳变现象. 而在同时段内，逐历元解算得到的实时滤波轨道结果则为连续变化的，且不存在大幅度跳变. 这也侧面验证了实时滤波轨道相较于超快速轨道具有更好的轨道连续性.
3.2 轨道精度分析
为了评定北斗实时滤波轨道的精度水平，这里分别统计了在收敛时段内(2021年年积日第107—113天)北斗超快速轨道、实时滤波轨道浮点解和固定解与参考轨道产品在A、C、R和三维(3D)方向上轨道互差的均方根(RMS)值. 图5与图6分别给出了MEO卫星和IGSO卫星的轨道比较结果. 相较于A和C，北斗卫星轨道R上受到了更多动力学模型的约束，因此具有更好的精度水平. 相较于超快轨道产品，实时滤波轨道浮点解和固定解在R和C上的精度水平有明显改善，且对于IGSO卫星改善幅度更大. 表3中按卫星类型统计了不同方案轨道各方向上所有卫星的平均RMS值. 对于北斗实时滤波轨道固定解，MEO卫星在C、A和R上的平均RMS值分别为5.5 cm、8.5 cm和4.8 cm，而IGSO卫星在C、A和R上的平均RMS值分别为10.5 cm、12.4 cm和9.9 cm. 受限于动力学模型精度及其本身具有的轨道特性，IGSO的轨道整体精度要低于MEO卫星. 而在3D方向上，MEO卫星和IGSO卫星在3D方向上的实时滤波轨道固定解平均RMS值分别为11.4 cm和19.2 cm，相较于超快速轨道产品分别提升了46%和68%，精度水平有了明显改善.
表 3 北斗超快速轨道产品以及实时滤波轨道浮点解和固定解与WUM产品轨道比较RMS均值统计结果cm 卫星
类型Ultra-rapid SRIF-浮点解 SRIF-固定解 A C R 3D A C R 3D A C R 3D MEO 17.9 7.6 8.3 21.3 11.9 6.5 4.7 14.4 8.5 5.5 4.8 11.4 IGSO 49.0 18.7 31.0 61.4 16.8 13.6 9.3 23.5 12.4 10.5 9.9 19.2 3.3 卫星激光测距(SLR)结果统计与分析
对上述北斗卫星轨道的三种方案进一步采用SLR数据进行检核[18]. 图7给出了实时滤波轨道固定解、超快轨道产品和WUM产品的SLR检核残差时序图. 对于IGSO卫星(C08、C10和C13)实时滤波轨道与WUM产品检核残差普遍分布较为集中，而超快速轨道检核残差则存在较大的跳变，这和前面所展示的超快速轨道产品中IGSO卫星具有显著的周期性跳变相吻合；对于MEO卫星，三种类型的轨道检核残差则分布相对接近，其中的实时滤波轨道检核残差的标准差要略优于超快速轨道. 表4则进一步对实时滤波轨道、超快轨道产品和WUM产品SLR检核残差的均值、标准差及均方根误差（RMSE）进行了统计. 相较于超快速轨道产品，北斗实时滤波轨道中的IGSO卫星的标准差有了明显的降低，侧面反映了其更好的连续性. 对于所有卫星，WUM产品具有最小的RMSE，精度最优；而实时滤波轨道的RMSE普遍小于超快速轨道产品，与表3中实时滤波轨道精度比超快速轨道产品更优的结果相一致.
表 4 北斗卫星轨道SLR检核统计结果cm 卫星编号 轨道类型 实时滤波轨道 超快速轨道 WUM产品 均值 标准差 RMSE 均值 标准差 RMSE 均值 标准差 RMSE C08 IGSO 3.2 5.4 6.3 −21.5 20.6 29.8 6.5 5.0 8.2 C10 IGSO −1.7 11.5 11.6 −10.4 9.8 14.3 −1.0 5.2 5.3 C13 IGSO 6.0 10.1 11.7 −18.3 19.2 26.5 6.6 2.5 7.1 C11 MEO 0.1 4.0 4.0 3.8 3.5 5.2 1.1 2.5 2.7 C20 MEO 6.5 3.4 7.3 6.8 3.0 7.4 5.4 1.8 5.7 C21 MEO 5.4 5.5 7.7 5.6 3.6 6.6 4.8 2.1 5.3 C29 MEO −4.6 5.7 7.3 −5.9 5.3 8.0 2.7 3.0 4.0 C30 MEO −1.8 5.9 6.1 −3.5 6.0 7.0 3.6 3.8 5.2 4. 总 结
本文采用SRIF方法对26颗北斗卫星进行了实时定轨解算，并对滤波轨道的时序特性与精度进行了全面评估. 轨道时序特性分析的结果表明，BDS MEO卫星能够在滤波开始20 h内完成收敛，IGSO卫星也能够在24 h左右达到收敛. 相比于超快速轨道，实时滤波轨道的连续性和稳定性更好，有效避免了弧段边界的轨道跳变问题. 进一步统计了滤波收敛后的轨道精度情况. 结果表明，滤波定轨模式下，BDS MEO和IGSO卫星能够分别取得11.4 cm和19.2 cm的三维轨道精度，相比于超快速产品，精度提升分别可达到46%和68%. 在此基础上，通过SLR对滤波轨道进行了外部检核，结果也表明，实时滤波轨道的精度要普遍优于超快速产品. 本文的研究对于提升北斗系统的实时高精度服务性能具有一定意义.
-
表 1 实时滤波轨道解算时的处理策略
相关内容 处理策略 观测值 MGEX观测数据；伪距和相位；非差IF组合 观测频率选用 B1I+B3I 观测处理间隔 300 s 截止高度角 7° 观测先验 伪距0.6 m；相位0.02周 对流层干延迟 采用Saastamonien模型和GMF的投影函数 电离层延迟 采用IF组合消除，不考虑高阶延迟 天线相位中心改正 采用igs14.atx天线文件改正 卫星天线相位缠绕 模型改正 地球潮汐影响 对其中固体潮汐，海洋潮汐和极移潮汐进行改正 相对论效应 模型改正 地球重力场 采用了12×12阶次的EGM08模型 N体引力 考虑月亮及太阳系内行星的影响；采用JPL DE405计算星体位置 地球潮汐 采用IERS2010模型和FES2004模型计算 相对论效应 采用IERS2010模型计算 太阳光压模型 采用ECOM2模型 卫星位置/速度参数 采用随机游走模型估计；模型噪声为10−12 卫星太阳光压参数 采用随机游走模型估计；模型噪声为10−12 测站坐标参数 采用IGS的SNX文件进行强约束 测站钟差/卫星钟差参数 采用白噪声模型估计 模糊度固定参数 采用常数模型进行估计 对流层参数 采用分段常数模型估计；时段间隔为2 h 系统间偏差参数 采用分段常数模型估计；时间间隔为24 h 地球自转定向参数 采用分段线性模型估计；时段间隔为24 h 模糊度固定策略 松约束固定所有双差基线 最短卫星共视时长 900 s 宽巷模糊度平滑处理时间间隔 30 s 双差基线长度阈值 3500 km 
下载: 导出CSV
表 3 北斗超快速轨道产品以及实时滤波轨道浮点解和固定解与WUM产品轨道比较RMS均值统计结果
cm 卫星
类型Ultra-rapid SRIF-浮点解 SRIF-固定解 A C R 3D A C R 3D A C R 3D MEO 17.9 7.6 8.3 21.3 11.9 6.5 4.7 14.4 8.5 5.5 4.8 11.4 IGSO 49.0 18.7 31.0 61.4 16.8 13.6 9.3 23.5 12.4 10.5 9.9 19.2
下载: 导出CSV
表 4 北斗卫星轨道SLR检核统计结果
cm 卫星编号 轨道类型 实时滤波轨道 超快速轨道 WUM产品 均值 标准差 RMSE 均值 标准差 RMSE 均值 标准差 RMSE C08 IGSO 3.2 5.4 6.3 −21.5 20.6 29.8 6.5 5.0 8.2 C10 IGSO −1.7 11.5 11.6 −10.4 9.8 14.3 −1.0 5.2 5.3 C13 IGSO 6.0 10.1 11.7 −18.3 19.2 26.5 6.6 2.5 7.1 C11 MEO 0.1 4.0 4.0 3.8 3.5 5.2 1.1 2.5 2.7 C20 MEO 6.5 3.4 7.3 6.8 3.0 7.4 5.4 1.8 5.7 C21 MEO 5.4 5.5 7.7 5.6 3.6 6.6 4.8 2.1 5.3 C29 MEO −4.6 5.7 7.3 −5.9 5.3 8.0 2.7 3.0 4.0 C30 MEO −1.8 5.9 6.1 −3.5 6.0 7.0 3.6 3.8 5.2
下载: 导出CSV
-
[1] 陈俊勇. 全球导航卫星系统进展及其对导航定位的改善[J]. 大地测量与地球动力学, 2009, 29(2): 1-3. DOI: 10.3969/j.issn.1671-5942.2009.02.001 [2] 杨元喜. 北斗卫星导航系统的进展, 贡献与挑战[J]. 测绘学报, 2010, 39(1): 1-6. [3] 李冬航. “北斗+”融合创新与“+北斗”时空应用[J]. 卫星应用, 2020(103): 32-36. [4] 徐黎, 袁运斌. 不同分析中心GNSS实时SSR产品研究与分析[J]. 导航定位学报, 2020(6): 71-81. DOI: 10.3969/j.issn.2095-4999.2020.06.011 [5] ZHANG Q, MOORE P, HANLEY J, et al. Auto-BAHN: Software for near real-time GPS orbit and clock computations[J]. Advances in space research, 2007, 39(10): 1531-1538. DOI: 10.1016/j.asr.2007.02.062
[6] LAURICHESSE D, CERRI L, BERTHIAS J-P, et al. Real time precise GPS constellation and clocks estimation by means of a Kalman filter[C]//Proceedings of the 26th International Technical Meeting of the Aatellite Division of the Institute of Navigation (ION GNSS+ 2013), 2013: 1155-1163.
[7] DAI X L, LOU Y D, DAI Z Q, et al. Real-time precise orbit determination for BDS satellites using the square root information filter[J]. GPS solutions, 2019, 23(2): 1-14. DOI: 10.1007/s10291-019-0827-1
[8] 匡开发. GNSS卫星实时精密定轨技术研究[D]. 武汉: 武汉大学, 2019. [9] MONTEBRUCK O, HUGENTOBLER U, DACH R, et al. Apparent clock variations of the Block IIF-1 (SVN62) GPS satellite[J]. GPS solutions, 2012, 16(3): 303-313. DOI: 10.1007/s10291-011-0232-x
[10] SAASTAMOINEN J. Contribution to the theory of atmospheric refraction[J]. Bull geodesique, 1972, 105(1): 279-298. DOI: 10.1007/BF02522083
[11] 赵齐乐. GPS导航星座及低轨卫星的精密定轨理论和软件研究[D]. 武汉: 武汉大学, 2004. [12] 刘林. 人造地球卫星轨道力学[M]. 北京: 高等教育出版社, 1992. [13] BIERMAN G J. The treatment of bias in the square-root information filter/smoother[J]. Journal of optimization theory and applications, 1975, 16(1): 165-178. DOI: 10.1007/BF00935630
[14] LI X X, HAN X J, LI X, et al. GREAT-UPD: An open-source software for uncalibrated phase delay estimation based on multi-GNSS and multi-frequency observations[J]. GPS solutions, 2021, 25(2): 1-9. DOI: 10.1007/s10291-020-01070-2
[15] BERTIGER W, DESAI S D, HAINES B, et al. Single receiver phase ambiguity resolution with GPS data[J]. Journal of geodesy, 2010, 84(5): 327-337. DOI: 10.1007/s00190-010-0371-9
[16] GE M, GENDT G, DICK G, et al. Improving carrier-phase ambiguity resolution in global GPS network solutions[J]. Journal of geodesy, 2005, 79(1): 103-110. DOI: 10.1007/s00190-005-0447-0
[17] 戴小蕾. 基于平方根信息滤波的GNSS导航卫星实时精密定轨理论与方法[D]. 武汉: 武汉大学, 2016. [18] ZHANG H F, LONG M L, YANG H F, et al. Overview of Satellite Laser Ranging for BeiDou Navigation Satellite System[J]. Aerospace china, 2020, 21(04): 31-41. DOI: 10.3969/j.issn.1671-0940.2020.04.004
-
期刊类型引用(2)
1. 张世龙，管育贤，龚晓颖. 不同太阳光照条件下BDS卫星轨道预报精度分析. 全球定位系统. 2025(01): 93-102 . 
本站查看
2. 杨乐涵，李岩林，魏娜. GNSS基准站调制季节性信号滤波估计方法. 全球定位系统. 2025(01): 78-85 . 本站查看
其他类型引用(0)
计量
- 文章访问数: 271
- HTML全文浏览量: 281
- PDF下载量: 57
- 被引次数: 2
