丁渃鵬,楊久東,張凌云
(華北理工大學 礦業(yè)工程學院,河北 唐山 063210)
將GPS高程測量測得的大地高轉(zhuǎn)換為符合國家高程系統(tǒng)的正常高,其關(guān)鍵在于如何提高高程異常的求解,提高的轉(zhuǎn)換方法通常有使用全球超高階地球重力場模型,EGM2008是由美國國家地理空間情報局(NGA)EGM(Earth Gravitational Model)團隊公開發(fā)布的全球重力場模型,該模型的空間分辨率為5′,該模型的階次完全至2 159(另外球諧系數(shù)的階擴展至2 190)[1]。XGM2019作為EGM2020的實驗版,是由慕尼黑技術(shù)大學天文和物理大地測量研究所在2019年發(fā)布[2],2020年由ICGEM(International Center for Global Gravity Field Models)發(fā)布了XGM2019的若干版本,其中XGM2019_2159是與EGM2008相同階次的全球重力場模型,XGM2019e的階次完全至5 399(另外球諧系數(shù)的階擴展至5 540),相比于EGM2008,XGM2019e在地面重力數(shù)據(jù)上提升了包含一個15′的重力異常網(wǎng)格,在陸地上添加了地形信息,在海洋上添加了高程異常數(shù)據(jù)[3]。
該項研究通過全球重力場模型EGM2008與XGM2019計算出擾動位,根據(jù)布隆斯公式計算出不同階次下的高程異常的長波部分[4],根據(jù)實測的GNSS/水準測量數(shù)據(jù)計算出高程異常值,用于計算對比出2種模型下高程異常的長波部分與GNSS/水準測量計算的高程異常值之間的差距,從而評估XGM2019相較于EGM2008在直接轉(zhuǎn)換為正常高精度上的提升。根據(jù)物理大地測量學關(guān)于高程異常的相關(guān)理論,對EGM2008與XGM2019使用"移去-擬合-恢復"的方法及使用顧及由于地形起伏變化引起的高程異常的短波部分,通過使用二次曲面法對高程異常的殘余部分進行擬合,并對2種重力場模型計算的結(jié)果進行對比。
在使用GPS進行定位時,所獲得的平面坐標和高程值分別是大地經(jīng)度,大地緯度和大地高,水準測量獲得的高程值為正常高。其中大地高HEH是地面上的點沿法線的方向延伸到參考橢球面上的距離,而正常高HNH是地面上的點沿鉛垂線的方向延伸到似大地水準面上的距離。由于我國的高程基準采用1985年國家高程基準,并且高程系統(tǒng)為正常高系統(tǒng),當已知高程異常時,可以通過式(1)將任一點的大地高轉(zhuǎn)為同一點的正常高。
ζ=HEH-HNH
(1)
根據(jù)物理大地測量學理論,高程異常ζ可以通過式(2)進行表達[1]
ζ=ζEGM+ζDEM+ζres
(2)
式中第一項是高程異常的長波部分,是通過全球重力模型計算擾動位,并根據(jù)布隆斯公式計算得出;ζDEM是由于地形起伏引起的高程異常的短波部分,用于描述20 km以下的短波特征[5];ζres是殘余的高程異常部分,也稱為局部似大地水準面中波特征,可以通過重力異常邊值問題計算得到。
"移去-擬合-恢復"法的主要思想是將第一項與第二項進行合并,并從已知點的ζ中進行移去,單獨剩下ζres,再通過適當?shù)臄?shù)學擬合方法,對ζres進行擬合,最后在擬合值上加上先前移去的ζEGM與ζDEM,此時就可以得到"移去-擬合-恢復"法后的高程異常值ζ,具體流程可見圖 1所示。

