XGM2019模型在似大地水準(zhǔn)面精化中的應(yīng)用

丁渃鵬,楊久東,張凌云

(華北理工大學(xué) 礦業(yè)工程學(xué)院,河北 唐山 063210)

引言

將GPS高程測(cè)量測(cè)得的大地高轉(zhuǎn)換為符合國(guó)家高程系統(tǒng)的正常高,其關(guān)鍵在于如何提高高程異常的求解,提高的轉(zhuǎn)換方法通常有使用全球超高階地球重力場(chǎng)模型,EGM2008是由美國(guó)國(guó)家地理空間情報(bào)局(NGA)EGM(Earth Gravitational Model)團(tuán)隊(duì)公開發(fā)布的全球重力場(chǎng)模型,該模型的空間分辨率為5′,該模型的階次完全至2 159(另外球諧系數(shù)的階擴(kuò)展至2 190)[1]。XGM2019作為EGM2020的實(shí)驗(yàn)版,是由慕尼黑技術(shù)大學(xué)天文和物理大地測(cè)量研究所在2019年發(fā)布[2],2020年由ICGEM(International Center for Global Gravity Field Models)發(fā)布了XGM2019的若干版本,其中XGM2019_2159是與EGM2008相同階次的全球重力場(chǎng)模型,XGM2019e的階次完全至5 399(另外球諧系數(shù)的階擴(kuò)展至5 540),相比于EGM2008,XGM2019e在地面重力數(shù)據(jù)上提升了包含一個(gè)15′的重力異常網(wǎng)格,在陸地上添加了地形信息,在海洋上添加了高程異常數(shù)據(jù)[3]。

該項(xiàng)研究通過全球重力場(chǎng)模型EGM2008與XGM2019計(jì)算出擾動(dòng)位,根據(jù)布隆斯公式計(jì)算出不同階次下的高程異常的長(zhǎng)波部分[4],根據(jù)實(shí)測(cè)的GNSS/水準(zhǔn)測(cè)量數(shù)據(jù)計(jì)算出高程異常值,用于計(jì)算對(duì)比出2種模型下高程異常的長(zhǎng)波部分與GNSS/水準(zhǔn)測(cè)量計(jì)算的高程異常值之間的差距,從而評(píng)估XGM2019相較于EGM2008在直接轉(zhuǎn)換為正常高精度上的提升。根據(jù)物理大地測(cè)量學(xué)關(guān)于高程異常的相關(guān)理論,對(duì)EGM2008與XGM2019使用"移去-擬合-恢復(fù)"的方法及使用顧及由于地形起伏變化引起的高程異常的短波部分,通過使用二次曲面法對(duì)高程異常的殘余部分進(jìn)行擬合,并對(duì)2種重力場(chǎng)模型計(jì)算的結(jié)果進(jìn)行對(duì)比。

1 相關(guān)原理

1.1 高程異常

在使用GPS進(jìn)行定位時(shí),所獲得的平面坐標(biāo)和高程值分別是大地經(jīng)度,大地緯度和大地高,水準(zhǔn)測(cè)量獲得的高程值為正常高。其中大地高HEH是地面上的點(diǎn)沿法線的方向延伸到參考橢球面上的距離,而正常高HNH是地面上的點(diǎn)沿鉛垂線的方向延伸到似大地水準(zhǔn)面上的距離。由于我國(guó)的高程基準(zhǔn)采用1985年國(guó)家高程基準(zhǔn),并且高程系統(tǒng)為正常高系統(tǒng),當(dāng)已知高程異常時(shí),可以通過式(1)將任一點(diǎn)的大地高轉(zhuǎn)為同一點(diǎn)的正常高。

ζ=HEH-HNH

(1)

1.2 "移去-擬合-恢復(fù)"法

根據(jù)物理大地測(cè)量學(xué)理論,高程異常ζ可以通過式(2)進(jìn)行表達(dá)[1]

ζ=ζEGM+ζDEM+ζres

(2)

式中第一項(xiàng)是高程異常的長(zhǎng)波部分,是通過全球重力模型計(jì)算擾動(dòng)位,并根據(jù)布隆斯公式計(jì)算得出;ζDEM是由于地形起伏引起的高程異常的短波部分,用于描述20 km以下的短波特征[5];ζres是殘余的高程異常部分,也稱為局部似大地水準(zhǔn)面中波特征,可以通過重力異常邊值問題計(jì)算得到。

"移去-擬合-恢復(fù)"法的主要思想是將第一項(xiàng)與第二項(xiàng)進(jìn)行合并,并從已知點(diǎn)的ζ中進(jìn)行移去,單獨(dú)剩下ζres,再通過適當(dāng)?shù)臄?shù)學(xué)擬合方法,對(duì)ζres進(jìn)行擬合,最后在擬合值上加上先前移去的ζEGM與ζDEM,此時(shí)就可以得到"移去-擬合-恢復(fù)"法后的高程異常值ζ,具體流程可見圖 1所示。

