基于EGM2008模型處理DEM數(shù)據(jù)的應(yīng)用研究

宋學(xué)山 柴志勇 王冰 馮迎旭 王雪嬌

中水北方勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，中國·天津 300222

EGM2008；DEM；高程擬合

1 引言

隨著測(cè)繪技術(shù)的發(fā)展，尤其是合成孔徑雷達(dá)干涉測(cè)量技術(shù)（簡稱INSAR）、衛(wèi)星遙感技術(shù)、機(jī)載激光雷達(dá)（LiDAR）測(cè)量技術(shù)等成為新興的空間測(cè)量技術(shù)，由此獲得了數(shù)字高程模型（簡稱DEM）、數(shù)字表面模型（DSM）以及數(shù)字正攝影像（DOM）等大批量數(shù)據(jù)產(chǎn)品，具有精度高、信息豐富、直觀逼真、現(xiàn)實(shí)性強(qiáng)等優(yōu)點(diǎn)，可作為背景控制信息評(píng)價(jià)其他數(shù)據(jù)的精度、現(xiàn)實(shí)性和完整性；可從中提取自然信息和人文信息，并派生出新的信息和產(chǎn)品，為地形圖的修測(cè)和更新提供良好的數(shù)據(jù)和更新手段；同時(shí)克服了傳統(tǒng)測(cè)繪技術(shù)的人力投入大、通視等問題的影響，能夠快速、高效地獲得空間三維地理信息數(shù)據(jù)，引起了測(cè)繪等相關(guān)行業(yè)的濃厚興趣并得到廣泛應(yīng)用。

實(shí)際工程項(xiàng)目中，特別是在沒有勘測(cè)資料的工程區(qū)域，為了迅速而有效地實(shí)施規(guī)劃設(shè)計(jì)，大批量的DEM 數(shù)據(jù)應(yīng)運(yùn)而生。但DEM 數(shù)據(jù)產(chǎn)品的坐標(biāo)系統(tǒng)主要為WGS84 坐標(biāo)系，高程基準(zhǔn)為大地高，基于WGS84 橢球參考面，而在實(shí)際工程設(shè)計(jì)中采用的高程基準(zhǔn)為正常高，基于似大地水準(zhǔn)面。如何實(shí)現(xiàn)以WGS 橢球參考面為基準(zhǔn)的大地高到以似大地水準(zhǔn)面為基準(zhǔn)的正常高的轉(zhuǎn)換成為限制DEM 數(shù)據(jù)應(yīng)用的常見問題。

本文提出利用最新的EGM2008 地球重力場(chǎng)模型，對(duì)已有的DEM 數(shù)據(jù)進(jìn)行高程數(shù)據(jù)轉(zhuǎn)換，以滿足實(shí)際生產(chǎn)的需要，為以后在無人能到達(dá)的區(qū)域、面積較大區(qū)域、山區(qū)或密林區(qū)域，尤其是前期資料缺乏無法實(shí)施水準(zhǔn)的海外項(xiàng)目奠定基礎(chǔ)，以期達(dá)到快速獲取高程，服務(wù)前期設(shè)計(jì)，減少成本，降低人力消耗的目標(biāo)。

2 基于EGM2008 模型的移去-擬合-恢復(fù)法確定區(qū)域似大地水準(zhǔn)面

2.1 地球重力場(chǎng)模型確定似大地水準(zhǔn)面的基本原理

重力場(chǎng)模型可以通過多種表達(dá)式來表示，其中使用最多的是球諧展開式。地球重力場(chǎng)模型實(shí)則為調(diào)和函數(shù)，該函數(shù)能夠不斷逼近地球質(zhì)體外部引力位，并且在無窮遠(yuǎn)處收斂，極值為零。一般情況下重力場(chǎng)模型可以展開成一個(gè)在理論上收斂的整階次球諧函數(shù)的無窮級(jí)數(shù)形式，而由級(jí)數(shù)的系數(shù)所組成的集合定義了一個(gè)與之相對(duì)應(yīng)的地球重力場(chǎng)模型。

使用現(xiàn)有的重力場(chǎng)模型可以很容易的獲取大地水準(zhǔn)面差距和重力異常值。利用全球重力場(chǎng)模型可以方便的求出任何一點(diǎn)的高程異常值，但這個(gè)值通常是含有很大的誤差，和真實(shí)值之間有一個(gè)差值，假設(shè)一定數(shù)量點(diǎn)的真實(shí)高程異常值是已知的，通過真實(shí)值與模型值就可以得到這些差值。對(duì)這些差值進(jìn)行數(shù)學(xué)擬合，用一個(gè)數(shù)學(xué)上的平面或者曲面去逼近這個(gè)值，在擬合出來的曲面上求得待求點(diǎn)的這部分差值，最后恢復(fù)模型計(jì)算出來的中長波部分的高程異常，就可以得到這些點(diǎn)的高程異常。這個(gè)過程被概括為移去擬合恢復(fù)法?；诖?，本文提出了從已知點(diǎn)的高程異常值中去掉高程異常長波部分，對(duì)于剩余的短波部分的高程異常值與中波部分整體看作為殘余高程異常值，對(duì)殘余高程異常值采用數(shù)學(xué)模型進(jìn)行擬合計(jì)算，得到整個(gè)測(cè)區(qū)的殘余高程異常函數(shù)計(jì)算式，再輸入任意一點(diǎn)的平面坐標(biāo)，得到該點(diǎn)的殘余高程異常值，再加上EGM2008 重力場(chǎng)模型計(jì)算的高程異常值的長波項(xiàng)，從而獲取未知點(diǎn)的精確高程異常值，從而實(shí)現(xiàn)了大地高到正常高的精準(zhǔn)轉(zhuǎn)化[1]。

