紀(jì)志剛,余 銳,宣 偉,郭劍琴
(1.天津市市政工程設(shè)計(jì)研究院,天津300051 2.武漢大學(xué) 測(cè)繪學(xué)院,湖北 武漢430079)
在許多大型線狀工程中,地球曲率引起的長度變形對(duì)項(xiàng)目影響顯著,需要采用高斯投影的方式對(duì)數(shù)據(jù)進(jìn)行處理。現(xiàn)有程序多基于克拉索夫斯基橢球或1975國際橢球[1,2],不利于解算。筆者將CGCS2000橢球參數(shù)代入原始方程推出新公式,并將公式編制成程序,利用計(jì)算機(jī)進(jìn)行解算,滿足了需求。
2000國家大地坐標(biāo)系的原點(diǎn)為包括海洋和大氣的整個(gè)地球的質(zhì)量中心。坐標(biāo)系的Z軸由原點(diǎn)指向歷元2000.0的地球參考極方向,X軸指向格林尼治參考子午線與地球赤道面(歷元2000.0)的交點(diǎn),Y軸與Z軸、X軸構(gòu)成右手正交坐標(biāo)系[3]。采用廣義相對(duì)論意義下的尺度,2000國家大地坐標(biāo)系采用的地球橢球?yàn)镃GCS2000橢球[4],其基本參數(shù)如表1所示。
表1 CGCS2000橢球基本參數(shù)表[5]
該程序的實(shí)現(xiàn)過程如圖1所示。
圖1 程序運(yùn)算流程圖
該程序基于VS2008平臺(tái)用C#語言編寫,主要模塊為6個(gè)函數(shù),分別為solve_x0()、conversion_degree()、conversion_Text()、Gauss_ZhengSuan()、Gauss_FanSuan()和Gauss_HuanDai()。
1)Double solve_x0(doubleB),該函數(shù)的功能是利用給定的參數(shù)B,去求解赤道到緯度B所對(duì)應(yīng)平行圈的子午線弧長:
式中,M為子午圈曲率半徑。將M按級(jí)數(shù)展開,并進(jìn)行積分得弧長計(jì)算公式為:
將CGCS2000橢球參數(shù)代入得各系數(shù)數(shù)值,如表2所示。
表2 弧長計(jì)算公式系數(shù)表
將參數(shù)代入式(2)得弧長計(jì)算公式為:
2)conversion_degree(),該函數(shù)的功能是將rad轉(zhuǎn)換為(°)。由于坐標(biāo)反算的結(jié)果為大地坐標(biāo),坐標(biāo)換帶時(shí)也需要將投影坐標(biāo)首先反算得到大地坐標(biāo)。
3)conversion_Text(),該函數(shù)的功能是將大地坐標(biāo)(只考慮經(jīng)緯度,不考慮高程)轉(zhuǎn)換為rad。程序中大地坐標(biāo)的輸入形式用下例說明:如大地緯度為95°23′18.3″,輸入形式為95.231 83。
4)Gauss_ZhengSuan(),該函數(shù)的功能是高斯正算,即將已知大地坐標(biāo)轉(zhuǎn)換為目標(biāo)帶的高斯投影坐標(biāo)。Gauss_ZhengSuan()所解算出來的大地坐標(biāo),其y坐標(biāo)在原來的基礎(chǔ)上增加了500 km?;贑GCS2000橢球元素,適合電算的高斯正算公式為[6]:
式中,
將CGCS2000坐標(biāo)系參數(shù)代入式(7),然后匯編入Gauss_ZhengSuan()函數(shù)中。
5)Gauss_FanSuan(),該函數(shù)的功能是將高斯投影坐標(biāo)反算為大地坐標(biāo)。其所接收的大地坐標(biāo)的y坐標(biāo)也是在原來的基礎(chǔ)上增加了500 km。Gauss_FanSuan()函數(shù)的核心為高斯反算公式,基于CGCS2000橢球元素的高斯反算公式為:
接下來求解Bf。用迭代法,首先令
式中,X為已知投影點(diǎn)X坐標(biāo),6 367 449.145 71為式(5)的第1個(gè)系數(shù)。每次按式(11)迭代計(jì)算,直至
6)Gauss_HuanDai(),該函數(shù)的功能是高斯換帶,即將已知帶投影坐標(biāo)換算為目標(biāo)帶坐標(biāo)。換帶類型不同,所得結(jié)果會(huì)有所不同。常規(guī)的分帶有6°帶和3°帶2種,函數(shù)根據(jù)分帶的不同,又細(xì)分為4種類型:6°帶換算為 6°帶、6°帶換算為 3°帶、3°帶換算為 3°帶、3°帶換算為6°帶[7]。其中6°帶和3°帶之間的相互轉(zhuǎn)換具有相對(duì)復(fù)雜性,下面重點(diǎn)介紹二者的相互轉(zhuǎn)換過程。
① 6°帶換算為3°帶。以圖2為例,首先在6°帶內(nèi)利用已知點(diǎn)A的坐標(biāo),采用Gauss_FanSuan()函數(shù)計(jì)算出A點(diǎn)的緯度B以及經(jīng)度與中央子午線的經(jīng)差l。若l<1.5°,在3°帶內(nèi),A點(diǎn)位于以117°為中央子午線的投影帶內(nèi),此時(shí)l不變,再利用Gauss_ZhengSuan()函數(shù)計(jì)算出A點(diǎn)在該帶內(nèi)的投影坐標(biāo);若l>1.5°,在3°帶內(nèi),A點(diǎn)位于以120°為中央子午線的投影帶內(nèi),此時(shí)l=6°-l,再利用Gauss_ZhengSuan()函數(shù)計(jì)算出A點(diǎn)在該帶內(nèi)的投影坐標(biāo)。
② 3°帶換算為6°帶。首先利用Gauss_FanSuan函數(shù)計(jì)算出A點(diǎn)的緯度B和經(jīng)度與中央子午線的經(jīng)差l。判斷A點(diǎn)所在3°帶中央子午線是否滿足L0=6n-3。若滿足,則A點(diǎn)坐標(biāo)不變;若不滿足,l=1.5°+l,再利用Gauss_ZhengSuan函數(shù)計(jì)算出A點(diǎn)在該帶內(nèi)的投影坐標(biāo)。程序?qū)嵗龍D如圖3、圖4所示。
圖3 程序?qū)嵗龍D(1)
圖4 程序?qū)嵗龍D(2)
為檢驗(yàn)轉(zhuǎn)換結(jié)果的準(zhǔn)確性,下面均采用Powercoor軟件轉(zhuǎn)換結(jié)果作參考值,同類型的轉(zhuǎn)換結(jié)果進(jìn)行比較,以獲取轉(zhuǎn)換精度。為了同軟件中的顯示格式相同,各表中經(jīng)緯度均采用如下表示形式:30°30′15.2″表示為30.301 52。以下表中初始數(shù)據(jù)均為模擬數(shù)據(jù)。
表3 正算結(jié)果分析表/m
坐標(biāo)分量中誤差計(jì)算公式分別為:
點(diǎn)位中誤差為:
對(duì)表3各點(diǎn)分別作中誤差分析,得出最大點(diǎn)位中誤差為0.011 49 m。
表3和表4中,原始大地坐標(biāo)經(jīng)過正反算后又得到新的大地坐標(biāo),將二者進(jìn)行比較,結(jié)果如表5所示。
表4 反算結(jié)果分析表
表5 新舊大地坐標(biāo)對(duì)照表
表3~表5中數(shù)據(jù)顯示,在X方向和B方向轉(zhuǎn)換結(jié)果較參考值均偏大,Y方向和L方向轉(zhuǎn)換結(jié)果較參考值均偏小,存在一定的系統(tǒng)誤差。
將坐標(biāo)(5 635 146.875 0,629 151.817 9)分別以6°→ 6°,6°→ 3°,3°→ 6°,3°→ 3°4 種方式進(jìn)行換帶計(jì)算,并將結(jié)果與參考值進(jìn)行比較,見表6。對(duì)表6結(jié)果同樣利用式(11)進(jìn)行點(diǎn)位中誤差分析,經(jīng)計(jì)算得各點(diǎn)轉(zhuǎn)換后最大中誤差為1.047 cm。
表6 換帶計(jì)算表
綜上所述,坐標(biāo)正算結(jié)果中誤差最大差值可達(dá)到1.149 cm,一般情況其精度在mm級(jí);反算得到的大地坐標(biāo)精度可以達(dá)到10-4(″)量級(jí);坐標(biāo)換帶精度一般也可以達(dá)到mm級(jí),基本滿足了工程建設(shè)的需要。
[1]林曉靜,張小紅.ITRF2005與CGCS2000坐標(biāo)轉(zhuǎn)換方法與精度分析[J].大地測(cè)量與地球動(dòng)力學(xué),2010,30(2):117-119
[2]唐運(yùn)海,何誠.坐標(biāo)轉(zhuǎn)換及換帶計(jì)算的研究與實(shí)驗(yàn)分析[J].測(cè)繪與空間信息,2011,34(2):1-5
[3]李征航,黃勁松.GPS測(cè)量與數(shù)據(jù)處理[M].武漢:武漢大學(xué)出版社,2005
[4]程鵬飛,文漢江.2000國家大地坐標(biāo)系橢球參數(shù)與GRS80和WGS84的比較[J].測(cè)繪學(xué)報(bào),2009,38(3):89-94
[5]孔祥元,郭際明,劉宗泉.大地測(cè)量學(xué)基礎(chǔ)[M].武漢:武漢大學(xué)出版社,2007
[6]劉飛,周琳琳.GPS大地坐標(biāo)向地方坐標(biāo)轉(zhuǎn)換的實(shí)用方法研究[J].華東師范大學(xué)學(xué)報(bào),2005(1):73-77
[7]胡圣武.對(duì)高斯投影分帶的研究[J].地理空間信息, 2012,10(1):54-56