應(yīng)用重力地質(zhì)法反演馬里亞納海溝地形

王永康周興華唐秋華王盼龍姜傳苓劉 豪張孝首

(1.自然資源部 第一海洋研究所,山東 青島266061;2.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院,山東 青島266590;3.南京大學(xué) 地理與海洋科學(xué)學(xué)院,江蘇 南京210023)

海底地形是海洋研究的重要內(nèi)容,高精度、大范圍的海底地形數(shù)據(jù)是研究海底地理特征、地殼結(jié)構(gòu)、海流、內(nèi)波等的有力資料。傳統(tǒng)的水深數(shù)據(jù)是通過船載聲吶測量獲取,盡管聲吶精度高,但是作業(yè)效率太低,不適合進(jìn)行大范圍的海底地形獲取。自從Seasat、Geosat等衛(wèi)星和Cryosat、ERS任務(wù)的成功實(shí)施,許多學(xué)者利用測高數(shù)據(jù),采用最小二乘配置法、逆Stokes法和逆Vening Meinesz法反演了高精度、高分辨率的海洋重力異常[1-5],并將測高重力異常應(yīng)用到海底地形的反演中[6-7]。重力地質(zhì)法是利用重力異常反演海洋水深的一種有效方法,最初應(yīng)用于冰川沉積物下的基巖高程測量[8-9],但是由于地下物質(zhì)的密度變化,該方法在陸地區(qū)域的使用受到了限制,海底洋殼和海水之間密度差變化較小,這使得重力地質(zhì)法十分適宜于利用衛(wèi)星測高重力異常反演海底地形的研究[10]。重力地質(zhì)法已經(jīng)成功應(yīng)用于許多海域的地形反演:Roman[11]反演了巴倫支海和格陵蘭島周圍海域的海底地形;Kim 和Yun[12]反演了韓國西海海底地形;Kim 等[13]反演了德雷克海峽的水深;胡敏章等[14]反演了皇帝海山的海底地形,并且指出參考深度對反演結(jié)果影響有限;歐陽明達(dá)等[10]、彭聰?shù)萚15]、李倩倩和鮑李峰[16]和Hsiao等[17]反演了中國南海的海底地形;范雕等[18]反演了馬里亞納海溝的海底地形。

利用充足的研究區(qū)域水深數(shù)據(jù)計(jì)算不同密度差常數(shù)下的水深模型,與檢核點(diǎn)對比得到最優(yōu)密度差常數(shù),在此密度差常數(shù)下利用重力地質(zhì)法對測高衛(wèi)星重力異常進(jìn)行海底地形反演,得到1′×1′分辨率的海底地形;通過和船測水深點(diǎn)、ETOPO1模型對比,評(píng)估了GGM 模型精度,分析了誤差原因。前人主要通過直接觀察來判斷GGM 模型比ETOPO1模型和直接格網(wǎng)化模型更加精細(xì),本文通過功率譜分析,檢測3個(gè)水深模型在不同波長下的能量大小,評(píng)價(jià)3個(gè)地形的精細(xì)程度,最后通過2條船測水深測線進(jìn)一步探究GGM 預(yù)測水深和實(shí)測水深的關(guān)系。

1 GGM 原理

重力異?？煞譃槎滩ㄖ亓Ξ惓：烷L波重力異常,短波重力異常是受當(dāng)?shù)氐匦纹鸱挠绊懚a(chǎn)生,長波重力異常是由地球深層的質(zhì)量變化引起的。則觀測到的海面重力異常[17]:

式中:Δgobs為自由空間重力異常,Δgshort為短波重力異常,Δglong為長波重力異常。

圖1是GGM 的原理圖,點(diǎn)j n(n=1,2,3,…,k)為水深控制點(diǎn),根據(jù)水深控制點(diǎn)來計(jì)算重力異常短波分量,可用布格板公式:

