基于球面重力反演苏拉威西地区莫霍面结构

褚伟， 徐亚*， 郝天珧,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.将本文的莫霍面反演结果与全球模型、前人的重力反演结果以及地震数据对比，证明本文的反演结果总体合理，能够反映苏拉威西地区的莫霍面变化特征.此外，本文的反演结果与研究区断裂体系具有较好的相关性，为该区域的深部结构研究提供了参考依据.

致谢感谢中国科学院大学张健教授，中国科学院地质与地球物理研究所黄松副研究员、南方舟博士，剑桥大学吕川川博士及中国石油大学(华东)杨庭威在成文过程中的讨论和建议.感谢两位审稿人的意见和建议，进一步提升了本文质量.