利用重力異常反演馬里亞納海溝海底地形

范 雕，李姍姍，孟書宇，邢志斌，馮進(jìn)凱

1.信息工程大學(xué)地理空間信息學(xué)院， 鄭州 450001 2.西安航天天繪數(shù)據(jù)技術(shù)有限公司, 西安 710054

0 引言

傳統(tǒng)的測深方法以船舶作為測量儀器的載體，以船舶定位技術(shù)和水聲技術(shù)作為測深技術(shù)的主體，對海深和海底障礙物進(jìn)行測繪。由于技術(shù)條件的限制，依然還有大面積的海域存在數(shù)據(jù)空白，導(dǎo)致人們對大洋底的地貌形態(tài)知之甚少，認(rèn)識還很貧乏[1]。衛(wèi)星測高技術(shù)在研究海洋表面及其內(nèi)部現(xiàn)象、海底構(gòu)造、全球重力場等方面具有無可比擬的優(yōu)勢。它可以大范圍、周期性、快速地觀測海洋上的各種現(xiàn)象及其變化，極大地提高人們對海洋認(rèn)識的廣度和深度，可以使以前幾個(gè)世紀(jì)的測量任務(wù)在幾個(gè)月內(nèi)完成。

如今，利用衛(wèi)星測高技術(shù)確定重力異常的理論和方法基本成熟[2]，通過測高數(shù)據(jù)獲取的重力異常與實(shí)測重力異常的差值均方根在幾個(gè)毫伽范圍內(nèi)[3]，使得依據(jù)重力數(shù)據(jù)與海深的關(guān)系反演海底地形成為可能[4-11]。目前，海底地形反演技術(shù)大多基于Paker[11]異常擾動位公式和Watts[12]的3個(gè)板塊地殼均衡模型，如Smith等[13]的方法。然而，經(jīng)典的海底地形反演方法通常需要引入先驗(yàn)?zāi)Ｐ?、考慮海底均衡狀況，且計(jì)算模型復(fù)雜，如Braitenberg等[14]在利用重力數(shù)據(jù)反演中國南海海底地形的過程中，考慮了地殼均衡、沉積層厚度、有效彈性厚度等因素的影響，過程較為復(fù)雜。重力地質(zhì)法(gravity-geologic method，GGM)與經(jīng)典方法相比，具有模型簡單、易于計(jì)算[15]的優(yōu)點(diǎn)。該方法最初被用于地球物理中反演沉積物基巖深度的研究，由于陸地上地質(zhì)層的密度變化大，使得該方法在陸地上的應(yīng)用效果不佳。而在海底地形的探索方面，洋殼密度和海水密度的差異變化較小，通常認(rèn)為海水內(nèi)部密度不隨深度變化[15]。從而，GGM在反演海洋深度方面有較高的可行性[16]。

以重力異常為數(shù)據(jù)源，依據(jù)GGM反演海底地形的關(guān)鍵在于密度差異常數(shù)的確定。歐陽明達(dá)等[17]利用GGM模型反演的海深值與檢核點(diǎn)海深的相關(guān)系數(shù)，以及反演海深值與檢核點(diǎn)海深差值的標(biāo)準(zhǔn)差，獲得該海域密度差異常數(shù)為1.32 g/cm3，進(jìn)而對中國南海海域進(jìn)行反演分析；胡敏章等[15]應(yīng)用GGM反演皇帝海山的海底地形時(shí)，密度差異常數(shù)的獲取也是依據(jù)檢核點(diǎn)處GGM海深與船測海深的相關(guān)系數(shù)以及海深差值的標(biāo)準(zhǔn)差得到的。本文通過重力異常向下延拓方法確定密度差異常數(shù)，進(jìn)而開展海底地形反演實(shí)驗(yàn)，并結(jié)合反演結(jié)果對研究區(qū)域的海底地形特征進(jìn)行分析。

1 原理與方法

1.1 海底地形反演基本原理

GGM的原理如圖1所示。海洋表面的測深點(diǎn)分布設(shè)為J1,J2, …,Jn，測深點(diǎn)對應(yīng)的船測海深為d1,d2,d3, …,dn，d為參考深度，對應(yīng)測量任務(wù)中最深點(diǎn)的測量結(jié)果，P′為最深測量點(diǎn)對應(yīng)的海底點(diǎn)。

