國際高程參考系統(tǒng)在珠峰地區(qū)的實(shí)現(xiàn)

蔣 濤，黨亞民，郭春喜，陳 斌，章傳銀

1. 中國測繪科學(xué)研究院，北京 100830; 2. 自然資源部大地測量數(shù)據(jù)處理中心，陜西 西安 710054; 3. 中國地質(zhì)調(diào)查局自然資源航空物探遙感中心，北京 100083

2019—2020年期間，我國完成2020珠峰高程測量，與尼泊爾合作確定了基于國際高程參考系統(tǒng)的最新珠峰雪面正高，并于2021年12月8日共同宣布。中國和尼泊爾有各自國家的法定高程基準(zhǔn)，為實(shí)現(xiàn)共同宣布唯一珠峰高程，雙方?jīng)Q定：根據(jù)2015年國際大地測量協(xié)會(huì)關(guān)于國際高程參考系統(tǒng)定義和實(shí)現(xiàn)的決議[1]，在珠峰地區(qū)實(shí)現(xiàn)國際高程參考系統(tǒng)，作為珠峰高程的起算基準(zhǔn)。因此，建立基于國際高程參考系統(tǒng)的高精度珠峰區(qū)域大地水準(zhǔn)面模型、精密確定珠峰峰頂大地水準(zhǔn)面差距，是2020珠峰高程測量的關(guān)鍵問題。

大地水準(zhǔn)面確定需要區(qū)域內(nèi)密集分布的重力數(shù)據(jù)，珠峰及周邊地區(qū)是世界最高海拔區(qū)域，地形地貌復(fù)雜險(xiǎn)峻，嚴(yán)寒缺氧、人跡罕至，大部分區(qū)域無法開展地面重力測量，特別是珠峰鄰近區(qū)域，地形和重力場變化劇烈，幾乎是地面重力數(shù)據(jù)空白，嚴(yán)重制約了珠峰區(qū)域大地水準(zhǔn)面模型的精度。為解決重力數(shù)據(jù)稀少/空白問題，此次珠峰測量首次在珠峰區(qū)域開展高精度航空重力測量，并首次采集珠峰峰頂?shù)孛嬷亓?shù)據(jù)，這為構(gòu)建高精度珠峰區(qū)域大地水準(zhǔn)面模型奠定了數(shù)據(jù)基礎(chǔ)。

珠峰及周邊地區(qū)具備地面和航空重力數(shù)據(jù)可供使用，兩類數(shù)據(jù)具有不同的空間分布、信號(hào)頻譜和誤差特性，揭示不同波段(或頻譜)的地球重力場信息。航空重力數(shù)據(jù)主要貢獻(xiàn)中短波重力場信息，地面重力數(shù)據(jù)則主要反映短波信號(hào)。如何實(shí)現(xiàn)航空和地面重力數(shù)據(jù)優(yōu)化融合處理是高精度重力大地水準(zhǔn)面確定的難點(diǎn)問題[2-8]。

目前聯(lián)合航空和地面重力數(shù)據(jù)解算大地水準(zhǔn)面的方法主要分為3類：

第一類方法是先將航空重力數(shù)據(jù)從飛行高度處向下延拓至地面或大地水準(zhǔn)面，與地面重力數(shù)據(jù)融合得到格網(wǎng)平均重力異常，再利用Stokes積分或最小二乘配置計(jì)算大地水準(zhǔn)面[5,9-12]。該方法是兩步法，需要對(duì)航空重力數(shù)據(jù)進(jìn)行兩次處理，每一次處理都會(huì)引入邊界效應(yīng)，造成數(shù)據(jù)資源浪費(fèi)，并且數(shù)據(jù)聯(lián)合時(shí)未按照每種類型數(shù)據(jù)的頻譜和誤差特性配權(quán)。

第二類方法是最小二乘配置。該方法可一步聯(lián)合航空和地面重力數(shù)據(jù)解算大地水準(zhǔn)面[10,13]，其優(yōu)點(diǎn)是能夠同時(shí)容納不同類型和空間分布的重力數(shù)據(jù)求解擾動(dòng)位及其泛函。最小二乘配置涉及高階協(xié)方差矩陣組建和大型矩陣求逆計(jì)算，對(duì)于大范圍區(qū)域，計(jì)算工作量龐大，這一缺陷制約了該方法的普遍應(yīng)用。

第三類方法是譜組合法[14-15]。該方法從重力場參量之間的解析關(guān)系出發(fā)，將重力場參量按階作譜分解，依據(jù)最小二乘原理確定各類數(shù)據(jù)的譜權(quán)，按譜權(quán)由各類數(shù)據(jù)通過大地測量積分解析關(guān)系求解重力場參量，可一步聯(lián)合航空和地面重力數(shù)據(jù)解算大地水準(zhǔn)面。譜組合引入譜權(quán)，考慮了不同類型數(shù)據(jù)之間頻譜重疊的情況，積分求解避免了大型矩陣求逆計(jì)算，計(jì)算效率高。譜組合法的關(guān)鍵在于合理確定各類重力數(shù)據(jù)的譜權(quán)，即準(zhǔn)確估計(jì)各類重力數(shù)據(jù)的誤差階方差[6,16]。

