姜韶,張坤先,匡志威
(1.四川省冶金地質(zhì)勘查局測(cè)繪工程大隊(duì),四川 成都 610212; 2.長(zhǎng)沙市規(guī)劃勘測(cè)設(shè)計(jì)研究院,湖南 長(zhǎng)沙 410007)
目前GNSS技術(shù)已經(jīng)基本取代了傳統(tǒng)大地測(cè)量技術(shù),在城市測(cè)量、工程測(cè)量中已經(jīng)成為主要的控制網(wǎng)測(cè)量技術(shù)。GNSS技術(shù)直接確定的是大地坐標(biāo)系坐標(biāo)(B、L、H),在實(shí)際應(yīng)用中需要通過(guò)高斯投影正算計(jì)算出平面直角坐標(biāo)系坐標(biāo)。根據(jù)高斯投影邊長(zhǎng)改化公式,離中央子午線的距離、高程均會(huì)引起長(zhǎng)度變形,因此在海拔較高的城市、地區(qū)以及線性工程的控制網(wǎng)測(cè)量中,均應(yīng)考慮選擇合適的中央子午線、相對(duì)參考橢球面來(lái)限制投影變形,精度要求應(yīng)滿足長(zhǎng)度變形值不應(yīng)大于 25 mm/km的要求。
隨著無(wú)人機(jī)航空攝影測(cè)量技術(shù)的廣泛應(yīng)用,利用無(wú)人機(jī)機(jī)載RTK實(shí)時(shí)獲取航攝相機(jī)曝光瞬間的高精度坐標(biāo)信息,加上無(wú)人機(jī)IMU慣性測(cè)量單元的高度集成化、高性能化、組合化發(fā)展,現(xiàn)階段無(wú)人機(jī)在小于1∶1 000比例尺的地形圖測(cè)量中可以完全實(shí)現(xiàn)免像控。但是POS數(shù)據(jù)中的坐標(biāo)一般都是地理坐標(biāo)和橢球高,在海拔較高城市、地區(qū)以及線性工程為了限制變形,仍需要在空三前將坐標(biāo)系統(tǒng)轉(zhuǎn)換到獨(dú)立坐標(biāo)系和1985國(guó)家高程基準(zhǔn)。少數(shù)的商業(yè)軟件雖然具備獨(dú)立坐標(biāo)系投影計(jì)算功能,但是對(duì)數(shù)據(jù)格式有固定要求,往往需要單獨(dú)將機(jī)載POS數(shù)據(jù)中的大地坐標(biāo)數(shù)據(jù)提取,轉(zhuǎn)換到獨(dú)立坐標(biāo)系后再替換到原POS數(shù)據(jù)中進(jìn)行空三處理,數(shù)據(jù)格式的往返和人工編輯替換比較費(fèi)事,不方便且容易出錯(cuò)。因此,經(jīng)過(guò)程序開(kāi)發(fā)實(shí)現(xiàn)高精度無(wú)人機(jī)機(jī)載POS數(shù)據(jù)的坐標(biāo)批量自動(dòng)轉(zhuǎn)換可有效減少工作環(huán)節(jié),有助于提高工作效率。
橢球膨脹法是保持參考橢球的定位、定向和扁率不變,對(duì)橢球進(jìn)行縮放,使得縮放后的參考橢球面與獨(dú)立坐標(biāo)系所選定的平面相切。因此,采用橢球膨脹法投影后得到的平面坐標(biāo)在數(shù)值上與國(guó)家參考橢球的橢球面上的平面坐標(biāo)接近,只是存在縮放關(guān)系,不會(huì)引起經(jīng)度(東坐標(biāo))的變化,只要選擇合適的投影面,就可以把高程所產(chǎn)生的變形降到最小,繼而忽略不計(jì)。因此在建立獨(dú)立坐標(biāo)系時(shí),通常采用橢球膨脹法。
橢球膨脹法關(guān)鍵是計(jì)算膨脹后的新橢球的長(zhǎng)半軸改正值dα,以及緯度改正值dB,具體的算法流程如圖1所示。
圖1 橢球膨脹法算法流程
1.2新橢球的長(zhǎng)半軸變化量計(jì)算公式
要滿足橢球膨脹法的縮放條件,要求項(xiàng)目區(qū)的平均曲率半徑R變化dR=H,其中H為區(qū)域的大地高。根據(jù)平均曲率半徑計(jì)算公式:
(1)
(2)
(3)
將式(2)、式(3)代入(1)求微分得到:
(4)
式中:M為子午圈曲率半徑,N為卯酉圈曲率半徑,B為大地緯度,e為橢球的第一偏心率,α為原橢球長(zhǎng)半軸,dα為縮放后橢球相對(duì)于原橢球的長(zhǎng)半軸改正數(shù),H為大地高。
1.3大地坐標(biāo)改正公式
從國(guó)家參考橢球到縮放后的參考橢球的坐標(biāo)由莫洛金斯基公式可得:
L=L0
(5)
(6)
上式中,L、B為縮放后的參考橢球?qū)?yīng)的大地坐標(biāo),L0、B0為原國(guó)家參考橢球?qū)?yīng)的大地坐標(biāo),H為大地高。
高斯投影正反算算法在很多文獻(xiàn)中均給出了詳細(xì)公式,在本文中將采用適用于不同橢球的高斯平面坐標(biāo)正算的實(shí)用公式[3],本算法可以在編程語(yǔ)言或者Excel表格中簡(jiǎn)單實(shí)現(xiàn)。
(7)
(8)
X=c0B-cosB(c1sinB+c2sin3B+c3sin5B)
(9)
式中,c0=Aa(1-e2)/ρ;c1=(C-B-D)a(1-e2)
(10)
基于本算法,在C#語(yǔ)言環(huán)境中開(kāi)發(fā)小程序“基于橢球膨脹法任意投影面高斯正反算工具”,小程序界面如圖2所示。選取四川西部某山區(qū)城市為例,投影面大地高為 2 600 m,投影中央子午線為102°,城市中心緯度為31°54′,現(xiàn)有GNSS控制點(diǎn)6個(gè)(CGCS2000坐標(biāo)系),利用該程序轉(zhuǎn)換到獨(dú)立坐標(biāo)系。
圖2 基于橢球膨脹法任意投影面高斯正反算軟件界面
為方便將無(wú)人機(jī)POS數(shù)據(jù)中的坐標(biāo)轉(zhuǎn)換到工程獨(dú)立坐標(biāo)系和1985國(guó)家高程基準(zhǔn),并按照規(guī)范要求進(jìn)行航片編號(hào)等,本文利用C#語(yǔ)言實(shí)現(xiàn)POS數(shù)據(jù)坐標(biāo)轉(zhuǎn)換和航攝文件整理功能,程序界面如圖3所示。其中平面投影轉(zhuǎn)換采用上文算法,由于無(wú)人機(jī)航飛區(qū)域一般跨度較小,一般僅僅幾平方公里,因此高程采用GNSS高程平面擬合即可。POS數(shù)據(jù)高程轉(zhuǎn)換目標(biāo)值(正常高)h由下式計(jì)算得到:
圖3 無(wú)人機(jī)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)點(diǎn)地理坐標(biāo);H為待轉(zhuǎn)POS坐標(biāo)點(diǎn)橢球高;B0、L0為測(cè)區(qū)中心地理坐標(biāo);a0、a1、a2為轉(zhuǎn)換系數(shù),由地面控制點(diǎn)計(jì)算得到。
(1)HGO、CosaGPS軟件投影計(jì)算
為驗(yàn)證本文實(shí)用算法投影后的長(zhǎng)度精度,采用相同的參數(shù)用CoordTool軟件(HGO軟件內(nèi)置系統(tǒng)工具)對(duì)以上六點(diǎn)進(jìn)行坐標(biāo)轉(zhuǎn)換。打開(kāi)CoordTool,“參數(shù)設(shè)置”界面上選擇“源橢球”“當(dāng)?shù)貦E球”均為“國(guó)家2000”;“投影”設(shè)置界面輸入中央子午線為“102:00:00.00000E”,“投影面高程”為“2600”,“平均緯度”為“031:54:00.00000N”,其他默認(rèn);“選項(xiàng)”界面“橢球變形方法”勾選“膨脹”;返回程序主界面,打開(kāi)源坐標(biāo)數(shù)據(jù)進(jìn)行高斯投影。
本文算法與HGO軟件計(jì)算結(jié)果如表1所示。
表1 本文算法與HGO軟件計(jì)算結(jié)果比較
采用相同的參數(shù)用在CosaGPS軟件里面設(shè)置好坐標(biāo)系統(tǒng)、中央子午線、測(cè)區(qū)平均緯度、投影面大地高等參數(shù)進(jìn)行坐標(biāo)轉(zhuǎn)換,具體設(shè)置如圖5所示。
圖5 CosaGPS軟件坐標(biāo)轉(zhuǎn)換設(shè)置界面
本文算法與CosaGPS軟件計(jì)算結(jié)果如表2所示。
表2 本文算法與CosaGPS軟件計(jì)算結(jié)果比較
由于與HGO、CosaGPS軟件在計(jì)算新橢球的長(zhǎng)半軸變化量及大地坐標(biāo)改正值采用的公式、參數(shù)精度存在差異,因而計(jì)算結(jié)果存在細(xì)微差別,這種差值在隨著投影面的高程增大而變大,投影點(diǎn)的緯度與中心緯度差值增大則變大。根據(jù)表1、表2可知,這種差值對(duì)邊長(zhǎng)影響較小,一般可以忽略不計(jì)。
(2)全站儀測(cè)距
為驗(yàn)證本算法的精度,采用全站儀對(duì)C1-C2、C3-C4、C5-C6邊長(zhǎng)進(jìn)行觀測(cè),經(jīng)過(guò)大氣改正的邊長(zhǎng)觀測(cè)值與本文計(jì)算結(jié)果反算的邊長(zhǎng)值對(duì)照如表3所示,可見(jiàn)均明顯優(yōu)于《工程測(cè)量標(biāo)準(zhǔn)》《城市測(cè)量規(guī)范》中對(duì)四等控制網(wǎng)最弱邊相對(duì)中誤差的要求。
表3 邊長(zhǎng)精度驗(yàn)算成果
(1)本文采用的獨(dú)立坐標(biāo)系高斯投影算法簡(jiǎn)單實(shí)用,適用于各種橢球,不僅可以通過(guò)計(jì)算機(jī)編程語(yǔ)言實(shí)現(xiàn),也可以利用Excel表格進(jìn)行投影計(jì)算,投影精度可以滿足一般城市地區(qū)、工程項(xiàng)目控制網(wǎng)測(cè)量要求。
(2)在進(jìn)行獨(dú)立坐標(biāo)系高斯投影時(shí),應(yīng)選擇合適的投影中央子午線、投影高程面、中心緯度。對(duì)于線性工程,比如高差較大的線路工程測(cè)量,應(yīng)進(jìn)行分帶分區(qū)投影,確保投影后的邊長(zhǎng)變形滿足工程精度要求。
(3)高精度無(wú)人機(jī)機(jī)載RTK技術(shù)在線路工程測(cè)量中已得到廣泛應(yīng)用,特別是帶狀地形圖測(cè)量基本由無(wú)人機(jī)航測(cè)技術(shù)替代,將POS數(shù)據(jù)中的高精度坐標(biāo)統(tǒng)一到工程獨(dú)立坐標(biāo)系中,有利于后期的空三加密處理工作,提高工作效率。
(4)本文算法有一定的局限性,單次投影的范圍不宜過(guò)大,特別是高差較大區(qū)域應(yīng)分區(qū)進(jìn)行投影計(jì)算,注意控制計(jì)算范圍。