基于平方根信息滤波的北斗卫星实时定轨方法

姚夏，郑鸿杰，袁勇强

(武汉大学测绘学院,武汉 430070)

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.7cm的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]：

BDS 卫星轨道状态参数由卫星位置、速度和太阳辐射压参数(SRP)构成.对于逐历元更新的北斗滤波定轨，需要实时构建下一历元的卫星轨道状态初值及前后历元间的状态转移矩阵.其求解可通过数值积分中的单步法和多步法进行：首先基于前一历元的初始轨道参数，利用龙格-库塔(Runge-Kutta)单步积分法获取足够的起算点；在此基础上，采用Adams 多步积分法完成后续弧段上积分点的计算[11-12].相应的轨道状态转移方程可表示为：

1.2 SRIF 方法

计算机进行数值计算过程中存在截断误差，导致传统卡尔曼滤波容易因数值计算误差而发散，因此扩展出了一些其他形式的平方根滤波算法.其中，SRIF因其具有数值稳定性和计算高效性，在导航及控制领域被广泛采用[13].相较于传统卡尔曼滤波需不断更新状态量及其协方差，SRIF 则是维护了基于平方根信息矩阵构成的信息方程.

SRIF 算法中包含时间更新和量测更新两个部分.假设ti时刻的先验信息方程和观测方程分别为和zi=Ai x+vi，式中：为先验的平方根信息矩阵；为方程噪声；zi为观测向量；vi为对应的观测噪声；Ai为系数矩阵.则相应量测更新过程可表示为：

式中：矩阵T为QR 分解得到的正交变换矩阵；构成了量测更新后的信息方程；ei为正交变换后的残余项.在北斗实时定轨中，根据历元间的时变特性，状态参数x可划分为确定性参数y和不确定性参数pi.则式(2)中的状态转移方程可进一步转换为

式中，w和D(w)为过程噪声及其协方差.结合量测更新后的信息方程和式(4)，可以得到时间更新的表达式为：

式中，参数pi可以直接消除以减少后续计算负担，而参数pi和y所对应的正交变换后的平方根信息矩阵则构成了时间更新后的信息方程.至此，SRIF 完成了当前历元的解算，通过循环上述过程不断更新相应的信息方程.

2 实验数据及处理策略

为充分测试基于SRIF 方法解算得到的北斗实时轨道精度，本文基于GREAT 软件[14]解算了2021 年年积日第106—113 天的地面测站的北斗观测数据.其中，选取了全球分布的125 个MGEX(Multi-GNSS Experiment)测站，具体分布如图1 所示.

图1 测站分布图

表1 给出了逐历元解算北斗滤波轨道时的相关模型及处理策略.由于处理卫星中同时包含BDS-2和BDS-3 卫星，相应的IF 组合观测值所采用的频率为B1I/B3I.对于北斗轨道状态参数，这里采用随机游走模型进行估计并设置相应的经验性模型噪声.而对于非保守力模型，在轨道状态参数估计时主要考虑了太阳光压模型，并采用了ECOM 五参数模型，其余非保守力模型则暂不考虑；对于其他保守力模型，则采用了常用的模型进行计算[15].另外，为最大程度降低实时模糊度错误固定对轨道精度的影响，采用了对所有模糊度双差基线进行松约束固定的策略[16]，并将固定约束仅作用在当前历元而非传递到后续历元.

表1 实时滤波轨道解算时的处理策略

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

图2 北斗卫星实时滤波轨道浮点解与WUM 产品轨道比较互差时序图

图3 北斗卫星实时滤波轨道固定解与WUM 产品轨道比较互差时序图

对于收敛后的北斗实时滤波轨道浮点解的时序结果，MEO 卫星在A、C 和R 上的变化幅度基本为±20 cm、±15 cm 和±10 cm，而IGSO 卫星在A、C 和R 上的变化幅度则基本为±25 cm.对于收敛后的固定解时序结果，其变化幅度略小于浮点解结果.同时由于前述的模糊度固定策略，固定解时序结果中包含了更多的锯齿状变化.而无论是浮点解和固定解，所有卫星在天边界处都包含了不同程度的跳变，其主要由参考轨道所有的天边界不连续性导致.进一步地，图4给出了超快轨道产品的预报弧段与WUM 产品互差的时序图.由于超快轨道分时段更新的特性，因此互差结果存在不同程度的周期性跳变，同时由于IGSO卫星轨道特性，其展现了比MEO 卫星更为显著的跳变现象.而在同时段内，逐历元解算得到的实时滤波轨道结果则为连续变化的，且不存在大幅度跳变.这也侧面验证了实时滤波轨道相较于超快速轨道具有更好的轨道连续性.

图4 北斗超快轨道产品与WUM 产品轨道比较互差时序图

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

图5 MEO 卫星超快速轨道产品以及实时滤波轨道浮点解和固定解与WUM 产品轨道比较RMS 统计图

图6 IGSO 卫星超快速轨道产品以及实时滤波轨道浮点解和固定解与WUM 产品轨道比较RMS 统计图

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

图7 北斗卫星轨道SLR 检核残差时序图

4 总结

本文采用SRIF 方法对26 颗北斗卫星进行了实时定轨解算，并对滤波轨道的时序特性与精度进行了全面评估.轨道时序特性分析的结果表明，BDS MEO 卫星能够在滤波开始20 h 内完成收敛，IGSO卫星也能够在24 h 左右达到收敛.相比于超快速轨道，实时滤波轨道的连续性和稳定性更好，有效避免了弧段边界的轨道跳变问题.进一步统计了滤波收敛后的轨道精度情况.结果表明，滤波定轨模式下，BDS MEO 和IGSO 卫星能够分别取得11.4 cm 和19.2 cm的三维轨道精度，相比于超快速产品，精度提升分别可达到46%和68%.在此基础上，通过SLR 对滤波轨道进行了外部检核，结果也表明，实时滤波轨道的精度要普遍优于超快速产品.本文的研究对于提升北斗系统的实时高精度服务性能具有一定意义.