2.2 基于EGM2008 模型確定似大地水準(zhǔn)面的基本原理

該思路綜合利用EGM2008 模型及數(shù)學(xué)擬合方法，無須實(shí)測(cè)重力數(shù)據(jù)，只須聯(lián)測(cè)測(cè)區(qū)少量水準(zhǔn)點(diǎn)，尤其是在地形起伏較大、水準(zhǔn)聯(lián)測(cè)點(diǎn)較少的地區(qū)，能方便地進(jìn)行GPS高程擬合。根據(jù)物理大地測(cè)量學(xué)理論，高程異常ζ 可以表示為：

式中，ζEGM是由 EGM2008模型求得的高程異常長波部分；ζres是殘余高程異常。

“移去-擬合-恢復(fù)”法思路為，從已知點(diǎn)的已知高程異常ζ 中移去ζEGM，獲得ζres 進(jìn)行模型擬合；基于擬合模型計(jì)算待求點(diǎn)的ζres，再加上ζEGM，就可得到待求點(diǎn)的高程異常ζ。

圖1 “移去-擬合-恢復(fù)”法基本流程圖

2.3 擬合優(yōu)化函數(shù)模型

利用高精度EGM2008 重力場(chǎng)模型數(shù)據(jù)，可以迅速獲得測(cè)區(qū)任意一點(diǎn)的重力高程異常。因此求取精準(zhǔn)高程異常值最關(guān)鍵的一步就是殘余高程異常值的求取。

任意一點(diǎn)的殘余高程異常無法準(zhǔn)確獲取，只能通過現(xiàn)有的GPS/水準(zhǔn)測(cè)量方式獲取一定數(shù)量的區(qū)域控制點(diǎn)。通過已知控制點(diǎn)可以獲取該區(qū)域的殘余高程，然后選擇合適的數(shù)學(xué)模型擬合該區(qū)域的優(yōu)化函數(shù)模型[2]。

本文選擇多項(xiàng)式曲面擬合法作為擬合優(yōu)化函數(shù)模型，其原理是：根據(jù)區(qū)域中已知點(diǎn)的平面坐標(biāo)X,Y(或大地坐標(biāo)B、L)和高程異常值?，用數(shù)學(xué)法擬合測(cè)區(qū)似大地水準(zhǔn)面，再內(nèi)插出待求點(diǎn)的高程異常，從而求得待求點(diǎn)的正常高。多項(xiàng)式曲面擬合的一般模型為：

式中a0、a1、a2、a3、a4、a5、a6、a7、a8、a9、為模型待求參數(shù)，?為已知點(diǎn)高程異常值。

當(dāng)已知控制點(diǎn)為n 時(shí)，其方程式可寫成：

即：

通過高斯消元求出模型參數(shù)A，然后代入待求點(diǎn)的平面坐標(biāo)計(jì)算出高程異常值，進(jìn)而計(jì)算出正常高Hr。

當(dāng)已知控制點(diǎn)數(shù)目較多時(shí)，采用最小二乘法計(jì)算待求參數(shù)A，其高程異常值?與已知控制點(diǎn)平面坐標(biāo)x、y 關(guān)系式可寫成：

根據(jù)最小二乘法原理，可計(jì)算得A=(X TX)-1XT?，進(jìn)而代入待求點(diǎn)平面坐標(biāo)，計(jì)算其高程異常值，然后求出其正常高。

3 平臺(tái)開發(fā)及案例分析

根據(jù)上面所述內(nèi)容，結(jié)合本文提出的基于EGM2008 模型的移去-擬合-恢復(fù)法求取高程異常值的核心思想，進(jìn)行了相關(guān)程序的編寫與計(jì)算平臺(tái)的開發(fā)，以人機(jī)相互的方式實(shí)現(xiàn)高程異常值的計(jì)算。開發(fā)的平臺(tái)名稱為“基于EGM2008模型處理DEM 數(shù)據(jù)的程序設(shè)計(jì)”[3]。

3.1 計(jì)算平臺(tái)的功能

