褚永海，阮文飛，柯寶貴，汪海洪

剩余地形模型及多面函數(shù)改進(jìn)基準(zhǔn)面的研究

褚永海1,2，阮文飛1，柯寶貴3，汪海洪1,2

（1. 武漢大學(xué) 測(cè)繪學(xué)院，武漢 430079；2. 武漢大學(xué) 地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430079；3. 中國(guó)測(cè)繪科學(xué)研究院 大地測(cè)量與地球動(dòng)力學(xué)研究所，北京 100036)

為了進(jìn)一步提高高程基準(zhǔn)面的精度，研究剩余地形模型（RTM）在近海區(qū)域高程基準(zhǔn)面中的影響：利用陸海統(tǒng)一巖石等效地形（RET）構(gòu)建近海區(qū)域RTM和殘余高程異常；并聯(lián)合多面函數(shù)擬合殘余高程異常改正面，實(shí)現(xiàn)近海區(qū)域高程基準(zhǔn)面模型精度的提高。實(shí)驗(yàn)結(jié)果表明：在近海區(qū)域，利用RTM高程異?？商岣咝录夹g(shù)改進(jìn)后的歐洲地球重力聯(lián)合模型（EIGEN-6C4）的精度約為3 mm；多面函數(shù)擬合后，可以將EIGEN-6C4模型的精度由9.0 cm提高到3.4 cm。

高程異常；剩余地形模型；巖石等效地形；多面函數(shù)；高程基準(zhǔn)

0 引言

在區(qū)域高程基準(zhǔn)面確定中，地形數(shù)據(jù)主要用于高頻信息的獲取。處理地形數(shù)據(jù)通常有4種策略：①考慮整個(gè)地形的影響，包括布格改正和地形改正；②只考慮局部地形影響；③對(duì)地形進(jìn)行均衡歸算，有艾里-海斯卡涅以及普拉特-海福德模型[1-2]；④使用剩余地形模型（residual terrain model, RTM）[3]。

RTM的主要思想是用1個(gè)真實(shí)地形減去1個(gè)參考地形得到剩余地形，然后再計(jì)算剩余地形的影響[3]。本質(zhì)上來(lái)說(shuō)，RTM將參考面以上的地形質(zhì)量移去，而參考面以下質(zhì)量虧損部分將被填充。在物理大地測(cè)量學(xué)中的短尺度重力建模[4]或高頻重力正演[5]中，經(jīng)常使用RTM。此外，還利用RTM來(lái)改善全球引力場(chǎng)模型（global geopotential model, GGM）計(jì)算的高程基準(zhǔn)面。例如文獻(xiàn)[6-7]聯(lián)合地球引力模型2008（Earth gravitational model 2008, EGM2008）和數(shù)字地面模型2006.0（digital terrain model 2006.0, DTM2006.0），來(lái)改善重力數(shù)據(jù)稀少的山地區(qū)域似大地水準(zhǔn)面，精度可以提高2 cm，這主要是考慮山區(qū)重力數(shù)據(jù)稀少，而且直接使用的高階GGM也沒(méi)有完全包含地球重力場(chǎng)的高頻信息，所產(chǎn)生的截?cái)嗾`差在山區(qū)對(duì)高程異常的影響可以達(dá)到10 cm量級(jí)。利用RTM可以對(duì)截?cái)嗾`差進(jìn)行估計(jì)并有效改善GGM[8]。這種思想也被國(guó)內(nèi)學(xué)者用于區(qū)域基準(zhǔn)面的改進(jìn)。例如：文獻(xiàn)[9]在華南地區(qū)選擇了1個(gè)帶狀和1個(gè)面狀區(qū)域，利用RTM模型來(lái)改進(jìn)EGM2008；文獻(xiàn)[10]在廣東省及周邊區(qū)域?qū)TM高程異常的精度進(jìn)行了數(shù)值計(jì)算與分析；文獻(xiàn)[11]將RTM用于改進(jìn)長(zhǎng)江南京段的高程控制。