式中:Δgs(j)為j點(diǎn)上的短波重力異常;G 為引力常數(shù)(6.67×10-8cm3·g-1·s-2);Δρ為海水和海底洋殼的密度差異常數(shù);E j為j點(diǎn)處的海深;D為參考深度,一般取控制點(diǎn)中最大水深。

控制點(diǎn)長波重力異常可以由觀測到的自由空間重力異常減去控制點(diǎn)處短波重力異常計(jì)算得出,再把控制點(diǎn)長波重力異常插值生成長波重力異常場,把長波重力異常場從觀測的自由空間重力異常場中移除,即可得到測區(qū)短波重力異常場。任意一點(diǎn)i的水深值E i為

式中:Δg is為i點(diǎn)的短波重力異常,由測區(qū)短波重力異常場通過適當(dāng)插值方法得到。

圖1 GGM 原理Fig.1 Principe of the gravity geology method(GGM)

圖2 控制點(diǎn)和檢核點(diǎn)分布Fig.2 The distribution of controlling points and checking points

2 海底地形反演

2.1 測區(qū)概況及數(shù)據(jù)準(zhǔn)備

馬里亞納海溝是太平洋板塊自東向西俯沖于菲律賓板塊之下形成的一條向東弧形凸出,近南北向延伸的深溝,北起硫黃列島、西南至雅浦島附近,全長2550 km,平均寬70 km,大部分水深8000 m 以上[19]。本文研究區(qū)域位于142°36′～147°18′E,23°～27°N,是馬里亞納海溝的一部分,該區(qū)域地形復(fù)雜,最深約9000 m,最淺約700 m,多海山和海溝,地形起伏大,非常適合用來檢驗(yàn)GGM 的反演效果。

船測水深來自NGDC(National Geophysical Data Center)的單波束海深測量數(shù)據(jù),因?yàn)榇瑴y水深存在粗差,數(shù)千米深海海域的NGDC 船測水深值只有數(shù)百米,經(jīng)過和ETOPO1模型對比,人工剔除粗差點(diǎn),最終選取6736個(gè)水深控制點(diǎn)和1683個(gè)檢核點(diǎn),船測水深控制點(diǎn)和檢核點(diǎn)分布均勻且不重復(fù),控制點(diǎn)約是檢核點(diǎn)5倍(圖2)。測區(qū)自由空間重力異常數(shù)據(jù)來自SIO(Scripps Institution of Oceanography)(圖3a),版本為v24.1,與船測重力對比,精度約2 mgal[20],適合用來進(jìn)行水深反演,船測控制點(diǎn)的重力異常由此內(nèi)插而得到。檢核模型為ETOPO1模型,ETOPO1是由NGDC 和NOAA(National Oceanic and Atmospheric Administration)聯(lián)合發(fā)布的全球地形高程模型,分辨率為1′×1′。

圖3 研究區(qū)域重力異常Fig.3 Gravity anomalies in the study area

2.2 密度差常數(shù)的確定

采用重力地質(zhì)法反演海底地形時(shí),海水和洋殼的密度差常數(shù)的確定是一個(gè)關(guān)鍵問題。密度差常數(shù)的精確確定,可使短波和長波重力異常從自由空間重力異常中分離出來,確保短波重力異常和水深保持良好的線性關(guān)系,由圖3可知,短波重力異常和水深有較強(qiáng)的線性關(guān)系,長波重力異常和水深的關(guān)系不明顯。如果密度差常數(shù)預(yù)測過小,預(yù)測海深將會(huì)偏大,反之,預(yù)測海深將會(huì)接近于參考深度[10]。全球洋殼平均密度為2.67 g/cm3,海水平均密度為1.03 g/cm3,因而全球平均密度差常數(shù)為1.64 g/cm3,但是局部地形不宜使用平均密度差常數(shù),應(yīng)該精確估計(jì)研究區(qū)域的密度差常數(shù)。密度差常數(shù)一般由2種方法來確定:一種是基于快速傅里葉變換的向下延拓法[7],適用于水深數(shù)據(jù)缺乏的區(qū)域;另一種是利用水深控制點(diǎn)計(jì)算不同密度差常數(shù)下的水深模型,再與檢核點(diǎn)進(jìn)行比較,得到不同密度差常數(shù)下的預(yù)測水深和實(shí)測水深的相關(guān)關(guān)系和標(biāo)準(zhǔn)差圖,從而得到最優(yōu)密度差參數(shù)[14],這種方法適應(yīng)于水深數(shù)據(jù)較多的區(qū)域。本研究水深控制點(diǎn)數(shù)量充足,所以采用第2種方法,由圖4可知,依據(jù)標(biāo)準(zhǔn)差較小而相關(guān)系數(shù)較大原則,密度差常數(shù)選取1.20 g/cm3時(shí),標(biāo)準(zhǔn)差為最小值152.9 m,相關(guān)系數(shù)為最大值0.9943。