將測深點(diǎn)上的重力異常Δg分為“長波段重力異?！焙汀岸滩ǘ沃亓Ξ惓！?個(gè)分量，即

Δg=Δglong+Δgshort。

(1)

(2)

式中：Δρ為洋殼與海水的密度之差；f為引力常數(shù)(通常取6.672×10-8cm3/(g·s2))。

圖1 GGM模型原理示意圖Fig.1 Principle of the GGM

GGM模型反演海底地形的步驟如下：

1)依據(jù)公式(2)得到測深點(diǎn)的Δgshort；

2)由Δgshort根據(jù)公式(1)得到測深點(diǎn)的Δglong；

5)根據(jù)式(2)，得到格網(wǎng)點(diǎn)的海深為

(3)

測深點(diǎn)的重力異常Δg由重力異常模型內(nèi)插得到。

1.2 向下延拓方法

如圖2所示，設(shè)海水的平均密度為ρsea，海面船測點(diǎn)I處的重力異常為Δg0，I點(diǎn)對應(yīng)平均海深面上的J點(diǎn)。將I點(diǎn)位置的重力異常Δg0向下延拓到平均海深面上，得到J點(diǎn)的重力異常ΔgJ。那么，根據(jù)文獻(xiàn)[17]中方法，Δρ可表示為

(4)

常用的向下延拓方法是頻率域內(nèi)的傅里葉變換法[19]。對海面重力異常Δg0進(jìn)行二維離散傅里葉變換：

(5)

式中：j2=-1；(x,y)為重力異常在時(shí)域的坐標(biāo)；Δg0為時(shí)域重力異常；G0為頻率域重力異常；(fx,fy)為重力異常在頻率域上的坐標(biāo)；M、N為波數(shù)[20]。通過式(5)，可將重力異常轉(zhuǎn)換到頻率域進(jìn)行處理。海面重力異常在頻率域延拓的公式為

(6)

式中：Gh為延拓面上的頻率域重力異常；Δh為延拓高度。Δh<0為向下延拓，Δh>0為向上延拓。

圖2 向下延拓法示意圖Fig.2 Downward continuation

向下延拓是誤差放大的過程，特別是高頻噪聲[21]。為去除噪聲的影響，需將Gh(fx,fy)進(jìn)行降噪處理。本文選取高斯低通濾波器FGau：

(7)

(8)

然后，對頻率域上的重力異常進(jìn)行傅里葉逆變換得到時(shí)域上的重力異常：

(9)

實(shí)際計(jì)算過程中，顧及計(jì)算的速度和效率，常采用快速傅里葉變換(FFT)方法進(jìn)行時(shí)頻轉(zhuǎn)換。

利用移去恢復(fù)技術(shù)，首先依據(jù)EGM2008重力場模型計(jì)算海面上的重力異常和延拓面上的重力異常。在實(shí)際操作過程中，由于延拓面上的重力異常涉及到地球內(nèi)部重力場，因此，考慮了海洋表面到延拓面(平均海深面)間海水層質(zhì)量的影響。海面實(shí)際重力異常減去海面模型重力異常得到海面剩余重力異常；然后，海面剩余重力異常向下延拓得到延拓面重力異常；最后，延拓面重力異常加上延拓面模型重力異常(顧及海洋表面與延拓面間海水層質(zhì)量的影響)得到最終的延拓面重力異常。由于延拓過程中涉及對重力異常進(jìn)行傅里葉變換，為消除邊緣效應(yīng)的影響，實(shí)際處理時(shí)將研究區(qū)域分別沿經(jīng)度和緯度方向向外擴(kuò)大1°范圍。而后，將結(jié)果按照研究區(qū)域大小進(jìn)行截取處理得到最終結(jié)果。

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

2.1 實(shí)驗(yàn)區(qū)域選擇和數(shù)據(jù)準(zhǔn)備

