基于球面重力反演蘇拉威西地區(qū)莫霍面結(jié)構(gòu)

褚偉， 徐亞*， 郝天珧,4

1 中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049 2 中國(guó)科學(xué)院地質(zhì)與地球物理研究所, 中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室， 北京 100029 3 中國(guó)科學(xué)院地球科學(xué)研究院， 北京 100029 4 海底科學(xué)重點(diǎn)實(shí)驗(yàn)室， 自然資源部第二海洋研究所， 杭州 310012

0 引言

在20世紀(jì)早期，莫霍面深度及結(jié)構(gòu)特征主要是通過(guò)天然地震的波速變化來(lái)探測(cè)，20世紀(jì)40年代之后，炸藥震源的使用提供了更多高分辨率的結(jié)果，1980年之后，使用基于地震的接收函數(shù)方法研究地殼結(jié)構(gòu)逐漸流行.地震波折射、接收函數(shù)和面波成像等方法為現(xiàn)今我們認(rèn)識(shí)地殼結(jié)構(gòu)和莫霍面深度做出了重要貢獻(xiàn).然而，在地震數(shù)據(jù)覆蓋較少的區(qū)域，使用重力異常數(shù)據(jù)確定莫霍面深度也是一個(gè)較為可靠的方法.Teng等(2013)探討中國(guó)地區(qū)莫霍面深度時(shí)，發(fā)現(xiàn)使用主動(dòng)源地震確定的莫霍面深度與重力反演或者面波成像的結(jié)果基本一致，差距在3～5 km之內(nèi)；Van der Meijde等(2013)使用衛(wèi)星重力數(shù)據(jù)確定南美地區(qū)的莫霍面深度時(shí)，發(fā)現(xiàn)除了活動(dòng)造山帶，重力反演結(jié)果與地震觀測(cè)基本一致；Aitken等(2013)在反演澳大利亞地區(qū)莫霍面深度時(shí)，測(cè)試了存在地震約束時(shí)重力反演的有效性，表明研究區(qū)的地震數(shù)據(jù)覆蓋中等且均勻時(shí)重力反演效果良好.

蘇拉威西地區(qū)位于太平洋板塊、菲律賓海板塊、歐亞板塊和印度洋—澳大利亞板塊的匯聚地帶，包括多條海溝和斷裂帶，地質(zhì)背景與地質(zhì)演化過(guò)程復(fù)雜.研究區(qū)大部分位于海域，目前地球物理觀測(cè)尚不充分，因此使用衛(wèi)星重力觀測(cè)數(shù)據(jù)反演莫霍面深度是目前較為可行的方法.反演得到的莫霍面結(jié)構(gòu)可為該區(qū)域的三維密度結(jié)構(gòu)反演提供約束，進(jìn)而為動(dòng)力學(xué)模擬提供邊界條件，最終為認(rèn)識(shí)該區(qū)域俯沖轉(zhuǎn)換邊緣、對(duì)向俯沖系統(tǒng)以及俯沖初始機(jī)制等獨(dú)特地質(zhì)現(xiàn)象提供依據(jù).

重力方法反演莫霍面深度時(shí)，對(duì)于形狀不規(guī)則的地下介質(zhì)，可將其分割為不同的幾何單元分別求取異常，再相加得到總異常.在球坐標(biāo)系中，可以按照經(jīng)緯度網(wǎng)格剖分模型(Li et al.，2011；Uieda et al.，2016)，基于球坐標(biāo)系的模型剖分與重力正演，能夠在區(qū)域研究中避免地球曲率帶來(lái)的重力正演誤差.本文使用球坐標(biāo)系下的重力方法反演研究區(qū)的莫霍面深度.首先結(jié)合頻譜分析及匹配濾波方法分析布格重力異常，并獲取反映莫霍面結(jié)構(gòu)的區(qū)域重力異常，同時(shí)利用頻譜分析給出莫霍面的參考深度，使用兩次隨機(jī)子抽樣交叉驗(yàn)證確定反演超參數(shù)，最終獲取合理的莫霍面深度.

1 地質(zhì)與地球物理概況

