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

范 雕，李姍姍，孟書宇，邢志斌，馮進凱

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

0 引言

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

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

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

1 原理與方法

1.1 海底地形反演基本原理

GGM的原理如圖1所示。海洋表面的測深點分布設為J1,J2, …,Jn，測深點對應的船測海深為d1,d2,d3, …,dn，d為參考深度，對應測量任務中最深點的測量結果，P′為最深測量點對應的海底點。

將測深點上的重力異常Δg分為“長波段重力異常”和“短波段重力異常”2個分量，即

Δg=Δglong+Δgshort。

(1)

(2)

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

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

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

1)依據公式(2)得到測深點的Δgshort；

2)由Δgshort根據公式(1)得到測深點的Δglong；

5)根據式(2)，得到格網點的海深為

(3)

測深點的重力異常Δg由重力異常模型內插得到。

1.2 向下延拓方法

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

(4)

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

(5)

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

(6)

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

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

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

(7)

(8)

然后，對頻率域上的重力異常進行傅里葉逆變換得到時域上的重力異常：

(9)

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

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

2 實驗結果與分析

2.1 實驗區域選擇和數據準備

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

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

黑色標志點為檢核點；藍色標志點為控制點；背景為ETOPO1模型。圖3 研究區船測點分布Fig.3 Distribution of ship datum in the study area

2.2 密度差異常數確定

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

2.3 反演結果比較

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

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

a.長波段重力異常；b.短波段重力異常。圖5 研究區重力異常分量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 研究區海深模型與相對誤差Fig.6 Bathymetry model and relative error in the study area

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

對比GGM海深模型、ETOPO1模型、模型1(圖6a1、b1、c1)可以看出：GGM海深模型對海底地形的反映比ETOPO1模型更加精細，可以有效填補無船測數據區域的海深空白；模型1在船測點密集的海域對海底地形的反映比較詳細，而在船測點稀疏的海域對海底地形的呈現不如GGM海深模型和ETOPO1模型詳細，如12°N—14°N、 138°E—140°E海域。

3種海深模型與檢核點差值的統計結果如表1所示。由表1看出，GGM海深模型的精度明顯優于另外2個模型，然而，GGM海深模型與檢核點差值統計結果中的最值項卻超過了千米量級。對數據進行再次檢查發現，出現該結果是因為在實驗初始階段有部分粗差未被發現。如本次出現差值最大值的檢核點處水深為143 m，而其周圍檢核點水深基本在1 200 m以上，說明該處測量值明顯出現錯誤，而在粗差剔除初始階段未被發現，進而影響到最終的統計結果。將該錯誤數據剔除后，GGM海深模型相對誤差最大值由3.053 1減小為1.647 9，相對誤差標準差由0.079 5減小為0.075 3。

表1研究區GGM海深模型、ETOPO1模型、模型1與檢核點差值統計結果

Table1StatisticalresultsofthedifferencebetweenGGMseadepthmodel,ETOPO1model,model1andcheckingpointsinthestudyarea

m

可以看出：若將帶有粗差的船測點作為控制點，會影響模型的建立；若將含有粗差的船測點作為檢核點，又將對精度統計結果產生誤導。基于此，在數據處理初始階段，對船測海深粗差的剔除和發現顯得尤為重要。

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

對比圖6a3、b3、c3發現，GGM海深模型的相對誤差為-0.01～0.01的占絕大多數，相對誤差為-0.05～0.05的比例占78.82%，而在ETOPO1模型和模型1中，相對誤差為-0.05～0.05所占比例分別為69.75%和77.83%；可見，GGM海深模型的相對誤差小值所占比例明顯多于另外2個模型。3種海深模型的相對誤差統計結果見表2。

表2研究區GGM海深模型、ETOPO1模型、模型1相對誤差統計結果

Table2StatisticalresultsofrelativeerrorofGGMseadepthmodel,ETOPO1modelandmodel1inthestudyarea

數據類型最大值最小值平均值標準差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

進一步比較GGM海深模型和ETOPO1模型，將GGM海深模型和ETOPO1模型作差，結果如圖7a所示。將圖7a與船測點分布圖(圖3)比較分析發現：模型差值較大區域為船測點匱乏或者過于稀疏的海域，如12°N—14°N、138°E—139°E海域。同時，在多海溝、多海山區域模型差值也普遍較大；另外，船測點分布的密集程度對模型差值的影響顯著，反映了GGM海深模型和ETOPO1模型在模型建立時源數據的差異，GGM海深模型和ETOPO1模型的測深基準面差異對結果也會產生影響，但不是主要影響。GGM海深模型與船測點格網化模型的差值受船測點的分布情況更加明顯(圖7b)。3個海深模型的差值統計結果見表3。

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

Table3ComparisonsbetweenGGMseadepthmodelandETOPO1model,model1inthestudyarea

數據類型海深/m最大值最小值平均值標準差GGM海深模型-ETOPO1模型2 454.20-4 443.8033.85418.31GGM海深模型-模型13 991.70-5 321.6079.36488.59

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

2.4 海底地形分析

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

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

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

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

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

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

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

3 結論

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

1)GGM利用衛星測高重力異常，可以高分辨率地反演海底地形，結果優于ETOPO1模型，優于按本文格網化方法將船測點直接格網化的結果，為大面積地研究地質構造、區域板塊和深度構造有積極作用。

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

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

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

向下延拓理論上本應在整個球面上進行，文章在向下延拓過程中，直接在局部范圍內將海洋表面的重力異常向下延拓，同時，忽略了海水質量的影響，對最終結果可能有影響，還需進一步論證。