圖1 "移去-擬合-恢復(fù)"法流程圖

1.3 高程異常的長(zhǎng)波項(xiàng)

由全球超高階重力場(chǎng)模型EGM2008,XGM2019_2159和XGM2019e分別計(jì)算擾動(dòng)位,配合布隆斯公式計(jì)算出高程異常的長(zhǎng)波項(xiàng),其中EGM2008與XGM2019_2159兩者的階次完全至2 159(另外球諧系數(shù)的階拓展至2 190,次為2 159)[1],而XGM2019e的階次完全至5 399(另外球諧系數(shù)的階拓展至5540,次為5 399)[3],并且以上模型均可以進(jìn)行截?cái)?。?jì)算高程異常的長(zhǎng)波項(xiàng)[1]可參考式(3)。

(3)

1.4 地形改正

由地形起伏變化引起的高程異常變化量δ可以通過式(4)進(jìn)行表達(dá)[6]

(4)

公式(4)中,G為萬有引力常數(shù);ρ為地球質(zhì)量密度;hi為流動(dòng)單元的高程;h為待求點(diǎn)高程;γ為待求點(diǎn)正常重力;l為待求點(diǎn)到流動(dòng)單元的距離;dxdy為劃分單元的面積。

高程異常短波項(xiàng) 的求解步驟如下:

(1)獲取研究區(qū)域的DEM數(shù)據(jù),目前DEM數(shù)據(jù)可以通過地理空間數(shù)據(jù)云中進(jìn)行下載,目前提供90 m分辨率SRTM數(shù)據(jù)和30 m分辨率的GDTM數(shù)據(jù)。

(2)根據(jù)研究區(qū)域范圍劃分流動(dòng)單元。

(3)根據(jù)公式(4)計(jì)算出每一個(gè)格網(wǎng)點(diǎn)由地形起伏變化引起的高程變化量。

(4)根據(jù)格網(wǎng)的雙線性內(nèi)插公式計(jì)算出待定點(diǎn)高程異常的短波項(xiàng),計(jì)算公式如式(5)[1]。

(5)

公式(5)中δa,δb,δc,δd為流動(dòng)單元4個(gè)格網(wǎng)點(diǎn)由地形起伏變化引起的高程異常變化量,其順序依次是左下角、左上角、右上角、右下角;l為流動(dòng)單元的邊長(zhǎng);Δx,Δy為待定點(diǎn)到流動(dòng)單元左下角格網(wǎng)點(diǎn)的坐標(biāo)差。

2 多項(xiàng)式曲面擬合

2.1 二次曲面擬合

在使用"移去-擬合-恢復(fù)"方法中,通過對(duì)殘余的高程異常進(jìn)行數(shù)學(xué)擬合的方法有很多,常用的數(shù)學(xué)擬合方法有,多項(xiàng)式曲線擬合,多面函數(shù)法和多項(xiàng)式曲面擬合,其中二次曲面擬合適用于測(cè)區(qū)比較平坦的情況,二次曲面擬合方程如式(6)所示:

ζres=a0+a1B+a2L+a3BL+a4B2+a5L2

(6)

公式(6)中ζres為殘余的高程異常的擬合值,a0、a1、a2、a3、a4、a5為多項(xiàng)式的擬合系數(shù),通過最小二乘求得,其中要求至少需要6個(gè)聯(lián)測(cè)的GNSS/水準(zhǔn)點(diǎn),由模型顯著性檢驗(yàn)的結(jié)果可知,系數(shù)只取到二次項(xiàng)[6]。

3 實(shí)驗(yàn)分析

3.1 數(shù)據(jù)介紹

實(shí)驗(yàn)數(shù)據(jù)選取于河南省平頂山市某實(shí)習(xí)基地,共采集了30個(gè)GPS靜態(tài)數(shù)據(jù)和其聯(lián)測(cè)的四等水準(zhǔn)測(cè)量數(shù)據(jù),其中大地高的最大值為137.09 m最小值為96.87 m,水準(zhǔn)測(cè)量的正常高最大值為157.12 m,最小值為116.94 m,具體點(diǎn)位分布如圖2和圖3所示。其中選取21個(gè)點(diǎn)作為"移去-擬合-恢復(fù)"當(dāng)中的擬合點(diǎn),用于對(duì)殘余高程異常部分進(jìn)行擬合,剩余9個(gè)點(diǎn)用于檢驗(yàn),利用內(nèi)、外符合精度來評(píng)價(jià)建立模型的質(zhì)量。在顧及地形起伏變化引起的高程異常變化量當(dāng)中,DEM數(shù)據(jù)選取的來自地理數(shù)據(jù)云當(dāng)中的SRTM90 m分辨率數(shù)據(jù),使用流動(dòng)單元?jiǎng)澐譃?的格網(wǎng)來求解由于地形起伏變化引起的高程異常變化量。

