雙線程滑動窗口多系統GNSS超快精密定軌實現及評估

金煒桐,蔚保國,盛傳貞,張京奎,武子謙,陳永昌

(1. 中國電子科技集團公司第五十四研究所,石家莊 050081;2. 衛星導航系統與裝備技術國家重點實驗室,石家莊 050081)

0 引言

作為位置服務產業全面跨越升級的重要手段,全球導航衛星系統(global navigation satellite system, GNSS)高精度定位在交通、公安及林業等領域已形成規模化應用并全面進入大眾市場。隨著我國數字經濟的高速發展,位置服務產業將升級到以時空服務為核心的綜合定位、導航與授時(positioning, navigation and timing, PNT)時空服務產業[1],對GNSS定位的需求也不僅僅局限于高精度,而是邁向“又快(快速收斂)又精(高精度)”定位,因此GNSS提供的時空基準也需要同步響應快速、高精度的GNSS定位需求。

作為GNSS衛星空間基準的信息來源,精密衛星軌道產品是GNSS衛星為用戶提供連續精確的PNT服務的先決條件和必要前提[2]。從2000年起,國際GNSS服務組織(international GNSS service, IGS)為滿足實時用戶需求開始提供超快速軌道產品(IGS ultra-rapid, IGU),后續如武漢大學(Wuhan University, WHU)、歐洲定軌中心(Center for Orbit Determina-tion in Europe, CODE)及德國地學研究中心(Deuts-ches GeoForschungsZentrum, GFZ)等各大分析中心也逐步開始提供相應GNSS超快速軌道產品。由于GNSS衛星在近地空間自由受力(除軌道機動、地影及偏航姿態切換等非平穩運動外),其軌道較為平滑,在較短弧段內進行預報能夠維持初始歷元的軌道確定精度,因此超快速軌道確定主要基于GNSS基準站的事后觀測數據,生成的軌道產品中包含事后觀測弧段和預報弧段兩部分。相較于濾波法實時定軌[3],目前在同時需要保證精度和穩定性的實時應用中廣泛使用的(如實時鐘差估計、軌道鐘差改正實時數據流生成等)仍是超快精密定軌中的預報部分。IGS官方網站給出了最近60個GPS周期間8個分析中心的超快軌道產品與IGS超快軌道產品的互差量級[4],其中5個分析中心的超快軌道產品與IGS超快軌道產品的互差穩定在5 cm以下,3個分析中心的超快軌道產品與IGS超快軌道產品的互差穩定在5～10 cm。在此期間,各大分析中心提供的超快軌道產品較為穩定,除法國圖盧茲空間大地測量團隊(Grou-pe de Recherche en Géodésie Spatiale, GRGS)分析中心在GPS 2240～GPS 2242周與IGS超快軌道產品互差超過10 cm外,其他分析中心與IGS超快軌道產品的互差一直保持毫米級別的波動水平。

盡管各大分析中心的GNSS超快速軌道產品目前可在互聯網自由獲取,但受限于產品自主性需求等因素,仍有必要進行GNSS超快精密定軌的獨立實現。本文面向下一代GNSS系統建設,提出了天地協同PNT網絡的概念架構,并進行了基礎設施、服務中心和自主終端的研發和建設部署,依托天地協同PNT網絡服務中心云平臺對GNSS四系統衛星進行了超快速定軌,主要貢獻如下:介紹了天地協同PNT網絡的提出背景、概念內涵以及其服務中心部署的超快速精密定軌軟件,給出了服務中心運行的GNSS超快速精密定軌理論模型及算法,以及雙線程滑動窗口的精密定軌策略,最后使用重疊弧段比較、與外部分析中心軌道產品比較以及衛星激光測距(satellite laser ranging, SLR)檢核3種手段,全面驗證了30天的GNSS四系統超快速精密定軌精度。

1 天地協同PNT網絡服務中心及超快速精密定軌軟件

本文進行的GNSS超快速精密定軌業務部署在天地協同PNT網絡服務中心,其中天地協同PNT網絡的概念參見文獻[5]。服務中心生成的相應軌道產品通過數據接收模塊接收的“Multi-GNSS實驗”(The Multi-GNSS Experiment, MGEX)測站小時觀測數據和廣播星歷,利用部署在云平臺上的超快速精密定軌軟件解算得到,所采用的精密定軌軟件架構如圖1所示。

圖1 超快速GNSS精密定軌軟件架構