圖4 密度差常數(shù)確定Fig.4 The determination of density contrast

2.3 GGM 水深模型建立

先根據(jù)船測水深控制點(diǎn),依照式(2)計(jì)算船測控制點(diǎn)處的短波重力異常;接著用控制點(diǎn)處的自由空間重力異常減去控制點(diǎn)處的短波重力異常,得到控制點(diǎn)處的長波重力異常,使用GMT(the Generic Mapping Tools)中張力樣條函數(shù)將控制點(diǎn)處長波重力異常插值得到1′×1′分辨率的測區(qū)長波重力異常場(圖3b);然后利用測區(qū)自由空間重力異常場減去測區(qū)長波重力異常場,得到測區(qū)短波重力異常場(圖3c);再依據(jù)式3便可得到測區(qū)任意一點(diǎn)的水深值,最終使用張力樣條函數(shù)插值生成測區(qū)1′×1′分辨率的海底地形圖(以下稱為GGM 模型),見圖5b。

圖5 研究區(qū)域水深模型Fig.5 Bathymetry model for the study area

3 結(jié)果分析

為了對GGM 反演的海底地形進(jìn)行精度分析,使用了ETOPO1模型、控制點(diǎn)直接插值得到的1′×1′分辨率的水深模型(以下稱之為模型I,見圖5c)、船測檢核點(diǎn)水深和兩條測線。

3.1 模型對比分析

首先對GGM 模型、ETOPO1模型和模型I進(jìn)行統(tǒng)計(jì)分析和相關(guān)性分析,模型統(tǒng)計(jì)和相關(guān)性數(shù)據(jù)見表1。根據(jù)統(tǒng)計(jì)信息可知GGM 模型、ETOPO1模型和模型I之間的最大值、最小值、平均值等差別很小,表明3個(gè)模型的相似度很高。根據(jù)相關(guān)性分析可知,GGM 模型和ETOPO1模型的相關(guān)系數(shù)為0.9883,GGM 模型和模型I的相關(guān)系數(shù)為0.9993,ETOPO1模型和模型I的相關(guān)系數(shù)為0.9882,3個(gè)模型具有很強(qiáng)的相關(guān)性。其中,GGM 模型和模型I的相關(guān)性更強(qiáng)的原因主要有2個(gè)方面:1)GGM 模型是在控制點(diǎn)的基礎(chǔ)上建立的;2)原因是控制點(diǎn)分布相對均勻,直接插值得出的模型也可以較好地反映出海底地形地貌,所以兩者相關(guān)性強(qiáng)。

利用檢核點(diǎn)對3個(gè)模型進(jìn)行精度檢驗(yàn),在檢核點(diǎn)處內(nèi)插得出的水深同檢核點(diǎn)水深的較差信息見表2。由表2可知,GGM 水深值和檢核點(diǎn)水深值較差的標(biāo)準(zhǔn)差和均方根分別為152.9 和153.0 m,相較于ETOPO1模型有著更高的精度。此外,模型I比ETOPO1模型精度高的原因,主要是控制點(diǎn)分布均勻,直接把水深點(diǎn)插值計(jì)算也可以獲取較為精確的海底地形。GGM 模型和模型I與檢核點(diǎn)校差的最大值均超過了2000 m,原因主要有3個(gè):存在粗差水深點(diǎn)未成功剔除,測高重力異常存在較大誤差和地形起伏劇烈導(dǎo)致反演精度下降。具體的原因需要根據(jù)具體的分析而確定。

