基于雙頻非差無電離層延遲組合觀測量的導航衛星(GPS)精密星歷計算算法與實踐

蘇貝

(哈爾濱市勘察測繪研究院，黑龍江 哈爾濱 150010)

1 引 言

精密定軌理論的研究主要開始于20 世紀50年代末，經過半個多世紀的發展已日趨完善，定軌精度隨著衛星跟蹤技術的發展，測量精度的提高以及軌道動力學模型的演變和精細化得到了明顯的提升［1］。

在這之中，IGS 扮演了重要的角色，它匯集了世界知名的數據分析中心，包括JPL(Jet Propulsion Laboratory，美國噴氣推進實驗室)、CODE(Center for Orbit Determination in Europe，歐洲定軌中心)、GFZ(Geo Forschungs Zentrum，德國地學研究中心)、MIT(Massachusetts Institute of Technology，美國麻省理工學院)、SIO(Scripps Institution of Oceanography，美國斯克里普斯海洋學研究所)、USNO(U．S．Naval Observatory，美國海軍天文臺)、NOAA(National Oceanic and Atmospheric Administration，美國國家海洋和大氣管理局)等，計算并發布精密星歷，為地球科學研究、多學科應用、商業開發及教育提供支撐。在此基礎上，各家數據分析中心形成了業界知名的衛星導航定位數據分析軟件，如JPL 的GIPSY，CODE 的BERNESE，MIT 的GAMIT 等。

2 基本算法

從解微分方程的角度來看，精密定軌就是將一個常微分方程初值問題轉化為邊值問題，由邊值條件來確定初值［2］。這些微分方程通常也需要應用數值方法來解算。

2.1 觀測模型

導航衛星精密定軌的基本觀測量是L1 和L2 載波相位非差無電離層延遲線性組合觀測量電離層延遲可抵消掉。

基本觀測方程如下:

式中，λ1、λ2為L1/L2 載 波 波 長(m);f1、f2為L1/L2載波頻率(Hz);φL1、φL2為L1/L2 載波相位觀測值(周);ρ 為站星幾何距離(m);dt、dT為接收機和衛星的時鐘鐘差(s);Vtron為信號對流層路徑延遲(m);NLC為無電離層延遲線性組合的載波相位偏差(m);△pcvs、△pcvr為衛星和接收機的天線相位中心變化(m);△rel為相對論效應改正(m);△phw為相位纏繞效應改正(m);εLC為無電離層延遲線性組合測量噪聲(m)。

式(1)里的站星幾何距離ρ 通過解光時方程得到的，該方程考慮了衛星運動和接收機鐘差導致的誤差，為:

式中，rs(t)為衛星的空間位置(m);r'r為觀測站的空間位置(m);U(t)為慣性系到地固系的轉換矩陣;為衛星和接收機的天線位置改正(m);為測站位移(m)。

地固坐標通過轉換矩陣可以轉換成慣性系坐標，考慮歲差章動及地球自轉，轉換矩陣為:

式中，N(t)、P(t)為歲差章動矩陣;Rx、Ry、Rz分別為繞x、y、z 軸的旋轉矩陣;GAST 為格林尼治真恒星時;GMST 為格林尼治平恒星時;xp、yp為地極偏移;UT1－UTC 為地球自轉角偏;△ψ、ε 為章動經線分量和傾角;r 為世界時和恒星時之比。

衛星位置rs須使用統一的參考時間t 來表達，通常是采用GPS 時。因信號傳播時間與接收機鐘差之和△t 通常不超過0.1 s，時間同步方程可簡化為式(8)，該方程僅考慮衛星速度和地球點質量效應。同時，式(9)為φLC相對于時刻t 衛星位置rs的偏導數。

式中，GME為地球引力系數;vs(t)為t 時刻衛星s的速度(m/s);Φrr(t，t0)為t0時刻到t 時刻的衛星狀態轉移矩陣。

應用映射函數NMF，對流層延遲Vtron表示如下:

式中，ZTD 為天頂總對流層延遲(m);ZHD 為天頂對流層靜力學延遲(m);GE、GN為對流層大氣水平梯度參數;Mdry、Mwet為映射函數的干、濕分量;Az、El 為衛星的方位角、高度角;P0為平均海水面處的大氣壓(hPa);φ、H 為測站的緯度、大地高(m)。

衛星天線相位中線變化表示如下:

式中，Esat→eci為星固坐標系到慣性系的轉換;△satao為以星固坐標表示的天線偏移;θ 為相對于地面站的天底角。

接收機天線相位中心變化表示如下:

式中，Elocal→ecef為測站坐標系到地固坐標系的轉換;△ecc為測站標志到天線參考點的距離(m);△apcr1、△apcr2為天線參考點到L1、L2 載波相位中心的偏差(m);△pcvr1、△pcvr2為L1、L2 載波相位中心偏差(m)。由地球潮汐引起的測站位移表示如下:

式中，△solid為地球固體潮引起的位移(m);△ocean為海洋負荷引起的位移(m);△polar為極潮引起的位移(m)。