圖2 點(diǎn)位分布圖

圖3 點(diǎn)位地形起伏圖

3.2 實(shí)驗(yàn)設(shè)計(jì)與結(jié)果

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

圖4 2種高程異常的長(zhǎng)波部分與高程異常的殘差圖

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

圖5 全球重力場(chǎng)模型在不同階數(shù)下均方根誤差/m

實(shí)驗(yàn)三:利用EGM2008,XGM2019_2159,XGM2019e進(jìn)行"移去-擬合-恢復(fù)"法來驗(yàn)證不同模型對(duì)似大地水準(zhǔn)面精化精度上的提升,選取數(shù)據(jù)當(dāng)中21個(gè)點(diǎn)作為二次多項(xiàng)式曲面擬合的擬合數(shù)據(jù)點(diǎn),采用9個(gè)點(diǎn)作為檢核點(diǎn),使用內(nèi)符合精度與外符合精度來評(píng)判建模水平,內(nèi)符合精度和外符合精度計(jì)算公式如式(7)和式(8)所示。

(7)

(8)

式中μ表示內(nèi)符合精度,v為已知點(diǎn)與擬合點(diǎn)的殘差值,n為參與擬合點(diǎn)的個(gè)數(shù)。RMS表示外符合精度,Δ為已知點(diǎn)與檢驗(yàn)點(diǎn)的殘差值,m為檢核點(diǎn)的個(gè)數(shù)。內(nèi)符合精度反映了利用已知點(diǎn)進(jìn)行建模的模型質(zhì)量,外符合精度反映了利用所建模型求解檢核點(diǎn)的整體質(zhì)量。利用"移去-擬合-恢復(fù)"法計(jì)算的內(nèi)、外符合精度計(jì)算結(jié)果如表 1所示。

表1 高程異常的殘余部分/m

2種全球重力場(chǎng)模型在不同階次下使用"移去-擬合-恢復(fù)"法中內(nèi)符合精度基本一致,在外符合精度上,XGM2019優(yōu)于EGM2008和XGM2019e。

實(shí)驗(yàn)四:利用EGM2008,XGM2019_2159,XGM2019e進(jìn)行"移去-擬合-恢復(fù)"法在顧及由于地形起伏變化引起的高程異常的短波部分,驗(yàn)證不同模型下似大地水準(zhǔn)面精化上的提升,內(nèi)符合精度與外符合精度的計(jì)算結(jié)果如表2所示。

表1 顧及地形改正的高程異常的殘余部分/m

XGM2019的內(nèi)符合精度略低于EGM2008與XGM2019e,但是其外符合精度均優(yōu)于EGM2008與XGM2019e,說明XGM2019的更適合研究區(qū)域的似大地水準(zhǔn)面的精化工作。

4 結(jié)論

(1)GM2020的實(shí)驗(yàn)版本XGM2019當(dāng)其階次在2 190階時(shí),相比EGM2008在高程異常長(zhǎng)波部分上有厘米級(jí)的提升,適用于精度要求不高時(shí)將大地高直接轉(zhuǎn)換為正常高使用,表明在研究區(qū)域內(nèi)使用XGM2019下的似大地水準(zhǔn)面相比EGM2008在符合程度上更優(yōu)的結(jié)論。但由于XGM2019為實(shí)驗(yàn)版本,所以在其他地方的適用性還有待研究。

(2)XGM2019當(dāng)其階次在2 190時(shí),其高程異常長(zhǎng)波項(xiàng)與高程異常的均方根誤差已經(jīng)達(dá)到厘米級(jí),隨著重力場(chǎng)模型的階次不斷提高,其精度有厘米級(jí)提升,但對(duì)于完全規(guī)格化締合函數(shù)的選擇就更為重要。在ICGEM中提供的完全規(guī)格化位系數(shù)高達(dá)百萬行,對(duì)于數(shù)據(jù)預(yù)處理也是不小的挑戰(zhàn)。

(3)"移去-擬合-恢復(fù)"方法中,XGM2019的內(nèi)外符合精度均優(yōu)于同階次的EGM2008模型與4 700階的XGM2019e模型,但在顧及地形起伏變化的"移去-擬合-恢復(fù)"方法中相比未進(jìn)行顧及地形起伏變化的"移去-擬合-恢復(fù)法"在外符合精度上的提升不大。