一種偽距單點定位的數學模型研究及程序實現

郭平，周適，段太生，李學仕，王靠省

(1.中鐵二局集團有限公司,四川 成都 610031,2.中鐵二院工程集團有限責任公司，四川 成都 610031)

0 引 言

單點定位計算測站坐標已是比較成熟的技術．很多商業軟件著重解決不同測站點的基線向量解算,采用偽距和載波相位觀測值參與基線解算,筆者查閱文獻資料,現存大量關于單點定位甚至精密單點定位(PPP)的論文,但其講述的數學模型不夠詳細,通過論文依然較難整理出能直接應用于程序設計的一整套數學公式．筆者對單點定位的數學模型進行研究,提出一種能正確計算單點定位的數學模型,利用偽距觀測值,通過C#編程實現計算測站單點定位坐標,并提出一些結論．

1 單點定位數學模型分析

單點定位可利用接收機實測獲得原始數據,并轉換為通用的RINEX格式．其中N文件為GPS衛星的廣播星歷文件,G文件為GLONASS衛星的廣播星歷文件,C文件為BDS衛星的廣播星歷文件,N、G、C文件也可整合形成為一個文件,即公共導航星歷P文件[1],記載了各衛星系統的星歷數據．通過星歷文件可以計算衛星某一時刻的空間直角坐標XYZ, 再通過O文件各歷元的偽距和相位觀測值可以計算測站點的坐標．本文主要研究利用O文件的偽距觀測值,在GPS衛星系統下,采用無電離層模型,計算測站單點定位的坐標．首先對利用一個歷元的觀測數據實現測站單點定位功能的數學模型進行詳細介紹,在此基礎上實現多歷元數據共同解算測站單點定位的坐標．

若在歷元時刻t,接收機接收到n顆GPS衛星,選取高度角最高的衛星作為參考衛星,記作衛星r,其余衛星記作衛星s．

對于參考衛星r,列出該衛星與地面測站點在L1頻率的偽距觀測方程[2],如式(1)所示．

(1)

(2)

式中：x0，y0，z0為測站初始坐標,初值可設為(0,0,0),或采用該歷元所有GPS衛星的重心坐標在地面的投影值,初值不同對最終測站點的坐標計算沒有任何影響,只是迭代次數和運算時間上的差異而已．xr，yr，zr為衛星r的空間直角坐標．

同理,可列出其余衛星s與地面測站點在L1頻率下的偽距觀測方程

(3)

需注意的是,1個歷元共接收n顆GPS衛星,除了參考衛星r外,其余衛星數量為n-1,式(3)可列出n-1個方程．

利用其余衛星和參考衛星進行星間作差,即用式(3)減去式(1),可消去式(1)和式(3)的公共項C·dtr,即接收機鐘差可消去,可得:

(4)

(5)

602)C·dTrs+(772-602)Trs+

(6)

將式(6)移項可得:

C·dTrs-(772-602)Trs

(7)

簡寫為

L=Bx+V.

(8)

對于1個歷元接收到n顆GPS衛星,可列出n-1個誤差方程,如式(8)所示．其中常數項矩陣L(n-1)×1,系數矩陣B(n-1)×3,未知數矩陣x3×1,誤差項矩陣V(n-1)×1．根據最小二乘原理,可計算未知數x,如式(9)所示[3]．

x=(BTPB)-1BTPL.

(9)

此時,測站坐標可按式(10)計算．

(10)

由于測站初始坐標為(0,0,0),因此一次計算的結果還不準確,需進行迭代計算．迭代結束計算的判別指標可令式(9)計算的改正數x<1 cm．即改正數小于1 cm時程序停止迭代循環計算,筆者采用編制的單點定位程序運行時發現,迭代計算一般在5次以內,即可得到改正數x<1 cm的結果,大部分歷元數據迭代次數在3～4次即可完成計算．

式(9)中的權陣P可通過方差陣D求逆運算得到．考慮到歷元的每顆衛星高度角不同,高度角越高的衛星其方差應越小[4],考慮到式(7)隨機誤差所帶的系數,將衛星不同頻率偽距觀測值的方差視作相同,采用式(11)確定其方差．

(11)

根據誤差傳播定律,不難推導方差陣D：