超快速GNSS精密定軌軟件由三部分核心模塊組成。數據接收下載模塊采用多線程方式對全球測站的小時觀測數據以及小時廣播星歷進行并行下載,并針對24 h弧段的觀測數據和廣播星歷進行合并存儲在相應路徑;GNSS數據處理模塊首先進行衛星廣播星歷動力學擬合,觀測數據鐘跳修復、粗差剔除和周跳探測,然后建立精密定軌法方程,結合最小二乘法對衛星軌道、光壓參數及衛星鐘差等多類型參數進行聯合估計,利用估計后的軌道參數和力模型參數作為初值進行GNSS衛星軌道預報,最后對參數估計中殘差較大的觀測值進行剔除,該過程迭代進行直至達到最大迭代次數或相應統計指標達到閾值范圍內停止迭代;GNSS產品生成模塊主要將軟件自定義的軌道文件格式轉換為標準sp3格式,并按照指定規則命名存儲在相應路徑。

2 GNSS衛星超快速精密定軌算法及解算策略

2.1 GNSS衛星超快速精密定軌算法及策略

用于GNSS衛星精密定軌的基本觀測方程可以表示為

(1)

通常采用無電離層組合消除一階電離層影響,線性化后的誤差方程如下

(2)

(3)

基于精密定軌基本觀測方程和誤差方程,本文采用的精密定軌策略如表1所示。

表1 GNSS衛星超快精密定軌模型配置策略

其中Galileo衛星各面板光學參數采用其官網發布的元數據(https://www.gsc-europa.eu/support-to-developers/galileo-satellite-metadata),BDS衛星各面板光學參數參見(http://www.beidou.gov.cn/yw/gfgg/201912/t20191209_19613.html)。鐘差基準選取主要有兩種方式:一是將衛星鐘差對齊到廣播星歷鐘差;另一種是選擇外接高精度原子鐘的地面測站鐘作為基準鐘。本文軟件采用第二種方式。為避免被選為基準鐘的測站沒有觀測數據的情況發生,在處理策略上選擇多個地面測站外接高精度氫原子鐘為基準鐘候選。確定基準鐘的原則為:可用觀測值比例(可用觀測值數量與觀測值總數的比值)大于85%,且可用觀測值數量最多。這里可用觀測值的具體含義是經過鐘跳修復、粗差剔除和周跳探測后實際用于參數估計過程的觀測值。

2.2 面向實時服務需求的解算策略

天地協同PNT網絡服務中心運行的GNSS超快速精密定軌最終需服務于有實時精密定位需求的用戶,因此針對后續的實時估鐘和軌道鐘差改正數播發等應用,需要部署GNSS超快速精密定軌的運行策略和后臺定時任務。目前武漢大學IGS數據中心、地殼動力學數據信息系統(the crustal dynamics data information system, CDDIS)等提供的小時觀測數據和各測站廣播星歷文件滯后約30～45 min,即每天第i(0≤i≤23) h 30～ 45 min期間,分析中心才會完全上傳完畢截止到第i(0≤i≤23) h之前的觀測數據。在接收數據完畢后,軟件運行一次完整定軌流程(指從對數據進行預處理至軌道產品生成完畢)所需時間約1.5 h,圖2給出了從2023年年積日151～157天100次完整定軌流程的運行時間統計。

圖2 2023年年積日151～157天100次完整定軌流程運行時間統計

從圖2中可以看出,一次完整定軌流程的運行時間在1.5 h上下波動(主要由于網絡波動以及服務中心其他運行程序的影響),大部分運行時長在1.5 h以內。因此在實際運行一次定軌過程結束后往往還未到第(i+2) h的45 min,仍需等待小時觀測數據和廣播星歷下載完畢才能進行下一次定軌流程,這樣一次定軌過程往往會滯后2.75 h。而對于后續實時應用,軌道產品的實際滯后時間更長,這是因為定軌處理流程較為耗時,在這1.5 h內并未有新的軌道生成,所使用的實際仍是上一次定軌得到的軌道產品。為降低后續實時應用所需軌道產品的滯后性,本文采用雙線程滑動窗口超快速精密定軌部署策略,具體如下:

1)每個小時的第30 min開啟數據下載線程,從武漢大學IGS數據中心和CDDIS下載小時觀測數據、各測站廣播星歷、差分碼偏差(differential code bias, DCB)文件,以及測站坐標周解文件和更新地球自轉參數(Earth rotation parameters, ERP)文件,具體下載流程如圖3所示。

圖3 數據下載進程具體算法流程

2)每個奇數小時和偶數小時的第30 min分別開啟單獨的定軌線程,每個定軌線程的具體流程如圖4所示。

圖4 定軌進程具體算法流程

