吳 勇，常燕敏

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

0 引言

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

1 原理與方法

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

ζ=ζEGM+ζres

(1)

高程異常長(zhǎng)波項(xiàng)(ζEGM)可以利用EGM2008模型計(jì)算求得，殘余高程異常值(ζres)需要利用測(cè)區(qū)內(nèi)已知點(diǎn)的殘余高程異常值通過(guò)數(shù)學(xué)函數(shù)擬合求得。

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

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

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

1.1 求解高程異常長(zhǎng)波項(xiàng)

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

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

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

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

1.2 求解殘余高程異常

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

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

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

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

(2)

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

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

ζres=a0+a1x+a2y

(3)

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

二次曲面擬合主要運(yùn)用在地形起伏不大，控制點(diǎn)較少的情況下。因?yàn)槎吻鏀M合的模型比較平滑，可以認(rèn)為在一定范圍高程異常是平緩的變化。二次曲面函數(shù)表達(dá)式為：

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

(4)

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

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

2.1 計(jì)算高程異常長(zhǎng)波項(xiàng)

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

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

按行循環(huán)讀取格網(wǎng)高程異常長(zhǎng)波數(shù)據(jù)，調(diào)用距離計(jì)算函數(shù)(get＿dist)計(jì)算每個(gè)格網(wǎng)點(diǎn)到計(jì)算點(diǎn)的距離，將距離計(jì)算結(jié)果添加到列表中，直到添加的距離個(gè)數(shù)等于K(搜索格附件網(wǎng)點(diǎn)個(gè)數(shù))值時(shí)，再將計(jì)算的結(jié)果與列表中的最大值進(jìn)行對(duì)比，如果小于列表中的最大值就將該值替換列表的最大值，實(shí)現(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計(jì)算函數(shù)擬合系數(shù)

主要調(diào)用 Python 的leastsq庫(kù)對(duì)給定的三維數(shù)據(jù)點(diǎn)進(jìn)行最小二乘擬合，需先定義擬合函數(shù)和錯(cuò)誤函數(shù)，利用numpy庫(kù)讀取已知點(diǎn)數(shù)據(jù)，實(shí)現(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 案例計(jì)算

3.1 高程異常長(zhǎng)波項(xiàng)精度檢驗(yàn)

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

表1 部分點(diǎn)高程異常長(zhǎng)波項(xiàng)計(jì)算對(duì)比表Tab.1 Comparison of long-wave component calculations for height anomalies at selected points

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

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

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

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

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

圖4 控制點(diǎn)和檢核點(diǎn)分布圖Fig.4 Distribution of control points and check points

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

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

4 結(jié)束語(yǔ)

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