利用Python編程實現基于EGM2008模型的高程異常計算*

吳 勇，常燕敏

(四川省地質礦產勘查開發(fā)局測繪隊，四川 成都 610017)

0 引言

全球導航衛(wèi)星系統(tǒng)(GNSS)以其定位快、精度高等優(yōu)勢在測繪、水利、交通、資源調查等領域被廣泛應用，特別是連續(xù)運行參考站(CORS)已經成為外業(yè)數據采集的主要測量手段[2]，但在高程方面，GNSS提供的不是正常高，而是大地高。將大地高轉換成正常高，關鍵是高程異常的求解。高程異常分解成長波項、短波部分和殘余部分，利用美國國家地理空間情報局(National Geospatial-intelligence Agency，NGA)研制的全球超高階地球重力場模型EGM2008，能夠精確表達高程異常的長波項，地形變化引起的短波部分可以利用剩余地形模型計算所得[3]，也可以利用數字高程模型(Digital Elevation Model，DEM)格網數據求解所得。本研究主要運用在地形起伏不大，滿足一定范圍內的測繪需要，因此不單獨考慮短波部分，將高程異常分解為高程異常長波項和殘余高程異常值。采用“移去-擬合-恢復”法的思想，利用Python編程實現高程異常長波項的計算以及利用測區(qū)內極少數的GNSS水準點的殘余高程異常值和一定的擬合函數計算出待求點的殘余高程異常，從而求得待求點的高程異常值。

1 原理與方法

本文選擇 “移去-擬合-恢復”法求解高程異常，高程異常值分解為高程異常長波項(ζEGM)和殘余高程異常值(ζres)。高程異常ζ可以表示為：

ζ=ζEGM+ζres

(1)

高程異常長波項(ζEGM)可以利用EGM2008模型計算求得，殘余高程異常值(ζres)需要利用測區(qū)內已知點的殘余高程異常值通過數學函數擬合求得。

根據 “移去-擬合-恢復”的思想，先利用EGM2008模型計算已知GNSS水準點的高程異常長波項(ζEGM)[2]，從已知點的高程異常(ζ)中移去已知點的ζEGM，獲得已知點的ζres，將已知點的ζres代入數學函數中計算出擬合系數。

再利用EGM2008模型計算待求點的高程異常長波項(ζEGM)，利用數學模型和擬合系數計算待求點的殘余高程異常值(ζres)。將待求點的高程異常長波項和殘余高程異常值相加，從而求得待求點的高程異常值(ζ)。具體處理流程見圖1。

圖1 高程異常計算流程圖Fig.1 Flow of height anomaly calculation

1.1 求解高程異常長波項

EGM2008重力場模型使用了最先進的算法與建模技術，以 PGM2007B 模型( PGM2007A 的變種模型)為依據，計算中采用了 GRACE 重力衛(wèi)星數據、全球5′×5′重力異常數據、TOPEX 衛(wèi)星測量得到的高程數據、地形數據[3]。目前 EGM2008 全球重力場模型是使用最廣泛的重力場模型，階次完全至2 159。目前可用通過NGA網站下載1′×1′分辨率的格網數據，是目前計算高程異常長波項最理想的模型。

利用EGM2008模型來計算高程異常長波項可以采用現有的軟件(Alltrans EGM2008 Calculator)計算，該軟件調用1′×1′或10′×10′分辨率的格網數據插值計算出高程異常長波項[3]。考慮到程序的實用性和計算數據格式統(tǒng)一性，況且1′×1′分辨率的格網數據量較大(約980 M)，不方便程序運行，因此本文編寫的程序沒有直接利用EGM2008模型求取高程異常長波項，而是利用Alltrans EGM2008 Calculator 1.2軟件導出計算區(qū)域內各格網點的高程異常長波項數據，從格網點的高程異常值長波項數據中通過移動二次曲面擬合求得待求點的高程異常長波項。

在Alltrans EGM2008 Calculator 1.2軟件的“Grid-area”選項輸入需要求取高程異常值區(qū)域的緯度和經度范圍。在“spacing”選項輸入格網間距，可以根據計算精度選擇不同的格網間距，點擊“Calc！”即可導出每個格網點的高程異常長波項數據文件。如圖2所示，左邊是格網數據導出界面，在該界面中輸入相應的參數，右邊則是導出后的格網數據，即是本程序計算時需要調用的數據。

圖2 導出格網數據的界面和導出后數據格式Fig.2 Interface for exporting grid data and exported data format