D=(774+604)

(12)

式中,r為高度角最高的參考衛星,其余n-1顆衛星的方差根據不同高度角求得．

通過以上分析可求得每個歷元計算的單點定位坐標及精度信息．綜合各歷元數據計算的測站單點定位的坐標可通過法方程疊加的方法實現．具體分析如下:首先計算得到的每個歷元的單點定位坐標,以前面兩個歷元的數據為例,未知數x可寫為

(13)

現需要用前面兩個歷元的數據共同計算坐標,其誤差方程式為

(14)

根據最小二乘原理,未知數矩陣x為

(15)

展開后可得式(16)

(16)

同理,若有k個歷元共同參與解算,將k-1個歷元計算的坐標作為初始值,通過第k-1個歷元計算得到的坐標和累加的(BTPB)-1和BTPL,計算以k個歷元數據共同計算的單點定位坐標:

(17)

采用k個歷元數據共同計算測站坐標的單位權中誤差σ0，

(18)

以k個歷元共同計算測站點坐標的點位精度以及表征衛星分布狀態的DOP值計算過程如下:首先計算空間直角坐標系下的協因素陣Qxx，可由式(19)[5]得出，

(19)

該協因素陣是在空間直角坐標系中給出的,為了便于估算測站的位置精度,常采用其在站心坐標系的表達形式．站心坐標系統中的協因素陣設為QB．

(20)

其中旋轉矩陣H為[6]

(21)

式中,B、L分別為測站點對應的經緯度．

測站N、E、U方向上的點位精度MN、ME、MU可由式(22)計算得出．

(22)

平面位置精度因子(HDOP)值、高程精度因子(VDOP)值、空間位置精度因子(PDOP)值可由式(23)計算得出．

(23)

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

(24)

RMSE的定義應是各歷元坐標分量減去坐標真值,而非減去平均值(真值可采用精密星歷計算的測站坐標,精度可達cm級,相對于廣播星歷m級的精度,可近似認為是真值)．若沒有真值的情況下,可簡化用平均值代替．

在計算坐標各分量RMSE值后,易于得到一種簡易判別粗差的數據處理方法．即通過計算得到的RMSE值,判別各歷元的坐標與平均值的吻合程度,若坐標分量較差超過2倍RMSE值,則不采用該歷元的數據進行坐標計算的結果．

2 計算過程中需要用到的坐標轉換模型

在進行單點定位計算時,需要用到的坐標轉換有四個坐標轉換模型．分別是:1)大地坐標BLH轉換為XYZ;2)空間直角坐標XYZ轉換為大地坐標BLH;3)空間直角坐標XYZ轉換為站心坐標NEU;4)站心坐標NEU轉換為空間直角坐標XYZ．

由BLH轉換為XYZ的數學公式為

(25)

由XYZ轉換成BLH的數學公式為[7]

(26)

式中,緯度B需要迭代計算[8]，可令初值

由某點的空間直角坐標XYZ和以測站點X0、Y0、Z0作為站心原點的站心坐標NEU的數學公式為

(27)

將式(27)方程兩邊乘以H-1,易求由站心坐標NEU轉換為空間直角坐標XYZ的公式：

(28)

3 根據廣播星歷計算GPS衛星坐標的程序設計

由GPS廣播星歷可計算衛星在WGS-84系統下的空間直角坐標,星歷文件每2小時發布一次,每顆衛星都有參考歷元瞬間的開普勒軌道6參數、反映攝動力影響的9參數以及參考時刻參數,共計16個星歷參數[9]．廣播星歷計算衛星的精度約為2 m,若需要衛星坐標達到厘米級精度,需下載精密星歷產品,本文采用GPS導航N文件(廣播星歷)計算衛星坐標．廣播星歷的數據格式定義及衛星坐標計算公式可參見相關文獻[2，6]．

在程序設計時,需建立各種類,采用面向對象的方法進行編程．首先應建立RINEX類用于讀取N文件和O文件,RINEX類中定義兩種方法,用于讀取N文件和O文件,并將讀取的數據進行存儲．為使用方便,還需要定義衛星類(Satellite),星歷類(Ephemeris),空間笛卡爾坐標類(Cartesian)．在星歷類(Ephemeris)中定義N文件中所有星歷參數的變量,在笛卡爾坐標類(Cartesian)中定義空間直角坐標XYZ的變量, 衛星類(Satellite)包含星歷、接收機鐘差改正等信息,在衛星類中編寫一個Cartesian類型的函數,通過給定的時間參數,計算得到Cartesian類型的衛星空間直角坐標XYZ．

