基于雙頻非差無電離層延遲組合觀測(cè)量的導(dǎo)航衛(wèi)星(GPS)精密星歷計(jì)算算法與實(shí)踐

蘇貝

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

1 引 言

精密定軌理論的研究主要開始于20 世紀(jì)50年代末，經(jīng)過半個(gè)多世紀(jì)的發(fā)展已日趨完善，定軌精度隨著衛(wèi)星跟蹤技術(shù)的發(fā)展，測(cè)量精度的提高以及軌道動(dòng)力學(xué)模型的演變和精細(xì)化得到了明顯的提升［1］。

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

2 基本算法

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

2.1 觀測(cè)模型

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

基本觀測(cè)方程如下:

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

式(1)里的站星幾何距離ρ 通過解光時(shí)方程得到的，該方程考慮了衛(wèi)星運(yùn)動(dòng)和接收機(jī)鐘差導(dǎo)致的誤差，為:

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

地固坐標(biāo)通過轉(zhuǎn)換矩陣可以轉(zhuǎn)換成慣性系坐標(biāo)，考慮歲差章動(dòng)及地球自轉(zhuǎn)，轉(zhuǎn)換矩陣為:

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

衛(wèi)星位置rs須使用統(tǒng)一的參考時(shí)間t 來表達(dá)，通常是采用GPS 時(shí)。因信號(hào)傳播時(shí)間與接收機(jī)鐘差之和△t 通常不超過0.1 s，時(shí)間同步方程可簡(jiǎn)化為式(8)，該方程僅考慮衛(wèi)星速度和地球點(diǎn)質(zhì)量效應(yīng)。同時(shí)，式(9)為φLC相對(duì)于時(shí)刻t 衛(wèi)星位置rs的偏導(dǎo)數(shù)。

式中，GME為地球引力系數(shù);vs(t)為t 時(shí)刻衛(wèi)星s的速度(m/s);Φrr(t，t0)為t0時(shí)刻到t 時(shí)刻的衛(wèi)星狀態(tài)轉(zhuǎn)移矩陣。

應(yīng)用映射函數(shù)NMF，對(duì)流層延遲Vtron表示如下:

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

衛(wèi)星天線相位中線變化表示如下:

式中，Esat→eci為星固坐標(biāo)系到慣性系的轉(zhuǎn)換;△satao為以星固坐標(biāo)表示的天線偏移;θ 為相對(duì)于地面站的天底角。

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

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

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

相對(duì)論效應(yīng)和相位纏繞效應(yīng)表示如下:

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

表1歸納了觀測(cè)模型中精密測(cè)量改正的相關(guān)信息。

表1 精密測(cè)量改正項(xiàng)

2.2 衛(wèi)星軌道模型

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

衛(wèi)星軌道方程如下式所示，為衛(wèi)星位置和速度的常微分方程。其中，狀態(tài)轉(zhuǎn)移矩陣采用方程數(shù)值積分方式解算。表2列出了計(jì)算時(shí)所采用的一些精密衛(wèi)星軌道模型。

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

表2 精密衛(wèi)星軌道模型

2.3 其他狀態(tài)轉(zhuǎn)移模型

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

2.4 參數(shù)估計(jì)

參數(shù)估計(jì)采用擴(kuò)展卡爾曼濾波的方法。下列方程組分別表示觀測(cè)更新和時(shí)序更新。除了前向?yàn)V波外，計(jì)算時(shí)還可采用后向?yàn)V波及其更平滑的方法。

對(duì)于GPS 衛(wèi)星定軌而言，卡爾曼濾波器所估計(jì)的狀態(tài)向量，包括衛(wèi)星位置和速度、衛(wèi)星時(shí)鐘、接收機(jī)時(shí)鐘、對(duì)流層參數(shù)。此外，如必要，衛(wèi)星的太陽輻射壓和地球自轉(zhuǎn)參數(shù)也能納入進(jìn)來。在估計(jì)中，衛(wèi)星和接收機(jī)的時(shí)鐘當(dāng)作相對(duì)于參考時(shí)刻而變化的時(shí)鐘來處理，并固定參考時(shí)刻為0。

輸入觀測(cè)向量zk后，tk時(shí)刻的狀態(tài)向量^xk通過卡爾曼濾波估計(jì)過程逐步得到。

3 算例

3.1 HRBCORS

哈爾濱雙星導(dǎo)航服務(wù)系統(tǒng)(GNSS)由哈爾濱市勘察測(cè)繪研究院于2007年建成，一期包含5 個(gè)連續(xù)運(yùn)行參考站點(diǎn)，分別是阿城、哈西、對(duì)青、方臺(tái)、雙城。每站配備一臺(tái)雙頻72 通道接收機(jī)，型號(hào)為NetR5，可同時(shí)接收GPS 和GLONASS 衛(wèi)星信號(hào)，是本次計(jì)算的數(shù)據(jù)來源，如圖1所示。

圖1 哈爾濱市連續(xù)運(yùn)行參考站分布圖

3.2 軌道精度

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

表3列出了衛(wèi)星軌道的估計(jì)精度，以IGS 綜合精密星歷為標(biāo)準(zhǔn)，給出了所有衛(wèi)星位置均方跟誤差的平均值。圖2比較了所有衛(wèi)星的軌道估計(jì)誤差。

表3 衛(wèi)星軌道精度

圖2 各衛(wèi)星的軌道估計(jì)誤差

4 結(jié) 語

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

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

［2］施闖，趙齊樂，李敏等．北斗衛(wèi)星導(dǎo)航系統(tǒng)的精密定軌與定位研究［J］．中國科學(xué):地球科學(xué)，2012，42:845～861．

［3］鄭作亞，黃城，盧秀山．星載GPS 精密定軌進(jìn)展及其數(shù)學(xué)模型［J］．大地測(cè)量與地球動(dòng)力學(xué)，2007，27:112～118．

［4］周建華，楊龍，徐波等．一種導(dǎo)航衛(wèi)星中長(zhǎng)期軌道預(yù)報(bào)方法［J］．測(cè)繪學(xué)報(bào)，2011，40(S):39～45．

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

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

［7］葉世榕．GPS 非差相位精密單點(diǎn)定位理論與實(shí)現(xiàn)［D］．武漢:武漢大學(xué)，2002．

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