表1 GGM 模型、ETOPO1模型和模型I的統(tǒng)計(jì)數(shù)據(jù)和相關(guān)系數(shù)Table 1 Statistics for GGM model,ETOPO1 model and Model I and their correlation coefficients

表2 GGM 模型、ETOPO1模型和模型I與檢核點(diǎn)比較Table 2 Comparison between checking points and models(GGM model,ETOPO1 model,and Model I)

從3個(gè)水深模型圖對比可知,ETOPO1模型中有一個(gè)區(qū)域的地形在GGM 模型和模型I上沒有清晰表現(xiàn)出來(圖5中的矩形區(qū)域A),經(jīng)過和圖2對比發(fā)現(xiàn),區(qū)域A 控制點(diǎn)分布稀疏,所以未能反演出來區(qū)域A 的精確地形,因此,GGM 在控制點(diǎn)分布密集海域的預(yù)測精度優(yōu)于控制點(diǎn)分布稀疏海域。把GGM 模型和ETOPO1模型進(jìn)行做差處理(圖6),從GGM 模型和ETOPO1模型差異圖可以清晰看出,在地形變化大的海溝和海山區(qū)域,模型之間差距很大;在地形平坦的區(qū)域,差距很小。

圖6 GGM 模型與ETOPO1模型差異Fig.6 The differences between GGM model and ETOPO1 model

3.2 功率譜密度分析

Kim 等[21]利用功率譜密度法計(jì)算了不同水深模型在不同波長下的能量大小,對比不同模型在不同波長下對地形的表現(xiàn)能力,能量越大,說明模型在此波長下對地形表現(xiàn)能力越強(qiáng)。因此,對GGM 模型、ETOPO1模型和模型I進(jìn)行功率譜密度分析,利用GMT 中g(shù)rdfft命令,計(jì)算3個(gè)模型在不同波長下的能量變化,功率譜密度計(jì)算方法:

式中:P為不同波長下的能量;PSD為功率譜密度(d B)。本文范圍為4.7°×4.0°,所以最大波長約為500 km,結(jié)果見圖7。圖7中以17 km(虛線處)為界分為:在17～500 km 的中長波范圍,3種模型表現(xiàn)能力相當(dāng);在小于17 km 的短波范圍,GGM 模型能量明顯大于模型I和ETOPO1模型,表明GGM 對于微小細(xì)致的短波地形的刻畫優(yōu)于ETOPO1模型和直接格網(wǎng)化結(jié)果。

圖7 GGM 模型、ETOPO1模型和模型Ⅰ功率譜密度比較Fig.7 The comparison in power spectral densities between GGM model,ETOPO model and Model I

為了進(jìn)一步證實(shí)功率譜密度分析的結(jié)果,把3種模型在檢核點(diǎn)處的水深值內(nèi)插出來并和檢核點(diǎn)水深值做差并取絕對值,統(tǒng)計(jì)信息如圖8所示。根據(jù)圖8可知,在較差0～100 m 范圍內(nèi),GGM 模型、ETOPO1模型和模型I所占比例分別是71.95%,59.77%和67.5%;在較差0～300 m 范圍內(nèi),GGM 模型、ETOPO1模型和模型I所占比例分別為94.76%,85.8%和92.81%;進(jìn)一步表明GGM 模型比ETOPO1模型和模型I精度更高,對短波長的水深表現(xiàn)能力更強(qiáng)。

3.3 測線剖面分析