1.1 地質(zhì)與構(gòu)造特征

研究區(qū)位于太平洋板塊、菲律賓海板塊、歐亞板塊和印度洋—澳大利亞板塊的匯聚地帶，北部為蘇拉威西海，南部為蘇拉威西島，西側(cè)為加里曼丹島，東側(cè)為馬魯古海.該區(qū)域有復(fù)雜的俯沖和斷裂體系，包括由蘇拉威西海向蘇拉威西島之下俯沖形成的北蘇拉威西海溝，背向俯沖的桑義赫海溝和哈馬黑拉海溝，以及Palu-Koro斷裂帶，Sula-Sorong斷裂帶，Mantano斷裂帶等(圖1).

圖1 蘇拉威西地區(qū)區(qū)域構(gòu)造簡(jiǎn)圖 部分構(gòu)造單元下文簡(jiǎn)稱(chēng)：科托巴托海溝(CTB)，桑義赫海溝(SGH)，北蘇拉威西海溝(NSW)，Palu-Koro斷裂帶(PK)，Sula-Sorong斷裂帶(SS)，古俯沖帶(PSZ)，Mantano斷裂帶(MTN)，UNA-UNA火山(UU)，俯沖轉(zhuǎn)換傳遞邊緣(TSF)，哈馬黑拉海溝(HMH)，馬魯古海(MLC).Fig.1 Map of the regional structure in Sulawesi areaSome tectonic units hereinafter referred to as：Cotabato trench (CTB), Sangihe trench (SGH), North Sulawesi trench (NSW), Palu-Koro fault (PK), Sula-Sorong fault (SS), Palaeo-subduction zone (PSZ), Mantano fault (MTN), UNA-UNA volcano (UU), Transition zone of subduction and fracture (TSF), Halmahera trench (HMH), Molucca Sea (MLC).

根據(jù)印度洋和印度尼西亞地區(qū)的古板塊重建研究(Hall，2012)，蘇拉威西海自始新世(約47 Ma)開(kāi)始擴(kuò)張，中始新世(約40 Ma)停止擴(kuò)張；中中新世(約15 Ma)蘇拉威西海開(kāi)始向蘇祿群島的西北邊緣俯沖，同時(shí)蘇祿海開(kāi)始弧后海底擴(kuò)張過(guò)程，晚中新世至早上新世(約12～5 Ma)，蘇拉威西—蘇祿俯沖帶邊俯沖、邊后撤，直至早上新世(5 Ma)俯沖帶停止俯沖，蘇祿海也停止了弧后海底擴(kuò)張的過(guò)程.這時(shí)，蘇拉威西海板塊開(kāi)始向南俯沖到蘇拉威西島北支下面，形成了北蘇拉威西海溝，同時(shí)發(fā)生了俯沖后撤.由于蘇拉威西海板塊的俯沖作用，蘇拉威西島北支以其東端為軸，開(kāi)始以4°/Ma的角速度順時(shí)針?lè)较蛐D(zhuǎn)(Zhang et al., 2021)，導(dǎo)致該區(qū)古近紀(jì)斷裂重新活動(dòng)并發(fā)展成為了左旋走滑的Palu-Koro斷裂帶.

1.2 莫霍面研究現(xiàn)狀

目前蘇拉威西地區(qū)的莫霍面深度研究資料以及地球物理研究成果較少，隨著近年來(lái)對(duì)該地區(qū)的地球物理調(diào)查逐漸增多，不同學(xué)者陸續(xù)從地震學(xué)、重力、磁法、地?zé)釋W(xué)等方面開(kāi)展對(duì)該區(qū)域深部結(jié)構(gòu)的研究.