相對論效應和相位纏繞效應表示如下:

式中，D、D'為衛星、接收機的有效偶極向量;x，y，z為本地接收機單位向量(東、北、高度);x'，y'，z'為衛星體坐標單位向量;k 為衛星至接收機的單位向量。

表1歸納了觀測模型中精密測量改正的相關信息。

表1 精密測量改正項

2.2 衛星軌道模型

一般而言，衛星軌道的高度決定其所受擾動的程度。就精密定軌而言，任何大于10－10m/s2量級的加速度都需要考慮進來。因此，計算GPS 衛星軌道時，模型必須考慮到至少8 階重力位影響及相關潮汐效應，還有日月引力、太陽輻射壓和相對論效應。

衛星軌道方程如下式所示，為衛星位置和速度的常微分方程。其中，狀態轉移矩陣采用方程數值積分方式解算。表2列出了計算時所采用的一些精密衛星軌道模型。

式中，ageop為重力位加速度改正(m/s2);a3rdbody為三體重力加速度(m/s2);asrp為太陽輻射壓加速度(m/s2);arel為由相對論效應產生的加速度(m/s2)。

表2 精密衛星軌道模型

2.3 其他狀態轉移模型

其他的狀態轉移模型如下:衛星鐘差采用一階高斯馬爾科夫;接收機鐘為白噪聲;對流層參數為隨機游走;地球自轉參數為隨機游走。在一個弧段內，載波相位偏差認為是一個固定值。如果在該弧段探測到了周跳，則將其重新初始化。

2.4 參數估計

參數估計采用擴展卡爾曼濾波的方法。下列方程組分別表示觀測更新和時序更新。除了前向濾波外，計算時還可采用后向濾波及其更平滑的方法。

對于GPS 衛星定軌而言，卡爾曼濾波器所估計的狀態向量，包括衛星位置和速度、衛星時鐘、接收機時鐘、對流層參數。此外，如必要，衛星的太陽輻射壓和地球自轉參數也能納入進來。在估計中，衛星和接收機的時鐘當作相對于參考時刻而變化的時鐘來處理，并固定參考時刻為0。

輸入觀測向量zk后，tk時刻的狀態向量^xk通過卡爾曼濾波估計過程逐步得到。

3 算例

3.1 HRBCORS

哈爾濱雙星導航服務系統(GNSS)由哈爾濱市勘察測繪研究院于2007年建成，一期包含5 個連續運行參考站點，分別是阿城、哈西、對青、方臺、雙城。每站配備一臺雙頻72 通道接收機，型號為NetR5，可同時接收GPS 和GLONASS 衛星信號，是本次計算的數據來源，如圖1所示。

圖1 哈爾濱市連續運行參考站分布圖

3.2 軌道精度

在上面敘述的情況下，以開源的GAMIT 軟件代碼為基礎，添加了非差相位觀測方程解算模塊，修改了部分參數估計算法，得到了GPS 衛星的定軌結果，然后與IGS 綜合精密星歷比較，得出結論。后者由IGS 組織應用全球范圍內的GPS 站網，對所有分析中心給出的結果加權平均得到，其中每家分析中心應用不同的精密分析軟件包進行數據分析，由此得到的軌道精度一般能夠達到5 cm以內。

表3列出了衛星軌道的估計精度，以IGS 綜合精密星歷為標準，給出了所有衛星位置均方跟誤差的平均值。圖2比較了所有衛星的軌道估計誤差。

表3 衛星軌道精度

圖2 各衛星的軌道估計誤差

4 結 語

本文敘述了一種GPS 衛星精密定軌算法，結合實測數據計算了衛星星歷，與IGS 綜合精密星歷相比較，軌道外符合精度略大于5 cm，即衛星位置的三維均方差約為5.2 cm。作者下一步可繼續進行精密鐘差估計的工作，為區域非差分用戶提供精密定位和授時服務。

［1］劉林．航天器軌道理論［M］．北京:國防工業出版社，2000．

［2］施闖，趙齊樂，李敏等．北斗衛星導航系統的精密定軌與定位研究［J］．中國科學:地球科學，2012，42:845～861．

［3］鄭作亞，黃城，盧秀山．星載GPS 精密定軌進展及其數學模型［J］．大地測量與地球動力學，2007，27:112～118．

［4］周建華，楊龍，徐波等．一種導航衛星中長期軌道預報方法［J］．測繪學報，2011，40(S):39～45．

［5］李敏．多模GNSS 融合精密定軌理論及其應用研究［D］．武漢:武漢大學，2011．

［6］韓保民．基于星載GPS 的低軌衛星幾何法定軌理論研究［D］．武漢:中科院測量與地球物理研究所，2003．

［7］葉世榕．GPS 非差相位精密單點定位理論與實現［D］．武漢:武漢大學，2002．

［8］Oliver M．etc．．Satellite Orbits:Models，Methods and Applications［M］．Springer－Verlag，Heidelberg，2000．