基于球面重力反演蘇拉威西地區莫霍面結構

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

1 中國科學院大學地球與行星科學學院， 北京 100049 2 中國科學院地質與地球物理研究所, 中國科學院油氣資源研究重點實驗室， 北京 100029 3 中國科學院地球科學研究院， 北京 100029 4 海底科學重點實驗室， 自然資源部第二海洋研究所， 杭州 310012

0 引言

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

蘇拉威西地區位于太平洋板塊、菲律賓海板塊、歐亞板塊和印度洋—澳大利亞板塊的匯聚地帶，包括多條海溝和斷裂帶，地質背景與地質演化過程復雜.研究區大部分位于海域，目前地球物理觀測尚不充分，因此使用衛星重力觀測數據反演莫霍面深度是目前較為可行的方法.反演得到的莫霍面結構可為該區域的三維密度結構反演提供約束，進而為動力學模擬提供邊界條件，最終為認識該區域俯沖轉換邊緣、對向俯沖系統以及俯沖初始機制等獨特地質現象提供依據.

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

1 地質與地球物理概況

1.1 地質與構造特征

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

圖1 蘇拉威西地區區域構造簡圖 部分構造單元下文簡稱：科托巴托海溝(CTB)，桑義赫海溝(SGH)，北蘇拉威西海溝(NSW)，Palu-Koro斷裂帶(PK)，Sula-Sorong斷裂帶(SS)，古俯沖帶(PSZ)，Mantano斷裂帶(MTN)，UNA-UNA火山(UU)，俯沖轉換傳遞邊緣(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).

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

1.2 莫霍面研究現狀

目前蘇拉威西地區的莫霍面深度研究資料以及地球物理研究成果較少，隨著近年來對該地區的地球物理調查逐漸增多，不同學者陸續從地震學、重力、磁法、地熱學等方面開展對該區域深部結構的研究.

全球模型是對區域開展整體認識的重要依據，Crust1.0全球模型(Laske et al., 2013)給出了研究區范圍1°×1°的莫霍面深度(圖2a)，莫霍面深度在10～33 km之間，平均深度為25 km.全球層析成像模型可以用于研究區深部地幔結構的研究，目前有較多的全球層析成像模型，例如MIT-P08(Li et al., 2008)、LLNL-G3Dv3(Simmons et al., 2012)，Fukao和Obayashi(2013)等給出的多個P波層析模型.此外，在東南亞和西太平洋地區，區域地震層析成像研究也給出了研究區地幔的三維成像結果(Hall and Spakman, 2015; Huang et al., 2015; Zenono et al., 2019)，但是層析成像結果在淺層的分辨率較低，無法為本文莫霍面深度反演提供有效的約束.

圖2 蘇拉威西地區莫霍面深度(a) 可用于研究區莫霍面深度對比的數據，其中藍色圓點為天然地震數據位置，紅色實線是雙船折射地震測線位置，灰色虛線是Crust1.0全球模型的莫霍面深度等值線; (b) 研究區莫霍面深度重力反演結果，數據來自郝天珧等(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).

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

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

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

1.3 重力異常特征

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

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

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

對比空間重力異常及布格重力異常，可見在Sula-Sorong斷裂帶北側以及馬魯古海的位置均呈現明顯的低值異常區，由于馬魯古海兩側是背向俯沖的桑義赫海溝與哈馬黑拉海溝，南側是走滑斷裂帶，因此布格重力異常低值可能與該區域復雜的地質構造背景有關.

1.4 重力異常場分離

重力異常是地下所有地質體密度差異引起重力響應的疊加場，通常淺層或小規模密度異常體與深層或大規模密度異常體引起的異常具有明顯的頻譜差異，通過頻譜分析可獲得對地下密度異常體分布的總體認識，并據此進行深-淺場源重力異常的分離.

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

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

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

2 莫霍面反演方法

2.1 球面重力正演方法

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

圖6 球棱柱示意圖Q為積分點，在球棱柱內部，位于地心坐標系(X,Y,Z).P為計算點，位于局部坐標系(x,y,z)，r,φ,λ分別是P點的半徑、緯度和經度.l為P和Q點在笛卡爾坐標系中的距離(據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).

位于球棱柱外的任意觀測點在球坐標系中的重力位和重力加速度的體積分公式如下：

(1)

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

(2)

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

2.2 莫霍面反演方法

莫霍面反演基于區域重力異常，該重力異常消除了地殼內部其他構造及密度體分布的影響，根據數據誤差和正則化項構建反演目標函數Γ(p)(Uieda and Barbosa，2017)，

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

+μpTRTRp,

(3)

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

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

(4)

(5)

Hk=2AkTAk+2μRTR,

(6)

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

(7)

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

A=2πGΔρI,

(8)

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

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

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

對最優化參數選取的交叉驗證實現過程如下：

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

第二次隨機子抽樣交叉驗證方法確定最優的密度差Δρ，具體過程為：

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

4)對應的均方根誤差最小的zref與Δρ為最優超參數，對應的莫霍面深度pl,m為最終反演結果.

3 蘇拉威西地區莫霍面反演與分析

3.1 反演參數選取

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

表1 不同超參數組合測試Table 1 Tests of different combination of hyperparameters

可以看到，反演算法對參考莫霍面深度敏感性很強.進一步比較上述不同超參數組合的反演結果，統計重力異常擬合數據、重力異常擬合數據誤差、反演莫霍面深度的平均值及變化幅度，統計結果如表2所示.重力異常擬合數據誤差的直方圖如圖9所示.綜合重力異常的擬合程度等統計數據，認為Test 2的反演結果最合理，Test 2反演得到的重力異常擬合數據總誤差在±15 mGal范圍內，平均值為0.04 mGal，標準差為1.80 mGal，能夠滿足反演精度，因此將該結果作為本文最終的反演結果.

表2 反演結果統計分析Table 2 Statistical analysis of the inversion results

3.2 反演結果分析

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

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

3.2.1 反演結果與Crust1.0全球模型對比

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

圖10 本文莫霍面深度反演結果與Crust1.0全球模型對比(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 反演結果與前人重力反演結果對比

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

圖11 本文莫霍面深度反演結果與郝天珧等(2014)對比(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 反演結果與地震數據對比

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

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

同時，我們使用了一個尚未發表的天然地震臺站數據對結果做進一步驗證，該天然地震數據來自臺站KINA，臺網YC，如圖2a中藍色圓點所示，經緯度位置為(6.06°N，116.56°E).對天然地震數據做接收函數反演，得到該點莫霍面深度為25 km(呂川川等，個人通訊)，本文反演得到該點莫霍面深度約為25.5 km，兩者結果十分接近.

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

4 蘇拉威西地區莫霍面結構特征

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

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

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

5 結論

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

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

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