全球模型是對(duì)區(qū)域開(kāi)展整體認(rèn)識(shí)的重要依據(jù)，Crust1.0全球模型(Laske et al., 2013)給出了研究區(qū)范圍1°×1°的莫霍面深度(圖2a)，莫霍面深度在10～33 km之間，平均深度為25 km.全球?qū)游龀上衲Ｐ涂梢杂糜谘芯繀^(qū)深部地幔結(jié)構(gòu)的研究，目前有較多的全球?qū)游龀上衲Ｐ?，例如MIT-P08(Li et al., 2008)、LLNL-G3Dv3(Simmons et al., 2012)，F(xiàn)ukao和Obayashi(2013)等給出的多個(gè)P波層析模型.此外，在東南亞和西太平洋地區(qū)，區(qū)域地震層析成像研究也給出了研究區(qū)地幔的三維成像結(jié)果(Hall and Spakman, 2015; Huang et al., 2015; Zenono et al., 2019)，但是層析成像結(jié)果在淺層的分辨率較低，無(wú)法為本文莫霍面深度反演提供有效的約束.

圖2 蘇拉威西地區(qū)莫霍面深度(a) 可用于研究區(qū)莫霍面深度對(duì)比的數(shù)據(jù)，其中藍(lán)色圓點(diǎn)為天然地震數(shù)據(jù)位置，紅色實(shí)線是雙船折射地震測(cè)線位置，灰色虛線是Crust1.0全球模型的莫霍面深度等值線; (b) 研究區(qū)莫霍面深度重力反演結(jié)果，數(shù)據(jù)來(lái)自郝天珧等(2014).Fig.2 Moho depth in Sulawesi area(a) The data can be used to compare the Moho depth in study area, among them the blue dot represents the natural earthquake data, the red solid lines represent the two-ship seismic refraction lines, and the gray dotted line is the Moho depth contour of the Crust1.0 global model. (b) Moho depth from gravity inversion in study area, data from Hao et al. (2014).

研究區(qū)雖有若干地震測(cè)線數(shù)據(jù)(Murauchi et al., 1973; Kopp et al., 1999)，但限于這些地震測(cè)線時(shí)間久遠(yuǎn)，長(zhǎng)度較短，深部界面不清晰，無(wú)法從整體上認(rèn)識(shí)研究區(qū)莫霍面的分布特征.Murauchi等(1973)發(fā)布的6條雙船折射地震剖面數(shù)據(jù)中，測(cè)線18、19、20位于蘇拉威西海，測(cè)線位置如圖2a所示，數(shù)據(jù)質(zhì)量相對(duì)較好，可以輔助驗(yàn)證本文的莫霍面反演結(jié)果.近年來(lái)也有學(xué)者在研究區(qū)開(kāi)展三維地震觀測(cè)(呂川川等，2019)以及天然地震研究，但研究結(jié)果尚未發(fā)表.

研究區(qū)的高精度地震數(shù)據(jù)較少，因此使用重力數(shù)據(jù)反演莫霍面深度是目前較可靠的方法.郝天珧等(2014)編繪的中國(guó)海陸1∶500萬(wàn)莫霍面深度圖利用重力數(shù)據(jù)和Parker界面反演法獲得了研究區(qū)北部的莫霍面結(jié)構(gòu)，如圖2b所示，給出的莫霍面深度在6.7～34.3 km之間，平均深度為19.8 km.Zhang等(2021)利用重力異常遺傳-有限單元方法反演莫霍面密度差，獲取莫霍面深度，反演結(jié)果在5.9～25.5 km之間，平均深度為18.4 km.

以上全球模型、地震數(shù)據(jù)以及莫霍面深度的重力反演結(jié)果可為本研究提供對(duì)比參考.

1.3 重力異常特征

根據(jù)全球衛(wèi)星測(cè)高重力異常(Sandwell et al., 2014)可以獲得蘇拉威西地區(qū)空間重力異常，如圖3a所示，數(shù)據(jù)精度為1′.該區(qū)域空間重力異常總體變化范圍大概為-250～400 mGal，海域以低值異常為主，陸地區(qū)域以中高異常為主.在海陸過(guò)渡帶出現(xiàn)明顯的重力異常正負(fù)交替，例如桑義赫海溝和Sula-Sorong斷裂帶處為明顯的負(fù)異常，在Sulawesi島北支為正異常.

圖3 蘇拉威西地區(qū)重力異常圖(a) 空間重力異常； (b) 布格重力異常.Fig.3 Map of the gravity anomaly in Sulawesi area(a) Free air gravity anomaly; (b) Bouguer gravity anomaly.

