基于EGM2008重力場(chǎng)高程異常值模型在礦山測(cè)量中的應(yīng)用

王 兵

（貴州地礦測(cè)繪院有限公司，貴州 貴陽(yáng) 550018）

隨著全球衛(wèi)星定位系統(tǒng)GNSS（Global Navigation Satellite System，全球?qū)Ш叫l(wèi)星系統(tǒng)）不斷發(fā)展和民用精度的不斷提高，在工程測(cè)繪和礦山測(cè)繪等領(lǐng)域的重要性日益明顯。目前，基于GNSS的RTK（Real Time Kinematic，實(shí)時(shí)動(dòng)態(tài)）實(shí)時(shí)動(dòng)態(tài)差分測(cè)量不僅在傳統(tǒng)地面測(cè)繪還是在航空攝影測(cè)量中都發(fā)揮著重要作用。而RTK實(shí)時(shí)動(dòng)態(tài)差分測(cè)量所得到高程系統(tǒng)為以參考橢球面為基準(zhǔn)的大地高系統(tǒng)，而我國(guó)在工程領(lǐng)域或采礦等領(lǐng)域采用的是以我國(guó)1985國(guó)家高程為基準(zhǔn)的似大地水準(zhǔn)面正常高系統(tǒng)，故RTK實(shí)時(shí)動(dòng)態(tài)差分測(cè)量所得到大地高不能直接使用，必須通過(guò)高程異常改正為正常高后才能使用，故建立全省高程異常值模型是必要的。建立全省高程異常值模型既利于在水準(zhǔn)點(diǎn)稀少的區(qū)域獲得高精度的正常高，又利于在測(cè)量點(diǎn)大地高向正常高批量轉(zhuǎn)換電算化時(shí)有比較精確的模型數(shù)據(jù)基礎(chǔ)，提高電算化精度和計(jì)算機(jī)運(yùn)行效率，對(duì)貴州范圍內(nèi)高程異常值變化進(jìn)行應(yīng)用研究和掌握貴州范圍內(nèi)高程異常值變化趨勢(shì)和值域分布，可為工程測(cè)繪或礦山測(cè)量等領(lǐng)域的正常高成果精度分析提供數(shù)據(jù)基礎(chǔ)。

1 EGM2008全球重力場(chǎng)概述（Earth Gravitational Model EGM2008）

EGM2008重力場(chǎng)模型是較EGM84、EGM96之后由美國(guó)國(guó)家地理空間情報(bào)局地球重力場(chǎng)研發(fā)小組(US National Geospatial-Intelligence Agency)發(fā)布的較新地球超高階重力場(chǎng)模型，該模型的階數(shù)達(dá)到了2159次，并且加球諧系數(shù)的階擴(kuò)展到2190次，在精度和分辨率上都有巨大的提升。該機(jī)構(gòu)發(fā)布EGM2008重力場(chǎng)模型的目的在于解算更加精確的WGS84坐標(biāo)系統(tǒng)所對(duì)應(yīng)的似大地水準(zhǔn)面[1]，如下圖1所示。

圖1 EGM2008似大地水準(zhǔn)面圖

2 EGM2008高程異常計(jì)算

采用EGM2008重力場(chǎng)模型計(jì)算參考橢球面與EGM2008似大地水準(zhǔn)面高程異常值，計(jì)算公式為：

式中：φ為地心緯度；λ為地心經(jīng)度；ρ為地心向經(jīng)；GM為地心引力常數(shù)；α為參考橢球長(zhǎng)半軸；γP為正常重力值；Cmm和Smm分別為完全規(guī)格化位系數(shù)；為規(guī)格化Legndra函數(shù)；N為EGM2008重力場(chǎng)模型展開(kāi)的最高階數(shù)[2]。

本次采用AllTrans EGM2008 Calculator軟件進(jìn)行網(wǎng)格計(jì)算，經(jīng)查詢貴州省級(jí)行政區(qū)域在WGS84坐標(biāo)系下范圍為東經(jīng)103.60°～109.59°，北緯24.62°～29.23°，以該區(qū)域?yàn)榫W(wǎng)格區(qū)域進(jìn)行采樣計(jì)算（如下圖2所示），采樣間隔為1′×1′，計(jì)算結(jié)果選擇Bi-Quadratic值，EGM文件采用Und_min1x1_egm2008_isw=82_WGS84_TideFree_SE文件1′×1′網(wǎng)格[3]。