圖1 "移去-擬合-恢復"法流程圖
由全球超高階重力場模型EGM2008,XGM2019_2159和XGM2019e分別計算擾動位,配合布隆斯公式計算出高程異常的長波項,其中EGM2008與XGM2019_2159兩者的階次完全至2 159(另外球諧系數(shù)的階拓展至2 190,次為2 159)[1],而XGM2019e的階次完全至5 399(另外球諧系數(shù)的階拓展至5540,次為5 399)[3],并且以上模型均可以進行截斷。計算高程異常的長波項[1]可參考式(3)。
(3)

由地形起伏變化引起的高程異常變化量δ可以通過式(4)進行表達[6]
(4)
公式(4)中,G為萬有引力常數(shù);ρ為地球質(zhì)量密度;hi為流動單元的高程;h為待求點高程;γ為待求點正常重力;l為待求點到流動單元的距離;dxdy為劃分單元的面積。
高程異常短波項 的求解步驟如下:
(1)獲取研究區(qū)域的DEM數(shù)據(jù),目前DEM數(shù)據(jù)可以通過地理空間數(shù)據(jù)云中進行下載,目前提供90 m分辨率SRTM數(shù)據(jù)和30 m分辨率的GDTM數(shù)據(jù)。
(2)根據(jù)研究區(qū)域范圍劃分流動單元。
(3)根據(jù)公式(4)計算出每一個格網(wǎng)點由地形起伏變化引起的高程變化量。
(4)根據(jù)格網(wǎng)的雙線性內(nèi)插公式計算出待定點高程異常的短波項,計算公式如式(5)[1]。
(5)
公式(5)中δa,δb,δc,δd為流動單元4個格網(wǎng)點由地形起伏變化引起的高程異常變化量,其順序依次是左下角、左上角、右上角、右下角;l為流動單元的邊長;Δx,Δy為待定點到流動單元左下角格網(wǎng)點的坐標差。
在使用"移去-擬合-恢復"方法中,通過對殘余的高程異常進行數(shù)學擬合的方法有很多,常用的數(shù)學擬合方法有,多項式曲線擬合,多面函數(shù)法和多項式曲面擬合,其中二次曲面擬合適用于測區(qū)比較平坦的情況,二次曲面擬合方程如式(6)所示:
ζres=a0+a1B+a2L+a3BL+a4B2+a5L2
(6)
公式(6)中ζres為殘余的高程異常的擬合值,a0、a1、a2、a3、a4、a5為多項式的擬合系數(shù),通過最小二乘求得,其中要求至少需要6個聯(lián)測的GNSS/水準點,由模型顯著性檢驗的結(jié)果可知,系數(shù)只取到二次項[6]。
實驗數(shù)據(jù)選取于河南省平頂山市某實習基地,共采集了30個GPS靜態(tài)數(shù)據(jù)和其聯(lián)測的四等水準測量數(shù)據(jù),其中大地高的最大值為137.09 m最小值為96.87 m,水準測量的正常高最大值為157.12 m,最小值為116.94 m,具體點位分布如圖2和圖3所示。其中選取21個點作為"移去-擬合-恢復"當中的擬合點,用于對殘余高程異常部分進行擬合,剩余9個點用于檢驗,利用內(nèi)、外符合精度來評價建立模型的質(zhì)量。在顧及地形起伏變化引起的高程異常變化量當中,DEM數(shù)據(jù)選取的來自地理數(shù)據(jù)云當中的SRTM90 m分辨率數(shù)據(jù),使用流動單元劃分為 的格網(wǎng)來求解由于地形起伏變化引起的高程異常變化量。

圖2 點位分布圖

圖3 點位地形起伏圖
實驗一:使用與EGM2008同階的全球重力場模型XGM2019_2159,利用公式(1)直接計算30個點位的高程異常的長波部分,并與GNSS/水準計算所得的高程異常值做差,從而確定由全球重力場模型計算出的高程異常長波部分與GNSS/水準計算的高程異常之間的差距。計算結(jié)果如圖 4所示,從圖4可看出,在使用XGM2019_2159計算下的殘差明顯小于EGM2008,其精度有近1倍提升,故在研究區(qū)域內(nèi)使用XGM2019模型下的似大地水準面的符合程度相比于EGM2008更優(yōu)。