本文利用譜組合方法聯(lián)合航空和地面重力數(shù)據(jù)建立珠峰區(qū)域重力似大地水準(zhǔn)面模型，并確定基于國際高程參考系統(tǒng)的珠峰峰頂大地水準(zhǔn)面差距，在珠峰地區(qū)實(shí)現(xiàn)國際高程參考系統(tǒng)。

1 數(shù)學(xué)模型

基于國際高程參考系統(tǒng)的珠峰峰頂大地水準(zhǔn)面差距N為

N=ζ+Δ+N0

(1)

式中，ζ為峰頂高程異常，由珠峰區(qū)域重力似大地水準(zhǔn)面模型插值計(jì)算得到；Δ為高程異常轉(zhuǎn)換為大地水準(zhǔn)面差距的改正項(xiàng)；N0為基于國際高程參考系統(tǒng)的大地水準(zhǔn)面差距零階項(xiàng)。

1.1 譜組合方法

珠峰地區(qū)重復(fù)覆蓋有地面和航空重力數(shù)據(jù)，可利用譜組合方法聯(lián)合地面和航空重力數(shù)據(jù)建立珠峰區(qū)域重力似大地水準(zhǔn)面模型。區(qū)域內(nèi)某點(diǎn)的高程異常可分解為地面和航空重力數(shù)據(jù)貢獻(xiàn)2個(gè)部分

ζGra=ζTer+ζAir

(2)

式中，ζGra為高程異常；ζTer為地面重力數(shù)據(jù)的高程異常貢獻(xiàn)；ζAir為航空重力數(shù)據(jù)的高程異常貢獻(xiàn)。

采用移去-計(jì)算-恢復(fù)技術(shù)，優(yōu)選地球重力場模型作為參考重力場，利用殘余地形模型(RTM)方法由高分辨率數(shù)字高程模型(DEM)計(jì)算高頻重力地形效應(yīng)[17]，地面和航空重力數(shù)據(jù)的高程異常貢獻(xiàn)可進(jìn)一步分解，則式(2)改寫為

(3)

(4)

式中，r為計(jì)算點(diǎn)地心距離；HP和H分別為計(jì)算點(diǎn)和流動(dòng)點(diǎn)的地形高度；δζHC表示高程異常的解析延拓改正項(xiàng)；KTer(ψ)為帶譜權(quán)的Stokes核函數(shù)，寫為

(5)

式中，ψ為計(jì)算點(diǎn)與流動(dòng)點(diǎn)的球面角距；NTer為地面重力異常格網(wǎng)對(duì)應(yīng)的最大展開階數(shù)；WTer(n)為第n階地面重力譜權(quán)。

(6)

式(4)中的高程異常解析延拓改正項(xiàng)δζHC為

(7)

式中，殘余重力異常的垂直梯度按照式(8)計(jì)算[18]

(8)

(9)

式中，KAir(ψ)為帶譜權(quán)的廣義Hotine核函數(shù)，寫為

(10)

式中，hM為平均飛行高度(大地高)；NAir為航空重力數(shù)據(jù)對(duì)應(yīng)的最大展開階數(shù)；WAir(n)為第n階航空重力譜權(quán)。

(11)

(12)

譜組合方法的關(guān)鍵在于合理確定各類重力數(shù)據(jù)的相對(duì)譜權(quán)，這需要各類重力數(shù)據(jù)的誤差信息。地面和航空重力數(shù)據(jù)互為獨(dú)立觀測，誤差階協(xié)方差可設(shè)為零，在確定譜權(quán)時(shí)只需要考慮各自的誤差階方差信息[14-15]。地面重力數(shù)據(jù)的譜權(quán)WTer(n)和航空重力數(shù)據(jù)的譜權(quán)WAir(n)按照式(13)—式(15)計(jì)算

(13)

(14)

(15)

式中，IEGM(n)、ITer(n)和IAir(n)分別為參考重力場模型、地面和航空重力數(shù)據(jù)的誤差階方差的倒數(shù)

(16)