需注意的是,一天24小時共發布12個星歷,可用一個變量表示星歷號,程序通過給定的時間判斷星歷號,若O文件的歷元時刻所對應的星歷號在N文件中沒有,則需要搜素與之相鄰的星歷,這樣才能獲取N文件中的星歷參數用于計算GPS衛星坐標．

4 根據O、N文件計算測站單點定位的程序設計

通過對N文件的讀取,獲取各GPS衛星的星歷參數,可計算任一時刻GPS衛星的空間直角坐標XYZ．通過對O文件的讀取,獲取各歷元在L1、L2頻率下的偽距觀測值P1、P2[10].

在進行程序設計時,應建立如下類:數據類(data)、 橢球類(Ellipsoid)、大地坐標類(Geodetic)、站心坐標類(Local-Coordinates)、坐標轉換類(Coordinate-Transform)、對流層改正類(Troposphere)、測站類(Station)．數據類(data)里定義如下變量:偽距觀測值、GPS衛星的PRN號(1～32),GPS衛星在每個歷元數據中的索引號．橢球類Ellipsoid定義橢球的參數變量,如橢球長半軸、第一和第二偏心率、扁率．Geodetic定義某點的大地坐標BLH．Local-Coordinates定義某點相對于測站點為原點的站心坐標NEU．坐標轉換類Coordinate-Transform中定義四種靜態類型的坐標轉換的方法,如前文所述．Troposphere定義對流層改正模型的方法,本文采用Saastamonien對流層改正模型進行對流層改正．測站類Station中實例化之前定義的變量類型,以及定義式(9)中的各矩陣數組,通過之前所述的單點定位數學模型,實現利用偽距觀測值進行單點定位計算,可計算測站接收機相位中心的三維坐標以及得到DOP值、點位精度等精度信息．最后根據接收機天線相位參考點ARP和天線垂高信息,可計算測站點的空間直角坐標X，Y，Z．程序按照如下流程圖進行設計,可實現單歷元單點定位計算, 在單歷元計算得到定位坐標后,根據前述數學模型,可計算多歷元數據共同作用下單點定位的三維坐標．

圖1 根據O、N文件計算測站單點定位程序流程圖

5 實測數據的坐標及精度分析

5.1 采用各歷元單獨計算的坐標和RMSE值分析

運用本文所述的數學模型編制單點定位程序,利用實測某靜態RINEX數據,采樣間隔10 s,共396個歷元,高度截止角為15°．程序自動計算采用從1～396個歷元每個歷元單獨計算的單點定位坐標．如圖2～4所示．

圖2 單獨每個歷元計算的X坐標

圖3 單獨每個歷元計算的Y坐標

圖4 單獨每個歷元計算的Z坐標

從圖1～3所示,單獨利用每個歷元的數據計算的單點定位坐標是離散不連續的,沒有規律可循,這是因為各歷元之間的數據本身是獨立的,各歷元的坐標分量變化幅度在16 m以內,如表1所示．

表1 單獨采用每個歷元數據計算的坐標統計表

由表1可見,X坐標變化范圍為10.81 m,Y坐標變化范圍為15.56 m,Z坐標變化范圍為11.47 m．各歷元之間坐標分量變化范圍在16 m以內,震蕩稍大．這是因為僅采用了GPS偽距觀測值,未使用相位觀測值數據．若數學模型能加上相位觀測值,則各歷元單獨計算的坐標值震蕩會更小．

將各歷元單獨計算的坐標取平均值后與O文件中的近似坐標比較,X坐標較差為0.84 m,Y坐標較差為-1.78 m,Z坐標較差為0.07 m,各坐標分量較差均在2 m范圍內,一方面,O文件中的近似坐標也并非準確坐標,是接收機自帶程序計算的單點定位坐標,算法未知,可能使用了GPS和GLONASS的偽距觀測值或者相位觀測值．而本文僅使用了GPS衛星系統的偽距觀測值．

