姜韶,張坤先,匡志威
(1.四川省冶金地質(zhì)勘查局測繪工程大隊,四川 成都 610212; 2.長沙市規(guī)劃勘測設(shè)計研究院,湖南 長沙 410007)
目前GNSS技術(shù)已經(jīng)基本取代了傳統(tǒng)大地測量技術(shù),在城市測量、工程測量中已經(jīng)成為主要的控制網(wǎng)測量技術(shù)。GNSS技術(shù)直接確定的是大地坐標(biāo)系坐標(biāo)(B、L、H),在實際應(yīng)用中需要通過高斯投影正算計算出平面直角坐標(biāo)系坐標(biāo)。根據(jù)高斯投影邊長改化公式,離中央子午線的距離、高程均會引起長度變形,因此在海拔較高的城市、地區(qū)以及線性工程的控制網(wǎng)測量中,均應(yīng)考慮選擇合適的中央子午線、相對參考橢球面來限制投影變形,精度要求應(yīng)滿足長度變形值不應(yīng)大于 25 mm/km的要求。
隨著無人機航空攝影測量技術(shù)的廣泛應(yīng)用,利用無人機機載RTK實時獲取航攝相機曝光瞬間的高精度坐標(biāo)信息,加上無人機IMU慣性測量單元的高度集成化、高性能化、組合化發(fā)展,現(xiàn)階段無人機在小于1∶1 000比例尺的地形圖測量中可以完全實現(xiàn)免像控。但是POS數(shù)據(jù)中的坐標(biāo)一般都是地理坐標(biāo)和橢球高,在海拔較高城市、地區(qū)以及線性工程為了限制變形,仍需要在空三前將坐標(biāo)系統(tǒng)轉(zhuǎn)換到獨立坐標(biāo)系和1985國家高程基準(zhǔn)。少數(shù)的商業(yè)軟件雖然具備獨立坐標(biāo)系投影計算功能,但是對數(shù)據(jù)格式有固定要求,往往需要單獨將機載POS數(shù)據(jù)中的大地坐標(biāo)數(shù)據(jù)提取,轉(zhuǎn)換到獨立坐標(biāo)系后再替換到原POS數(shù)據(jù)中進行空三處理,數(shù)據(jù)格式的往返和人工編輯替換比較費事,不方便且容易出錯。因此,經(jīng)過程序開發(fā)實現(xiàn)高精度無人機機載POS數(shù)據(jù)的坐標(biāo)批量自動轉(zhuǎn)換可有效減少工作環(huán)節(jié),有助于提高工作效率。
橢球膨脹法是保持參考橢球的定位、定向和扁率不變,對橢球進行縮放,使得縮放后的參考橢球面與獨立坐標(biāo)系所選定的平面相切。因此,采用橢球膨脹法投影后得到的平面坐標(biāo)在數(shù)值上與國家參考橢球的橢球面上的平面坐標(biāo)接近,只是存在縮放關(guān)系,不會引起經(jīng)度(東坐標(biāo))的變化,只要選擇合適的投影面,就可以把高程所產(chǎn)生的變形降到最小,繼而忽略不計。因此在建立獨立坐標(biāo)系時,通常采用橢球膨脹法。
橢球膨脹法關(guān)鍵是計算膨脹后的新橢球的長半軸改正值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大地坐標(biāo)改正公式
從國家參考橢球到縮放后的參考橢球的坐標(biāo)由莫洛金斯基公式可得:
L=L0
(5)
(6)
上式中,L、B為縮放后的參考橢球?qū)?yīng)的大地坐標(biāo),L0、B0為原國家參考橢球?qū)?yīng)的大地坐標(biāo),H為大地高。
高斯投影正反算算法在很多文獻中均給出了詳細(xì)公式,在本文中將采用適用于不同橢球的高斯平面坐標(biāo)正算的實用公式[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)
基于本算法,在C#語言環(huán)境中開發(fā)小程序“基于橢球膨脹法任意投影面高斯正反算工具”,小程序界面如圖2所示。選取四川西部某山區(qū)城市為例,投影面大地高為 2 600 m,投影中央子午線為102°,城市中心緯度為31°54′,現(xiàn)有GNSS控制點6個(CGCS2000坐標(biāo)系),利用該程序轉(zhuǎn)換到獨立坐標(biāo)系。

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

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