利用全球地形數(shù)據(jù)(Smith and Sandwell, 1997)對(duì)該地區(qū)空間重力異常進(jìn)行布格改正(Fullea et al., 2008)，改正半徑取166.7 km，地形改正密度取2.67 g·cm-3，海水密度取1.03 g·cm-3.經(jīng)過(guò)上述改正獲得研究區(qū)布格重力異常，如圖3b所示，其總體變化范圍大概為-100～400 mGal，在蘇拉威西海呈現(xiàn)明顯的高異常圈閉，在陸區(qū)如蘇拉威西島和加里曼丹島為明顯的低異常.

對(duì)比空間重力異常及布格重力異常，可見(jiàn)在Sula-Sorong斷裂帶北側(cè)以及馬魯古海的位置均呈現(xiàn)明顯的低值異常區(qū)，由于馬魯古海兩側(cè)是背向俯沖的桑義赫海溝與哈馬黑拉海溝，南側(cè)是走滑斷裂帶，因此布格重力異常低值可能與該區(qū)域復(fù)雜的地質(zhì)構(gòu)造背景有關(guān).

1.4 重力異常場(chǎng)分離

重力異常是地下所有地質(zhì)體密度差異引起重力響應(yīng)的疊加場(chǎng)，通常淺層或小規(guī)模密度異常體與深層或大規(guī)模密度異常體引起的異常具有明顯的頻譜差異，通過(guò)頻譜分析可獲得對(duì)地下密度異常體分布的總體認(rèn)識(shí)，并據(jù)此進(jìn)行深-淺場(chǎng)源重力異常的分離.

莫霍面是地下具有明顯密度差異的全局性界面，其重力異常通常對(duì)應(yīng)頻譜中的低頻部分.在反演莫霍面之前，需要分離出反映莫霍面特征的重力異常響應(yīng).本文結(jié)合頻譜分析及匹配濾波方法，分析布格重力異常并提取出反映莫霍面結(jié)構(gòu)的區(qū)域重力異常.如圖4所示，蘇拉威西地區(qū)布格重力異常的頻譜明顯具有兩段性，根據(jù)頻譜計(jì)算平均地質(zhì)體埋深的基本方法，可以計(jì)算得到對(duì)應(yīng)頻譜低頻段的構(gòu)造面平均深度為29.6 km，可基本代表該地區(qū)莫霍面的平均深度，并作為后續(xù)確定莫霍面反演參數(shù)的主要依據(jù).

圖4 蘇拉威西地區(qū)布格重力異常頻譜圖Fig.4 Spectrum of the Bouguer gravity anomaly in Sulawesi area

結(jié)合頻譜分析結(jié)果，利用匹配濾波方法(王家林等，1991)提取布格重力異常頻譜的低頻成分，獲得研究區(qū)區(qū)域重力異常(圖5).區(qū)域重力異?？傮w變化較為平緩，與該地區(qū)總體的海陸構(gòu)造具有較為明顯的鏡像特征，在海區(qū)以區(qū)域性的高值異常為主，在陸區(qū)以低值異常為主.

2 莫霍面反演方法

2.1 球面重力正演方法

正演是反演工作的基礎(chǔ)，在大尺度構(gòu)造區(qū)必須考慮地球曲率的影響，因此本文基于球坐標(biāo)系對(duì)地下模型空間進(jìn)行剖分和計(jì)算，避免地球曲率在傳統(tǒng)重力異常正演時(shí)的近似和改正.重力正演時(shí)，使用自適應(yīng)離散算法將模型剖分為球棱柱(Li et al.，2011).球棱柱定義在地心球坐標(biāo)系中，由兩條經(jīng)線、兩條緯線和兩個(gè)同心圓內(nèi)面組成，如圖6所示.

