剩余地形模型及多面函數改進基準面的研究

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

剩余地形模型及多面函數改進基準面的研究

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

（1. 武漢大學 測繪學院，武漢 430079；2. 武漢大學 地球空間環(huán)境與大地測量教育部重點實驗室，武漢 430079；3. 中國測繪科學研究院 大地測量與地球動力學研究所，北京 100036)

為了進一步提高高程基準面的精度，研究剩余地形模型（RTM）在近海區(qū)域高程基準面中的影響：利用陸海統(tǒng)一巖石等效地形（RET）構建近海區(qū)域RTM和殘余高程異常；并聯(lián)合多面函數擬合殘余高程異常改正面，實現近海區(qū)域高程基準面模型精度的提高。實驗結果表明：在近海區(qū)域，利用RTM高程異常可提高新技術改進后的歐洲地球重力聯(lián)合模型（EIGEN-6C4）的精度約為3 mm；多面函數擬合后，可以將EIGEN-6C4模型的精度由9.0 cm提高到3.4 cm。

高程異常；剩余地形模型；巖石等效地形；多面函數；高程基準

0 引言

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

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

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

1 原理方法

1.1 剩余地形建模

1.2 高程異常

式中：G為萬有引力常數；ρ為長方體的密度；為長方體頂點坐標；l為距離，用坐標表示成。長方體頂點E的坐標可以用表示；頂點K的坐標用表示，其他點都可以用點E和點K的坐標表示，具體位置關系參見圖1。

根據布隆斯（Bruns）公式，每個長方體地形引力位對高程異常的貢獻[15]表示成

1.3 殘余高程異常及多面函數擬合

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

2 實驗與結果分析

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

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

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

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

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

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

此外，將GPS/水準點共計152個分為2組：第1組作為多面函數擬合計算點（如圖3所示），共計88個；第2組為檢核點（如圖4所示），共計64個點，用來對擬合改正面進行獨立檢核。

圖3 擬合GPS水準點分布

圖4 檢核GPS水準點分布

2.2 積分半徑對RTM高程異常的影響

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

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

圖5 RTM高程異常與積分半徑的關系

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

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

圖6 RTM高程異常分布

圖7 多面函數法擬合殘余高程異常

圖8 擬合改正后的高程基準面

2.4 精度檢核

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

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

表3 EIGEN-6C4/RTM高程異常與GPS/水準高程異常比較

續(xù)表3

3 結束語

剩余地形模型在高海拔地形起伏地區(qū)，特別是在無實測重力地區(qū)，能夠有效地改善全球重力場模型短波分量因模型截斷的影響[7,9-11]；在近海地區(qū)，剩余地形模型對高程異常的改進可提高的精度約為6.5 mm[8]。為了分析剩余地形模型在中國近海的應用情況，選取青島及周邊區(qū)域進行計算分析，包括數值積分半徑的選取、剩余地形模型的貢獻、擬合改正后的精度情況等，數值結果表明：

剩余地形模型對重力場模型短波分量的改善，不僅可以用于高程異常的改善，也可以用于重力、垂線偏差的改善[8]，這對于無實測重力數據地區(qū)重力場模型的構建，具有重要的作用。近幾十年來，RTM技術已經被廣泛用于大地測量領域，但以目前的方式使用時，還有一些問題需要解決：例如，計算點位于較低的山谷時，為了保證引力位的調和性，需要進行調和改正；當數值積分半徑較大時，需要顧及地球曲率的影響。另外，在殘余地形擬合改正面中，海域因無GPS/水準約束出現扭曲，這可以將海面地形作為“虛擬”GPS/水準點，以保證將陸地高程傳遞到海域范圍?？傊绻赗TM建模中進行調和改正和地球曲率改正，在擬合曲面時對海洋進行約束，可以進一步提高程基準面的精度，并且實現陸海高程基準的統(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高程轉換方法[J].測繪學報, 2012, 41(1): 25-32.

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

[11] 翟長治, 姚宜斌, 岳順.基于EGM2008和剩余地形模型的區(qū)域似大地水準面精化方法[J]. 大地測量與地球動力學, 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ū)域與全球高程基準差異的確定[J].測繪學報, 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

褚永海，阮文飛，柯寶貴，等. 剩余地形模型及多面函數改進基準面的研究[J]. 導航定位學報, 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

國家重點研發(fā)計劃項目（2016YFB0501702）；國家自然科學基金項目（41774010, 41974016）。

褚永海（1976—），男，貴州普定人，博士，副教授，碩士研究生導師，研究方向為衛(wèi)星測高技術應用、高程基準確定等。