星間頻間偽距定位模型的精度分析及程序實現*

龔櫟澎，楊久東，李 韌，梁 鵬

(華北理工大學 礦業工程學院，河北 唐山 063210)

偽距定位是通過利用GNSS接收機的偽距觀測數據計算測站坐標的定位方法，現在已經發展成為比較成熟的技術，主要用于導航定位和差分定位中初始大概位置的定位[1]。隨著精密單點定位(PPP)的逐漸深入研究，偽距單點定位方法及理論作為其基礎方法及定位理論，對偽距單點定位方法的深入研究是極其重要的[2-3]。偽距定位已經不僅僅局限于傳統的偽距定位，通過研究的不斷深入，已經衍生出了多種偽距定位解算模型[4-5]。研究發現，采用最小二乘迭代法對偽距進行解算,其迭代次數少但定位精度較低[6-7]，所以采用卡爾曼( Kalman) 濾波法解算可使得精度有所提高，但迭代次數過多[8]；很多學者研究使用最小二乘法和Kalman濾波法組合使用的方法進行偽距解算，減少迭代次數，提高定位精度[9]。偽距單點定位方法的精度一般為米級，組合偽距觀測值定位精度可高于單一偽距觀測值精度[10]，GPS靜態精密單點定位模型，點位精度可達毫米級[11]。BDS+GPS雙系統組合偽距定位模型的定位精度可高于單一系統偽距定位精度[12-13]。

偽距單點定位方法及理論的研究雖然很多[14]，但關于偽距定位的數學模型的表述多不夠詳細，筆者通過對偽距定位的進一步研究，整理出加權最小二乘偽距定位和星間頻間偽距定位模型，以Python3.8為工具進行編程實現，并利用實測數據驗證該數學模型的可行性與正確性以及程序的穩定性，同時將兩種偽距定位方法的定位結果進行比較，分析兩種偽距定位方法的定位效果及精度。

1 偽距定位模型

筆者通過對偽距單點定位模型的研究，利用接收機中的GPS衛星系統的觀測文件o文件中的偽距觀測值和星歷文件P文件，采用無電離層模型較短基線測量的實驗方法，分別計算兩種偽距單點定位模型中的測站坐標。

1.1 加權最小二乘偽距單點定位模型

若在歷元時刻t，接收機接收到n顆GPS衛星，選取高度角最高的衛星r為參考衛星，其余的衛星記作衛星s，對于參考衛星r列出該衛星與地面測站點在L1頻率的偽距觀測方程[15-16]，如式(1)所示：

(1)

(2)

(3)

可得式(4)：

(4)

對未知參數進行權陣為高度角定權模型的加權最小二乘求解，求得坐標結果如式(5)所示：

(5)

將計算出的偽距定位結果與(x0,y0,z0)進行比較，從而判斷各分量位置差值的大小，如果各分量差值大于0.01 m，則令x0=Xi，返回衛星位置繼續進行迭代計算，如果小于0.01 m，則退出迭代計算，輸出此歷元的定位結果，繼續下一歷元的坐標計算。

由于單獨每個歷元的數據都可解算一個測站點坐標，各歷元數據可視作是獨立觀測值，可通過均方根誤差(RMSE)指標判別各歷元解算坐標與平均值的變化幅度，可通過式(6)求得坐標各分量的RMSE值和坐標的RMSE值。

(6)

1.2 星間頻間偽距定位模型

文獻中根據加權最小二乘傳統偽距定位模型(星間單差偽距單點定位模型)提出了一種利用同一衛星系統中的衛星間及頻率間的雙差偽距定位模型[17]。由偽距定位原理可知，將式(2)～式(4)進行變形整理后可得式(7)：

(7)

根據參考衛星和其他衛星的觀測方程進行星間作差，可消去式(7a)和式(7b)的公共項cVtR,如式(8)所示：

(8)

GPS衛星系統信號播發頻率L1是調制基本頻率10.23 MHz×154獲得的頻率，L2是調制基本頻率10.23 MHz×120獲得的頻率，如式(9)所示：

(9)

(10)

依據公式進行最小二乘求解，因為星間頻間偽距定位方法是依據L1、L2頻率間作差進行計算定位，但是衛星的頻率發射是不穩定的，存在一定的頻率波動范圍，在后期解算中需要進行一定的系數調整后再計算出測站的偽距定位坐標。在進行最小二乘解算中的權陣采用同一歷元下衛星高度角定權模型，權陣為方差陣的逆，其方差公式如式(11)所示：

(11)

式中,σ為某顆衛星方差;θ為衛星高度角。根據誤差傳播定律可推導出式(12)：

D=(774+604)

(12)

式中,參考星為衛星高度角最大的衛星，其n-1余顆衛星的方差依據其高度角求得。再計算出權陣后通過式(4)～式(6)求出定位坐標結果以及均方根誤差。

2 偽距定位的程序設計

在利用Python語言進行兩種偽距定位模型的程序設計時，如圖1所示，建立了不同的類來實現偽距定位坐標的計算，主要包括文件讀取類、坐標轉換類、時間統一類、測站類、偽距定位模型計算類、精度評價類等六大類。