針對(duì)近海區(qū)域區(qū)域高程基準(zhǔn)面確定中RTM的貢獻(xiàn)及影響，本文利用30″和15″分辨率陸地高程和海洋深度模型（SRTM30_PLUS和SRTM15_PLUS）[12]、巖石等效球諧地形模型2012(rock-equivalent terrain 2012, RET2012)[13]、歐洲地球重力場(chǎng)聯(lián)合模型（European improved gravity model of the earth by new techniques, EIGEN-6C4）[14]、全球定位系統(tǒng)（global positioning system，GPS）及水準(zhǔn)測(cè)量數(shù)據(jù)，分析RTM對(duì)EIGEN-6C4模型的改善情況。同時(shí)，采用多面函數(shù)法對(duì)殘余高程異常進(jìn)行擬合，再聯(lián)合GGM高程異常、RTM高程異常，恢復(fù)得到近海區(qū)域高程基準(zhǔn)面，并用獨(dú)立GPS/水準(zhǔn)測(cè)量數(shù)據(jù)進(jìn)行檢核。

1 原理方法

1.1 剩余地形建模

1.2 高程異常

式中：G為萬(wàn)有引力常數(shù)；ρ為長(zhǎng)方體的密度；為長(zhǎng)方體頂點(diǎn)坐標(biāo)；l為距離，用坐標(biāo)表示成。長(zhǎng)方體頂點(diǎn)E的坐標(biāo)可以用表示；頂點(diǎn)K的坐標(biāo)用表示，其他點(diǎn)都可以用點(diǎn)E和點(diǎn)K的坐標(biāo)表示，具體位置關(guān)系參見(jiàn)圖1。

根據(jù)布隆斯（Bruns）公式，每個(gè)長(zhǎng)方體地形引力位對(duì)高程異常的貢獻(xiàn)[15]表示成

1.3 殘余高程異常及多面函數(shù)擬合

在重力場(chǎng)參量的擬合計(jì)算中，廣泛使用多面函數(shù)，它的基本思想是，任何1個(gè)規(guī)則或不規(guī)則的連續(xù)曲面均可以由若干簡(jiǎn)單面來(lái)疊加逼近。具體做法是：在每個(gè)數(shù)據(jù)點(diǎn)上建立1個(gè)曲面，然后將各個(gè)旋轉(zhuǎn)曲面按一定比例疊加成一張整體的連續(xù)曲面，使之嚴(yán)格地通過(guò)各個(gè)數(shù)據(jù)點(diǎn)。其函數(shù)模型[17]可以表示成

2 實(shí)驗(yàn)與結(jié)果分析

2.1 研究區(qū)域和使用數(shù)據(jù)

當(dāng)前高精度、高分辨的全球陸地地形數(shù)據(jù)主要有3″分辨率的SRTM地形數(shù)據(jù)，海域數(shù)據(jù)除了少數(shù)船測(cè)數(shù)據(jù)之外，主要是依據(jù)衛(wèi)星測(cè)高反演海域重力異常而獲得海底地形[18]。將陸地SRTM和海洋深度模型融合后，國(guó)外先后發(fā)布了包含陸地地形及海底地形，分辨率為30″的SRTM30_PLUS模型及分辨率為15″的SRTM15_PLUS模型。其中SRTM15_PLUS（如圖2所示）在研究范圍（115°E~125°E，31°N~41°N）的統(tǒng)計(jì)信息如表1所示，最低點(diǎn)的高度-128.2 m，位于（122°13′30″E，38°50′15″N），最高點(diǎn)的高度2770.0 m，位于（115°02′30″E，39°56′30″N）。利用海域部分基于巖石等效地形的概念[8,19]，利用海水進(jìn)行壓縮，構(gòu)建密度均勻（巖石密度2670 kg·m-3）的巖石等效地形，相應(yīng)的數(shù)值結(jié)果列于表1。

表1 地形模型高度值統(tǒng)計(jì) 單位：m

圖2 SRTM15_PLUS陸地地形/海洋深度模型

基于SRTM30_PLUS，利用球諧分析確定了全球地形及引力位的球諧模型2012和2014[13,20]，其中包含巖石等效球諧地形模型（RET2012和RET2014）。RET2012模型的階數(shù)為2160階，相當(dāng)于5′分辨率，與EGM2008、EIGEN_6C4階次相同。RET2012按15″分辨率計(jì)算的球諧地形，在研究區(qū)域內(nèi)最大值高度約為2234.1 m，最小值為-67.6 m，平均高度為57.2 m（見(jiàn)表1）。