圖6 球棱柱示意圖Q為積分點(diǎn)，在球棱柱內(nèi)部，位于地心坐標(biāo)系(X,Y,Z).P為計(jì)算點(diǎn)，位于局部坐標(biāo)系(x,y,z)，r,φ,λ分別是P點(diǎn)的半徑、緯度和經(jīng)度.l為P和Q點(diǎn)在笛卡爾坐標(biāo)系中的距離(據(jù)Uieda et al., 2016).Fig.6 View of a tesseroidThe integration point Q inside the tesseroid, a geocentric coordinate system (X, Y, Z), the computation P and it′s local coordinate system (x,y,z); r,φ,λ are the radius, latitude, and longitude of point P, respectively, and l is the Cartesian distance between P and Q (After Uieda et al.,2016).

位于球棱柱外的任意觀測(cè)點(diǎn)在球坐標(biāo)系中的重力位和重力加速度的體積分公式如下：

(1)

其中α∈{x,y,z}，ρ為密度，重力常數(shù)G=6.674×10-11m3/(kg·s)，

(2)

計(jì)算中依據(jù)觀測(cè)點(diǎn)到源點(diǎn)的距離自適應(yīng)劃分網(wǎng)格間距，采用高斯-勒讓德數(shù)值積分公式進(jìn)行重力正演計(jì)算，從而在滿足計(jì)算精度的條件下提高計(jì)算速度(Li et al., 2011).

2.2 莫霍面反演方法

莫霍面反演基于區(qū)域重力異常，該重力異常消除了地殼內(nèi)部其他構(gòu)造及密度體分布的影響，根據(jù)數(shù)據(jù)誤差和正則化項(xiàng)構(gòu)建反演目標(biāo)函數(shù)Γ(p)(Uieda and Barbosa，2017)，

Γ(p)=[do-d(p)]T[do-d(p)]

+μpTRTRp,

(3)

其中p是莫霍面深度的M維參數(shù)向量，do為N維重力異常觀測(cè)數(shù)據(jù)，d(p)為重力異常預(yù)測(cè)數(shù)據(jù)，μ是正則化參數(shù)，R是L×M的有限差分矩陣.

反演過(guò)程中，采用高斯-牛頓迭代法最小化目標(biāo)函數(shù)，其中求解參數(shù)擾動(dòng)向量Δp0非常重要.擾動(dòng)向量用于更新參數(shù)向量：p1=p0+Δp0，其中p0為初值.迭代過(guò)程中不斷重復(fù)更新參數(shù)向量，直到目標(biāo)函數(shù)Γ(p)達(dá)到最小值.對(duì)于高斯-牛頓迭代方法，第k次迭代的參數(shù)擾動(dòng)向量Δpk通過(guò)求解線性方程組得到：

(4)

(5)

Hk=2AkTAk+2μRTR,

(6)

[AkTAk+μRTR]Δpk=AkT[do-d(pk)]-μRTRpk.

(7)

由于雅可比矩陣以及方程組的計(jì)算和存儲(chǔ)，使得上述方法計(jì)算量很大.為了解決這一問(wèn)題，參考Bott (1960)的方法，用布格板估計(jì)代替雅可比矩陣(Uieda and Barbosa，2017)，新的近似矩陣定義為：

A=2πGΔρI,

(8)

其中I是單位陣.這種方法不需要計(jì)算和存儲(chǔ)雅可比矩陣，可以顯著提高計(jì)算速度.

反演過(guò)程中影響反演結(jié)果但是不能直接估計(jì)的參數(shù)稱(chēng)超參數(shù)，包括正則化因子、莫霍面密度差、參考莫霍面深度，確定這3個(gè)參數(shù)對(duì)于反演結(jié)果非常重要.Uieda和 Barbosa(2017)提出使用隨機(jī)子抽樣交叉驗(yàn)證估計(jì)正則化因子，根據(jù)地震約束點(diǎn)的深度信息估計(jì)密度差和參考深度，最終反演得到莫霍面深度.蘇拉威西地區(qū)地震觀測(cè)資料少，缺少莫霍面深度的可靠約束，因此本文使用頻譜分析給出莫霍面參考深度，采用兩次隨機(jī)子抽樣交叉驗(yàn)證分別確定正則化參數(shù)和最優(yōu)的密度差，反演的總體技術(shù)路線如圖7所示.

