郭平，周適，段太生，李學(xué)仕，王靠省

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

0 引 言

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

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

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

若在歷元時刻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為測站初始坐標(biāo),初值可設(shè)為(0,0,0),或采用該歷元所有GPS衛(wèi)星的重心坐標(biāo)在地面的投影值,初值不同對最終測站點的坐標(biāo)計算沒有任何影響,只是迭代次數(shù)和運算時間上的差異而已．xr，yr，zr為衛(wèi)星r的空間直角坐標(biāo)．

同理,可列出其余衛(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)星進(jìn)行星間作差,即用式(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)

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

(10)

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

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

(11)

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

D=(774+604)

(12)

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

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

(13)

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

(14)

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

(15)

展開后可得式(16)

(16)

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

(17)

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

(18)

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

(19)

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

(20)

其中旋轉(zhuǎn)矩陣H為[6]

(21)

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

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

(22)

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

(23)

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

(24)

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

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

2 計算過程中需要用到的坐標(biāo)轉(zhuǎn)換模型

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

由BLH轉(zhuǎn)換為XYZ的數(shù)學(xué)公式為

(25)

由XYZ轉(zhuǎn)換成BLH的數(shù)學(xué)公式為[7]

(26)

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

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

(27)

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

(28)

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

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

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

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

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

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

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

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

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

5.1 采用各歷元單獨計算的坐標(biāo)和RMSE值分析

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

圖2 單獨每個歷元計算的X坐標(biāo)

圖3 單獨每個歷元計算的Y坐標(biāo)

圖4 單獨每個歷元計算的Z坐標(biāo)

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

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

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

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

5.2 采用多歷元共同計算的坐標(biāo)和精度分析

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

圖5 多歷元共同解算的X坐標(biāo)

圖6 多歷元共同解算的Y坐標(biāo)

圖7 多歷元共同解算的Z坐標(biāo)

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

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

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

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

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

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

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

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

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

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

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

6 結(jié)束語

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