板塊之間的俯沖、后退、擠壓，使得俯沖帶上發(fā)育了大量構(gòu)造地形。在西太平洋海域，太平洋板塊自東向西俯沖于菲律賓板塊之下形成了馬里亞納海溝[22]。海溝呈現(xiàn)近南北延伸、向東弧向凸出的地形走勢。北起硫磺列島，南至雅浦島附近，全長大約2 550 km，寬70 km左右。馬里亞納海溝因其在地質(zhì)構(gòu)造、區(qū)域板塊和深度構(gòu)造研究方面的積極作用，成為學(xué)者們研究的熱點(diǎn)海域。本文選取“挑戰(zhàn)者深淵”(因英國測量船“挑戰(zhàn)者號”首次在該海域測量而得名[23])附近6°×8°(8°N—14°N、138°E—146°E)海域范圍為研究區(qū)域。使用的地形模型為2008年8月由美國地球物理數(shù)據(jù)中心(NGDC)和美國海洋和大氣管理局(NOAA)發(fā)布的分辨率為1′×1′的ETOPO1模型[24]，該模型利用陸地和海洋實(shí)測數(shù)據(jù)，并融合部分其他地形模型生成，以前的ETOPO2v2[25]和ETOPO5[26]全球地形模型現(xiàn)在基本不再使用。重力異常數(shù)據(jù)來自于丹麥科技大學(xué)(technical university of denmark)空間實(shí)驗(yàn)室(DTU space)發(fā)布的分辨率為1′×1′的DTU10重力場模型[27]。船測海深數(shù)據(jù)來源于NOAA發(fā)布的多波束海深測量數(shù)據(jù)(船測數(shù)據(jù)下載網(wǎng)址：http://www.ngdc.noaa.gov/)。

首先依據(jù)3σ法則對研究海域的實(shí)測數(shù)據(jù)進(jìn)行初步剔除，篩選得到75 110個(gè)船測點(diǎn)數(shù)據(jù)。均勻選取實(shí)測數(shù)據(jù)中14 278個(gè)數(shù)據(jù)作為檢核點(diǎn)數(shù)據(jù)，大約占船測數(shù)據(jù)總數(shù)的五分之一；余下60 832個(gè)數(shù)據(jù)作為控制點(diǎn)數(shù)據(jù)。船測控制點(diǎn)和檢核點(diǎn)分布情況如圖3所示。

黑色標(biāo)志點(diǎn)為檢核點(diǎn)；藍(lán)色標(biāo)志點(diǎn)為控制點(diǎn)；背景為ETOPO1模型。圖3 研究區(qū)船測點(diǎn)分布Fig.3 Distribution of ship datum in the study area

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

研究海域海面重力異常如圖4a所示(以DTU10模型值作為海面重力異常)。通過計(jì)算，得到研究區(qū)域平均海深為3 751.5 m，將海面重力異常采用FFT法向下延拓到平均海深面上，高斯濾波參數(shù)k取10。濾波后的重力異常經(jīng)傅里葉逆變換得到延拓面(平均海深面)上的重力異常，如圖4b所示。分別對海面重力異常和延拓面上的重力異常取均值[7]，海水密度為1.03 g/cm3，依式(4)求得Δρ=2.32 g/cm3。

2.3 反演結(jié)果比較

依據(jù)GGM模型計(jì)算流程，采用四次曲面的局部多項(xiàng)式插值法內(nèi)插得出控制點(diǎn)的重力異常，與式(2)計(jì)算的控制點(diǎn)短波段重力異常作差，使用GMT(generic mapping tools)中的張力樣條函數(shù)進(jìn)行格網(wǎng)化處理，得到1′×1′分辨率的長波段重力異常(圖5a)；根據(jù)式(1)得到1′×1′分辨率的短波段重力異常(圖5b)；由式(3)最終得到1′×1′分辨率的GGM海深模型(圖6a1)；將GGM海深模型內(nèi)插到檢核點(diǎn)，并統(tǒng)計(jì)相對誤差(模型與檢核點(diǎn)海深的差值與檢核點(diǎn)海深之比)，結(jié)果如圖6a2、圖6a3所示。同理，對ETOPO1模型(圖6b1)內(nèi)插到檢核點(diǎn)，其相對誤差統(tǒng)計(jì)結(jié)果分別如圖6b2、圖6b3所示。使用GMT中的張力樣條函數(shù)對海深控制點(diǎn)進(jìn)行網(wǎng)格化處理，生成1′×1′的格網(wǎng)化海深模型(以下稱為模型1)，結(jié)果如圖6c1所示，其相對誤差統(tǒng)計(jì)結(jié)果分別如圖6c2、圖6c3所示。

a．海面重力異常；b. 平均海深面重力異常。圖4 研究區(qū)重力異常Fig.4 Gravity anomaly in the study area

a.長波段重力異常；b.短波段重力異常。圖5 研究區(qū)重力異常分量Fig.5 Component of gravity anomaly in the study area

a1.GGM海深模型；a2.GGM海深模型相對誤差分布;a3.GGM海深模型相對誤差;b1.ETOPO1模型;b2.ETOPO1模型相對誤差分布;b3.ETOPO1模型相對誤差;c1.模型1; c2.模型1相對誤差分布;c3.模型1相對誤差。圖6 研究區(qū)海深模型與相對誤差Fig.6 Bathymetry model and relative error in the study area

將海深模型(圖6a1)與圖5對比可以看出，海深與短波段重力異常(圖5b)之間存在明顯的相關(guān)性，而與長波段重力異常(圖5a)的相關(guān)性不明顯；說明利用短波段重力異常反演海底地形在理論上具有可行性。

對比GGM海深模型、ETOPO1模型、模型1(圖6a1、b1、c1)可以看出：GGM海深模型對海底地形的反映比ETOPO1模型更加精細(xì)，可以有效填補(bǔ)無船測數(shù)據(jù)區(qū)域的海深空白；模型1在船測點(diǎn)密集的海域?qū)５椎匦蔚姆从潮容^詳細(xì)，而在船測點(diǎn)稀疏的海域?qū)５椎匦蔚某尸F(xiàn)不如GGM海深模型和ETOPO1模型詳細(xì)，如12°N—14°N、 138°E—140°E海域。