由式(13)—式(16)可知，譜權(quán)確定問題轉(zhuǎn)化為各類重力數(shù)據(jù)的誤差階方差估計(jì)問題。本文采用一種直接從測量數(shù)據(jù)本身出發(fā)估計(jì)其誤差和誤差階方差的方法，可稱為數(shù)據(jù)驅(qū)動(dòng)法，其基本思想是：在低階部分(如NCut=200階以下)，以參考重力場模型作為基準(zhǔn)，將地面和航空重力數(shù)據(jù)與參考重力場模型進(jìn)行對(duì)比，利用球諧分析進(jìn)行譜分解估計(jì)地面和航空重力數(shù)據(jù)的長波誤差及低階誤差階方差；在中、高階部分(如NCut=200階以上)，使用舍一交叉驗(yàn)證方法(leave-out-one cross validation,LOOCV)，直接由地面和航空重力數(shù)據(jù)本身估計(jì)其中、短波誤差和中、高階誤差階方差。數(shù)據(jù)驅(qū)動(dòng)法的優(yōu)勢在于不需要重力數(shù)據(jù)的先驗(yàn)誤差信息，獲取的誤差階方差和譜權(quán)與實(shí)際重力數(shù)據(jù)達(dá)到最佳符合。關(guān)于數(shù)據(jù)驅(qū)動(dòng)的誤差階方差估計(jì)方法，詳細(xì)論述見文獻(xiàn)[6]。

1.2 高程異?！蟮厮疁?zhǔn)面差距轉(zhuǎn)換改正項(xiàng)

因珠峰地區(qū)屬于特大山區(qū)，地形和重力場起伏大，應(yīng)采用顧及地形質(zhì)量影響的嚴(yán)密公式(式(17))計(jì)算高程異?！蟮厮疁?zhǔn)面差距轉(zhuǎn)換改正項(xiàng)[19]

(17)

(18)

1.3 基于國際高程參考系統(tǒng)的大地水準(zhǔn)面差距零階項(xiàng)

(19)

式中，GM=3.986 004 415×1014m3s-2，為地心引力常數(shù)；GM0=3.986 005×1014m3s-2，為GRS80橢球的地心引力常數(shù)；W0=62 636 853.4 m2s-2，為國際高程參考系統(tǒng)定義的重力位值；U0=62 636 860.850 m2s-2，為GRS80橢球的正常重力位值；r為大地水準(zhǔn)面上相應(yīng)點(diǎn)的地心距離；γ為橢球面上相應(yīng)點(diǎn)的正常重力值。

2 數(shù)據(jù)處理與分析

2.1 珠峰區(qū)域和數(shù)據(jù)情況

建立珠峰地區(qū)的高精度大地水準(zhǔn)面模型是精確測定珠峰高程的關(guān)鍵所在，這需要密集、均勻分布的高精度重力數(shù)據(jù)。在珠峰及周邊區(qū)域(25°N—32°N、83°E—91°E)已有8022點(diǎn)歷史地面重力數(shù)據(jù)，但點(diǎn)位分布很不均勻，特別是在珠峰及鄰近區(qū)域地面重力點(diǎn)十分稀少。2020年5月1日—2020年6月11日期間，在珠峰地區(qū)成功開展航空重力測量，有效解決了該區(qū)域的重力數(shù)據(jù)稀少/空白問題。

航空重力測量使用航空地質(zhì)一號(hào)(空中國王350ER型飛機(jī))，從拉薩貢嘎機(jī)場起降，地面架設(shè)3個(gè)GNSS基準(zhǔn)站，同機(jī)搭載GT-2A型航空重力儀和DGA-01型國產(chǎn)航空重力儀，平均飛行速度441.7 km/h，平均飛行高度10 249 m(大地高)，共采集了東西向數(shù)據(jù)測線39條，南北向檢驗(yàn)測線9條，形成264個(gè)交叉點(diǎn)，測線間距5 km，在珠峰鄰近區(qū)域數(shù)據(jù)測線間距加密為2.5 km，測線總長度5 635.2 km，覆蓋面積1.27萬km2，數(shù)據(jù)采樣率2 Hz，共獲取83 803個(gè)數(shù)據(jù)點(diǎn)。選取同架次、同測線觀測數(shù)據(jù)進(jìn)行比較，GT-2A型和DGA-01型航空重力儀的內(nèi)符合精度達(dá)0.34 mGal，具有良好的一致性，最終采用GT-2A型航空重力儀的數(shù)據(jù)。航空重力數(shù)據(jù)處理采用GT-2A隨機(jī)軟件，經(jīng)100 s卡爾曼測線濾波處理和交叉點(diǎn)平差(調(diào)平)后，測線網(wǎng)交叉點(diǎn)差值RMS為1.1 mGal。利用GNSS差分動(dòng)態(tài)定位測定空中數(shù)據(jù)點(diǎn)的大地坐標(biāo)，測線觀測值為重力擾動(dòng)。利用CG-5型相對(duì)重力儀，采用往返重復(fù)觀測方式與國家重力基本點(diǎn)聯(lián)測，獲取飛機(jī)停機(jī)坪上航空重力儀正下方點(diǎn)位的絕對(duì)重力值，測量精度±7 μGal。珠峰地區(qū)航空重力測線網(wǎng)如圖1所示，其中黑色線表示航空重力測線，紅色點(diǎn)表示GNSS水準(zhǔn)點(diǎn)，紫色矩形框表示似大地水準(zhǔn)面模型范圍，底圖采用SRTM地形數(shù)據(jù)繪制。

