GPS偽距單點定位程序實現若干問題＊

李鶴峰，黨亞民，王世進，3，王霞迎

（1.中國測繪科學研究院大地測量與地球動力研究所，北京100830；2.山東科技大學測繪科學與工程學院，山東 青島266510；3.遼寧工程技術大學測繪與地理科學學院，遼寧 阜新123000）

0 引 言

在全球衛星導航系統（GNSS）的數據處理中，偽距單點定位是經常涉及的問題，并廣泛應用于實時位置導航和為RTK提供實時衛星位置和鐘差等信息［1－3]。作為GNSS最具代表性的全球定位系統（GPS），國內外學者有非常多的研究，偽距單點定位技術已經相當成熟。但作者在研究中發現，定位程序的編寫在相當多的文獻中僅一筆略過，對學習GPS編程中遇到問題的讀者造成很大的困擾，浪費大量研究的時間去摸索程序的編寫。選取GPS系統，基于Visual C＋＋平臺，編寫偽距單點定位程序，詳細研究偽距單點定位的原理、算法及程序實現過程，重點就程序編寫過程中的關鍵部分和易于出錯之處總結了詳細的解決思路，通過算例計算分析，驗證了問題解決的正確性和程序編寫的合理性。

1 偽距單點定位

1.1 定位原理

GPS偽距定位是通過空間距離的后方交會來實現，如圖1所示。理論上通過測算三顆衛星到用戶的距離，組成三元方程組，從而計算出用戶的三維坐標。實際應用中由于接收機鐘差難以預先準確確定，通常將其看作一個未知參數，與用戶三維坐標一并求解。這時需要第四顆衛星參與解算，組成四元方程組，進而精確計算出用戶的三維坐標［2］。

圖1 偽距定位原理

設t在時刻用戶測站點到S1，S2，S3，S4四顆衛星的偽距觀測值為ρj，j＝1，2，3，4，通過衛星導航電文解算，譯出該時刻四顆衛星的三維地心坐標（Xj，Yj，Zj），j＝1，2，3，4.用上述空間距離后方交會的定位原理，計算用戶位置的三維地心坐標（X，Y，Z），偽距觀測方程為

式中：假設偽距觀測量ρj已經通過星歷中的電離層、對流層改正和衛星鐘差改正；Rj為接收機天線相位中心到衛星天線相位中心的幾何距離；c為光速；δtk為接收機鐘差。

1.2 數學解算模型

式（1）為方便理解偽距定位原理，假設衛星鐘差、電離層和對流層延遲均已改正過。在定位解算中，考慮到這些誤差改正，偽距觀測方程式（1）改寫為

實際解算中，可以從觀測文件提取接收機的近似坐標 （X0，Y0，Z0），令 （δx，δy，δz）為接收機坐標的改正值，將式（2）在近似坐標 （X0，Y0，Z0）處用泰勒級數展開，取一次微小項，得線性化偽距觀測方程：

式中：

這些都可以根據已知值計算出結果，衛星鐘差δtj可以根據δtj＝af0＋af1（tc－toe）＋af2（tc－toe）2進行修正（af0，af1，af2，toe從衛星廣播星歷中讀取，tc為衛星發射信號時刻，通過計算求得），電離層延遲可以通過雙頻觀測消除，對流層延遲可以通過模型進行改正，所以，令常數項部分－δρjtrop＝Lj，四個未知數，需要四個觀測方程，將式（3）寫為矩陣形式：

式中：

根據最小二乘法求得改正數δX ＝ （ATA）－1ATL.求出δX＝［δxδyδzcδtk］T后，即可按

求出接收機位置坐標 （Xk，Yk，Zk）.

2 程序實現中問題及解決思路

根據上述GPS單點定位原理，基于Visual C＋＋平臺，編寫單點定位程序，程序設計流程如圖2所示。

由于在GPS單點定位程序編寫過程中，需要用到不少函數的調用及矩陣循環計算，數值分析，具體程序代碼較長且繁瑣，這里就程序實現中的關鍵點也是易于出錯點給予詳細的解決思路。

1）偽距電離層改正

偽距的電離層改正一般采用雙頻觀測消除。GPS信號中的C1，P1，P2測距碼調制在L1，L2載波上 （L1－C1，L1－P1，L2－P2），但是N 文件中用于電離層改正的偽距觀測數據類型（C1，P1，P2），并不是在每個歷元中出現的情況相同，為了達到最好的改正效果，就要有選擇的去利用，程序實現時要分情況考慮，優先選擇偽距觀測值P1，P2組合為

式中：P為雙頻改正后的偽距值；f1，f2為L1，L2載波頻率。沒有P1情況下使用C1，P2組合，如果沒有P2或同時沒有C1，P1的情況下，就不能使用雙頻觀測對電離層進行消除，需要另行考慮在單頻情況下的電離層模型改正。

圖2 GPS偽距單點定位程序設計流程圖

2）參考星歷選取