3種海深模型與檢核點(diǎn)差值的統(tǒng)計(jì)結(jié)果如表1所示。由表1看出，GGM海深模型的精度明顯優(yōu)于另外2個(gè)模型，然而，GGM海深模型與檢核點(diǎn)差值統(tǒng)計(jì)結(jié)果中的最值項(xiàng)卻超過了千米量級。對數(shù)據(jù)進(jìn)行再次檢查發(fā)現(xiàn)，出現(xiàn)該結(jié)果是因?yàn)樵趯?shí)驗(yàn)初始階段有部分粗差未被發(fā)現(xiàn)。如本次出現(xiàn)差值最大值的檢核點(diǎn)處水深為143 m，而其周圍檢核點(diǎn)水深基本在1 200 m以上，說明該處測量值明顯出現(xiàn)錯(cuò)誤，而在粗差剔除初始階段未被發(fā)現(xiàn)，進(jìn)而影響到最終的統(tǒng)計(jì)結(jié)果。將該錯(cuò)誤數(shù)據(jù)剔除后，GGM海深模型相對誤差最大值由3.053 1減小為1.647 9，相對誤差標(biāo)準(zhǔn)差由0.079 5減小為0.075 3。

表1研究區(qū)GGM海深模型、ETOPO1模型、模型1與檢核點(diǎn)差值統(tǒng)計(jì)結(jié)果

Table1StatisticalresultsofthedifferencebetweenGGMseadepthmodel,ETOPO1model,model1andcheckingpointsinthestudyarea

m

可以看出：若將帶有粗差的船測點(diǎn)作為控制點(diǎn)，會影響模型的建立；若將含有粗差的船測點(diǎn)作為檢核點(diǎn)，又將對精度統(tǒng)計(jì)結(jié)果產(chǎn)生誤導(dǎo)?；诖?，在數(shù)據(jù)處理初始階段，對船測海深粗差的剔除和發(fā)現(xiàn)顯得尤為重要。