對于下載線程,在下載過程結束后會在指定目錄下生成下載成功或失敗的標志文件;對于定軌線程,由于網速在實際運行時存在波動現象,該策略在線程的起始處每1 min檢測一次下載標志文件是否存在,若超過20 min還未檢測到下載成功標志則表示本次數據下載失敗,繼續進行后續的數據預處理流程。對于數據預處理,首先對各個小時觀測文件和各個測站廣播星歷文件分別進行融合,生成24 h 弧段長度的各測站觀測文件和全球廣播星歷文件;然后將全球廣播星歷視為觀測值對各GNSS衛星進行動力學擬合,剔除不健康衛星;接著對各測站觀測數據進行鐘跳修復、粗差剔除和周跳探測,根據處理結果剔除數據質量較差的測站,更新測站列表,至此定軌預處理過程完畢。

對于定軌核心進程,本文分GPS+GLONASS(GR)、GPS+Galileo(GE)和GPS+BDS(GC)3個線程并行執行定軌流程,既可以提升定軌整體運行效率,也能夠保證在其中一個線程失敗時不影響其他線程;最后將3個線程所確定的初始參數進行合并,執行軌道積分,生成精密軌道產品,拷貝到指定產品文件目錄。

此外,為了同時滿足定軌精度與處理速度要求,從現有的多系統測站中,挑選出了100個測站用于精密定軌。測站選取時,根據緯度將地球表面劃分成格網,盡可能保證每個格網點都有用于定軌的測站,并盡量均勻分布,其分布圖如圖5所示。

圖5 GNSS超快精密定軌測站分布圖

3 定軌精度驗證及分析

本章采用重疊弧段軌道比較、與外部參考軌道比較和SLR檢核3種方式對2023年年積日151～180天的定軌結果進行精度驗證。其中重疊弧段軌道比較是指將相鄰兩次定軌弧段重疊部分的軌道進行比較,是比較軌道內符合精度的常用方式,可驗證算法、模型和軟件實現的自洽性;與外部參考軌道比較和SLR檢核旨在比較軌道的外符合精度,前者是指將同時段定軌結果與分析中心的軌道產品進行比較,本文選擇武漢大學分析中心提供的單日最終事后精密軌道產品進行比較;后者則是固定解算的軌道,計算SLR觀測量的理論值并將其與SLR實際觀測值做差比較[15]。SLR檢核對軌道徑向偏差較為敏感,目前除GPS衛星外,BDS、GLONASS、Galileo等不同時期的衛星均搭載有SLR反射棱鏡,這也為相應的軌道產品SLR檢核提供了條件[2]。

進行軌道精度驗證的常用指標為所比較的2個軌道序列之差的均方根(root mean square, RMS)誤差,具體公式如下

(4)

其中,dRMS表示同時段內2個軌道的三維RMS誤差;dA、dC和dR分別表示同時段內2個軌道在沿跡、法向和徑向的RMS誤差。dA、dC和dR的計算公式具體如下

(5)

式中,n表示相同時段內參與軌道比較的歷元個數;ΔpiA、ΔpiC和ΔpiR分別表示相同時段內的第i個歷元,參與比較的2個軌道在沿跡、法向和徑向之差。此外,考慮到本文所用軌道框架與各個分析中心軌道的不一致性,軌道比較時采用了Helmert 7參數轉換,以便消除由于框架原點、尺度和旋轉所引起的系統性差異[16]。

3.1 重疊弧段軌道比較

本文采用相鄰兩次的定軌結果進行重疊弧段比較,分事后重疊弧段比較和預報重疊弧段比較兩種方式。對于奇數/偶數小時30 min定時運行的定軌線程而言,滑動窗口均為2 h,兩種比較方式的示意圖如圖6所示。

圖6 事后和預報重疊弧段示意圖

圖6中,T-24和T表示某一次生成的軌道產品的初始歷元和結束歷元。考慮到實時應用不會使用滯后時間過長的超快軌道產品,因此本文選擇的預報重疊弧段為6 h。

圖7和圖8分別顯示了2023年年積日151～180天各GNSS衛星事后和預報重疊弧段的平均軌道誤差。兩圖從上到下分別表示GPS衛星、GLONASS衛星、Galileo衛星、BDS2衛星和BDS3衛星,橫軸表示衛星PRN號,縱軸表示重疊弧段軌道誤差。Along、Cross、Radial和3D-direction分別表示沿跡方向、法向、徑向和三維誤差。表2給出了GNSS四系統衛星的平均重疊弧段軌道誤差統計情況,其中G、R、E以及C2I、C2M和C3M分別表示GPS、GLONASS、Galileo以及BDS2 IGSO、BDS2 MEO以及BDS3 MEO衛星。

表2 各GNSS衛星平均重疊弧段軌道誤差統計表