選取SRTM15_PLUS巖石等效地形作為真實(shí)地形模型，參考地形使用RET2012巖石等效球諧地形，2者相減得到RTM模型。RTM模型在研究區(qū)域最大值高度約為851.4 m，最小值高度約為-684.2 m，平均值高度約為-0.3 m（見(jiàn)表1）。

此外，將GPS/水準(zhǔn)點(diǎn)共計(jì)152個(gè)分為2組：第1組作為多面函數(shù)擬合計(jì)算點(diǎn)（如圖3所示），共計(jì)88個(gè)；第2組為檢核點(diǎn)（如圖4所示），共計(jì)64個(gè)點(diǎn)，用來(lái)對(duì)擬合改正面進(jìn)行獨(dú)立檢核。

圖3 擬合GPS水準(zhǔn)點(diǎn)分布

圖4 檢核GPS水準(zhǔn)點(diǎn)分布

2.2 積分半徑對(duì)RTM高程異常的影響

為了解長(zhǎng)方體模型積分半徑的選取情況，選取4個(gè)點(diǎn)、、和進(jìn)行測(cè)試（如圖2所示）。4個(gè)點(diǎn)的海拔高程分別約為382.0、103.6、38.3和78.6 m。其中點(diǎn)位于山區(qū)，點(diǎn)位于二山中間，和點(diǎn)均在近海附近。此外，4個(gè)點(diǎn)的剩余地形分別約為-77.3、-8.4、25.3和-6.8 m。從剩余地形來(lái)看，除了點(diǎn)高于參考面之外，其他3個(gè)點(diǎn)都低于參考面。不分內(nèi)區(qū)和外區(qū)，數(shù)據(jù)分辨率取7.5″，積分半徑取0.1、0.2、0.5、1、2、3、4、5、10、15、…、300 km后分別進(jìn)行計(jì)算。

為了進(jìn)一步說(shuō)明積分半徑的影響，以積分半徑40、60、…、300 km計(jì)算了152個(gè)GPS水準(zhǔn)點(diǎn)上的RTM高程異常，其統(tǒng)計(jì)結(jié)果列于表2。從極值可以看出，RTM對(duì)高程異常的影響為-4.6~1.7 cm。由于構(gòu)建剩余地形的初衷之一是剩余量的平均值趨近于零，因此積分半徑取240 km比較合適。

圖5 RTM高程異常與積分半徑的關(guān)系

表2 不同積分半徑下RTM高程異常統(tǒng)計(jì)

2.3 RTM高程異常及殘余高程異常擬合

圖6 RTM高程異常分布

圖7 多面函數(shù)法擬合殘余高程異常

圖8 擬合改正后的高程基準(zhǔn)面

2.4 精度檢核

重力模型高程異常顧及剩余地形影響，并利用GPS/水準(zhǔn)擬合改正得到的區(qū)域高程基準(zhǔn)面，除了海域無(wú)GPS/水準(zhǔn)控制存在局部扭曲之外，其他區(qū)域的精度可以用另一組64個(gè)GPS/水準(zhǔn)數(shù)據(jù)（如圖4所示）進(jìn)行外部檢核。為了比較RTM高程異常改正、殘余地形擬合改正對(duì)高程基準(zhǔn)面的貢獻(xiàn)，進(jìn)行了以下精度估計(jì)。

在顧及RTM高程異常的基礎(chǔ)上，如果再顧及殘余高程異常改正，即利用GPS/水準(zhǔn)數(shù)據(jù)檢核擬合后的基準(zhǔn)面（見(jiàn)圖8），高程異常差值統(tǒng)計(jì)列于表3中的最后3行。其中，由于擬合殘余高程異常改正面使用的是88個(gè)擬合點(diǎn)，多面函數(shù)法擬合的曲面通過(guò)擬合點(diǎn)，因此，這組GPS/水準(zhǔn)點(diǎn)檢核統(tǒng)計(jì)數(shù)值均為零，即內(nèi)符合檢核無(wú)誤差。而使用未參與擬合的64個(gè)GPS/水準(zhǔn)數(shù)據(jù)檢核，最大差值為0.0666 m，最小值為-0.1051 m，平均差異為-0.0015 m，標(biāo)準(zhǔn)差為0.0341 m，均方差為0.0339 m。說(shuō)明了外符合檢核精度約為3.4 cm。如果使用全部152個(gè)點(diǎn)檢核，精度約為2.2 cm。

表3 EIGEN-6C4/RTM高程異常與GPS/水準(zhǔn)高程異常比較

