基于平方根信息濾波的北斗衛星實時定軌方法

姚夏，鄭鴻杰，袁勇強

(武漢大學測繪學院,武漢 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 對濾波軌道進行了外部檢核，結果也表明，實時濾波軌道的精度要普遍優于超快速產品.本文的研究對于提升北斗系統的實時高精度服務性能具有一定意義.