圖7 2023年年積日151～180天各GNSS衛星事后重疊弧段的平均軌道誤差

圖8 2023年年積日151～180天各GNSS衛星預報重疊弧段的平均軌道誤差

對于四系統MEO衛星,無論是事后還是預報重疊弧段,徑向軌道誤差均穩定在1～3 cm,但預報弧段在沿跡和法向兩個方向的軌道誤差發散較快,導致預報重疊弧段的三維軌道誤差大于事后重疊弧段。其中事后重疊弧段的四系統MEO衛星三維軌道誤差穩定在2～3 cm,預報重疊弧段的四系統MEO衛星三維軌道誤差則穩定在8～10 cm。BDS2 IGSO衛星觀測幾何不同于MEO衛星,其整體重疊弧段軌道誤差相比于MEO衛星較大:其事后和預報重疊弧段一維軌道誤差最小的方向分別約為3 cm和11 cm,三維軌道誤差分別約為6 cm和22 cm。

3.2 與外部參考軌道比較

本節選擇將定軌結果與武漢大學分析中心的最終事后精密軌道產品對比。圖9和圖10分別顯示了2023年年積日151～180天各GNSS衛星的事后和預報弧段與同時段分析中心軌道產品對比的平均軌道誤差。兩圖從上到下分別表示GPS衛星、GLONASS衛星、Galileo衛星、BDS2衛星和BDS3衛星,橫軸表示衛星PRN號,縱軸表示重疊弧段軌道誤差。Along、Cross、Radial和3D-direction分別表示沿跡方向、法向、徑向和三維誤差。表3給出了GNSS四系統衛星的平均重疊弧段軌道誤差統計情況,其中G、R、E以及C2I、C2M和C3M分別表示GPS、GLONASS、Galileo和BDS2 IGSO、BDS2 MEO以及BDS3 MEO衛星。

表3 各GNSS衛星與外部軌道產品對比軌道誤差統計表

圖9 2023年年積日151～180天各GNSS衛星與外部軌道產品比較的事后弧段平均軌道誤差

圖10 2023年年積日151～180天各GNSS衛星與外部軌道產品比較的預報弧段平均軌道誤差

可以看出,當與外部軌道產品進行比較時,各系統GNSS衛星的整體定軌精度較重疊弧段方式略低。對于事后弧段,GPS和Galileo衛星與外部軌道產品對比的徑向平均軌道誤差最小,穩定在1～2 cm之間;三維軌道誤差也整體最小,穩定在3～4 cm之間。特別地,G04、G11、G14、G18和G23這幾顆BLOCK IIIA衛星沿跡方向軌道誤差較其他衛星量級略大,這是因為兩者采用的光壓模型和相應軌道參數估計策略方式不同。武漢大學分析中心的軌道產品針對GPS BLOCK IIIA衛星采用7參數ECOM模型,在參數估計方面引入了沿跡方向的常數參數吸收誤差,而GPS III衛星的箱體是長方體,傳統的5參數ECOM光壓模型適用程度低。其他分析中心的軌道產品也采用了較新光壓模型和相應的經驗加速度參數估計策略,如CODE和GFZ均采用7參數模型,且在沿跡、法向和徑向3個方向引入了偽隨機參數并加入先驗約束,本軟件后續將針對光壓模型和經驗加速度估計方面進行升級改進。GLONASS衛星的平均徑向軌道誤差穩定在3～4 cm之間,但沿跡和法向軌道誤差均在6～8 cm之間,導致整體的平均軌道誤差達到10 cm左右量級;BDS2 MEO和BDS3 MEO衛星平均徑向軌道誤差均穩定在2～3 cm之間,三維軌道誤差整體穩定在5～7 cm水平。BDS2 IGSO衛星觀測構型不同于MEO衛星,其一維最小軌道誤差整體穩定在8～9 cm 之間,三維軌道誤差位于15～18 cm量級。

對于預報弧段,四系統MEO衛星的徑向平均軌道誤差相比于事后弧段增加了約1～2 cm量級。其中GPS和Galileo衛星與外部軌道產品對比的徑向平均軌道誤差穩定在2～3 cm之間,GLONASS衛星與外部軌道產品對比的徑向平均軌道誤差穩定在4～5 cm 之間,BDS2 MEO和BDS3 MEO衛星徑向平均軌道誤差均穩定在3～5 cm之間。但軌道預報在沿跡和法向發散較快,導致三維軌道誤差量級增大。其中GPS衛星三維軌道誤差最小,穩定在6～7 cm水平;Galileo衛星次之,其整體三維軌道誤差穩定在9 cm 水平;GLONASS、BDS2 MEO和BDS MEO衛星的整體平均三維軌道誤差則穩定在10～15 cm水平。BDS2 IGSO衛星觀測構型不同于MEO衛星,其預報弧段的一維最小軌道誤差整體穩定在15 cm水平,三維軌道誤差位于35 cm水平。