3 定位坐標解算及精度分析

運用本文所述的兩種偽距定位模型，利用2019年11月在河北省唐山市曹妃甸生態城控制網靜態觀測數據，選用C401、C410和GPS1共3個控制點，C401與C410間距5 493 m，C401與GPS1間距1 950 m，C410與GPS1間距4 348 m。同步觀測時間為100 min，靜態接收機型號為海星達IRTK2型雙頻接收機，且經過質檢鑒定，靜態數據文件為RINEX3.02文件，采樣間隔1 s，高度截止角為10°，因觀測時間較長，故選取觀測歷元中的300個歷元參與模型計算。程序自動計算采用從1～300個歷元計算兩個偽距定位模型中的定位坐標，同時使用中海達HGO基線解算軟件進行控制網的解算并進行嚴密二維約束平差后，最弱點中誤差的最大值為1.19 mm，符合精度要求，相對偽距定位模型的米級定位結果，可以認為HGO解算的二維約束平差結果為真值，并將兩種偽距定位模型的定位結果進行比較，如圖(2)～圖(10)所示。

圖(2)～圖(4)為測站C401的XYZ坐標值，圖(5)～圖(7)為測站C410的XYZ坐標值，圖(8)～圖(10)為測站GPS1的XYZ坐標值。從圖中分析可知，三個測站的加權最小二乘偽距定位模型和星間頻間偽距定位模型的每個歷元解算的定位坐標都出現了不同程度的波動，符合偽距定位的解算規律。在兩種定位方法模型解算的300個歷元坐標中，加權最小二乘方法的圖線波動較小，圖形表現較為平穩，最大波動較差不足5 m，但是均在精密坐標圖線下方，反觀星間頻間偽距定位模型解算的300歷元定位坐標，波動較大，大部分圖線波動在精密坐標圖線的上方與下方，其平均值與精密坐標圖線更為接近。

圖1 程序設計流程圖

圖2 測站C401的X坐標值

圖3 測站C401的Y坐標值

圖4 測站C401的Z坐標值

圖5 測站C410的X坐標值

圖6 測站C410的Y坐標值

圖7 測站C410的Z坐標值

圖8 測站GPS1的X坐標值

圖 9 測站GPS1的Y坐標值

圖10 測站GPS1的Z坐標值

為了更加直觀地分析兩種偽距定位方法在3個測站的定位精度，繪制表1，其中，ΔX、ΔY、ΔZ為300歷元的解算坐標平均值與精確坐標的差值，DX、DY、DZ為300歷元解算坐標的變化范圍。可以看出，3個測站在經過300個歷元解算后，對比兩種偽距定位方法的坐標平均值與精密平差的坐標的較差值發現，加權最小二乘偽距定位的X坐標精度要優于星間頻間偽距定位模型，但是Y、Z的平均值坐標解算精度均是星間頻間偽距定位的定位精度較高，加權最小二乘偽距定位的最小較差為0.863 7 m，最大較差為4.781 2 m，星間頻間偽距定位的最小較差為0.451 1 m，最大較差為2.526 4 m，但在解算中的坐標變化范圍波動較大，最大波動達到14.403 7 m，而加權最小二乘偽距定位最大波動僅為4.507 7 m。

表1 測站定位坐標誤差及坐標變化范圍/m

RMSE值反映觀測值偏離真實有效數據值的波動幅度情況，這一值越小誤差就越低，精度越高，偽距定位越準確。通過表2分析可以得出，用兩種偽距定位方法在3個測站進行300個歷元的解算，3個測站的RMSE值表現一致，均是加權最小二乘偽距定位的X坐標解算精度高于星間頻間偽距定位，而星間頻間偽距定位在3個測站的RMSEY、RMSEZ、RMSEXYZ的值均低于加權最小二乘偽距定位。綜合分析，星間頻間偽距定位方法的定位精度要高于加權最小二乘偽距定位，但是X坐標解算精度低于加權最小二乘偽距定位，主要原因可能在于選擇的解算模型不夠精確，在X定位計算中發生偏移，導致定位精度不高。

表2 測站RMSE值

4 結 論

本文對兩種偽距定位的數學模型進行分析和研究，用Python程序實現了加權最小二乘偽距定位模型、星間頻間偽距定位模型的計算，并利用靜態實測數據進行分析，測站點的偽距定位精度在米級，加權最小二乘偽距定位方法計算的XYZ坐標誤差在5 m內，星間頻間偽距定位方法計算的XYZ坐標誤差在3 m內，星間頻間偽距定位的精度高于加權最小二乘偽距定位。經過大量的實測數據的計算檢驗，筆者編制的偽距定位程序是可靠的，其數學模型是可行的。由于實驗數據區域較小，屬于短基線解算，對流層及電離層誤差影響較小，且僅采用了GPS單系統偽距觀測值進行偽距計算，未加入其它衛星系統進行聯合解算。若在后期加入其它衛星系統進行有效組合偽距定位，能進一步提高偽距定位精度及穩定性。