如圖6a2、b2、c2所示，在13°N—14°N、144°E—146°E，3個(gè)模型的海深差值相對誤差普遍很大。探究其原因，對比分析不同海域的相對誤差分布情況發(fā)現(xiàn)：13°N—14°N、144°E—146°E船測控制點(diǎn)密集(圖3)，海底地形變化劇烈，相對誤差大；同樣船測控制點(diǎn)密集、海底地形變化不劇烈的11°N、142°E周圍海域的相對誤差小得多；9°N—10°N、143°E—144°E海域，船測控制點(diǎn)分布均勻，海底地形變化平緩，相對誤差也很小；船測控制點(diǎn)分布均勻，海底地形變化劇烈海域(如8°N—10°N、139°E—140°E附近)，相對誤差較大。從而，說明海底地形變化的劇烈程度是影響相對誤差的主要原因。

對比圖6a3、b3、c3發(fā)現(xiàn)，GGM海深模型的相對誤差為-0.01～0.01的占絕大多數(shù)，相對誤差為-0.05～0.05的比例占78.82%，而在ETOPO1模型和模型1中，相對誤差為-0.05～0.05所占比例分別為69.75%和77.83%；可見，GGM海深模型的相對誤差小值所占比例明顯多于另外2個(gè)模型。3種海深模型的相對誤差統(tǒng)計(jì)結(jié)果見表2。

表2研究區(qū)GGM海深模型、ETOPO1模型、模型1相對誤差統(tǒng)計(jì)結(jié)果

Table2StatisticalresultsofrelativeerrorofGGMseadepthmodel,ETOPO1modelandmodel1inthestudyarea

數(shù)據(jù)類型最大值最小值平均值標(biāo)準(zhǔn)差GGM海深模型1.647 9-0.682 20.000 20.075 3ETOPO1模型1.893 5-0.644 2-0.006 80.104 9模型1 1.598 7-0.750 3-0.001 70.077 9

進(jìn)一步比較GGM海深模型和ETOPO1模型，將GGM海深模型和ETOPO1模型作差，結(jié)果如圖7a所示。將圖7a與船測點(diǎn)分布圖(圖3)比較分析發(fā)現(xiàn)：模型差值較大區(qū)域?yàn)榇瑴y點(diǎn)匱乏或者過于稀疏的海域，如12°N—14°N、138°E—139°E海域。同時(shí)，在多海溝、多海山區(qū)域模型差值也普遍較大；另外，船測點(diǎn)分布的密集程度對模型差值的影響顯著，反映了GGM海深模型和ETOPO1模型在模型建立時(shí)源數(shù)據(jù)的差異，GGM海深模型和ETOPO1模型的測深基準(zhǔn)面差異對結(jié)果也會產(chǎn)生影響，但不是主要影響。GGM海深模型與船測點(diǎn)格網(wǎng)化模型的差值受船測點(diǎn)的分布情況更加明顯(圖7b)。3個(gè)海深模型的差值統(tǒng)計(jì)結(jié)果見表3。

表3研究區(qū)GGM海深模型與ETOPO1模型、模型1比較

Table3ComparisonsbetweenGGMseadepthmodelandETOPO1model,model1inthestudyarea

數(shù)據(jù)類型海深/m最大值最小值平均值標(biāo)準(zhǔn)差GGM海深模型-ETOPO1模型2 454.20-4 443.8033.85418.31GGM海深模型-模型13 991.70-5 321.6079.36488.59

圖7和表3結(jié)果顯示，GGM海深模型與ETOPO1更為接近。那么建立綜合模型時(shí)，應(yīng)該考慮船測點(diǎn)的分布情況；船測點(diǎn)分布均勻的海域，用GGM海深模型；在船測點(diǎn)匱乏海域，考慮到ETOPO1模型源數(shù)據(jù)的充分性和來源多樣性，選擇ETOPO1模型。

2.4 海底地形分析

依據(jù)GGM反演的海深模型結(jié)果如圖8a所示，從圖中能清晰地看出本文實(shí)驗(yàn)的馬里亞納海溝海域?yàn)闁|西走勢。海溝的平均深度大于8 000 m。在11°N、142°E附近出現(xiàn)10 000 m以上的深度值。“挑戰(zhàn)者深淵”就位于該海域，選擇最深點(diǎn)周圍海域進(jìn)行分析(11°00′N—11°36′N、141°36′E—143°00′E)，海深結(jié)果如圖8b所示。