圖2 AllTrans EGM2008 Calculator軟件進(jìn)行網(wǎng)格計(jì)算

經(jīng)過(guò)該軟件計(jì)算得到1′×1′網(wǎng)格的EGM2008高程異常值，輸出的文件格式為文本文件數(shù)據(jù)交換格式，可以用文本文件或CSV格式EXCEL文件打開(kāi)查看和編輯。

3 起算面偏差值測(cè)定和偏差改正

由于EGM2008模型起算面為全球高程基準(zhǔn)所定義的起算面，與我國(guó)通過(guò)在青島驗(yàn)潮站測(cè)得的平均海水面所使用的起算面有一定的差值，即起算面偏差Δζ[4]。起算面偏差可以通過(guò)平均分布在全省的多個(gè)水準(zhǔn)點(diǎn)測(cè)得，采用在已知1985高程基準(zhǔn)高程值的水準(zhǔn)點(diǎn)上架設(shè)RTK設(shè)備測(cè)得該水準(zhǔn)點(diǎn)的大地高，設(shè)所測(cè)得的大地高為H大地高，水準(zhǔn)點(diǎn)處EGM2008高程異常值為ξEGM2008，已知水準(zhǔn)點(diǎn)1985高程基準(zhǔn)高程值為H85高，則有：

設(shè)AllTrans EGM2008 Calculator軟件進(jìn)行網(wǎng)格計(jì)算得到1′×1′網(wǎng)格的EGM2008高程異常值為m×n常數(shù)矩陣

采用全省多個(gè)水準(zhǔn)點(diǎn)聯(lián)測(cè)求得起算面偏差改正數(shù)m×n常數(shù)矩陣

則改正后大地高至1985國(guó)家高程基準(zhǔn)高程異常值m×n矩陣為：

4 高程異常模型建立

根據(jù)公式3雖然可得到可用于RTK實(shí)時(shí)動(dòng)態(tài)差分測(cè)量所得的大地高向1985國(guó)家高程基準(zhǔn)轉(zhuǎn)換的高程異常值，但是為m×n高程異常值矩陣，采樣間隔為1′×1′，矩陣在實(shí)際使用中并不利于計(jì)算機(jī)查詢，且在行列號(hào)相鄰處存在截?cái)嗾`差[5]。

為解決這一問(wèn)題，筆者采用不規(guī)則三角網(wǎng)的數(shù)字高程模型（簡(jiǎn)稱(chēng)TIN）根據(jù)矩陣中記錄的高程異常值集合將區(qū)域劃分為若干相連的三角面網(wǎng)絡(luò)，形成一系列互不交叉、互不重疊的連接在一起的三角形來(lái)表示高程異常表面，形成一個(gè)三維空間的高程異常線性模型（如下圖3所示）。用該TIN模型線性插值法生成加密的數(shù)字柵格模型，可有減少截?cái)嗾`差對(duì)查詢結(jié)果的影響，并且利于計(jì)算機(jī)查詢和二次開(kāi)發(fā)[6]。我國(guó)使用的測(cè)繪成果坐標(biāo)系統(tǒng)目前有大地坐標(biāo)系統(tǒng)，有高斯克呂格投影三度分帶和六度分帶，由于采用高斯克呂格3度分帶在貴州范圍存在跨帶，而六度分帶變形過(guò)大，故模型采用大地坐標(biāo)系統(tǒng)經(jīng)緯度投影，可實(shí)現(xiàn)在全省內(nèi)坐標(biāo)統(tǒng)一，使用者在查詢時(shí)須將平面坐標(biāo)轉(zhuǎn)換為大地坐標(biāo)后根據(jù)坐標(biāo)點(diǎn)位查詢相應(yīng)的高程異常值。若采用電算化方式編程計(jì)算時(shí)須將坐標(biāo)轉(zhuǎn)換為數(shù)字柵格模型所在行列號(hào)之后進(jìn)行查詢，計(jì)算機(jī)返回的行列號(hào)值即為查詢坐標(biāo)處的高程異常值[7]。

圖3 貴州范圍高程異常數(shù)字柵格模型

5 實(shí)例應(yīng)用