圖7 反演算法流程圖Fig.7 Flow chart of inversion algorithm

對(duì)最優(yōu)化參數(shù)選取的交叉驗(yàn)證實(shí)現(xiàn)過(guò)程如下：

圖8 交叉驗(yàn)證的訓(xùn)練集(空心圓)和檢驗(yàn)集(實(shí)心圓)Fig.8 Sketch of a data grid separated into the training (open circles) and testing (black dots) data sets in hold-out cross-validation

第二次隨機(jī)子抽樣交叉驗(yàn)證方法確定最優(yōu)的密度差Δρ，具體過(guò)程為：

1)給定密度差Δρ的變化范圍，用zref、Δρ的不同組合對(duì)訓(xùn)練集做反演，得到莫霍面深度pl,m；

4)對(duì)應(yīng)的均方根誤差最小的zref與Δρ為最優(yōu)超參數(shù)，對(duì)應(yīng)的莫霍面深度pl,m為最終反演結(jié)果.

3 蘇拉威西地區(qū)莫霍面反演與分析

3.1 反演參數(shù)選取

反演莫霍面時(shí)需要選擇合理的超參數(shù)，以獲得良好的重力異常擬合結(jié)果和合理的莫霍面深度.在蘇拉威西地區(qū)莫霍面反演中，為分析不同參數(shù)對(duì)反演結(jié)果的影響，選取了正則化因子、莫霍面密度差、莫霍面平均參考深度這3個(gè)超參數(shù)的不同組合進(jìn)行計(jì)算分析，參數(shù)測(cè)試組合見(jiàn)表1.在用匹配濾波方法分離重力異常時(shí)，根據(jù)頻譜分析給出莫霍面平均參考深度為29.6 km，因此莫霍面平均參考深度選擇在29 km附近.將參考莫霍面深度與正則化因子、莫霍面密度差進(jìn)行不同的組合，以測(cè)試方法的穩(wěn)定性及反演結(jié)果的合理性.反演過(guò)程中給定的正則化因子的變化范圍為10-10～10-2，莫霍面密度差變化范圍為0.3～0.5 g·cm-3，步長(zhǎng)為0.02 g·cm-3.結(jié)合不同的參考莫霍面深度，測(cè)試分為四組進(jìn)行(表1(a))，每組獲取的最優(yōu)超參數(shù)組合見(jiàn)表1(b).

表1 不同超參數(shù)組合測(cè)試Table 1 Tests of different combination of hyperparameters

可以看到，反演算法對(duì)參考莫霍面深度敏感性很強(qiáng).進(jìn)一步比較上述不同超參數(shù)組合的反演結(jié)果，統(tǒng)計(jì)重力異常擬合數(shù)據(jù)、重力異常擬合數(shù)據(jù)誤差、反演莫霍面深度的平均值及變化幅度，統(tǒng)計(jì)結(jié)果如表2所示.重力異常擬合數(shù)據(jù)誤差的直方圖如圖9所示.綜合重力異常的擬合程度等統(tǒng)計(jì)數(shù)據(jù)，認(rèn)為T(mén)est 2的反演結(jié)果最合理，Test 2反演得到的重力異常擬合數(shù)據(jù)總誤差在±15 mGal范圍內(nèi)，平均值為0.04 mGal，標(biāo)準(zhǔn)差為1.80 mGal，能夠滿足反演精度，因此將該結(jié)果作為本文最終的反演結(jié)果.

表2 反演結(jié)果統(tǒng)計(jì)分析Table 2 Statistical analysis of the inversion results

3.2 反演結(jié)果分析

將最終反演得到的莫霍面深度與Crust1.0全球模型、前人的重力反演研究以及收集到的地震數(shù)據(jù)作比較，驗(yàn)證本文結(jié)果的可靠性.

圖9 重力異常擬合數(shù)據(jù)誤差的直方圖Fig.9 Histogram of fitted gravity anomaly residuals

3.2.1 反演結(jié)果與Crust1.0全球模型對(duì)比