圖4 2種高程異常的長波部分與高程異常的殘差圖
實驗二:利用ICGEM提供的不同階次的XGM2019模型,驗證階次提升的情況下,對高程異常殘差值的影響,研究選用EGM2008(階次為2190),XGM2019_2159(階次為2190)和XGM2019e(階次為4 700),并用均方根誤差作為評判依據(jù)。計算結(jié)果如圖5所示,當階次達到4 700時,其均方根誤差是明顯的優(yōu)于前兩項的,故使用的更高階次的全球重力場模型進行高程異常長波部分的計算更能提升似大地水準面的精化精度。

圖5 全球重力場模型在不同階數(shù)下均方根誤差/m
實驗三:利用EGM2008,XGM2019_2159,XGM2019e進行"移去-擬合-恢復"法來驗證不同模型對似大地水準面精化精度上的提升,選取數(shù)據(jù)當中21個點作為二次多項式曲面擬合的擬合數(shù)據(jù)點,采用9個點作為檢核點,使用內(nèi)符合精度與外符合精度來評判建模水平,內(nèi)符合精度和外符合精度計算公式如式(7)和式(8)所示。
(7)
(8)
式中μ表示內(nèi)符合精度,v為已知點與擬合點的殘差值,n為參與擬合點的個數(shù)。RMS表示外符合精度,Δ為已知點與檢驗點的殘差值,m為檢核點的個數(shù)。內(nèi)符合精度反映了利用已知點進行建模的模型質(zhì)量,外符合精度反映了利用所建模型求解檢核點的整體質(zhì)量。利用"移去-擬合-恢復"法計算的內(nèi)、外符合精度計算結(jié)果如表 1所示。

表1 高程異常的殘余部分/m
2種全球重力場模型在不同階次下使用"移去-擬合-恢復"法中內(nèi)符合精度基本一致,在外符合精度上,XGM2019優(yōu)于EGM2008和XGM2019e。
實驗四:利用EGM2008,XGM2019_2159,XGM2019e進行"移去-擬合-恢復"法在顧及由于地形起伏變化引起的高程異常的短波部分,驗證不同模型下似大地水準面精化上的提升,內(nèi)符合精度與外符合精度的計算結(jié)果如表2所示。

表1 顧及地形改正的高程異常的殘余部分/m
XGM2019的內(nèi)符合精度略低于EGM2008與XGM2019e,但是其外符合精度均優(yōu)于EGM2008與XGM2019e,說明XGM2019的更適合研究區(qū)域的似大地水準面的精化工作。
(1)GM2020的實驗版本XGM2019當其階次在2 190階時,相比EGM2008在高程異常長波部分上有厘米級的提升,適用于精度要求不高時將大地高直接轉(zhuǎn)換為正常高使用,表明在研究區(qū)域內(nèi)使用XGM2019下的似大地水準面相比EGM2008在符合程度上更優(yōu)的結(jié)論。但由于XGM2019為實驗版本,所以在其他地方的適用性還有待研究。
(2)XGM2019當其階次在2 190時,其高程異常長波項與高程異常的均方根誤差已經(jīng)達到厘米級,隨著重力場模型的階次不斷提高,其精度有厘米級提升,但對于完全規(guī)格化締合函數(shù)的選擇就更為重要。在ICGEM中提供的完全規(guī)格化位系數(shù)高達百萬行,對于數(shù)據(jù)預處理也是不小的挑戰(zhàn)。
(3)"移去-擬合-恢復"方法中,XGM2019的內(nèi)外符合精度均優(yōu)于同階次的EGM2008模型與4 700階的XGM2019e模型,但在顧及地形起伏變化的"移去-擬合-恢復"方法中相比未進行顧及地形起伏變化的"移去-擬合-恢復法"在外符合精度上的提升不大。