吳 勇，常燕敏

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

0 引言

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

1 原理與方法

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

ζ=ζEGM+ζres

(1)

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

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

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

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

1.1 求解高程異常長波項

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

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

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

圖2 導出格網(wǎng)數(shù)據(jù)的界面和導出后數(shù)據(jù)格式Fig.2 Interface for exporting grid data and exported data format

1.2 求解殘余高程異常

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

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

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

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

(2)

式中：a0，a1，…，am為待定系數(shù)。

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

ζres=a0+a1x+a2y

(3)

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

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

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

(4)

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

2 編程實現(xiàn)高程異常值計算

2.1 計算高程異常長波項

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

2.1.1搜索最近格網(wǎng)點

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

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計算函數(shù)擬合系數(shù)

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

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 高程異常長波項精度檢驗

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

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

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

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

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

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

表2 部分已知點數(shù)據(jù)計算結(jié)果Tab.2 Data calculated results for some known points

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

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

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

4 結(jié)束語

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