將本文莫霍面反演結(jié)果減去Crust1.0全球模型給出的莫霍面深度，繪制差值的平面圖與直方圖，如圖10所示.由于Crust1.0的精度較低，因此只做整體性參考.Crust1.0的莫霍面平均深度約25 km，本文反演結(jié)果平均深度約20 km，從圖10b也可以看到，反演得到的莫霍面深度整體比Crust1.0模型淺大約5 km，存在一個(gè)整體偏差.而在馬魯古海地區(qū)，顯示本文反演得到的莫霍面則較深.

圖10 本文莫霍面深度反演結(jié)果與Crust1.0全球模型對(duì)比(a) 莫霍面深度之差； (b) 莫霍面深度之差直方圖.Fig.10 Comparison of the Moho depth in this paper and Crust1.0 model(a) The Moho depth difference; (b) Histogram of the difference of Moho depth.

3.2.2 反演結(jié)果與前人重力反演結(jié)果對(duì)比

將本文反演的莫霍面深度結(jié)果減去郝天珧等(2014)的重力反演結(jié)果，繪制差值的平面圖與直方圖，如圖11所示.大部分地區(qū)的莫霍面之差在-2～3 km之間，最大的差異出現(xiàn)在哈馬黑拉海溝東側(cè)，顯示本文反演結(jié)果偏小，而在馬魯古海北部，顯示本文反演莫霍面較深.

圖11 本文莫霍面深度反演結(jié)果與郝天珧等(2014)對(duì)比(a) 莫霍面深度之差； (b) 莫霍面深度之差直方圖.Fig.11 Comparison of the Moho depth in this paper and Hao et al. (2014)(a) The Moho depth difference; (b) Histogram of the difference of Moho depth.

3.2.3 反演結(jié)果與地震數(shù)據(jù)對(duì)比

此外，本文還收集到了3條雙船折射地震剖面數(shù)據(jù)，剖面位置如圖2a中紅色實(shí)線所示.這3條剖面給出的莫霍面深度如圖12中紅色實(shí)線所示，藍(lán)色虛線為本文給出的對(duì)應(yīng)位置的莫霍面深度反演結(jié)果，綠色粗虛線為郝天珧等(2014)給出的對(duì)應(yīng)位置的莫霍面深度反演結(jié)果.可以看到雙船折射剖面18、19、20給出的莫霍面深度范圍分別為12.8～13.6 km、 12.0～13.6 km、15.0～15.3 km.郝天珧等(2014)的結(jié)果整體淺于雙船折射剖面結(jié)果，本文的反演結(jié)果則介于兩者之間.由圖2a可知，這3條地震剖面位于蘇拉威西海，海域的莫霍面深度在10 km附近是一個(gè)較為合理的值.雙船折射地震數(shù)據(jù)來(lái)自于20世紀(jì)70年代，因此數(shù)據(jù)采集精度可能是差異較大的主要原因.

圖12 本文莫霍面深度反演結(jié)果與雙船折射地震數(shù)據(jù)對(duì)比Fig.12 Comparison of the Moho depth in this paper and two-ship seismic refraction lines

同時(shí)，我們使用了一個(gè)尚未發(fā)表的天然地震臺(tái)站數(shù)據(jù)對(duì)結(jié)果做進(jìn)一步驗(yàn)證，該天然地震數(shù)據(jù)來(lái)自臺(tái)站KINA，臺(tái)網(wǎng)YC，如圖2a中藍(lán)色圓點(diǎn)所示，經(jīng)緯度位置為(6.06°N，116.56°E).對(duì)天然地震數(shù)據(jù)做接收函數(shù)反演，得到該點(diǎn)莫霍面深度為25 km(呂川川等，個(gè)人通訊)，本文反演得到該點(diǎn)莫霍面深度約為25.5 km，兩者結(jié)果十分接近.

圖13 蘇拉威西地區(qū)莫霍面深度與斷裂體系劃分Fig.13 The Moho depth and fault system subdivision in Sulawesi area

4 蘇拉威西地區(qū)莫霍面結(jié)構(gòu)特征