注：藍(lán)色三角形代表珠峰。圖1 珠峰地區(qū)航空重力測線網(wǎng)與GNSS水準(zhǔn)并置點(diǎn)Fig.1 Airborne gravity survey lines and GNSS leveling points in the region of Mount Qomolangma

利用國產(chǎn)Z400型相對(duì)重力儀，在世界上首次獲取了珠峰峰頂重力觀測值。2005年珠峰高程測量時(shí)重力測量推進(jìn)到海拔高度7790 m，通過推算得到峰頂重力值[24-26]。此次獲取的峰頂實(shí)測重力值，有助于提升高程異?！蟮厮疁?zhǔn)面差距轉(zhuǎn)換改正項(xiàng)的計(jì)算精度[19]。此外，在珠峰鄰近區(qū)域拓展了4條新路線，結(jié)合水準(zhǔn)路線和登山路線，使用CG-6型相對(duì)重力儀共新測了210點(diǎn)地面重力數(shù)據(jù)，重力值精度優(yōu)于±39.5 μGal，結(jié)合該區(qū)域已有8022點(diǎn)地面重力數(shù)據(jù)，共有8232點(diǎn)地面重力數(shù)據(jù)可供使用。所有地面重力數(shù)據(jù)均統(tǒng)一到2000國家重力基本網(wǎng)，在進(jìn)行粗差剔除時(shí)，沒有發(fā)現(xiàn)粗差點(diǎn)。珠峰及周邊區(qū)域地面重力數(shù)據(jù)空間分布如圖2所示。

圖2 珠峰及周邊區(qū)域地面重力數(shù)據(jù)Fig.2 Terrestrial gravity data in the region of Mount Qomolangma

高分辨率DEM數(shù)據(jù)采用3″×3″ SRTM數(shù)據(jù)[27]，覆蓋范圍為25°N—32°N、83°E—91°E，基于SRTM數(shù)據(jù)的珠峰及周邊區(qū)域地形起伏如圖3所示。

注：藍(lán)色三角形代表珠峰。圖3 珠峰及周邊區(qū)域SRTM地形起伏Fig.3 SRTM topography in the region of Mount Qomolangma

設(shè)立由61個(gè)點(diǎn)組成的珠峰局部GNSS控制網(wǎng)，每點(diǎn)開展GNSS觀測1～2個(gè)時(shí)段，時(shí)段長度約8～14 h，各點(diǎn)大地高的平均精度為3.5 mm。這61個(gè)點(diǎn)是GNSS水準(zhǔn)并置點(diǎn)，結(jié)合水準(zhǔn)測量獲取各點(diǎn)的正常高，相對(duì)于起算點(diǎn)(國家一等水準(zhǔn)點(diǎn)“日喀則基巖點(diǎn)北”)的高程中誤差為0.1～1.84 mm，由此得到61個(gè)GNSS水準(zhǔn)點(diǎn)的高程異常，用于珠峰地區(qū)重力似大地水準(zhǔn)面模型精度檢核。61個(gè)GNSS水準(zhǔn)并置點(diǎn)空間分布如圖1所示。

2.2 計(jì)算方案

利用譜組合方法聯(lián)合航空重力擾動(dòng)和地面重力異常數(shù)據(jù)計(jì)算珠峰地區(qū)重力似大地水準(zhǔn)面，進(jìn)而確定基于國際高程參考系統(tǒng)的峰頂大地水準(zhǔn)面差距，計(jì)算方案如圖4所示。重力似大地水準(zhǔn)面模型覆蓋范圍為27.75°N—28.9°N、86.4°E—87.7°E，空間分辨率1′×1′。計(jì)算中采用Molodensky解析延拓法及基于地球重力場模型和RTM的移去-計(jì)算-恢復(fù)技術(shù)。對(duì)比分析地球重力場模型EGM2008[28]、EIGEN6-C4[29]和XGM2019[30]，優(yōu)選合適的參考重力場模型及其截?cái)嚯A數(shù)。利用3″×3″ SRTM數(shù)據(jù)計(jì)算地面重力異常高頻地形效應(yīng)。重力數(shù)據(jù)格網(wǎng)化采用反距離加權(quán)插值方法，分別得到1′×1′殘余航空重力擾動(dòng)格網(wǎng)和2′×2′殘余地面重力異常格網(wǎng)。以61個(gè)GNSS水準(zhǔn)并置點(diǎn)的高程異常為基準(zhǔn)值，對(duì)重力似大地水準(zhǔn)面模型進(jìn)行精度測試分析，優(yōu)選合適的Stokes和Hotine積分球冠半徑。利用珠峰地區(qū)重力似大地水準(zhǔn)面模型，顧及高差改正內(nèi)插計(jì)算得到峰頂高程異常，根據(jù)嚴(yán)密公式(式(17))施加高程異?！蟮厮疁?zhǔn)面差距轉(zhuǎn)換改正，結(jié)合基于國際高程參考系統(tǒng)的大地水準(zhǔn)面零階項(xiàng)(式(19))，最終確定基于國際高程參考系統(tǒng)的珠峰峰頂大地水準(zhǔn)面差距。

