利用測高重力異常建立南中國海海底地形模型

歐陽明達 孫中苗 翟振和 劉曉剛

1 信息工程大學地理空間信息學院，鄭州市科學大道62號，450001

2 西安測繪研究所，西安市雁塔路中段1號，710054

3 地理信息工程國家重點實驗室，西安市雁塔路中段1號，710054

4 西安測繪信息技術總站，西安市西影路36號，710054

海洋衛星雷達測高技術的出現，使海底地形測量技術得到飛速發展［1－3］。利用測高重力異常反演海底地形主要包括兩類方法：解析法和統計法［2，4］。解析法基于Parker理論和板塊模型，利用響應函數反演海底地形［5－9］。統計法最早由Tscherning提出，將重力數據進行統計迭代，改善了已有的水深模型［10］。此后，翟國君等［11］也使用類似方法開展了地形反演研究。2008年，美國地球物理中心（NGDC）發布1′分辨率的ETOPO1海深模型［12］。2011年，丹麥國家空間中心（DNSC）發布1′分辨率的重力數值模型DTU10GRAV，其殘差的均值為0.39mGal，標準差為3.82mGal，最大值為36.89mGal［13］。本文利用DTU10GRAV 海洋重力異常模型和ETOPO1海底地形模型，采用解析法反演南中國海1′×1′的海底地形，并利用ETOPO1海底地形模型和實測檢核點水深對其精度進行評價。

1 海深反演的基本原理

Parker將引力位引入到頻率域，給出地下物質界面起伏所產生的重力異常的頻率域表達式［14］：

式中，G為地球引力常數；ρc為海底洋殼密度；ρw為海水密度；k＝（kx，ky）＝（1／λx，1／λy），其中，（kx，ky）和（λx，λy）為x和y方向的頻率和波長；F為傅里葉變換；d為平均海深；，為圓頻率。

當海底地形負荷的波長小于地殼的彈性厚度時，式（1）中的第一項起主要作用，它僅取決于平均海深。此時，重力異常和海底地形之間是線性關系［15］：

G（k）和H（k）分別為重力異常和水深的傅里葉變換，Z1（k）為響應函數。

Z1（k）中的2πG（ρc－ρw）是布格常數，exp（－2πkd是向上延拓（從海底到海面）。將上式逆變換，可以通過衛星測高重力異常計算海底地形。在頻率域，該模型如下：

上式為未考慮補償模型的響應函數關系式。當海底地形負荷波長大于海底洋殼彈性厚度時，需要顧及Airy補償模型。這時的響應函數需要顧及兩個層面，一個為平均海深，另一個為補償深度，響應函數表達式為［15］：

Tc為平均地殼厚度，其他參數含義同上。顧及撓曲補償模型的響應函數表達式為［16］：

式中，ρi為撓曲內填補物質的密度，ρ2、ρ3為上、下地殼層密度，ρm為地幔層密度，T2、T3為上、下地殼層厚度。φ（k）表達式如下：

D為巖石圈的撓曲剛度，它與巖石圈的有效彈性厚度Te有關：

υ為巖石圈的泊松比，E為彈性模量。

從式（3）～（5）可以看出，重力異常與海水深度的響應函數和地幔密度、海底洋殼密度、地殼密度、地殼厚度、巖石圈撓曲剛度和有效彈性厚度相關。

2 響應函數

頻率域內重力異常和海底地形之間的關系通過響應函數Z（k）建立，它具有線性、各向同性等性質。利用式（3）～（5）的逆變換計算響應函數理論值，其所采用的地球物理參數見表1。圖1 分別為未補償模型、Airy補償模型、撓曲補償模型的響應函數分別受平均海水深度d、平均地殼厚度Tc和巖石圈有效彈性厚度Te影響所發生的變化。可以看出，重力異常和海底地形的關系是非線性的，但是在波段10～110km 范圍內呈近似的線性關系。在小于110km 的波長范圍，d對未補償響應函數模型的影響是不同的。對Airy補償和撓曲補償響應函數模型而言，其形態基本相同。在大于110km 的波長范圍，響應函數不僅依賴于d，Tc、Te都會對其造成很大影響。因而，當采用響應函數建立兩者關系，實現海底地形反演時，獲知精確的地球物理參數十分必要，這些地球物理參數可以參閱相關文獻或根據諸如Crust2.0等全球地殼模型等獲取。

