一種偽距單點定位的數(shù)學模型研究及程序?qū)崿F(xiàn)

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

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

0 引 言

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

1 單點定位數(shù)學模型分析

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

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

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

(1)

(2)

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

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

(3)

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

利用其余衛(wèi)星和參考衛(wèi)星進行星間作差,即用式(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衛(wèi)星,可列出n-1個誤差方程,如式(8)所示．其中常數(shù)項矩陣L(n-1)×1,系數(shù)矩陣B(n-1)×3,未知數(shù)矩陣x3×1,誤差項矩陣V(n-1)×1．根據(jù)最小二乘原理,可計算未知數(shù)x,如式(9)所示[3]．

x=(BTPB)-1BTPL.

(9)

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

(10)

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

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

(11)

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

D=(774+604)

(12)

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

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

(13)

現(xiàn)需要用前面兩個歷元的數(shù)據(jù)共同計算坐標,其誤差方程式為

(14)

根據(jù)最小二乘原理,未知數(shù)矩陣x為

(15)

展開后可得式(16)

(16)

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

(17)

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

(18)

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

(19)

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

(20)

其中旋轉矩陣H為[6]

(21)

式中,B、L分別為測站點對應的經(jīng)緯度．

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

(22)

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

(23)

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

(24)

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

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

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

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

由BLH轉換為XYZ的數(shù)學公式為

(25)

由XYZ轉換成BLH的數(shù)學公式為[7]

(26)

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

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

(27)

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

(28)

3 根據(jù)廣播星歷計算GPS衛(wèi)星坐標的程序設計

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

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

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

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

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

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

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

5 實測數(shù)據(jù)的坐標及精度分析

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

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

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

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

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

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

表1 單獨采用每個歷元數(shù)據(jù)計算的坐標統(tǒng)計表

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結束語

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