本程序計(jì)算平臺(tái)，利用EGM2008 模型及優(yōu)化函數(shù)擬合，可以實(shí)現(xiàn)任何一個(gè)區(qū)域內(nèi)DEM 數(shù)據(jù)大地高到正常高的轉(zhuǎn)換，在此基礎(chǔ)上，在測(cè)區(qū)附近有一定數(shù)量已知點(diǎn)的前提下，可以求取區(qū)域內(nèi)高程異常值數(shù)學(xué)模型，從而進(jìn)一步精確求取DEM數(shù)據(jù)的正常高，進(jìn)而滿足實(shí)際工程項(xiàng)目需求。

圖2 程序設(shè)計(jì)界面

3.2 案例分析

現(xiàn)以開發(fā)的“基于EGM2008 模型處理DEM 數(shù)據(jù)的程序設(shè)計(jì)”平臺(tái)為計(jì)算手段，以中國福建省某岸線規(guī)劃類項(xiàng)目為測(cè)試案例，將DEM 點(diǎn)云數(shù)據(jù)進(jìn)行了大地高到正常高的轉(zhuǎn)換，并將比對(duì)數(shù)據(jù)進(jìn)行了分析與研究[4]。

本案例測(cè)區(qū)地處閩西南上古生代覆蓋層低山丘陵地貌，本次選取了8 個(gè)已知控制點(diǎn)作為起算數(shù)據(jù)，控制點(diǎn)分布為面狀分布，且基本能夠均勻覆蓋測(cè)區(qū)，控制點(diǎn)示意圖如圖3所示?？刂泣c(diǎn)高程值采用的四等水準(zhǔn)連測(cè)獲取，控制點(diǎn)成果如表1所示。

本項(xiàng)目采用機(jī)載LiDAR 技術(shù)獲取了測(cè)區(qū)地面基于WGS84 橢球大地高的點(diǎn)云數(shù)據(jù)，現(xiàn)進(jìn)行如下計(jì)算：

圖3 控制點(diǎn)位置示意圖

表1 控制點(diǎn)成果表

方案一：將DEM 數(shù)據(jù)直接通過EGM2008 模型計(jì)算器獲取EGM2008 高程異常擬合正常高；

方案二：選取測(cè)區(qū)內(nèi)所有控制點(diǎn)（8 個(gè)）作為已知點(diǎn)，將DEM 數(shù)據(jù)通過本科研項(xiàng)目開發(fā)的計(jì)算平臺(tái)進(jìn)行處理獲取正常高；

方案三：選取測(cè)區(qū)中心一個(gè)控制點(diǎn)（CT03）作為已知點(diǎn)，基于EGM2008 地球重力模型進(jìn)行單點(diǎn)校正，通過計(jì)算平臺(tái)將DEM 點(diǎn)云數(shù)據(jù)由大地高轉(zhuǎn)換成正常高；

將上述三種方案的計(jì)算后的DEM 數(shù)據(jù)與DEM 處理后的水準(zhǔn)高程值進(jìn)行對(duì)比，對(duì)比結(jié)果如下：

表2 各計(jì)算方案精度指標(biāo)

圖4 擬合成果部分展示圖

由計(jì)算成果可以看出，直接用EGM2008 模型擬合計(jì)算的DEM 點(diǎn)云數(shù)據(jù)的正常高誤差值較大，擬合中誤差為35.56cm；采用本次科研項(xiàng)目基于EGM2008 模型的移去-擬合-恢復(fù)法為核心思想而開發(fā)的計(jì)算平臺(tái)擬合精度較高，采用測(cè)區(qū)全部8 個(gè)控制點(diǎn)時(shí)，擬合中誤差為3.35cm，采用單個(gè)控制點(diǎn)擬合中誤差為3.85cm。

4 結(jié)語

通過不同地域、不同地形條件下的案例測(cè)試，可以得出無論是在平原、丘陵還是高山區(qū)，基于EGM2008 模型的移去-擬合-恢復(fù)法的擬合計(jì)算方法可以實(shí)現(xiàn)DEM 數(shù)據(jù)的由大地高向正常高的高精度轉(zhuǎn)換，且轉(zhuǎn)換精度能夠達(dá)到10cm 以內(nèi)。

在一些由于自然環(huán)境惡劣而造成作業(yè)困難的項(xiàng)目地區(qū)，尤其是一些國際項(xiàng)目中，地處偏遠(yuǎn)的高山區(qū)，已知水準(zhǔn)點(diǎn)極少，高程控制難以布設(shè)，加之工期非常緊張，而又亟待快速獲取基礎(chǔ)地形圖資料。在無法獲得足夠的控制點(diǎn)來進(jìn)行擬合或轉(zhuǎn)換參數(shù)求取的情況下，可采用基于EGM2008 模型的移去-擬合-恢復(fù)法的擬合計(jì)算方法來快速實(shí)現(xiàn)測(cè)區(qū)DEM 數(shù)據(jù)的大地高向正常高的轉(zhuǎn)換，從而提高項(xiàng)目生產(chǎn)效益。