歐陽明達(dá) 孫中苗 翟振和 劉曉剛

1 信息工程大學(xué)地理空間信息學(xué)院，鄭州市科學(xué)大道62號(hào)，450001

2 西安測(cè)繪研究所，西安市雁塔路中段1號(hào)，710054

3 地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安市雁塔路中段1號(hào)，710054

4 西安測(cè)繪信息技術(shù)總站，西安市西影路36號(hào)，710054

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

1 海深反演的基本原理

Parker將引力位引入到頻率域，給出地下物質(zhì)界面起伏所產(chǎn)生的重力異常的頻率域表達(dá)式［14］：

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

當(dāng)海底地形負(fù)荷的波長(zhǎng)小于地殼的彈性厚度時(shí)，式（1）中的第一項(xiàng)起主要作用，它僅取決于平均海深。此時(shí)，重力異常和海底地形之間是線性關(guān)系［15］：

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

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

上式為未考慮補(bǔ)償模型的響應(yīng)函數(shù)關(guān)系式。當(dāng)海底地形負(fù)荷波長(zhǎng)大于海底洋殼彈性厚度時(shí)，需要顧及Airy補(bǔ)償模型。這時(shí)的響應(yīng)函數(shù)需要顧及兩個(gè)層面，一個(gè)為平均海深，另一個(gè)為補(bǔ)償深度，響應(yīng)函數(shù)表達(dá)式為［15］：

Tc為平均地殼厚度，其他參數(shù)含義同上。顧及撓曲補(bǔ)償模型的響應(yīng)函數(shù)表達(dá)式為［16］：

式中，ρi為撓曲內(nèi)填補(bǔ)物質(zhì)的密度，ρ2、ρ3為上、下地殼層密度，ρm為地幔層密度，T2、T3為上、下地殼層厚度。φ（k）表達(dá)式如下：

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

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

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

2 響應(yīng)函數(shù)

頻率域內(nèi)重力異常和海底地形之間的關(guān)系通過響應(yīng)函數(shù)Z（k）建立，它具有線性、各向同性等性質(zhì)。利用式（3）～（5）的逆變換計(jì)算響應(yīng)函數(shù)理論值，其所采用的地球物理參數(shù)見表1。圖1 分別為未補(bǔ)償模型、Airy補(bǔ)償模型、撓曲補(bǔ)償模型的響應(yīng)函數(shù)分別受平均海水深度d、平均地殼厚度Tc和巖石圈有效彈性厚度Te影響所發(fā)生的變化?？梢钥闯觯亓Ξ惓：秃５椎匦蔚年P(guān)系是非線性的，但是在波段10～110km 范圍內(nèi)呈近似的線性關(guān)系。在小于110km 的波長(zhǎng)范圍，d對(duì)未補(bǔ)償響應(yīng)函數(shù)模型的影響是不同的。對(duì)Airy補(bǔ)償和撓曲補(bǔ)償響應(yīng)函數(shù)模型而言，其形態(tài)基本相同。在大于110km 的波長(zhǎng)范圍，響應(yīng)函數(shù)不僅依賴于d，Tc、Te都會(huì)對(duì)其造成很大影響。因而，當(dāng)采用響應(yīng)函數(shù)建立兩者關(guān)系，實(shí)現(xiàn)海底地形反演時(shí)，獲知精確的地球物理參數(shù)十分必要，這些地球物理參數(shù)可以參閱相關(guān)文獻(xiàn)或根據(jù)諸如Crust2.0等全球地殼模型等獲取。

表1 響應(yīng)函數(shù)中假設(shè)的地球物理參數(shù)Tab.1 Summary of parameters assumed in the gravity modeling and bathymetry prediction

圖1 響應(yīng)函數(shù)值Fig.1 Value of response function

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

3.1 數(shù)據(jù)準(zhǔn)備和前期處理

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

3.2 處理流程

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

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

圖2 船測(cè)航跡分布Fig.2 Distribution of shipborne tracks

2）對(duì)該海域的ETOPO1海深模型也采用該

3）將殘余重力異常通過未補(bǔ)償響應(yīng)函數(shù)模型進(jìn)行反演計(jì)算，得到殘余海深。該殘余海深和參考海深模型之和即為南中國(guó)海海深模型。

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

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

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

3.3 精度評(píng)價(jià)

用反演模型插值計(jì)算得到檢核點(diǎn)處水深，并與實(shí)測(cè)值比較獲得其殘差。殘差主要集中在0～200m 范圍，部分殘差過大的予以剔除。處理后的殘差和相對(duì)精度指標(biāo)見表2。表2同時(shí)示出了ETOPO1模型在該海域的精度統(tǒng)計(jì)信息?？梢钥闯?，本文反演模型的相對(duì)精度保持在10%左右，精度水平略低于ETOPO1模型。

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

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

表2 殘差、相對(duì)精度結(jié)果統(tǒng)計(jì)Tab.2 Statistics of the differences and relative precision

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

表3 不同水深區(qū)間的殘差、相對(duì)精度結(jié)果統(tǒng)計(jì)Tab.3 Statistics of the differences and relative precision in different intervals of deep

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

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

4 結(jié) 語

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

［1］黃謨濤，翟國(guó)君，管崢，等.海洋重力場(chǎng)測(cè)定及其應(yīng)用［M］.北京：測(cè)繪出版社，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］王勇，許厚澤，詹金剛.中國(guó)海及其鄰近海域高分辨率海底地形［J］.科學(xué)通報(bào)，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，F(xiàn)orsyth 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］黃謨濤，翟國(guó)君.利用衛(wèi)星測(cè)高資料反演海底地形研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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］李大煒，李建成，豐海，等.利用衛(wèi)星測(cè)高數(shù)據(jù)反演中國(guó)海及鄰近海域海底地形［J］.大地測(cè)量與地球動(dòng)力學(xué)，2009，29（4）：70－73（Li Dawei，Li Jiancheng，F(xiàn)eng 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］方劍，張赤軍.中國(guó)海及鄰近海域2′×2′海底地形［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，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］羅佳.由測(cè)高數(shù)據(jù)和海洋重力資料聯(lián)合反演海底地形［D］.武漢：武漢大學(xué)，2000（Luo Jia.Bathymetry Inversion Combining Satellite Altimetry and Marine Gravity Data［D］.Wuhan：Wuhan University，2000

［8］聶琳娟，吳云孫，金濤勇，等.基于海水質(zhì)量虧損引起的重力異常反演南海海底地形［J］.大地測(cè)量與地球動(dòng)力學(xué)，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，F(xiàn)oesberg 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