從圖8b可以看出，大于10 000 m的等深線圈東西兩端較窄，中間寬，形似紡錘。在萬米等深線圈內(nèi)能清楚地看見3個(gè)洼地，從東到西面積逐漸減小。東部洼地呈長條狀。洼地的中軸方向和海溝的中軸方向基本一致。但是，洼地的細(xì)部輪廓特征描述不夠細(xì)致。

a.GGM海深模型-ETOPO1模型；b.GGM海深模型-模型1。圖7 研究區(qū)模型較差Fig.7 Difference of the model in the study area

a. 依據(jù)GGM反演的馬里亞納海溝海深模型；b.“挑戰(zhàn)者深淵”附近海底深度。圖8 馬里亞納海溝Fig.8 Mariana trench

在圖8a中直線所在位置作海溝的剖面圖，結(jié)果如圖9所示。從圖9中可看到，4 000 m深的海溝南北兩側(cè)不同程度地發(fā)育了海山。海溝南側(cè)在海深小于5 000 m的淺區(qū)坡度較平緩，為2°～5°；海深大于5 000 m以下坡度明顯增大，為10°～15°，直抵海溝底部。海溝北測在整個(gè)下降階段，坡度很大，為10°～15°；在北緯11°45′附近出現(xiàn)一個(gè)緩沖地帶。南北兩側(cè)不對稱。

圖9 海溝剖面圖Fig.9 Topography of the trench

依據(jù)反演的海底地形，能夠大概看出海溝的基本形態(tài)。但是，對于海溝南北坡的細(xì)節(jié)反映不夠明顯。從而，依據(jù)模型反演的海底地形只能反映某種大范圍的海底地形，對海底地形細(xì)節(jié)的反映不夠清晰。

3 結(jié)論

文章選取西太平洋板塊俯沖菲律賓板塊而形成的馬里亞納海溝所在海域作為實(shí)驗(yàn)區(qū)域，分析比較 了GGM海深模型、ETOPO1模型以及直接將船測點(diǎn)海深數(shù)據(jù)格網(wǎng)化3種海深模型之間的差異，并對所選海域最深點(diǎn)附近的“挑戰(zhàn)者深淵”兩側(cè)地貌進(jìn)行了研究，得到結(jié)論如下：

1)GGM利用衛(wèi)星測高重力異常，可以高分辨率地反演海底地形，結(jié)果優(yōu)于ETOPO1模型，優(yōu)于按本文格網(wǎng)化方法將船測點(diǎn)直接格網(wǎng)化的結(jié)果，為大面積地研究地質(zhì)構(gòu)造、區(qū)域板塊和深度構(gòu)造有積極作用。

2)相對誤差與海深的關(guān)系不大，受海底地形變化影響劇烈。海底地形變化平緩，船測點(diǎn)分布均勻，相對誤差小。

3)GGM海深模型和ETOPO1模型的差值，以及GGM海深模型與船測點(diǎn)海深格網(wǎng)化的差值受船測點(diǎn)的分布影響明顯。船測點(diǎn)匱乏或者過于稀疏的海域模型差值較大。

4)繪制“挑戰(zhàn)者深淵”兩側(cè)南北坡剖面圖后發(fā)現(xiàn)，海溝兩側(cè)的地貌有明顯差異。在海溝南側(cè)海深小于5 000 m的淺海部分坡度平緩(2°～5°)，而在水深大于5 000 m部分坡度明顯增大(10°～15°)；海溝北側(cè)在整個(gè)下降階段坡度很大(10°～15°)，11°45′ N附近出現(xiàn)一個(gè)緩沖地帶。依據(jù)模型反演的海底地形，只能大范圍地反演某些海底環(huán)境，對海底的細(xì)節(jié)描述依然不夠清晰。

向下延拓理論上本應(yīng)在整個(gè)球面上進(jìn)行，文章在向下延拓過程中，直接在局部范圍內(nèi)將海洋表面的重力異常向下延拓，同時(shí)，忽略了海水質(zhì)量的影響，對最終結(jié)果可能有影響，還需進(jìn)一步論證。