表1 響應函數中假設的地球物理參數Tab.1 Summary of parameters assumed in the gravity modeling and bathymetry prediction

圖1 響應函數值Fig.1 Value of response function

3 實驗結果與分析

3.1 數據準備和前期處理

研究范圍位于南中國海112°～119°E、12°～20°N。作為檢核的船測水深數據來自美國地球物理中心（NGDC）網站，共計10 529 個，如圖2所示。其背景為ETOPO1海深模型，黃色點為檢核點，紅色和綠色點為用作檢核的Cruise＿01 和Cruise＿02船測軌跡。測高重力異常數據（分辨率為1′）來自丹麥科技大學空間中心，如圖3所示。

3.2 處理流程

顧及Airy補償和撓曲補償的響應函數模型涉及物理參數眾多，過程復雜，將另文給出。以下僅采用較簡單的未補償響應函數模型對海底地形在15～110km 的中高頻分量進行反演。根據文獻［2］，密度差異常數Δρ確定為1 640kg／m3，采用移去－恢復法進行海底地形反演。

1）采用110km 的高斯函數對初始重力異常進行低通濾波，獲得參考重力異常。將參考重力異常從初始異常值中移去，得到殘余重力異常。高斯函數進行濾波，獲得參考海深模型，其均值即為平均海深d。

圖2 船測航跡分布Fig.2 Distribution of shipborne tracks

2）對該海域的ETOPO1海深模型也采用該

3）將殘余重力異常通過未補償響應函數模型進行反演計算，得到殘余海深。該殘余海深和參考海深模型之和即為南中國海海深模型。

值得注意的是，由于響應函數中存在exp（2πkd）項，該項為向下延拓，會帶來很大的高頻噪聲。根據文獻［5］，選用波長為15km 的高斯濾波函數進行抑制，得到最終的海深模型。由圖4可以看出，大部分海域內兩模型符合較好，但在海島附近以及多海山地區兩者差異較大，說明反演模型在該部分海域內的精度不甚理想。

圖3 測高自由空間重力異常Fig.3 Altimetry－derived gravity anomalies

圖4 南中國海海底地形模型及其比較Fig.4 Bathymetry models in the South Sea and differences of this two

3.3 精度評價

用反演模型插值計算得到檢核點處水深，并與實測值比較獲得其殘差。殘差主要集中在0～200m 范圍，部分殘差過大的予以剔除。處理后的殘差和相對精度指標見表2。表2同時示出了ETOPO1模型在該海域的精度統計信息。可以看出，本文反演模型的相對精度保持在10%左右，精度水平略低于ETOPO1模型。

圖5示出了殘差、相對精度與水深和重力異常的關系。在深水區，殘差和相對精度的數值較小，較大數值大多出現在－2 500～0 m 的水域。殘差、相對精度與重力異常的相關性不太明顯，但總的來說，深水海域、低異常值區間的反演精度明顯好于淺水海域以及高重力異常值區間。

表3給出了殘差、相對精度在不同水深層的量化結果。在－3 000m 以下海域的大多數殘差值位于0～170 m 區間，相對精度保持在5%以內；在－3 000m 以上海域，殘差值和相對精度出現大的突變。結合圖4可知，該水深層以上多分布有海山和大陸架，量化顯示結果與§3.2結論一致。

表2 殘差、相對精度結果統計Tab.2 Statistics of the differences and relative precision

圖5 殘差、相對精度與水深和重力異常的關系Fig.5 The relationship between the differences，relative precision and the bathymetry and gravity anomaly

表3 不同水深區間的殘差、相對精度結果統計Tab.3 Statistics of the differences and relative precision in different intervals of deep

圖6為Cruise＿01和Cruise＿02船測航跡測深剖面比較，實線為實測水深，點線為反演模型插值水深，兩者走勢基本一致，但在海山處插值水深往往低于實測水深，且出現較大的殘差量級。反演海深模型過低估計了海山峰值，不能被水下潛器用于導航。

圖6 Cruise＿01和Cruise＿02航跡比較Fig.6 Comparisons for bathymetry along the cruise＿01and cruise＿02

4 結 語

本文采用解析法反演了南中國海在中短波范圍的1′×1′海底地形，反演模型的相對精度為10%左右，略低于ETOPO1模型，在多海山地區以及大陸架覆蓋海域精度稍差，且存在過低估計海山峰值的問題。在后續研究中，可以嘗試通過最小二乘配置方法對其作進一步改善。