圖4 珠峰地區(qū)重力似大地水準(zhǔn)面和峰頂大地水準(zhǔn)面差距計(jì)算方案Fig. 4 Computation scheme for gravimetric quasigeoid in the region of Mount Qomolangma and geoid undulation at the summit

2.3 航空和地面重力數(shù)據(jù)的譜權(quán)確定

由式(5)和式(10)可知，在譜組合中引入各類重力數(shù)據(jù)的譜權(quán)，實(shí)質(zhì)上是按照階數(shù)對(duì)Stokes、Hotine積分核函數(shù)進(jìn)行修改，實(shí)現(xiàn)對(duì)航空和地面重力數(shù)據(jù)對(duì)應(yīng)高程異常貢獻(xiàn)的調(diào)節(jié)，以獲得最優(yōu)的似大地水準(zhǔn)面結(jié)果。根據(jù)式(13)—式(16)，譜權(quán)確定有賴于各類重力數(shù)據(jù)誤差階方差的準(zhǔn)確估計(jì)，采用基于數(shù)據(jù)驅(qū)動(dòng)的誤差和誤差階方差估計(jì)方法[6]，在低階部分(2～200階)，以參考重力場模型(2～200階)作為基準(zhǔn)，分別將航空重力擾動(dòng)、地面重力異常與參考重力場模型值作差得到殘余航空重力擾動(dòng)和殘余地面重力異常，并插值形成格網(wǎng)。在珠峰區(qū)域以外的格網(wǎng)點(diǎn)，采用參考重力場模型值(≥201階)填充，形成全球范圍的殘余空中重力擾動(dòng)格網(wǎng)和殘余地面重力異常格網(wǎng)，進(jìn)行球諧分析得到表征重力數(shù)據(jù)長波誤差的200階球諧系數(shù)，再分別計(jì)算得到航空和地面重力數(shù)據(jù)在2～200階的誤差階方差。在中、高階部分(≥201階)，采用LOOCV方法，欲估計(jì)某點(diǎn)的重力誤差，先舍棄該待求點(diǎn)，利用該點(diǎn)周圍多個(gè)點(diǎn)的重力觀測值內(nèi)插得到待求點(diǎn)的重力內(nèi)插值，將內(nèi)插值與觀測值之差作為重力數(shù)據(jù)的中、短波誤差估值。分別獲得所有航空和地面重力數(shù)據(jù)點(diǎn)的誤差估值后，插值形成5′×5′的誤差格網(wǎng)，由誤差格網(wǎng)推算誤差協(xié)方差函數(shù)，進(jìn)而分別計(jì)算出航空和地面重力數(shù)據(jù)在201～2160階的誤差階方差。至此，筆者獲取了航空和地面重力數(shù)據(jù)在2～2160階的誤差階方差，即可按照式(13)—式(16)確定相應(yīng)的譜權(quán)。

圖5為珠峰區(qū)域航空和地面重力數(shù)據(jù)的譜權(quán)。由圖5可以看出，航空重力數(shù)據(jù)的貢獻(xiàn)譜段集中在140～700階之間，航空重力譜權(quán)大約從120階開始快速上升，在200階左右達(dá)到頂點(diǎn)，最大譜權(quán)達(dá)0.8，隨后逐漸下降。200階以上航空和地面重力譜權(quán)之和為1，兩者是此降彼升的關(guān)系，航空重力數(shù)據(jù)的貢獻(xiàn)越小，地面重力數(shù)據(jù)的貢獻(xiàn)則越大，1000階以上航空重力數(shù)據(jù)的貢獻(xiàn)基本可忽略不計(jì)，譜組合主要是來自地面重力數(shù)據(jù)的貢獻(xiàn)。

圖5 航空與地面重力數(shù)據(jù)的譜權(quán)Fig.5 Spectral weights of airborne and terrestrial gravity data

2.4 參考重力場模型和球冠積分半徑的測試分析