圖4 HGO內(nèi)置CoordTool軟件界面
h=H-ζ
(11)
其中:
ζ=a0+a1×(B-B0)+a2×(L-L0)
(12)
式中:B、L為待轉(zhuǎn)POS坐標(biāo)點地理坐標(biāo);H為待轉(zhuǎn)POS坐標(biāo)點橢球高;B0、L0為測區(qū)中心地理坐標(biāo);a0、a1、a2為轉(zhuǎn)換系數(shù),由地面控制點計算得到。
(1)HGO、CosaGPS軟件投影計算
為驗證本文實用算法投影后的長度精度,采用相同的參數(shù)用CoordTool軟件(HGO軟件內(nèi)置系統(tǒng)工具)對以上六點進行坐標(biāo)轉(zhuǎn)換。打開CoordTool,“參數(shù)設(shè)置”界面上選擇“源橢球”“當(dāng)?shù)貦E球”均為“國家2000”;“投影”設(shè)置界面輸入中央子午線為“102:00:00.00000E”,“投影面高程”為“2600”,“平均緯度”為“031:54:00.00000N”,其他默認(rèn);“選項”界面“橢球變形方法”勾選“膨脹”;返回程序主界面,打開源坐標(biāo)數(shù)據(jù)進行高斯投影。
本文算法與HGO軟件計算結(jié)果如表1所示。

表1 本文算法與HGO軟件計算結(jié)果比較
采用相同的參數(shù)用在CosaGPS軟件里面設(shè)置好坐標(biāo)系統(tǒng)、中央子午線、測區(qū)平均緯度、投影面大地高等參數(shù)進行坐標(biāo)轉(zhuǎn)換,具體設(shè)置如圖5所示。

圖5 CosaGPS軟件坐標(biāo)轉(zhuǎn)換設(shè)置界面
本文算法與CosaGPS軟件計算結(jié)果如表2所示。

表2 本文算法與CosaGPS軟件計算結(jié)果比較
由于與HGO、CosaGPS軟件在計算新橢球的長半軸變化量及大地坐標(biāo)改正值采用的公式、參數(shù)精度存在差異,因而計算結(jié)果存在細(xì)微差別,這種差值在隨著投影面的高程增大而變大,投影點的緯度與中心緯度差值增大則變大。根據(jù)表1、表2可知,這種差值對邊長影響較小,一般可以忽略不計。
(2)全站儀測距
為驗證本算法的精度,采用全站儀對C1-C2、C3-C4、C5-C6邊長進行觀測,經(jīng)過大氣改正的邊長觀測值與本文計算結(jié)果反算的邊長值對照如表3所示,可見均明顯優(yōu)于《工程測量標(biāo)準(zhǔn)》《城市測量規(guī)范》中對四等控制網(wǎng)最弱邊相對中誤差的要求。

表3 邊長精度驗算成果
(1)本文采用的獨立坐標(biāo)系高斯投影算法簡單實用,適用于各種橢球,不僅可以通過計算機編程語言實現(xiàn),也可以利用Excel表格進行投影計算,投影精度可以滿足一般城市地區(qū)、工程項目控制網(wǎng)測量要求。
(2)在進行獨立坐標(biāo)系高斯投影時,應(yīng)選擇合適的投影中央子午線、投影高程面、中心緯度。對于線性工程,比如高差較大的線路工程測量,應(yīng)進行分帶分區(qū)投影,確保投影后的邊長變形滿足工程精度要求。
(3)高精度無人機機載RTK技術(shù)在線路工程測量中已得到廣泛應(yīng)用,特別是帶狀地形圖測量基本由無人機航測技術(shù)替代,將POS數(shù)據(jù)中的高精度坐標(biāo)統(tǒng)一到工程獨立坐標(biāo)系中,有利于后期的空三加密處理工作,提高工作效率。
(4)本文算法有一定的局限性,單次投影的范圍不宜過大,特別是高差較大區(qū)域應(yīng)分區(qū)進行投影計算,注意控制計算范圍。