獨立坐標系高斯投影實用算法程序?qū)崿F(xiàn)及應用

姜韶,張坤先,匡志威

(1.四川省冶金地質(zhì)勘查局測繪工程大隊,四川 成都 610212; 2.長沙市規(guī)劃勘測設計研究院,湖南 長沙 410007)

0 引 言

目前GNSS技術已經(jīng)基本取代了傳統(tǒng)大地測量技術,在城市測量、工程測量中已經(jīng)成為主要的控制網(wǎng)測量技術。GNSS技術直接確定的是大地坐標系坐標(B、L、H),在實際應用中需要通過高斯投影正算計算出平面直角坐標系坐標。根據(jù)高斯投影邊長改化公式,離中央子午線的距離、高程均會引起長度變形,因此在海拔較高的城市、地區(qū)以及線性工程的控制網(wǎng)測量中,均應考慮選擇合適的中央子午線、相對參考橢球面來限制投影變形,精度要求應滿足長度變形值不應大于 25 mm/km的要求。

隨著無人機航空攝影測量技術的廣泛應用,利用無人機機載RTK實時獲取航攝相機曝光瞬間的高精度坐標信息,加上無人機IMU慣性測量單元的高度集成化、高性能化、組合化發(fā)展,現(xiàn)階段無人機在小于1∶1 000比例尺的地形圖測量中可以完全實現(xiàn)免像控。但是POS數(shù)據(jù)中的坐標一般都是地理坐標和橢球高,在海拔較高城市、地區(qū)以及線性工程為了限制變形,仍需要在空三前將坐標系統(tǒng)轉(zhuǎn)換到獨立坐標系和1985國家高程基準。少數(shù)的商業(yè)軟件雖然具備獨立坐標系投影計算功能,但是對數(shù)據(jù)格式有固定要求,往往需要單獨將機載POS數(shù)據(jù)中的大地坐標數(shù)據(jù)提取,轉(zhuǎn)換到獨立坐標系后再替換到原POS數(shù)據(jù)中進行空三處理,數(shù)據(jù)格式的往返和人工編輯替換比較費事,不方便且容易出錯。因此,經(jīng)過程序開發(fā)實現(xiàn)高精度無人機機載POS數(shù)據(jù)的坐標批量自動轉(zhuǎn)換可有效減少工作環(huán)節(jié),有助于提高工作效率。

1 橢球膨脹法

橢球膨脹法是保持參考橢球的定位、定向和扁率不變,對橢球進行縮放,使得縮放后的參考橢球面與獨立坐標系所選定的平面相切。因此,采用橢球膨脹法投影后得到的平面坐標在數(shù)值上與國家參考橢球的橢球面上的平面坐標接近,只是存在縮放關系,不會引起經(jīng)度(東坐標)的變化,只要選擇合適的投影面,就可以把高程所產(chǎn)生的變形降到最小,繼而忽略不計。因此在建立獨立坐標系時,通常采用橢球膨脹法。

1.1 算法流程

橢球膨脹法關鍵是計算膨脹后的新橢球的長半軸改正值dα,以及緯度改正值dB,具體的算法流程如圖1所示。

圖1 橢球膨脹法算法流程

1.2新橢球的長半軸變化量計算公式

要滿足橢球膨脹法的縮放條件,要求項目區(qū)的平均曲率半徑R變化dR=H,其中H為區(qū)域的大地高。根據(jù)平均曲率半徑計算公式:

(1)

(2)

(3)

將式(2)、式(3)代入(1)求微分得到:

(4)

式中:M為子午圈曲率半徑,N為卯酉圈曲率半徑,B為大地緯度,e為橢球的第一偏心率,α為原橢球長半軸,dα為縮放后橢球相對于原橢球的長半軸改正數(shù),H為大地高。

1.3大地坐標改正公式

從國家參考橢球到縮放后的參考橢球的坐標由莫洛金斯基公式可得:

L=L0

(5)

(6)

上式中,L、B為縮放后的參考橢球?qū)拇蟮刈鴺?L0、B0為原國家參考橢球?qū)拇蟮刈鴺?H為大地高。

2 高斯投影通用算法

高斯投影正反算算法在很多文獻中均給出了詳細公式,在本文中將采用適用于不同橢球的高斯平面坐標正算的實用公式[3],本算法可以在編程語言或者Excel表格中簡單實現(xiàn)。

(7)

(8)

X=c0B-cosB(c1sinB+c2sin3B+c3sin5B)

(9)

式中,c0=Aa(1-e2)/ρ;c1=(C-B-D)a(1-e2)

(10)

3 算法編程實現(xiàn)及精度驗證

3.1 橢球膨脹法任意投影面高斯投影算法實現(xiàn)