由于航空和地面重力數(shù)據(jù)覆蓋范圍有限，需利用參考重力場模型作為區(qū)域以外重力場信息的貢獻(xiàn)，具體實(shí)現(xiàn)采用移去-計(jì)算-恢復(fù)技術(shù)。選取不同的地球重力場模型作為參考場、同一個(gè)模型選取不同的截?cái)嚯A數(shù)，會(huì)得出不同的重力似大地水準(zhǔn)面結(jié)果，有必要對(duì)不同參考重力場模型和模型截?cái)嚯A數(shù)進(jìn)行測試分析。選擇地球重力場模型EGM2008、EIGEN6-C4和XGM2019作為參考重力場的備選，3個(gè)模型均展開至2190階，利用珠峰地區(qū)61點(diǎn)GNSS水準(zhǔn)高程異常數(shù)據(jù)對(duì)模型進(jìn)行精度評(píng)價(jià)。表1列出了由3個(gè)重力場模型計(jì)算的高程異常與61點(diǎn)GNSS水準(zhǔn)高程異常之間的差值統(tǒng)計(jì)信息。由表1可以看出，在珠峰地區(qū)，EIGEN-6C4模型精度最高，XGM2019模型次之，EGM2008模型最差。EIGEN-6C4模型和XGM2019模型分別由不同機(jī)構(gòu)研制，采用的數(shù)據(jù)和方法不盡相同，兩者最大的區(qū)別是XGM2019模型利用DEM正演建模技術(shù)推算高頻重力異常(719階以上)[30]，這可能是XGM2019模型精度稍低于EIGEN-6C4模型的原因。因此，選定EIGEN-6C4模型作為重力似大地水準(zhǔn)面計(jì)算的參考重力場模型。

表1 地球重力場模型高程異常與GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)Tab.1 Statistics of the differences between Earth gravity model based and GNSS leveling measured height anomalies m