本文以我單位承擔(dān)的位于貴州省范圍內(nèi)一礦山傾斜攝影測(cè)量為例，受業(yè)主方委托，我單位采用傾斜攝影測(cè)量方式向業(yè)主提供5cm分辨率的礦山地表OSGB格式三維模型，要求坐標(biāo)系統(tǒng)為2000國(guó)家大地坐標(biāo)系，高程系統(tǒng)為1985國(guó)家高程基準(zhǔn)。礦山為已開(kāi)采的砂石礦山，占地面積約20hm2。

首先，我單位工作人員在礦山范圍內(nèi)均勻布設(shè)像控點(diǎn)5個(gè)，采用購(gòu)買(mǎi)千尋CORS位置服務(wù)賬號(hào)獲得2000國(guó)家大地坐標(biāo)系和大地高實(shí)時(shí)動(dòng)態(tài)差分信息[4]。在像控點(diǎn)上正確架設(shè)RTK接收機(jī)，在固定解狀態(tài)下平滑10次采集獲得像控點(diǎn)2000國(guó)家大地坐標(biāo)和大地高。采用計(jì)算機(jī)編程將平面坐標(biāo)轉(zhuǎn)換為大地坐標(biāo)，在貴州范圍高程異常數(shù)字柵格模型中采用計(jì)算機(jī)編程獲得每個(gè)像控點(diǎn)處高程異常值，并將RTK獲得的大地高改正為1985國(guó)家高程基準(zhǔn)正常高[8]，具體信息見(jiàn)表1所示。

表1 像控點(diǎn)坐標(biāo)及高程表 單位：m

為驗(yàn)證以上通過(guò)貴州范圍高程異常數(shù)字柵格模型轉(zhuǎn)換正常高的準(zhǔn)確性，我單位在省測(cè)繪資料檔案館收集到距測(cè)區(qū)約2km處的已知點(diǎn)YZ01，采用RTK在已知點(diǎn)上校正高程后實(shí)測(cè)得到5個(gè)像控點(diǎn)正常高，并對(duì)其進(jìn)行進(jìn)精度比較分析，結(jié)果如表2所示。

表2 高程精度對(duì)比表 單位：m

經(jīng)過(guò)實(shí)測(cè)與采用貴州范圍高程異常數(shù)字柵格模型轉(zhuǎn)換得到的高程對(duì)比發(fā)現(xiàn)，該測(cè)區(qū)高程最大誤差約為2cm，最小誤差約為0.5cm，滿足業(yè)主方傾斜攝影測(cè)量的精度要求。經(jīng)計(jì)算轉(zhuǎn)換高程中誤差為±1.2cm，完全滿足四等水準(zhǔn)高程精度要求，證明貴州范圍高程異常數(shù)字柵格模型在該區(qū)域內(nèi)準(zhǔn)確可靠。

6 結(jié)語(yǔ)

本文通過(guò)具體的流程詳述了應(yīng)用EGM2008重力場(chǎng)模型計(jì)算參考橢球面與EGM2008似大地水準(zhǔn)面的高程異常值，對(duì)該高程異常值進(jìn)行起算面偏差改正，建立貴州范圍高程異常數(shù)字柵格模型。對(duì)采用三角網(wǎng)插值法建立貴州范圍高程異常數(shù)字柵格模型的算法公式進(jìn)行圖解和推導(dǎo)，并根據(jù)我單位承擔(dān)的位于貴州省范圍內(nèi)一礦山傾斜攝影測(cè)量為實(shí)例驗(yàn)證貴州范圍高程異常數(shù)字柵格模型精度和可靠性，本方法適用于各類(lèi)采用1985國(guó)家高程基準(zhǔn)的工程測(cè)繪和礦山測(cè)量等圖根測(cè)量和碎部測(cè)量工作中，采用RTK實(shí)時(shí)動(dòng)態(tài)差分測(cè)量進(jìn)行高程異常改正精度可滿足于四等水準(zhǔn)測(cè)量精度要求。該模型可直接用于相關(guān)測(cè)繪地理信息系統(tǒng)軟件開(kāi)發(fā)，可用于監(jiān)測(cè)定位傳感器大地高轉(zhuǎn)正常高自動(dòng)化程序開(kāi)發(fā)等，對(duì)類(lèi)似測(cè)繪以及軟件開(kāi)發(fā)工作具有重要的參考價(jià)值。