1.2 求解殘余高程異常

待求點殘余高程的求解主要分成兩步進行，先利用已知點的殘余高程異常值計算擬合函數的系數，再將待求點的數據代入擬合函數中計算出待求點的殘余高程異常值。

1)根據式(1)可以計算出已知GNSS水準點的殘余高程異常值。即殘余高程異常值由高程異常值(ζ)減去高程異常長波部分(ζEGM)，而高程異常值(ζ)由已知的大地高減去正常高。再將已知點坐標x，y和殘余高程異常值代入擬合函數計算出擬合系數?？梢愿鶕y區(qū)的地形起伏、已知點分布和測區(qū)形狀選擇適當的擬合函數進行擬合，常用的數學函數有曲面函數、平面函數、線性函數，各種函數都有一定的局限性和使用范圍，將在后面詳細敘述。

2)將待求點坐標x，y代入上一步選擇的擬合函數中，即可求出待求點的殘余高程異常值。因本次測試數據所在區(qū)域高程起伏不大，不是帶狀分布，因此本文只利用了二次曲面擬合函數進行驗證。在實際使用時可根據求解范圍、區(qū)域地形起伏和已知點數據分布情況選擇不同的數學函數進行擬合。

多項式曲線擬合主要是運用在線路上的高程控制點擬合，控制點呈狹長的帶狀分布，可看作一條曲線，因此選用一個多項式作為差值函數。殘余高程異常值ζres與坐標xi(或yi，i=1，2，…，n)之間函數關系，可用多項式曲線擬合表示為：

(2)

式中：a0，a1，…，am為待定系數。

平面擬合法主要運用于小區(qū)域且地勢較為平坦的區(qū)域，類似于利用平面區(qū)代替局部的似大地水準面。平面擬合表達式為：

ζres=a0+a1x+a2y

(3)

式中：x、y為點的坐標；a0、a1、a2為模型參數。

二次曲面擬合主要運用在地形起伏不大，控制點較少的情況下。因為二次曲面擬合的模型比較平滑，可以認為在一定范圍高程異常是平緩的變化。二次曲面函數表達式為：

ζres=a0+a1x+a2y+a3xy+a4x2+a5y2

(4)

式中：x、y為點的坐標；a0、a1、a2a3、a4、a5為模型參數。

2 編程實現高程異常值計算

2.1 計算高程異常長波項

讀取計算點的經緯和緯度，以計算點為中心在格網高程異常長波數據文件中搜索出距該點最近的12個格網點的數據，建立一個擬合曲面，即將這12個點的數據代入式(4)中求解出擬合系數a0、a1、a2、a3、a4、a5。再將計算點的經度、緯度和擬合系數代入式(4)中，求解出待計算點的高程異常長波項。

2.1.1搜索最近格網點

按行循環(huán)讀取格網高程異常長波數據，調用距離計算函數(get＿dist)計算每個格網點到計算點的距離，將距離計算結果添加到列表中，直到添加的距離個數等于K(搜索格附件網點個數)值時，再將計算的結果與列表中的最大值進行對比，如果小于列表中的最大值就將該值替換列表的最大值，實現代碼如下：

for xx in B＿L＿z：

x=float(xx[0])

y=float(xx[1])

dist＿1=get＿dist(b，l，x，y)

if len(x＿y＿z) < k：

dist＿all.append(dist＿1)

x＿y＿z.append(xx)

else：

if dist＿1 < max(dist＿all)：

site=dist＿all.index(max(dist＿all))

dist＿all[site]，x＿y＿z[site]=dist＿1，xx

2.1.2計算函數擬合系數

主要調用 Python 的leastsq庫對給定的三維數據點進行最小二乘擬合，需先定義擬合函數和錯誤函數，利用numpy庫讀取已知點數據，實現代碼如下：

def func(p，x，y)：

a，b，c，d，e，f=p

return a*x**2+b*y**2+c*x*y+d*x+e*y+f

def error(p，x，y，z)：

return func(p，x，y)-z

b=np.loadtxt(file＿path，delimiter=′，′)

Xi=b[：，0]

Yi=b[：，1]

Zi=b[：，2]

Para=leastsq(error，p0，args=(Xi，Yi，Zi))

a，b，c，d，e，f=Para[0]

3 案例計算

3.1 高程異常長波項精度檢驗