［1］黃謨濤，翟國君，管崢，等.海洋重力場測定及其應用［M］.北京：測繪出版社，2005（Huang Motao，Zhai Guojun，Guan Zheng，et al.The Marine Gravity Field Determination and Its Application［M］.Beijing：Surveying and Mapping Press，2005）

［2］王勇，許厚澤，詹金剛.中國海及其鄰近海域高分辨率海底地形［J］.科學通報，2001，46（11）：956－960（Wang Yong，Xu Houze，Zhan Jingang.High Resolution Sea Floor Topography of China Sea and Its Adjacent Areas［J］.Chinese Science Bulletin，2001，46（11）：956－960）

［3］Neumann G A，Forsyth D W，Sandwell D.Comparison of Marine Gravity from Shipboard and High－density Satellite Altimetry Along the Mid－Atlantic Ridge［J］.J Geophys Res，1993，20：1 639－1 642

［4］黃謨濤，翟國君.利用衛星測高資料反演海底地形研究［J］.武漢大學學報：信息科學版，2002，27（2）：133－137（Huang Motao，Zhai Guojun.The Recovery of Bathymetry from Altimeter Data［J］.Geomatics and Information Science of Wuhan University，2002，27（2）：133－137）

［5］李大煒，李建成，豐海，等.利用衛星測高數據反演中國海及鄰近海域海底地形［J］.大地測量與地球動力學，2009，29（4）：70－73（Li Dawei，Li Jiancheng，Feng Hai，et al.Research on Inversion of Seafloor Topography of China Sea and Its Adjacent Areas from Satellite Altimetric Data［J］.Journal of Geodesy and Geodynamics，2009，29（4）：70－73）

［6］方劍，張赤軍.中國海及鄰近海域2′×2′海底地形［J］.武漢大學學報：信息科學版，2003，28（1）：38－40（Fang Jian，Zhang Chijun.2′×2′Sea Floor Bathymetry Prediction of China Sea and Its Vicinity［J］.Geomatics and Information Science of Wuhan University，2003，28（1）：38－40

［7］羅佳.由測高數據和海洋重力資料聯合反演海底地形［D］.武漢：武漢大學，2000（Luo Jia.Bathymetry Inversion Combining Satellite Altimetry and Marine Gravity Data［D］.Wuhan：Wuhan University，2000

［8］聶琳娟，吳云孫，金濤勇，等.基于海水質量虧損引起的重力異常反演南海海底地形［J］.大地測量與地球動力學，2012，32（1）：43－46（Nie Linjuan，Wu Yunsun，Jin Taoyong，et al.Inversion of Submarine Topography of South China Sea by Using Gravity Anomaly Caused by Mass Deficiency［J］.Journal of Geodesy and Geodynamics，2012，32（1）：43－46）

［9］Watts A B，Sandwell D T，Smith W H F，et al.Global Gravity，Bathymetry，and the Distribution of Submarine Volcanism Through Space and Time［J］.J Geophys Res，2006，111：B08408

［10］Tscherning C C，Knudsen P，Foesberg R.First Experiments with Improvement of Depth Information Using Gravity Anomalies in the Mwditerranean Sea［R］.University of Thessaloniki，1994

［11］Zhai Guojun，Huang Motao，Ouyang Yongzhong，et al.A Modified Method for Recovering Bathymetry from Altimeter Data［J］.International Association of Geodesy Symposia，2003，126：84－88

［12］Amante C，Eakins B W.ETOPO1 1Arc－minute Global Relief Model：Procedures，Data Sources and Analysis［R］.NOAA Technical Memorandum，2009

［13］Baltazar O，Knudsen A P.The DNSC08GRA Global Marine Gravity Field from Double Retracked Satellite Altimetry［J］.Journal of Geodesy，2010，84：191－199

［14］Parker R L.The Rapid Calculation of Potential Anonmalies［J］.J Geophys Res，1972，31：447－455

［15］Watts A B.An Analysis of Isostasy in the World’s Oceans：1 Hawaiian Emperor Seament Chain［J］.J Geophys Res，1978，83：5 989－6 004

［16］Watts A B.Isostasy and Flexure of the Lithosphere［M］.New York：Cambridge University Press，2001