GPS星歷每2h播發一次，衛星發射信號時刻與星歷播發時刻相同的概率極小，所以計算衛星位置時，要選取最近的星歷作為衛星發射時刻的參考去計算衛星坐標。程序實現時，由于衛星信號傳播的時間非常短（0.07s左右）［2]，可以用觀測 O 文件中的時間t0與廣播星歷N文件中的時間tn求差值δt＝fabs（t0－tn）.當δt≤3 600s時，該時刻tn的廣播星歷即可作為參考用于計算觀測時刻tn對應的衛星坐標。

3）衛星發射信號時刻歸化

GPS時間系統采用原子時秒長作為時間基準，起算原點為1980年1月6日（星期日）0時0分0秒的協調世界時（UTC）.廣播星歷中播發的參考時刻toe是去除GPS整周數后不滿一周總秒長，GPST整周數的總秒長，而參與解算衛星位置的時間tk（歸化時間）是發射時刻tc相對于參考時刻toe的秒長，即tk－tc－toe.編寫程序時，要考慮到周的開始和結束，計算出tk后要對其進行如下判斷調整，如tk＞302400，那么tk＝tk－604800；如tk＜－302400，那么tk＝tk＋604800，302 400為半周的秒長，一周共604 800s.

4）偏近點角迭代計算

根據衛星軌道的偏近點角公式Ek＝Mk＋esinEk（平近點角Mk，衛星軌道偏心率e都已知），由于e非常小，sinEk＜1，故esinEk是一個微小值。可用程序通過迭代法實現Ek的快速求解，令Ek＝Mk連續回代求新的Ek，直到合適的精度停止迭代。截取C＋＋程序片段如下：

10－12為經驗值，取該值可以達到需要的收斂效果

Ek＝Ek1；／／迭代出偏近點真值

……

5）真近點角象限的判斷

真近點角計算式：

在程序的編寫過程中是最容易出錯的地方，其主要原因是用反正切函數沒有考慮fk所在的象限，若要通過程序正確計算fk，需要聯合考慮cos fk＝與值的符號，以確定fk所在的象限。截取C＋＋程序片段如下：

6）地球自轉改正

式中：ω為地球自轉角速度；t為信號傳播時間。

3 算例分析

本算例選用2012年1月5日河北某已知基準點的實驗數據，觀測采用雙頻GPS接收機，時間為7h41min 10s－8h51min 10s，共1h10min的GPS觀測數據，接收機采樣率為10s采集一次數據。根據編寫的C＋＋單點定位程序，計算接收機坐標，并將程序計算值與基準站已知坐標值對比求差，各分量的差值如圖3所示。

從圖3可以看出，通過1h10min的GPS連續觀測數據，用上述解決思路編寫的偽距單點定位程序解算基準點坐標，將解算值與已知基準點坐標值求差，結果X分量的差值在9m以內，Y分量的差值大都在15m以內，Z分量的差值在5m以內，表1示出了定位結果各分量殘差中誤差RMS統計分析表（坐標數據考慮到涉密問題僅取小數點前4位）。

表1 定位結果統計分析（單位：m）

由于主要著重解決程序編寫中易于出現的上述問題，所以程序在編寫過程中沒有考慮多路徑效應，接收機天線相位中心改正以及相對論效應等。但從定位結果圖和統計分析表可以看出，定位結果在5m以內，完全可以驗證GPS偽距單點定位解算結果的正確性，從而驗證了上述程序編寫中易于出錯點問題解決的合理性。

4 結 論

充分考慮各種誤差影響（電離層、對流層延遲，衛星、接收機鐘差，地球自轉改正，多路徑、接收機天線相位中心改正等），GPS單點定位精度一般在10m以內。雖然目前精密單點定位靜態已實現mm～cm級精度，動態精度也達到cm～dm級［7－8]，但是，由于GPS精密星歷在11天后才能得到，對實時導航定位沒有太大意義。利用實時GPS廣播星歷進行偽距單點定位，可以廣泛開展實時動態定位，對實現船只、飛機和車輛等各種運動目標的導航定位，監控管理具有重要意義。因此GPS偽距定位的研究不會過時，單點定位的原理和程序的編寫就尤為重要，更是學習和研究GNSS的基礎。

［1]徐紹銓，張華海，楊志強，等.GPS測量原理及應用［M].2版，武漢：武漢大學出版社，2003.

［2]黨亞民，秘金鐘，成英燕.全球導航衛星系統原理與應用［M].北京：測繪出版社，2007.

［3]廖 華.GPS偽距單點定位算法的綜合比較［J].測繪科學，2011，36（1）：20－21.

［4]朱 勇，方源敏，劉建忠.基于雙頻P碼的GPS偽距單點定位研究及精度分析［J].海洋測繪，2008，28（7）：15－18.

［5]何曉薇，牟其鋒，向淑蘭.GPS信號的電離層延遲誤差及改正［J].中國民航飛行學院學報，2008，19（1）：16－18.

［6]龔佑興.GPS單點定位研究［D].長沙：中南大學，2004.

［7]張金寶，王少閔.GPS精密單點定位程序設計與實現［J].城市勘察，2009（5）：60－62.

［8]CAI Changsheng，GAO Yang.Precise point positioning using combined GPS and GLONASS observations［J].Journal of Global Positioning Systems，2007，6（1）：13－22.