再測試參考重力場模型的截?cái)嚯A數(shù)，表2為EIGEN-6C4模型截取至不同階數(shù)時(shí)、聯(lián)合航空和地面重力數(shù)據(jù)解算的重力似大地水準(zhǔn)面與61點(diǎn)GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)信息。結(jié)果表明，參考重力場模型EIGEN-6C4截取至1080階時(shí)的結(jié)果是最優(yōu)的，精度達(dá)3.8 cm。若截取階數(shù)偏高，會(huì)削弱區(qū)域內(nèi)實(shí)測重力數(shù)據(jù)的貢獻(xiàn)，導(dǎo)致重力似大地水準(zhǔn)面難以達(dá)到最佳精度。若截取階數(shù)偏低，區(qū)域外重力場截?cái)嗾`差增大，導(dǎo)致結(jié)果精度偏低。

表2 重力似大地水準(zhǔn)面與GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)Tab.2 Statistics of the differences between gravimetric geoid and GNSS leveling measured height anomalies m

根據(jù)式(4)和式(9)可知，優(yōu)選合適的Stokes和Hotine積分球冠半徑，對(duì)獲取最佳的重力似大地水準(zhǔn)面模型也是至關(guān)重要的。表3為選取不同球冠積分半徑時(shí)、聯(lián)合航空和地面重力數(shù)據(jù)解算的重力似大地水準(zhǔn)面與61點(diǎn)GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)信息。結(jié)果顯示，當(dāng)球冠積分半徑取1°，重力似大地水準(zhǔn)面模型精度最高。

表3 重力似大地水準(zhǔn)面與GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)Tab.3 Statistics of the differences between gravimetric geoid and GNSS leveling measured height anomalies m

2.5 重力似大地水準(zhǔn)面模型精度評(píng)估

根據(jù)以上重力似大地水準(zhǔn)面模型與GNSS水準(zhǔn)高程異常的比較分析結(jié)果，選定參考重力場模型為EIGEN-6C4，并截取至1080階，Stokes和Hotine積分球冠半徑選定為1°。利用譜組合方法聯(lián)合航空和地面重力數(shù)據(jù)解算得到重力似大地水準(zhǔn)面模型(圖6)，經(jīng)61點(diǎn)GNSS水準(zhǔn)高程異常檢核，模型精度達(dá)到3.8 cm。圖7為該模型與61點(diǎn)GNSS水準(zhǔn)高程異常的差值(已扣除平均值)。

注：藍(lán)色三角形代表珠峰。圖6 聯(lián)合航空和地面重力數(shù)據(jù)構(gòu)建的重力似大地水準(zhǔn)面模型Fig.6 Gravimetric quasigeoid model computed from the combination of airborne and terrestrial gravity data

注：藍(lán)色三角形代表珠峰。圖7 重力似大地水準(zhǔn)面模型與GNSS水準(zhǔn)高程異常的差值Fig.7 Differences between gravimetric quasigeoid model and GNSS leveling measured height anomalies

將譜組合方法稱為方案1，本文還采用了另外兩種計(jì)算方案進(jìn)行對(duì)比分析。方案2是利用EIGEN-6C4模型和地面重力異常數(shù)據(jù)，不使用航空重力數(shù)據(jù)，采用基于移去-計(jì)算-恢復(fù)的Molodensky方法計(jì)算重力似大地水準(zhǔn)面[18,26]。在Molodensky級(jí)數(shù)解公式中以地形改正項(xiàng)代替Molodensky級(jí)數(shù)解重力一階改正項(xiàng)。利用移去-恢復(fù)技術(shù)由均衡重力異常計(jì)算得到格網(wǎng)平均重力異常，均衡改正計(jì)算選擇愛黎-海斯卡寧模型，均衡深度選34 km。

方案3是利用基于快速傅里葉變換(FFT)的泊松積分法將航空重力異常數(shù)據(jù)向下延拓至地面，與地面重力異常數(shù)據(jù)融合形成格網(wǎng)平均重力異常后，結(jié)合EIGEN-6C4模型，采用基于移去-計(jì)算-恢復(fù)的Molodensky方法計(jì)算重力似大地水準(zhǔn)面。

利用珠峰地區(qū)61點(diǎn)GNSS水準(zhǔn)高程異常數(shù)據(jù)對(duì)重力似大地水準(zhǔn)面模型進(jìn)行精度評(píng)估，基于3種方案的重力似大地水準(zhǔn)面模型相對(duì)于GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)信息見表4。結(jié)果表明，與單獨(dú)利用地面重力數(shù)據(jù)相比，加入航空重力數(shù)據(jù)能夠顯著提高重力似大地水準(zhǔn)面模型的精度，方案1和方案3的提升幅度分別為51.3%和38.5%。方案1的結(jié)果要優(yōu)于方案3，這主要?dú)w功于譜組合方法不需要對(duì)航空重力數(shù)據(jù)進(jìn)行顯式的向下延拓處理，而是一步聯(lián)合航空和地面重力數(shù)據(jù)計(jì)算似大地水準(zhǔn)面。方案3則是兩步法，需要將航空重力數(shù)據(jù)向下延拓、再與地面重力數(shù)據(jù)融合后進(jìn)行積分計(jì)算，航空重力數(shù)據(jù)的利用率顯然低于方案1。

表4中的平均差值代表珠峰區(qū)域重力似大地水準(zhǔn)面(不含零階項(xiàng))與我國1985國家高程基準(zhǔn)的垂直偏差。根據(jù)筆者和相關(guān)學(xué)者的研究結(jié)果，重力似大地水準(zhǔn)面(不含零階項(xiàng))與1985國家高程基準(zhǔn)之間的垂直偏差介于28～35 cm[31-34]。加入航空重力數(shù)據(jù)后，平均差值分別從16.9 cm增大為27 cm(方案1)和31.4 cm(方案3)，體現(xiàn)了航空重力數(shù)據(jù)對(duì)重力似大地水準(zhǔn)面的貢獻(xiàn)，使之更為接近真實(shí)值。

表4 重力似大地水準(zhǔn)面模型與GNSS水準(zhǔn)高程異常的差值統(tǒng)計(jì)Tab.4 Statistics of the differences between gravimetric geoid and GNSS leveling measured height anomalies m

2.6 基于國際高程參考系統(tǒng)的峰頂大地水準(zhǔn)面差距確定

2.6.1 顧及高差改正的峰頂高程異常插值計(jì)算

珠峰峰頂高程異?？捎芍亓λ拼蟮厮疁?zhǔn)面模型經(jīng)樣條插值或雙線性插值得到。珠峰地區(qū)地形起伏劇烈，由模型內(nèi)插計(jì)算峰頂高程異常不同于其他地區(qū)，應(yīng)考慮高程變化對(duì)內(nèi)插高程異常的影響，并進(jìn)行相應(yīng)的高差改正。

似大地水準(zhǔn)面模型的格網(wǎng)計(jì)算點(diǎn)是由計(jì)算者選定、用于表示地表的特定DEM確定。設(shè)在似大地水準(zhǔn)面模型計(jì)算時(shí)用1′×1′ DEM代表地表格網(wǎng)計(jì)算點(diǎn)，若某地面點(diǎn)P的高程為h，由1′×1′ DEM內(nèi)插得到該點(diǎn)在DEM面上對(duì)應(yīng)點(diǎn)P0的高程為h0，由對(duì)應(yīng)的1′×1′似大地水準(zhǔn)面模型內(nèi)插得到點(diǎn)P在DEM面上對(duì)應(yīng)點(diǎn)P0的高程異常為ζ0，則P點(diǎn)的實(shí)際高程異常ζ為

(20)

計(jì)算式(20)時(shí)，將點(diǎn)P0至點(diǎn)P的高差分為N=100個(gè)等分高差段Δh，δgi和γi為每個(gè)高差段中點(diǎn)Pi的重力擾動(dòng)和正常重力，δgi由式(21)計(jì)算

δgi=gP+Fi-γi

(21)

式中，gP為點(diǎn)P的重力值；Fi為點(diǎn)P沿垂線到點(diǎn)Pi的空間改正。

珠峰地區(qū)重力似大地水準(zhǔn)面計(jì)算時(shí)采用1′×1′ DEM作為地表計(jì)算點(diǎn)格網(wǎng)。利用峰頂點(diǎn)實(shí)際正高及其在1′×1′ DEM地形面上投影點(diǎn)的正高，結(jié)合峰頂實(shí)測重力值，根據(jù)式(21)計(jì)算得到峰頂內(nèi)插高程異常的高差改正Δζ=-0.05 m。

為檢核高差改正計(jì)算的正確性，根據(jù)珠峰峰頂精確三維坐標(biāo)，利用譜組合方法聯(lián)合航空和地面重力數(shù)據(jù)，直接單點(diǎn)積分計(jì)算得到峰頂高程異常。對(duì)比峰頂高程異常的單點(diǎn)積分結(jié)果與顧及高差改正的模型內(nèi)插結(jié)果，兩者僅相差0.9 cm，這證明了高差改正的必要性和正確性。

2.6.2 基于國際高程參考系統(tǒng)的峰頂大地水準(zhǔn)面差距計(jì)算

利用方案1構(gòu)建的重力似大地水準(zhǔn)面模型，經(jīng)樣條插值和高差改正得到峰頂高程異常ζ。根據(jù)式(17)和式(18)，利用珠峰峰頂實(shí)測重力值和3″×3″ SRTM數(shù)據(jù)計(jì)算得到峰頂高程異常轉(zhuǎn)換為大地水準(zhǔn)面差距的改正項(xiàng)Δ=-1.302 m。根據(jù)式(19)，采用國際高程參考系統(tǒng)定義的重力位值，選定GRS80參考橢球，計(jì)算得到大地水準(zhǔn)面差距零階項(xiàng)N0=-0.177 m，最終按照式(1)計(jì)算得到基于國際高程參考系統(tǒng)的峰頂大地水準(zhǔn)面差距。采用基于方案3的重力似大地水準(zhǔn)面模型計(jì)算得到峰頂大地水準(zhǔn)面差距，與方案1所得峰頂大地水準(zhǔn)面差距結(jié)果僅相差3.6 cm。

3 結(jié) 論

2020珠峰高程測量，成功獲取珠峰峰頂及周邊區(qū)域1.27萬km2的航空重力數(shù)據(jù)，測線網(wǎng)平均精度達(dá)1.1 mGal，有效填補(bǔ)了珠峰地區(qū)重力數(shù)據(jù)空白。聯(lián)合航空和地面重力數(shù)據(jù)，建立了珠峰區(qū)域1′×1′高精度重力似大地水準(zhǔn)面模型，基于國際高程參考系統(tǒng)定義的重力位值W0和GRS80參考橢球，確定了國際高程參考系統(tǒng)中的珠峰峰頂大地水準(zhǔn)面差距，為中國和尼泊爾合作確定基于國際高程參考系統(tǒng)的珠峰正高提供了高精度基準(zhǔn)。

基于譜組合理論和數(shù)據(jù)驅(qū)動(dòng)的譜權(quán)確定方法，聯(lián)合航空和地面重力數(shù)據(jù)解算的重力似大地水準(zhǔn)面模型，經(jīng)61個(gè)GNSS水準(zhǔn)并置點(diǎn)的高程異常檢核，模型精度達(dá)3.8 cm。加入航空重力數(shù)據(jù)能有效填充地面重力數(shù)據(jù)空白，譜組合中航空重力數(shù)據(jù)的貢獻(xiàn)譜段集中在140～700階。與單獨(dú)利用地面重力數(shù)據(jù)相比，航空重力數(shù)據(jù)使得重力似大地水準(zhǔn)面模型精度提升了51.3%(3.8 cm VS 7.8 cm)。

珠峰地區(qū)地形極端崎嶇、重力場變化劇烈，利用重力似大地水準(zhǔn)面模型計(jì)算珠峰峰頂大地水準(zhǔn)面差距時(shí)，有必要采用顧及高差改正的內(nèi)插方法和顧及地形質(zhì)量影響的高程異?！蟮厮疁?zhǔn)面差距轉(zhuǎn)換改正公式，結(jié)合此次采集的峰頂?shù)孛嬷亓?shù)據(jù)，以確保峰頂大地水準(zhǔn)面差距的計(jì)算精度。

國際大地測量協(xié)會(huì)2015年發(fā)布了國際高程參考系統(tǒng)的定義，并于2019年提出了建立國際高程參考框架(IHRF)的遠(yuǎn)景目標(biāo)[35]。2020珠峰高程測量成功在珠峰地區(qū)實(shí)現(xiàn)了國際高程參考系統(tǒng)，可為建立國際高程參考框架提供重要參考?，F(xiàn)階段，我國及全球范圍仍存在較大范圍地面重力數(shù)據(jù)空白區(qū)域，已有重力數(shù)據(jù)的精度和現(xiàn)勢性也亟須提升、更新。開展大規(guī)模航空重力測量，獲取均勻分布的高精度重力數(shù)據(jù)，進(jìn)一步發(fā)展多種類型重力數(shù)據(jù)融合處理的理論和方法，對(duì)于厘米級(jí)精度重力大地水準(zhǔn)面模型構(gòu)建和全球高程基準(zhǔn)統(tǒng)一具有重要價(jià)值。

致謝：特此向?yàn)?020珠峰高程測量作出貢獻(xiàn)的專家學(xué)者、測繪隊(duì)員和登山隊(duì)員表示衷心的感謝。