為了進(jìn)一步探究GGM 預(yù)測水深和實(shí)測水深的關(guān)系,本文選擇line1和line2兩條測線剖面進(jìn)行分析,測線貫穿整個(gè)測區(qū),位置見圖9。2條測線既經(jīng)過地勢平坦地區(qū),也經(jīng)過地勢起伏較大的海山和海溝,適合對GGM 模型進(jìn)行評(píng)價(jià)。圖10為水深剖面,經(jīng)過對比可知,在海底地形平坦的地區(qū),GGM 反演的水深和實(shí)測水深較差很小;當(dāng)海底地形開始出現(xiàn)起伏后,GGM 預(yù)測的水深和實(shí)測水深較差便出現(xiàn)較劇烈的變化,進(jìn)一步表明GGM 在地形起伏較小的區(qū)域優(yōu)于地形起伏較大區(qū)域的反演效果。GGM 反演水深和測線水深最大的較差出現(xiàn)在地形起伏劇烈的區(qū)域,最大差值約為1500 m,所以本文反演的GGM 模型與檢核點(diǎn)最大差值2000多m 的原因是地形起伏劇烈導(dǎo)致反演精度下降。

圖9 測線分布Fig.9 Distribution of checking lines

圖10 GGM 模型與檢查線水深剖面Fig.10 Comparison of profiles between GGM model and checking lines

4 結(jié) 語

利用重力地質(zhì)法反演了馬里亞納海溝的海底地形,經(jīng)過模型對比、功率譜分析和測線剖面分析,表明重力地質(zhì)法適合進(jìn)行大范圍的海底地形反演;重力地質(zhì)法模型在短波長地形地貌的刻畫優(yōu)于ETOPO1模型和控制點(diǎn)直接格網(wǎng)化模型;重力地質(zhì)法模型精度主要受到密度差常數(shù)、控制點(diǎn)分布和海底復(fù)雜程度的影響,精確的密度差常數(shù)、控制點(diǎn)的均勻分布和平坦的地形均可以提高重力地質(zhì)法反演精度。所以建議在建立水深模型時(shí)候,若控制點(diǎn)分布密集,采用GGM 方法;若控制點(diǎn)分布稀疏,采用ETOPO1模型。

本研究中重力異常數(shù)據(jù)與船測重力數(shù)據(jù)存在約2 mgal的誤差,若研究區(qū)域有船測重力數(shù)據(jù),可先把船測重力數(shù)據(jù)和衛(wèi)星測高重力數(shù)據(jù)融合,得到高精度融合重力場,然后再應(yīng)用于GGM 水深反演,這個(gè)方法有待于下一步實(shí)驗(yàn)驗(yàn)證。

參考文獻(xiàn)(References):

[1]SANDWELL D T,SMITH W H F.Marine gravity anomaly from Geosat and ERS-1 satellite altimetry[J].Journal of Geophysical Research Solid Earth,1997,102(B5):10039-10054.

[2]HWANG C,KAO E C,PARSONS B.Global derivation of marine gravity anomalies from Seasat,Geosat,ERS-1 and TOPEX/POSEIDON altimeter data[J].Geophysical Journal International,1998,134(2):449-459.

[3]ANDERSEN O B,KNUDSEN P.Global marine gravity field from the ERS-1 and Geosat geodetic mission altimetry[J].Journal of Geophysical Research Oceans,1998,103(C4):8129-8137.

[4]KNUDSEN P,ANDERSEN O B.Improved recovery of the global marine gravity field from the GEOSAT and the ERS-1 geodetic mission altimetry[J].Gravity,Geoid and Marine Geodesy,1997,117:429-436.

[5]HUANG M T,ZHAI G J,GUAN Z,et al.On the recovery of gravity anomalies from altimeter data[J].Acta Geodaetica et Cartograpica Sinica,2001(2):179-184.黃謨濤,翟國君,管錚,等.利用衛(wèi)星測高數(shù)據(jù)反演海洋重力異常研究[J].測繪學(xué)報(bào),2001(2):179-184.