基于本算法,在C#語言環(huán)境中開發(fā)小程序“基于橢球膨脹法任意投影面高斯正反算工具”,小程序界面如圖2所示。選取四川西部某山區(qū)城市為例,投影面大地高為 2 600 m,投影中央子午線為102°,城市中心緯度為31°54′,現(xiàn)有GNSS控制點6個(CGCS2000坐標系),利用該程序轉(zhuǎn)換到獨立坐標系。

圖2 基于橢球膨脹法任意投影面高斯正反算軟件界面

3.2 POS數(shù)據(jù)坐標自動轉(zhuǎn)換功能實現(xiàn)

為方便將無人機POS數(shù)據(jù)中的坐標轉(zhuǎn)換到工程獨立坐標系和1985國家高程基準,并按照規(guī)范要求進行航片編號等,本文利用C#語言實現(xiàn)POS數(shù)據(jù)坐標轉(zhuǎn)換和航攝文件整理功能,程序界面如圖3所示。其中平面投影轉(zhuǎn)換采用上文算法,由于無人機航飛區(qū)域一般跨度較小,一般僅僅幾平方公里,因此高程采用GNSS高程平面擬合即可。POS數(shù)據(jù)高程轉(zhuǎn)換目標值(正常高)h由下式計算得到:

圖3 無人機POS數(shù)據(jù)坐標轉(zhuǎn)換及航攝文件整理工具界面

圖4 HGO內(nèi)置CoordTool軟件界面

h=H-ζ

(11)

其中:

ζ=a0+a1×(B-B0)+a2×(L-L0)

(12)

式中:B、L為待轉(zhuǎn)POS坐標點地理坐標;H為待轉(zhuǎn)POS坐標點橢球高;B0、L0為測區(qū)中心地理坐標;a0、a1、a2為轉(zhuǎn)換系數(shù),由地面控制點計算得到。

3.3 橢球膨脹法任意投影面高斯投影算法精度驗證

(1)HGO、CosaGPS軟件投影計算

為驗證本文實用算法投影后的長度精度,采用相同的參數(shù)用CoordTool軟件(HGO軟件內(nèi)置系統(tǒng)工具)對以上六點進行坐標轉(zhuǎn)換。打開CoordTool,“參數(shù)設置”界面上選擇“源橢球”“當?shù)貦E球”均為“國家2000”;“投影”設置界面輸入中央子午線為“102:00:00.00000E”,“投影面高程”為“2600”,“平均緯度”為“031:54:00.00000N”,其他默認;“選項”界面“橢球變形方法”勾選“膨脹”;返回程序主界面,打開源坐標數(shù)據(jù)進行高斯投影。

本文算法與HGO軟件計算結(jié)果如表1所示。

表1 本文算法與HGO軟件計算結(jié)果比較

采用相同的參數(shù)用在CosaGPS軟件里面設置好坐標系統(tǒng)、中央子午線、測區(qū)平均緯度、投影面大地高等參數(shù)進行坐標轉(zhuǎn)換,具體設置如圖5所示。

圖5 CosaGPS軟件坐標轉(zhuǎn)換設置界面

本文算法與CosaGPS軟件計算結(jié)果如表2所示。

表2 本文算法與CosaGPS軟件計算結(jié)果比較

由于與HGO、CosaGPS軟件在計算新橢球的長半軸變化量及大地坐標改正值采用的公式、參數(shù)精度存在差異,因而計算結(jié)果存在細微差別,這種差值在隨著投影面的高程增大而變大,投影點的緯度與中心緯度差值增大則變大。根據(jù)表1、表2可知,這種差值對邊長影響較小,一般可以忽略不計。

(2)全站儀測距

為驗證本算法的精度,采用全站儀對C1-C2、C3-C4、C5-C6邊長進行觀測,經(jīng)過大氣改正的邊長觀測值與本文計算結(jié)果反算的邊長值對照如表3所示,可見均明顯優(yōu)于《工程測量標準》《城市測量規(guī)范》中對四等控制網(wǎng)最弱邊相對中誤差的要求。

表3 邊長精度驗算成果

4 結(jié) 論

(1)本文采用的獨立坐標系高斯投影算法簡單實用,適用于各種橢球,不僅可以通過計算機編程語言實現(xiàn),也可以利用Excel表格進行投影計算,投影精度可以滿足一般城市地區(qū)、工程項目控制網(wǎng)測量要求。

(2)在進行獨立坐標系高斯投影時,應選擇合適的投影中央子午線、投影高程面、中心緯度。對于線性工程,比如高差較大的線路工程測量,應進行分帶分區(qū)投影,確保投影后的邊長變形滿足工程精度要求。

(3)高精度無人機機載RTK技術在線路工程測量中已得到廣泛應用,特別是帶狀地形圖測量基本由無人機航測技術替代,將POS數(shù)據(jù)中的高精度坐標統(tǒng)一到工程獨立坐標系中,有利于后期的空三加密處理工作,提高工作效率。

(4)本文算法有一定的局限性,單次投影的范圍不宜過大,特別是高差較大區(qū)域應分區(qū)進行投影計算,注意控制計算范圍。