本文反演獲得的莫霍面深度最終結(jié)果如圖13所示，平均深度為20.0 km，深度變化范圍為9.2～33.3 km.總體上，海區(qū)莫霍面淺，大約為10～20 km；陸區(qū)莫霍面深，大約為25～33 km.蘇祿海、蘇拉威西海以及班達(dá)海的莫霍面深度最淺，整體小于20 km.反演得到的莫霍面深度與輸入重力異常的形狀相似，在蘇拉威西海和班達(dá)海成圈閉形態(tài).蘇拉威西島中部，加里曼丹島北部，馬魯古海南部的莫霍面深度大于30 km.

在馬魯古海地區(qū)，本文反演結(jié)果與前人結(jié)果存在明顯差異，原因分析如下：首先，Crust1.0模型精度低，因此不能對(duì)這一小區(qū)域精細(xì)成像；其次，郝天珧等(2014)的反演目標(biāo)研究區(qū)較大，本文研究區(qū)位于其研究區(qū)的東南邊界位置，二者結(jié)果存在的差異是在合理范圍內(nèi)的.

本文在馬魯古海地區(qū)反演得到的莫霍面較深，存在差異的原因猜測(cè)是由于15 Ma至今，背向俯沖的桑義赫海溝與哈馬黑拉海溝沿著左旋走滑的Sula-Sorong斷裂帶不斷向西移動(dòng)，同時(shí)發(fā)生俯沖后撤，在兩個(gè)板塊不斷向下俯沖的過(guò)程中，將馬魯古海帶入更深的位置，導(dǎo)致該區(qū)域的莫霍面更深.對(duì)比本文的莫霍面反演結(jié)果與研究區(qū)斷裂體系劃分結(jié)果(褚偉，2021)，可以看到，超殼斷裂位置往往對(duì)應(yīng)著莫霍面的顯著變化，桑義赫海溝、哈馬黑拉海溝與Sula-Sorong斷裂帶處存在明顯的莫霍面深度變化，本文的莫霍面反演結(jié)果與該區(qū)域的深部構(gòu)造具有很好的相關(guān)性.

5 結(jié)論

(1)基于球坐標(biāo)系的重力正、反演能夠考慮地球曲率的影響，適用于大尺度深部結(jié)構(gòu)研究.本文基于重力數(shù)據(jù)和球棱柱剖分模型，構(gòu)建了以高斯-牛頓迭代法反演莫霍面深度的方法.在缺少地震數(shù)據(jù)等深部有效約束的情況下，本文提出結(jié)合重力異常頻譜分析估計(jì)參考莫霍面平面深度，同時(shí)使用兩次交叉驗(yàn)證選取反演參數(shù)，優(yōu)化反演結(jié)果.

(2)將球坐標(biāo)系的莫霍面反演方法應(yīng)用到蘇拉威西地區(qū)，在實(shí)際應(yīng)用中測(cè)試反演方法，優(yōu)化反演參數(shù)并評(píng)價(jià)反演結(jié)果.本文的反演結(jié)果表明，蘇拉威西地區(qū)的莫霍面平均深度為20.0 km，深度變化范圍為9.2～33.3 km.總體上，海區(qū)莫霍面淺，約10～20 km，陸區(qū)莫霍面深，約25～33 km.將本文的莫霍面反演結(jié)果與全球模型、前人的重力反演結(jié)果以及地震數(shù)據(jù)對(duì)比，證明本文的反演結(jié)果總體合理，能夠反映蘇拉威西地區(qū)的莫霍面變化特征.此外，本文的反演結(jié)果與研究區(qū)斷裂體系具有較好的相關(guān)性，為該區(qū)域的深部結(jié)構(gòu)研究提供了參考依據(jù).

致謝感謝中國(guó)科學(xué)院大學(xué)張健教授，中國(guó)科學(xué)院地質(zhì)與地球物理研究所黃松副研究員、南方舟博士，劍橋大學(xué)呂川川博士及中國(guó)石油大學(xué)(華東)楊庭威在成文過(guò)程中的討論和建議.感謝兩位審稿人的意見(jiàn)和建議，進(jìn)一步提升了本文質(zhì)量.