5.2 采用多歷元共同計算的坐標和精度分析

筆者根據前文所述的數學模型,采用法方程疊加的原理編制單點定位程序可計算多歷元數據共同解算的坐標．從第2個歷元開始,使用1～2個歷元的數據共同解算,第3個歷元則采用1～3個歷元共同解算,直到第396個歷元完成1～396個歷元的數據共同解算．圖5～7為多歷元數據共同解算的XYZ坐標變化圖．

圖5 多歷元共同解算的X坐標

圖6 多歷元共同解算的Y坐標

圖7 多歷元共同解算的Z坐標

由圖4～6可知,采用多歷元數據共同解算的單點定位坐標XYZ,其坐標值隨著參與解算的歷元數據增加,相鄰多歷元的坐標變化可認為是連續而非離散的,這是由式(17)中的法方程疊加決定的．相鄰的多歷元之間(如1～25和1～26)解算的坐標是很接近的,坐標變化不到0.5 m,但綜合不同歷元解算后的坐標,坐標震蕩變化比較顯著．這主要有以下三個原因:1)個別歷元解算的精度較低影響了整體解算的精度,可通過粗差判別予以剔除;2)數學模型和個別未知參數估計有待進一步地優化和研究；3)本次解算只采用了GPS單系統的偽距觀測值,未加入相位觀測值及其它衛星系統的觀測值,導致坐標變化范圍較大．若采用多衛星系統,加入各衛星系統的相位觀測值,采用多歷元數據計算的坐標值變化范圍將大大縮小．

在不考慮粗差影響的情況下,采用所有歷元數據(1～396)解算的坐標是最終結果．其值與O文件中的近似坐標較差比較如表2所示．

表2 所有歷元共同解算的坐標與O文件中近似坐標較差比較

由表2可知,所有歷元數據共同解算的坐標與O文件近似坐標各分量較差在3 m以內,一方面O文件的近似坐標也不是準確的數值,只能作為一個參考．如果用所有歷元共同解算的坐標與前文中單獨各歷元計算的坐標平均值比較,如表3所示,坐標各分量較差在3.5 m以內．

表3 共同解算與單獨各歷元解算的坐標平均值比較

在計算坐標的同時,點位精度Mx、My、Mz和DOP值也一并計算得到．如圖8～9所示．

圖8 多歷元共同解算的點位精度Mx，My，Mz

圖9 多歷元共同解算的HDOP，VDOP，PDOP

采用所有歷元數據(1～396)解算的點位精度和DOP值如表4所示．

表4 所有歷元(1～396)解算的坐標和點位精度、DOP值m

由圖8～9的點位精度和DOP值變化圖可以看出,隨著參與解算的歷元數據增加,點位精度Mx，My，Mz的數值在降低,DOP值也降低,這符合法方程疊加的特點,參與解算的歷元數據越多,精度越高．采用所有歷元解算的點位精度在0.5 m以內,DOP值小于1,說明計算的結果較可靠．

6 結束語

本文對偽距單點定位的數學模型進行分析和研究,用C#程序實現偽距單點定位的計算,并采用實測數據作為分析,測站點單點定位的坐標精度在米級,本文中的偽距單點定位模型計算結果與O文件的近似坐標較為接近．利用雙頻消除電離層模型并非完全消除電離層影響[11],還存在高階項的電離層誤差,另外對流層延遲、接收機鐘差、粗差探測等仍需要通過大量研究來進一步減小誤差．筆者編制的單點定位程序計算的各歷元單點定位的坐標,通過大量實測數據的演算,該數學模型是可行的．由于僅采用GPS單系統偽距觀測值進行單點定位計算,未加入相位觀測值及其它衛星系統,各歷元計算的坐標仍有10余米的變化,震蕩變化比較顯著,還需要對數學模型進一步研究和優化．采用法方程疊加的原理計算多歷元數據共同解算的單點定位坐標, 其計算結果與O文件的近似坐標較為接近,點位精度和DOP值可一并計算得到．若在計算多歷元數據時,通過粗差探測的方法,將個別誤差較大的歷元觀測值剔除,并加入相位觀測值,加入其它衛星系統的觀測數據,能進一步提高最終解算坐標的精度．