夏 飛 羅文彬
(1湖北省交通規(guī)劃設計院, 湖北 武漢 430051;2中鐵第一勘察設計院集團有限公司, 陜西 西安 710043)
GAM IT/GLOBK軟件中的坐標參考系及其相互轉換
夏 飛1羅文彬2
(1湖北省交通規(guī)劃設計院, 湖北 武漢 430051;2中鐵第一勘察設計院集團有限公司, 陜西 西安 710043)
系統(tǒng)介紹了GAM IT/GLOBK軟件中使用的協(xié)議天球與協(xié)議地球參考系、ITRF參考框架、WGS84坐標系、IGS精密星歷參考框架、球面坐標系、站心地平坐標系、大地坐標系、空間直角坐標系以及NEU坐標系及其相互轉換關系。
GAM IT/GLOBK;坐標系統(tǒng)轉換;ITRF;參考框架
1.1 協(xié)議天球和地球參考系的定義
國際協(xié)議天球參考系(ICRS)往往用來描述天體和衛(wèi)星運動,天體和衛(wèi)星的星歷通常都在此系統(tǒng)中表示,技術人員可以利用GAMIT軟件來計算衛(wèi)星的運動方程,以及在實際運動過程中利用ARC的軌道坐標來進行計算,因此我們通常使用的 J2000協(xié)議來進行慣性坐標計算。這種計算方式的定義:在實際過程中可以利用地球質心來進行坐標原點的計算,以及選用2000年1月1日質心力學時(TDB)為標準歷元,利用瞬時來計算X軸和Y軸來計算,這樣X軸、Y軸Z軸就形成一個左右手結構體系,ICRS和瞬時天球坐標系之間可以相互轉換。
國際協(xié)議地球參考系(ITRS)的定義就是在進行地球力場的計算時,以及對地球表面點的定位計算的時候常常使用,這種計算方式的定義:通常需要以地球質心作為坐標原點,并控制z軸的方向為協(xié)議地極方向(CTP),控制好X軸指向格林尼治子午圈與協(xié)議地球赤道的交點,這樣就就可以構成Y軸X、Z軸構成左右手體系。ITRS與瞬時地球參考系可以通過極移改正進行相互轉換。
1.2 協(xié)議天球和地球參考系間的相互轉換
ICRS和ITRS之間的轉換需借助于瞬時真天球參考系與瞬時真地球參考系的坐標軸指向相同來實現(xiàn)。瞬時真天球參考系與協(xié)議天球參考系的差別在于地球自轉歲差和章動引起的坐標軸指向不同,轉換公式如下[4,5]:
式中,P為歲差旋轉矩陣,N為章動旋轉矩陣。AZ,Aθ和Aξ為歲差參數(shù),、Δ、ψΔ分別為黃赤交角、交角章動和黃經(jīng)章動。目前一般采用IAU1976年歲差常數(shù)和1980章動模型。
協(xié)議地球參考系與瞬時真地球參考系的差異是由于極移引起的,將瞬時真地球參考系轉換為協(xié)議地球參考系時,IERS的慣例是先繞 軸反時針旋轉 Py ,
再繞 軸旋轉 Px ,轉換公式為[4,5]:
式中,M 為極移旋轉矩陣。
因此,協(xié)議天球參考系與協(xié)議地球參考系的轉換關系為:
式中,E為地球自轉矩陣,GAST為真春分點的格林尼治時角。
2.1 ITRF框架相互變換
ITRF框架屬于一種動態(tài)性的地球參開框架,其定義就是利用框架的時間演變基準的明確定義、框架定向、尺度以及遠點等來實現(xiàn)的。在不同的時期,框架之間的四個基準分量定義也就存在著一定的差異,從而造成框架之間存在著較小的系統(tǒng)性差異,而這些差異通常都可以利用七個參數(shù)來進行表示,而且不同的框架之間能夠通過坐標系之間的相似變換來進行轉換,其轉換公式為:
式中 1T、 2T 、 3T 為平移量, 1R 、 2R 、 3R 為旋轉量,D為就是尺度改正因子,這就是ITRFxx框架到ITRFyy框架轉換的七個轉換參數(shù)。其中,任一參數(shù)P都會在指定時刻t的值都會與基準歷元的參數(shù) 0( )P t 相等加上基準歷元 0t 并轉換到歷元的變化量就會有:
P( t) = P( t0) + P˙( t - t0) (公式9)
這樣利用(8)式和(9)式就可以完成不同參考框架到指定歷元t的坐標轉換。
假設現(xiàn)在需要將 ITRFxx框架 t1歷元下的基準點坐標轉換到ITRFyy框架 t2 歷元下的,一般有兩種方法:(1)先轉框架,后轉歷元:首先將ITRFxx框架 t1歷元的基準點坐標轉換到ITRFyy框架 t1歷
在同一框架下不同歷元間轉換時,如果已知基準點的速度為 V I T RFxx,則基準點的坐標計算可按下式進行:
以上討論的是已知速度場的框架點間的轉換關系,對于不知道速度場的非框架點來說,可以采用目前國際上推薦使用的NNR-NUVEL-1A板塊運動模型進行近似計算。每個板塊的角速度分量可從地球物理模型計算得到[5],因此測站的速度可表示為:
式中 P i為測站所處的板塊,Ω為以年為單位的角速度矢量(單位為rad/a),
2.2 ITRF框架與IGS變換
GAMIT/GLOBK高精度數(shù)據(jù)處理一般采用國際GNSS服務中心提供的IGS精密星歷。IGS精密星歷選用的是ITRF參考框架,IGS采用的ITRF框架情況如表3-1所示。ITRF97框架之前IGS使用與ITRF相同的參考框架和參考歷元,ITRF97框架后,IGS開始使用自己的ITRF實現(xiàn),以保持一致性,IGS實現(xiàn)的框架與ITRF的差異在1cm精度范圍內。
表3 -1 IGS產(chǎn)品對應的ITRF框架[5]
高精度GNSS測量基準的統(tǒng)一,包括基準點坐標基準統(tǒng)一和基線解算時衛(wèi)星星歷基準統(tǒng)一兩方面的內容。假設某一GNSS測量中使用的IGS精密星歷所在的參考框架為ITRFxx,參考歷元為 t1;基準點所在的參考框架為ITRFyy,參考歷元為 t0 。采用“先轉框架,后轉歷元”的方法將基準點所在參考框架和參考歷元與 IGS精密星歷歸化到同一參考框架和參考歷元下,即此時再對歸化后兼容性較好的基準點坐標施加強約束,這樣建立的地心獨立坐標系屬于IGS精密星歷所采用的ITRF參考框架,實現(xiàn)了測量基準的統(tǒng)一。
2.3 ITRF框架與WGS84變換
WGS84坐標系屬于協(xié)議地球參考系,是美國GPS廣播星歷和美國國防制圖局NGS精密星歷的參考基準,最初是基于Transit衛(wèi)星多普勒數(shù)據(jù)建立的用于GPS廣播星歷的地球參考系,后來主要是基于GPS觀測數(shù)據(jù)實現(xiàn)。其定義是:原點為地球質心,Z軸指向BIH1984.0定義的協(xié)議地極CTP方向,X軸指向BIH1984.0零子午面與CTP對應的赤道的交點,Y軸與Z軸和X軸構成右手系。WGS84參考框架是由一組全球分布的GPS跟蹤站的坐標來具體實現(xiàn)的。WGS84系統(tǒng)經(jīng)過三次精化后,目前與ITRF框架的站坐標差異在1cm以內,在厘米級精度內可認為二者是同一參考框架[5]。一般來說,WGS84(1150GPS周)實用上被認為等同于ITRF2000;WGS84(873 GPS周)實用上被認為等同于ITRF94;WGS84(730 GPS周)實用上被認為等同于ITRF92。
3.1 球面坐標系與空間直角坐標系變換
GAMIT基線解算前需配置站坐標L文件,即輸入各測站點的先驗坐標,包括球面坐標(GAMIT的傳統(tǒng)格式)和指定歷元下的空間直角坐標和速度場(GLOBK的apr文件的格式)兩種形式。利用gapr_to_l程序可將GLOBK的apr文件轉換為指定歷元的站坐標L文件,ITRF框架坐標的apr文件可從MIT的ftp目錄獲得。ITRF參考框架下的球坐標與空間直角坐標的轉換公式如下:
式中, 為地心緯度, 為地心經(jīng)度,r為地心向徑。
使用球面坐標系能夠簡化球面坐標與空間的直角坐標之間的轉換,并且不需要大地坐標拔秧的迭代運算,降低了運算的難度和復雜程度,地心經(jīng)度與大地精度一樣,并且其經(jīng)度與緯度也比較接近,其地心向徑的變化幾乎被認識是大地高的變化。
3.2 站心地平坐標系與空間直角坐標系變換
GAMIT基線解算結果輸出文件為基線約束解q文件(詳細基線解)和o文件(簡略基線解)以及基線松弛解 h文件,主要包括基線解算過程參數(shù)和基線結果及其精度信息。在q文件和o文件中,基線解算結果各個分量及其方差協(xié)方差陣是以空間直角坐標系和站心地平坐標系兩種形式給出的。站心地平坐標系P-NEU定義為:以測站點P為原點,以P點的法線為U軸,指向天頂為正,以子午線方向為N軸,指北為正,E軸垂直于P點的大地子午面,向東為正,構成一個左手參考系。ITRF參考框架下的站心地平坐標系與空間直角坐標系的轉換公式如下:
式中,L為大地經(jīng)度,B為大地緯度。
需要說明的是,q文件和o文件中分別給出了XYZ和NEU形式基線分量的協(xié)方差陣,而NEU基線分量的中誤差是由XYZ
基線分量的轉換關系以及協(xié)方差陣并以誤差傳播定律為依據(jù)來進行計算的。通常都是將站心地平坐標中的基線NEU分量中的誤差當做基線高程以及水平方向的誤差。但這一觀念只能夠在基線相對較短時可這樣認為,若是基線較長就應對其基線NEU分量精度以及測站點中的精度之間的差別進行充分的考慮。
GLOBK的輸出文件一般為*.prt和*.org,在給出ITRF參考框架下的空間直角坐標的同時還給出了新的NEU坐標,它與站心地平坐標系定義的NEU有關系但又有不同,這種坐標系類似于平面坐標,屬于圓錐投影[6]。為了區(qū)分我們加用(G)表示GLOBK軟件中給出的NEU坐標系,記為 ( )N G 、 ( )E G 、 ( )U G 。NEU( G )坐標系統(tǒng)定義為:
N( G ):WGS84 橢球長半軸 a 與測站緯度之積,顯然它是一段弧長,北緯為正。
式中B是以弧度為單位的測站大地緯度。
E( G ):以二萬分之一弧度為最小度量單位,測站所在處最靠近的那條平行圈到起始子午線的平行圈弧長。
式中,L是以弧度為單位的測站大地經(jīng)度。r0為余緯為θ0時的緯圈半徑,其計算公式為:
式中,θ0定義為最接近二萬分之一弧度的余緯,其計算公式為:
式中, I nt(?)表示取整運算。
因此,如果知道了 GLOBK 軟件計算得到的 N EU( G )坐標, 我們尚不能嚴格反算出該測站的 XYZ 坐標,原因是無法得到精確的大地緯度B(因為有了取整運算)。另外從 N ( G )、 E ( G )、 U ( G )的計算過程可以看出,空間一組測站的 N ( G )、 E ( G )、 U ( G )并不構成一個統(tǒng)一的空間直角坐標系,因此也無法用兩個測站的 N ( G )、 E ( G )、U ( G )坐標差通過平方和開方的方法得到這兩個測站間的空間直線距離,這是它不同于一般真正的NEU坐標系的地方。N EU( G)基線分量的中誤差可由XYZ直角坐標系與大地坐標系的轉換關系以及大地坐標系與 N EU( G )坐標系的轉換關系根據(jù)誤差傳播定律計算得到。
希望本文內容對GAMIT/GLOBK學習者了解高精度GNSS數(shù)據(jù)處理的方法原理以及分析解算結果的精度有所幫助,從而更加合理地利用GNSS數(shù)據(jù)處理結果開展科研和生產(chǎn)工作。
[1]張捍衛(wèi)等. 天球參考系的基本理論和方法研究進展[J]. 測繪科學, 2005,30(2): 110-113
[2]孔祥元, 郭際明, 劉宗泉. 大地測量學基礎[M ]. 武漢:武漢大學出版社, 2005
U491.2
A
1007-6344(2017)09-0265-03