3.3 SLR檢核

本節對2023年年積日151～180天的超快定軌結果利用事后弧段和預報弧段兩種方式進行SLR殘差檢核,事后弧段和預報弧段的長度分別為24 h和6 h。在進行SLR檢核時,若殘差大于50 cm則被視為粗差加以剔除[17]。為保證預報弧段能夠覆蓋當天所有SLR觀測值,所參與計算的軌道在當天滑動窗口生成的所有軌道產品中進行了遍歷,直至能夠覆蓋當天所有SLR觀測值為止。在2023年年積日151～180天期間內,僅BDS和Galileo系統可以下載SLR觀測值,因此本節僅對這兩個系統的衛星進行SLR檢核。圖11和圖12分別顯示了事后弧段Galileo和BDS衛星的SLR殘差。其中為防止在一張圖中表示所有Galileo衛星造成混淆,分為了6個子圖進行展示。圖12中左側和右側分別表示BDS2和BDS3衛星的SLR檢核殘差。兩圖的橫軸表示時間,縱軸表示SLR殘差量級。圖13和圖14分別顯示了預報弧段Galileo和BDS衛星的SLR殘差,兩圖的具體排布與圖11和圖12相同,橫軸表示時間,縱軸表示SLR殘差量級。表4給出了事后弧段和預報弧段Galileo和BDS衛星的SLR殘差統計,表中均值和標準差均針對弧段內各自系統內所有衛星進行統計。

表4 BDS和Galileo衛星平均SLR殘差統計表

圖11 事后弧段2023年年積日151～180天Galileo衛星SLR檢核殘差

圖12 事后弧段2023年年積日151～180天BDS衛星SLR檢核殘差

圖13 預報弧段2023年年積日151～180天Galileo衛星SLR檢核殘差

圖14 預報弧段2023年年積日151～180天BDS衛星SLR檢核殘差

由于SLR檢核對軌道徑向誤差更敏感,鑒于軌道預報在徑向的發散程度不大,因此事后弧段和預報弧段的SLR殘差量級基本一致。其中,Galileo衛星SLR殘差均值穩定在1 cm水平,標準差穩定在2～3 cm水平;BDS衛星SLR殘差均值基本穩定在1～3 cm量級,標準差穩定在3～6 cm水平。

4 結論

本文基于天地協同PNT網絡服務中心云平臺實現了GNSS四系統超快速精密定軌,并對定軌結果進行了精度評價,介紹了天地協同PNT網絡概念內涵及服務中心超快速精密定軌軟件架構,提出了一種雙線程滑動窗口的超快速精密定軌策略,利用重疊弧段比較、與外部軌道產品比較以及SLR檢核3種方式對2023年年積日151～180天的超快精密定軌結果進行統計分析,得出如下結論:

1)對于軌道內符合精度,GPS、GLONASS、Galileo衛星以及BDS2和BDS3 MEO衛星事后和預報重疊弧段徑向軌道RMS誤差均穩定在1～3 cm,事后重疊弧段三維軌道RMS誤差穩定在2～3 cm,預報重疊弧段三維軌道RMS誤差穩定在8～10 cm;BDS2 IGSO衛星事后和預報重疊弧段最小一維軌道RMS誤差基本穩定在3～6 cm,三維軌道誤差分別為11 cm和22 cm水平。

2)對于軌道外符合精度,與武漢大學分析中心的最終事后精密軌道產品相比,GPS、GLONASS、Galileo衛星以及BDS2和BDS3 MEO衛星在事后弧段和預報弧段的徑向軌道RMS誤差整體在2～5 cm水平。GPS衛星三維軌道誤差最小,Galileo次之,然后依次是BDS3 MEO、BDS2 MEO和GLONASS衛星。BDS2 IGSO衛星觀測構型與MEO衛星不同,其軌道誤差大于MEO衛星;Galileo衛星SLR殘差均值穩定在1 cm水平,標準差穩定在2～3 cm水平;BDS衛星SLR殘差均值基本穩定在1～3 cm量級,標準差穩定在3～6 cm水平。

3)無論事后弧段還是預報弧段、進行軌道內符合還是外符合精度統計,MEO衛星沿跡和法向軌道誤差相對于徑向發散較快,其中徑向軌道誤差穩定在2～5 cm水平,能夠滿足后續厘米級實時定位應用對空間基準的精度需求。