續(xù)表3

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

剩余地形模型在高海拔地形起伏地區(qū)，特別是在無(wú)實(shí)測(cè)重力地區(qū)，能夠有效地改善全球重力場(chǎng)模型短波分量因模型截?cái)嗟挠绊慬7,9-11]；在近海地區(qū)，剩余地形模型對(duì)高程異常的改進(jìn)可提高的精度約為6.5 mm[8]。為了分析剩余地形模型在中國(guó)近海的應(yīng)用情況，選取青島及周邊區(qū)域進(jìn)行計(jì)算分析，包括數(shù)值積分半徑的選取、剩余地形模型的貢獻(xiàn)、擬合改正后的精度情況等，數(shù)值結(jié)果表明：

剩余地形模型對(duì)重力場(chǎng)模型短波分量的改善，不僅可以用于高程異常的改善，也可以用于重力、垂線偏差的改善[8]，這對(duì)于無(wú)實(shí)測(cè)重力數(shù)據(jù)地區(qū)重力場(chǎng)模型的構(gòu)建，具有重要的作用。近幾十年來(lái)，RTM技術(shù)已經(jīng)被廣泛用于大地測(cè)量領(lǐng)域，但以目前的方式使用時(shí)，還有一些問(wèn)題需要解決：例如，計(jì)算點(diǎn)位于較低的山谷時(shí)，為了保證引力位的調(diào)和性，需要進(jìn)行調(diào)和改正；當(dāng)數(shù)值積分半徑較大時(shí)，需要顧及地球曲率的影響。另外，在殘余地形擬合改正面中，海域因無(wú)GPS/水準(zhǔn)約束出現(xiàn)扭曲，這可以將海面地形作為“虛擬”GPS/水準(zhǔn)點(diǎn)，以保證將陸地高程傳遞到海域范圍。總之，如果在RTM建模中進(jìn)行調(diào)和改正和地球曲率改正，在擬合曲面時(shí)對(duì)海洋進(jìn)行約束，可以進(jìn)一步提高程基準(zhǔn)面的精度，并且實(shí)現(xiàn)陸海高程基準(zhǔn)的統(tǒng)一。

[1] HEISKANEN W A, MORITZ H. Physical Geodesy[M]. New York: Freeman W. H, 1967: 364.

[2] HOFMANN-WELLENHOF B, MORITZ H. Physical Geodesy[M]. Austria: Springer-Verlag Wien, 2006: 403.

[3] FORSBERG R. A study of terrain corrections, density anomalies and geophysical inversion methods in gravity field modelling, 355[EB/OL].[2019-12-28].https://apps.dtic.mil/dtic/tr/fulltext/u2/a150788.pdf.

[4] HIRT C, BUCHA B, YANG M, et al. A numerical study of residual terrain modelling (RTM) techniques and the harmonic correction using ultra-high-degree spectral gravity modelling[J]. Journal of Geodesy, 2019, 93(9): 1469-1486.

[5] REXER M, HIRT C, BUCHA B, et al. Solution to the spectral filter problem of residual terrain modelling (RTM)[J]. Journal of Geodesy, 2018, 92(6): 675-690.

[6] PAVLIS N K, HOLMES S A, KENYON S C, et al. The development and evaluation of the Earth gravitational model 2008 (EGM2008)[J]. Journal of Geophysical Research, 2012, 117(B04): 1-38.

[7] HIRT C, FEATHERSTONE W E, MARTI U. Combining EGM2008 and SRTM/DTM2006.0 residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data[J]. Journal of Geodesy, 2010, 88(9): 557-567.

[8] HIRT C. RTM gravity forward-modeling using topography/bathymetry data to improve high-degree global geopotential models in the coastal zone[J]. Marine Geodesy, 2013, 36(2): 1-20.

[9] 張興福, 劉成. 綜合EGM2008模型和SRTM/DTM2006.0剩余地形模型的GPS高程轉(zhuǎn)換方法[J].測(cè)繪學(xué)報(bào), 2012, 41(1): 25-32.

[10] 張永毅, 張興福, 周波陽(yáng), 等.剩余地形模型高程異常計(jì)算的積分法及精度分析[J].大地測(cè)量與地球動(dòng)力學(xué), 2016, 36(9): 770-774.