[6]SMITH W H F,SANDWELL D T.Bathymetric prediction from dense satellite altimetry and sparse shipboard bathymetry[J].Journal of Geophysical Research Solid Earth,1994,99(B11):21803-21824.

[7]HWANG C.A bathymetric model for the south china sea from satellite altimetry and depth data[J].Marine Geodesy,1999,22(1):15.

[8]IBRAHIM A,HINZE W J.Mapping buried bedrock topography with gravity[J].Ground Water,2010,10(3):18-23.

[9]NAGARAJAN R.Gravity-geologic investigation of buried bedrock topography in northwestern Ohio[D].Columbus,U.S.:Ohio State University,1994.

[10]OUYANG M D,SUN Z M,ZHAI Z H.Predicting bathymetry in South China Sea using the gravity-geologic method[J].Chinese Journal of Geophysics,2014,57(9):2756-2765.歐陽明達(dá),孫中苗,翟振和.基于重力地質(zhì)法的南中國海海底地形反演[J].地球物理學(xué)報(bào),2014,57(9):2756-2765.

[11]ROMAN D R.An integrated geophysical investigation of Greenland’s tectonic history[D].Columbus,U.S.:Ohio State University,1999.

[12]KIM K B,YUN H S.Satellite-derived bathymetry prediction in shallow waters using the gravity-geologic method:a case study in the West Sea of Korea[J].KSCE Journal of Civil Engineering,2018,22(7):2560-2568.

[13]KIM J W,FRESE R R B V,LEE B Y,et al.Altimetry-derived gravity predictions of bathymetry by the gravity-geologic method[J].Pure and Applied Geophysics,2011,168(5):815-826.

[14]HU M Z,LI J C,JIN T Y.Bathymetry inversion with gravity-geologic method in emperor seamount[J].Geomatics and Information Science of Wuhan University,2012,37(5):610-612,629.胡敏章,李建成,金濤勇.應(yīng)用重力地質(zhì)方法反演皇帝海山的海底地形[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2012,37(5):610-612,629.

[15]PENG C,WANG P,ZHOU X H,et al.Bathymetry inversion based on gravity-geologic method in South China Sea[J].Hydrographic Surveying and Charting,2018,38(2):8-11.彭聰,王鵬,周興華,等.基于重力地質(zhì)法反演南海海底地形[J].海洋測繪,2018,38(2):8-11.

[16]LI Q Q,BAO L F.Predicting submarine topography of the South China Sea from altimetry gravity field with high precision[J].Hydrographic Surveying and Charting,2016,36(2):1-5.李倩倩,鮑李峰.高精度測高重力場反演南海海底地形[J].海洋測繪,2016,36(2):1-5.

[17]HSIAO Y S,HWANG C,CHENG Y S,et al.High-resolution depth and coastline over major atolls of South China Sea from satellite altimetry and imagery[J].Remote Sensing of Environment,2016,176:69-83.

[18]FAN D,LI S S,MENG S Y,et al.Inversion of Mariana Trench seabed terrain using gravity anomalies[J].Journal of Jilin University(Earth Science Edition),2018,48(5):199-208.范雕,李?yuàn)檴?孟書宇,等.利用重力異常反演馬里亞納海溝海底地形[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2018,48(5):199-208.

[19]LIU F L,YANG S X,DENG X G,et al.Depth detection in the Challenger Deep of Mariana Trench[J].Hydrographic Surveying and Charting,2013,33(5):49-52.劉方蘭,楊勝雄,鄧希光,等.馬里亞納海溝“挑戰(zhàn)者深淵”最深點(diǎn)水深探測[J].海洋測繪,2013,33(5):49-52

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

[21]KIM K B,HSIAO Y S,KIM J W,et al.Bathymetry enhancement by altimetry-derived gravity anomalies in the East Sea(Sea of Japan)[J].Marine Geophysical Researches,2010,31(4):285-298.