為驗證本軟件高程異常長波項計算的可靠性，本次測試數據采用某市1′×1′的網格點坐標。該市最高點海拔為4 984.1 m，最低點海拔為308 m。共測試30 691個坐標點，測試前利用Alltrans EGM2008 Calculator軟件分別導出0.5′×0.5′和1′×1′的格網點高程異常長波數據，自編軟件分別調用這兩種分辨率的格網數據對待求點進行計算，求得高程異常長波項。再利用Alltrans EGM2008 Calculator軟件中批量轉換獲得的高程異常值長波項與本軟件計算所得的高程異常長波項進行對比。表1列舉了部分點高程異常長波項不同軟件計算的結果。

表1 部分點高程異常長波項計算對比表Tab.1 Comparison of long-wave component calculations for height anomalies at selected points

從表1中可以看出：采用1′×1′的格網數據計算的數據與Alltrans EGM2008 Calculator軟件計算所得的差值最大為7.12 cm，中誤差為1.85 cm。而采用0.5′×0.5′的格網數據計算的值與Alltrans EGM2008 Calculator1.2軟件計算所得的差值最大為4.65 cm，中誤差為0.98 cm。中誤差均小于4.0 cm，滿足區(qū)域似大地水準面精化基本技術規(guī)定中的城市級別的精度要求，因此本文提出的計算方法能滿足高程異常的長波項計算。本文也對差值在不同區(qū)間段的數據進行統(tǒng)計，如圖3所示。根據統(tǒng)計分析得知，采用0.5′×0.5′的格網數據計算的結果明顯優(yōu)于采用1′×1′的格網數據計算的結果，因此在實際生產運用中建議采用0.5′×0.5′的格網數據來計算高程異常的長波項，若只是滿足等外水準的精度要求，也可以采用1′×1′的格網數據來計算。

圖3 長波項差值在不同區(qū)間的統(tǒng)計圖Fig.3 Statistical chart of long-wave component differences in different range of values

3.2 以某市各區(qū)縣的控制點數據計算為例

某市總面積為5 911 km2，呈西北至東南的蠶形，西北為山區(qū)，中部為平原，東南是低山丘陵。考慮到西北山區(qū)沒有收集到控制點，因此只利用了中部和西南地區(qū)共13個C級GPS點(同時具有正常高數據和大地高數據)作為已知點，控制區(qū)域約4 745 km2。點位分布如圖4所示，已知點數據分布比較均勻，基本能覆蓋整個計算區(qū)域。數據中最大高程為763.502 m，最小高程為398.234 m?？紤]到本計算區(qū)域地勢起伏不大，已知點數據也大于6個，因此計算采用二次曲面擬合函數進行計算，根據式(1)計算得出這13個已知點的殘余高程異常最大為22.16 cm，最小為5.11 cm。利用該已知點數據計算二次曲面擬合函數的5個系數。表2列舉了部分已知點數據的計算成果。

表2 部分已知點數據計算結果Tab.2 Data calculated results for some known points

圖4 控制點和檢核點分布圖Fig.4 Distribution of control points and check points

利用實測的182個三等GPS水準點(同時具有正常高數據和大地高數據)作為檢核數據，檢核數據中最大高程為718.583 m，最小高程為329.472 m。將檢核數據的高程異常視為真值。利用本程序計算的高程異常值與真值對比，最大差值為5.86 cm，最小差值為0.007 cm，中誤差為2.56 cm，滿足城市測量規(guī)范中的衛(wèi)星定位高程控制測量中的相關要求。圖5列出了差值分布在不同區(qū)間的個數，其中差值小于3.0 cm占比80%，證明本文的計算方法和自編程序計算的結果比較可靠。

圖5 高程異常值的差值在不同區(qū)間的統(tǒng)計圖Fig.5 Statistical chart of height anomaly differences in different range of values

4 結束語

根據“移去-擬合-恢復”的計算流程，本文提出了基于EGM2008模型，利用Alltrans EGM2008 Calculator軟件導出高程異常格網數據，再利用自編程序利用高程異常格網數據計算出待求點的高程異常值數據，然后利用少量已知點的殘余高程異常通過一定的擬合函數計算出待求點的殘余高程[4]，根據式(1)即可求得待求點的高程異常。通過與某市的實測數據對比，結果表明本文的計算方法和程序能滿足大部分測繪工程的需求。本文未考慮地形因素的影響，如計算區(qū)域的地形起伏較大還應考慮地形因素的影響，后期將進一步完善程序，將地形影響因素也納入程序中。