[11] 翟長(zhǎng)治, 姚宜斌, 岳順.基于EGM2008和剩余地形模型的區(qū)域似大地水準(zhǔn)面精化方法[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2015, 35(6): 941-944.

[12] BECKER J J, SANDWELL D T, SMITH W H F, et al. Global bathymetry and elevation data at 30 arc seconds resolution: SRTM30_PLUS[J]. Marine Geodesy, 2009, 32(4): 355-371.

[13] HIRT C, KUHN M. Evaluation of high-degree series expansions of the topographic potential to higher-order powers[J]. Journal of Geophysical Research: Solid Earth, 2012, 17(B12): 1-12.

[14] F?RSTE C, BRUINSMA S L, ABRIKOSOV O, et al. EIGEN-6C4: the latest combined global gravity field model including GOCE data up to degree and order 2190 of GFZ Potsdam and GRGS Toulouse[EB/OL]. [2019-08-15]. http: //icgem. gfz-potsdam. de/Foerste-et-al-EIGEN-6C4. pdf.

[15] HIRT C. Prediction of vertical deflections from high-degree spherical harmonic synthesis and residual terrain model data[J]. Journal of Geodesy, 2010, 84(3): 179-190.

[16] NAGY D, PAPP G, BENEDEK J. The gravitational potential and its derivatives for the prism[J]. Journal of Geodesy, 2000, 74(7/8): 552-560.

[17] HARDY R L. Multiquadric equations of topography and other irregular surfaces[J]. Journal of Geophysical Research, 1971, 76(8): 1905-1915.

[18] SANDWELL D T, MüLLER R D, SMITH W H F, et al. New global marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure[J]. Science, 2014, 346(6205): 65-67.

[19] KUHN M, HIRT C. Topographic gravitational potential up to second-order derivatives: an examination of approximation errors caused by rock-equivalent topography (RET)[J]. Journal of Geodesy, 2016, 90(9): 883-902.

[20] HIRT C, REXER M. Earth2014: 1 arc-min shape, topography, bedrock and ice-sheet models-available as gridded data and degree-10, 800 spherical harmonics[J]. International Journal of Applied Earth Observation and Geoinformation, 2015, 39: 103-112.

[21] 李建成, 褚永海, 徐新禹.區(qū)域與全球高程基準(zhǔn)差異的確定[J].測(cè)繪學(xué)報(bào), 2017, 46(10): 1262-1273.

Research on improved datum surface of residual terrain model and polyhedral function

CHU Yonghai1,2, Van Phi NGUYEN1, KE Baogui3, WANG Haihong1,2

（1. Schoolof Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;2. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan 430079, China; 3. Institute of Geodesy and Geodynamic, Chinese Academy of Surveying and Mapping, Beijing 100036, China）

In order to further improve the accuracy of height datum, the paper studied on the effect of RTM in the determination of local height datum in offshore area: the land and sea unified rock-equivalent topography was used to establish the offshore area RTM and residual height anomaly, and the corrected surface of residual height anomaly was fitted by combining a polyhedral function to improve the accuracy of height datum surface model in offshore area. Experimental result showed that: in the offshore area, RTM height anomaly could help improve the accuracy of EIGEN-6C4 by about 3 mm; moreover, the accuracy of the EIGEN-6C4 model could be improved from 9.0 cm to 3.4 cm after the polyhedral function fitting.

height anomaly; residual terrain model(RTM); rock-equivalent terrain(RET); polyhedral function; height datum

P228

A

2095-4999(2020)03-0007-08

褚永海，阮文飛，柯寶貴，等. 剩余地形模型及多面函數(shù)改進(jìn)基準(zhǔn)面的研究[J]. 導(dǎo)航定位學(xué)報(bào), 2020,8(3): 7-14.（CHU Yonghai, Van Phi NGUYEN, KE Baogui, et al. Research on improved datum surface of residual terrain model and polyhedral function[J]. Journal of Navigation and Positioning, 2020, 8(3): 7-14.）

10.16547/j.cnki.10-1096.20200302.

2020-01-08

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2016YFB0501702）；國(guó)家自然科學(xué)基金項(xiàng)目（41774010, 41974016）。

褚永海（1976—），男，貴州普定人，博士，副教授，碩士研究生導(dǎo)師，研究方向?yàn)樾l(wèi)星測(cè)高技術(shù)應(yīng)用、高程基準(zhǔn)確定等。