-
地壳厚度是地质体的重要参数，不同成因的地质单元具有不同的地壳厚度特征，如正常洋壳平均厚度约为6 km，年轻正常的陆壳厚度平均32~35 km（Mutter and Mutter，1993; Christensen and Mooney，1995）。俯冲带地壳厚度是研究物质循环和计算俯冲通量的重要参数（Hacker，2008; David Völker et al.，2015）。马尼拉俯冲带是南海向菲律宾海板块发生俯冲消亡的场所，是距离中国大陆最近的活跃地震带与火山喷发带，了解该区域地壳厚度变化对研究这一区域的地壳形变、地震-火山灾害发生机制，以及量化边缘海的物质能量循环贡献至关重要。
-
地震方法是准确确定海洋地质体Moho深度的有效方法，如多道反射地震和海底地震仪（OBS）探测，因此被广泛用于区域地壳结构研究。前人在南海马尼拉俯冲带、洋陆转换带及花东海盆等区域开展了一系列高精度地震探测，获取较为准确的二维、三维岩石圈结构（McIntosh et al.，2013；Lester et al.，2014；Eakin et al.，2015; Wan Xiaoli et al.，2019；Fan Chaoyan et al.，2019）。但由于地震方法覆盖范围非常有限，并且在增生楔及岛弧处，强烈的变形以及地壳厚度在较窄空间上的急剧变化，深部的莫霍面反射（Eakin et al.，2014）会变得模糊甚至消失，无法获得准确的地壳厚度。利用重力数据反演地壳厚度，是有效获取大尺度莫霍面深度的常用手段（Wang Tingting et al.，2011；Bai Yongliang et al.，2014）。吴招才等（2017）、于传海等（2017）、Huang Liang et al.（2023）等曾对马尼拉俯冲带邻近区域开展了重力数据反演研究，取得了一定的成果。吴招才等（2017）将校正后的布格重力异常与地震剖面上的莫霍面深度值进行线性拟合之后利用Parker-Oldenburg重力反演法反复迭代，得出南海海盆的莫霍面深度为8~14 km，地壳厚度为3~9 km。随着地震剖面数据的增加和处理手段的进步，Huang Liang et al.（2023）在对重力异常进行了密度扰动和热扰动等精细校正后，以地震剖面数据为约束利用Parker-Oldenburg重力反演方法，反演得到南海区域的莫霍面深度。于传海等（2017）不仅获得了莫霍面深度，还利用沉积层重力异常改正后的地壳布格重力异常及其总水平导数模，来突出南海深部不同尺度的密度结构和莫霍面起伏特征。
-
已有大量研究利用神经网络进行重力反演以推测密度界面的深度。例如，何孟兴（2021）将遗传算法优化的BP神经网络应用于南海莫霍界面的反演，其结果与地震测线数据具有较高的一致性。He Mengxing et al.（2021）开发了一种基于卷积神经网络（CNN）直接从重力数据反演基底深度的新方法，并成功应用于实际的地质研究。Andrea et al.（2023）采用前馈神经网络的深度学习方法估算了潜在油田的基底深度，并将其应用于内华达州 Yucca Flat 沉积盆地，结果与基于传统重力反演方法对该地区的解释相符。Roy et al.（2024）使用生成对抗网络（cGAN）从观测到的重力异常数据中反演莫霍面地形起伏，通过与传统反演方法对比，结果表明该方法能够快速且准确地估算出高分辨率的莫霍面特征。
-
上述研究成果对研究密度界面特别是认识区域莫霍面深度有重要帮助，但也存在明显的争议，主要表现为不同的数据处理方式得到的区域地壳厚度存在很大的差异。特别是在俯冲增生带和洋陆转换带等构造变形复杂的区域，重力反演得到的地壳厚度往往与地震结果不一致。BP神经网络的方法在重力反演中具有广泛的应用潜力，BP神经网络具有强大的非线性映射能力和广泛的适用性。其次，我们在研究中采用了双隐藏层的结构，这在处理复杂的非线性问题上展现了更强的建模能力。BP神经网络作为一种经典的神经网络算法，经过了长时间的发展和优化，其算法成熟且稳定。相比之下，最新的网络架构可能尚处于研究和发展阶段，其稳定性和可靠性尚未得到充分验证。我们也在研究中针对其潜在的不足进行了相应的调整和优化，以确保其适用于本研究的问题。
-
本次研究在全球高精度水深和重力数据的基础上，以海底地震仪（ocean bottom seismometer，OBS）剖面结构作为约束条件，利用BP神经网络的方法反演马尼拉海沟及其邻近区域莫霍面深度和地壳厚度，探讨马尼拉俯冲带及其邻近区域构造作用和地壳厚度变化趋势。
-
1 地质背景
-
马尼拉海沟位于南海海盆与西菲律宾海盆之间，是菲律宾海板块与东亚板块的分界线。它从南部的民都洛岛火山带延伸至北部的台湾岛碰撞带（Sibuet and Hsu，1997），呈现向西略 “凸”的弧形（陈志豪等，2009）。南海东部次海盆沿马尼拉海沟向东俯冲，在海沟东侧形成了较厚的俯冲增生楔（李家彪等，2004）。Fan Jianke et al.（2015）应用层析成像走时反演技术对马尼拉海沟下的南海俯冲板块进行了追踪，提出洋脊俯冲和板块撕裂可能开始于约8 Ma。中国台湾、美国和法国联合开展的海洋TAIGER实验在横跨马尼拉海沟布设了T1、T2海底地震仪和多道地震剖面，揭示了在19°N~21.5°N区域存在一个构造体制过渡带，南部是南海洋壳向菲律宾西部俯冲主导，北部为中国大陆地壳向吕宋火山弧的俯冲和最终的碰撞（Eakin et al.，2014）。Ku and Hsu（2009）分析吕宋岛和台湾岛之间的12条反射地震剖面，结果表明，马尼拉海沟北部的基底一般向东和向南倾斜，海沟南部沉积物充填体积大于北部。
-
花东海盆位于马尼拉海沟和吕宋火山弧以东，是一个四边形的小海盆，西以吕宋岛弧为界、东以加瓜脊为界、向南与菲律宾活动带和吕宋岛相连。花东海盆、加瓜脊和西菲律宾盆地均沿琉球海沟向北俯冲于欧亚板块之下（Huang Chiyue et al.，2019）。花东海盆洋壳的年龄目前还存在很多争议（Hilde and Lee，1984；Deschamps et al.，2000；Qian Shengping et al.，2021），但越来越多的地球化学证据指示该海盆为中生代海盆（Deschamps et al.，2000；Qian Shengping et al.，2021）。基于海上广角地震资料，Eakin et al.（2015）建立了花东海盆和加瓜脊四个东西向剖面的速度结构，提出加瓜脊下面有12~18 km的地壳厚度。此外，速度结构还显示出地壳厚度具有显著的不对称性，这可能由西菲律宾海盆洋壳在花东海盆之下向西俯冲造成。
-
2 数据与方法
-
2.1 数据来源与数据处理
-
重力反演是利用地球的重力场变化来推断地下密度分布的方法。通过测量地表上的重力场数据，建立重力场与地下密度分布的数学关系，进而逆推地下结构。在计算Moho面深度时，只需提取与莫霍起伏相关的残余重力异常。为此需要移除海水、沉积物等因素的重力效应来计算布格重力异常。然后，我们利用布格重力异常向上延拓的方法提取与莫霍起伏相关的残余重力异常。本文的研究区域范围为117°E~125°E，12°N~24°N，涉及到的数据包括水深、自由空间重力异常、沉积物厚度和莫霍面深度。水深数据来源于GEBCO（the General Bathymetric Chart of the Oceans）最新公开的水深数据（http://www.gebco.net）（图1），精度是15弧秒；自由空间重力异常为Sandwell and Smith（2013）的1弧分精度数据（图2a）；沉积物厚度数据来源于美国地球物理数据中心（National Geophysical Data Center，NGDC），并在此基础上加入已发表的地震数据进行校正（图2b）；为了获取多道反射地震数据中的沉积物厚度，我们采用如下时深转换公式（Li Chunfeng et al.，2015）：
-
图1 马尼拉俯冲带区域水深和地震剖面分布图
-
Fig.1 Bathymetric and seismic profile distribution in the Manila subduction zone study area
-
红线为本文讨论中用于结果对比的海底地震仪（ocean bottom seismometer，OBS）剖面，灰线为未用于对比的OBS剖面，黄线为多道地震剖面，红色三角为Hung et al.（2021）布设的部分海底地震仪站位；其中未用于对比的OBS剖面作为神经网络的训练集，而用于结果对比的OBS剖面作为验证集
-
The red line is the ocean bottom seismometer (OBS) profiles used for the comparison of the results discussed in this paper, the gray line is the OBS profiles not used for the comparison, the yellow line is the multichannel seismic profiles, and the red triangles show some of the ocean bottom seismograph sites deployed by Hung et al. (2021) ; in the absence of OBS profiles for comparison, these are employed as the training set for the neural network; conversely, the OBS profiles used for result comparison are utilised as the validation set
-
图2 马尼拉俯冲带区域自由空间重力异常（a）、沉积物厚度（b）、布格重力异常（c）及向上延拓20 km的布格重力异常（d）
-
Fig.2 Free-air gravity anomaly (a) , sediment thickness (b) , Bouguer gravity anomaly (c) , Bouguer gravity anomaly after 20 km of upward continuation (d) in the Manila subduction zone region
-
其中，z为海底之下沉积物厚度，单位为m，t为双程旅行时，单位为ms。
-
Wessel et al.（2019）采用GMT软件的gravfft计算获得布格重力异常，假设地壳和海水的密度分别为2700 kg/m3和1035 kg/m3，计算海底之上的海水的影响（Parker公式计算到4阶），沉积物重力异常使用于传海等（2017）计算结果。考虑到花东海盆的年龄存在争议，我们并没有采用校正热结构影响的重力效应的方式来获取Moho面引发的重力异常（吴招才等，2017；Huang Liang et al.，2023），而是对布格重力异常进行频率域向上延拓的方式分离出Moho面引起的异常（韩波等，2008；何慧优等，2019）。本次研究将布格重力异常向上延拓多个高度进行分析，结果显示向上延拓20 km后既可以清晰识别加瓜海脊和马尼拉俯冲带等区域，又不至于损失很多细节信息。延拓后的重力异常如图2d所示。莫霍面深度数据通过对已发表的OBS剖面上数据化获得（表1）。此外，我们还收集整理了研究区先前的Moho深度反演数据进行对比（Huang Liang et al.，2023）。文中各反演参数与地壳厚度相关系数如表1所示。
-
2.2 地壳厚度反演
-
为更好地反映反演结果，我们采用PSO-BP神经网络进行区域地壳厚度反演并与Parker-Oldenburg法反演结果进行对比。
-
2.2.1 基于Parker-Oldenburg法的地壳厚度反演
-
Oldenburg（1974）根据Parker（1973）正演公式，推导出广泛应用于重磁界面反演的Parker-Oldenburg法，公式如下：
-
式中，ρ是上下两种媒介密度差，单位g/cm3，G是万有引力常数，z0是界面平均深度，单位km，是相对于z0界面起伏的高差，k为波数， F[]为傅里叶变换。
-
Gómez-Ortiz et al.（2005）提出一个基于此方法的程序：Matlab:3DINVER.M，用于与重力异常有关的密度界面的反演。反演基本参数涉及壳幔平均深度和壳幔密度差等。研究区域地壳厚度变化非常明显，不仅仅涉及到南海海盆，还包括吕宋岛弧、花东海盆和西菲律宾海盆等，最薄处5~6 km，最厚处达30 km以上（Eakin et al.，2014; Fan Chaoyan et al.，2019）。为了提高精度并减少Parker-Oldenburg迭代反演法存在的误差，本次研究对全区多个存在OBS剖面且可以看到Moho面的区域设置反演参数。为解决边界效应对反演结果的影响，先将区域扩边，在测量区域外围添加一个扩展区域，将原始测量区域扩大，选取扩展区域长度为目标区域的三分之一，之后利用裁剪出目标区域，得到区域的莫霍面起伏特征。
-
2.2.2 基于PSO-BP神经网络的地壳厚度反演
-
反向传播（back propagation，BP）神经网络是一种多层映射网络，该网络通过误差向后传播的方式进行学习（鲁娟娟等，2006）。数据从输入层输入，向前传播进入隐藏层后，传到输出层，得到预测结果。若输出层的输出与理论输出值存在误差，误差将向后传播，权重和偏差也会随之改变（Sibi et al.，2013）。BP神经网络具有实现任何复杂非线性映射的功能，无需建立模型或了解内部过程。为表明地壳厚度反演的复杂非线性映射关系，本文选择双隐藏层BP神经网络，其模型示意图如图3所示。
-
图3 BP神经网络双隐藏层模型示意图
-
Fig.3 Schematic diagram of the double hidden layer model of BP neural network
-
粒子群优化算法（particle swarm optimization， PSO）是一种基于群体的优化技术，其寻优过程是通过不断迭代更新粒子的位置和速度获得最优解，其更新公式如下：
-
其中，是粒子i经过第t次迭代时在第d维空间的速度，则是对应的粒子的位置，w是惯性权重，c1、c2为学习因子，也称加速常数（acceleration constant），r1、r2为[0，1]范围内的随机数，pibest是第i个粒子当前迭代的最优位置，pgbest是整个粒子群当前迭代的最优位置。
-
传统的BP神经网络收敛速度慢且易陷入局部极值点，因此本文提出基于粒子群优化的双隐藏层的BP算法，优化后的PSO-BP可充分利用PSO算法的全局寻优能力和BP算法的局部最优调整能力，使得神经网络能够更有效地学习输入数据中的模式和规律，改善传统BP神经网络的不足。PSO算法用于优化BP神经网络的权值和阈值。每个粒子表示一个解，通过不断迭代更新粒子的位置和速度，寻找最优解。反向传播算法则在每个训练周期结束后，根据误差反向调整网络的权值和阈值，以减小误差。结合上述优点以提高算法寻求最优解的效率和准确性。
-
本文拟采用基于PSO-BP神经网络的方法进行区域地壳厚度预测。PSO-BP网络实现了从输入到输出特征映射功能。具体的模型构建过程如下：
-
（1）输入数据选取。选取水深、自由空间重力异常、沉积物厚度、延拓计算得到的布格重力异常等作为输入数据，这些特征变量均是对输出具有关键影响的因素。
-
（2）数据预处理。对输入数据进行归一化处理，将输入和输出数据映射到[-1，1]之间，消除不同值之间过大的差异，减少网络收敛的难度，从而提高模型的性能和稳定性。归一化过程由式5得出：
-
其中，xnor为归一化后的数据，x为原始数据，xmin和xmax分别是x中的最小值和最大值。
-
（3）BP模型的构建。本文将BP神经网络模型作为基本网络框架，包括一个输入层，两个隐藏层和一个输出层，设置双隐藏层的目的是让其模型更适应复杂的非线性问题，提高预测的准确度。
-
（4）PSO-BP模型的训练。网络训练及预测的过程均在BP神经网络模型中完成，将研究区域内不同经纬度间隔的水深、自由空间重力异常、沉积物厚度和延拓计算得到的布格重力异常等数据校正到相同的经纬度，建立模型数据集。其中，如图1所示，灰线为BP神经网络的训练集，红线为BP神经网络的测试集，用于后续与地震测线结果以及传统重力反演结果的对比和讨论。以已建立的模型数据集作为训练集，训练神经网络，共5249组数据。在训练过程中，采用PSO算法对BP神经网络的权重（Wmi，Wij，Wjk）和阈值（bmi，bij，bjk）进行优化，使神经网络能够更好地反映输入特征与地壳厚度之间的非线性映射关系。
-
（5）反演目标函数。本文将均方差函数（mean-square error，MSE）作为目标函数，用来衡量模型的预测结果与真实值之间的差异，其值越小，说明预测结果越接近真实值。
-
（6）模型评估。为了评估神经网络模型的性能，本文采用均方根误差（root mean square error，RMSE）和平均绝对百分误差（mean absolute percentage error，MAPE）来衡量预测值与真实值之间的偏差，计算公式如下：
-
式中，yi表示真实值，f（xi）为预测值，n代表测试样本的数量。
-
（7）全区地壳厚度预测。经过多次训练，得到误差最小的网络模型，进行全区地壳厚度的预测，需输入全区的水深、自由空间重力异常和沉积物厚度。
-
综上，基于PSO-BP算法的地壳厚度预测模型流程图如图4所示。首先依据地壳厚度的关键影响因素建立数据集，对输入数据进行预处理，通过双隐藏层的结构，以及PSO算法对BP神经网络的初始权重和阈值的优化，使神经网络能够更好地反映输入特征与地壳厚度之间的非线性映射关系，通过提取最优的权重和阈值得到预测的参数，与真实地壳厚度进行比较并对进行评估。最后，为了获得全区的地壳厚度，利用训练好的PSO-BP模型，以全区的水深、自由空间重力异常、沉积物厚度、延拓20 km的布格重力异常作为输入，得到预测的地壳厚度，并与使用传统方法Parker-Oldenburg进行地壳厚度反演的值进行比较，评估其方法的性能。
-
3 结果与分析
-
将剖面的水深、自由空间重力异常、沉积物厚度和延拓得到的布格重力异常作为输入地壳厚度作为输出，构建训练集，对建立的PSO-BP模型进行训练，评估模型性能。整个输入量是一个4×5249 的矩阵。此矩阵中的延拓得到的布格重力异常也同样作为传统地壳厚度反演方法Parker-Oldenburg的输入。
-
图4 PSO-BP神经网络反演区域地壳厚度流程图
-
Fig.4 Flow chart of the PSO-BP neural network inversion of the regional crustal thickness
-
3.1 BP神经网络超参数选取
-
BP神经网络的输入层包含4个节点，输出层则设有1个节点，其数量依据输入与输出的维度具体确定。此外，在训练开始之前，还需设定一系列超参数，包括隐藏层的节点数量、学习率以及最大迭代次数等。部分超参数可依据经验公式得到，其中隐藏层节点个数和传递函数对网络性能影响较大，且难以确定。以下分析隐含层节点个数和传递函数对网络性能的影响。
-
3.1.1 双隐藏层节点数选取
-
隐藏层和输出层的传递函数使用默认函数分别为正切S型tansig函数和线性purelin函数（Haykin，1999）。学习率为0.01，选择迭代次数300次，模型的验证误差趋于稳定，目标值为10-7，训练函数采用Levenberg-Marquardt算法。采用PSO算法对BP进行优化，其中种群个体数目设置为30，迭代次数设置为50，惯性权重的取值为0.2，通过不断地迭代更新权重和阈值，使预测值与真实的地壳厚度值之间的差异达到最小。通过选取不同双隐藏层节点数的组合，训练PSO-BP神经网络模型，计算出预测的RMSE和MAPE。同时为了研究双隐藏层在模型预测上的优势，选取单隐藏层节点为5~30，5为间隔的模型进行对照。
-
图5为不同隐藏层节点数的地壳厚度预测结果，可以看出。随着单隐藏层节点数的不断增加，RMSE和MAPE总体上均呈现减小的趋势，增加隐藏层的层数时，RMSE和MAPE有所下降，当双隐藏层节点数为[20，20]时，误差均达到最小，当节点数增加，误差有所增大，因此在本研究中选取双隐藏层节点数为[20，20]。
-
3.1.2 传递函数的选取
-
传递函数对训练效果的影响取决于训练数据的非线性映射关系，不同传递函数的选择会对BP神经网络的预测精度带来较大的影响。常见的传递函数有线性purelin传递函数、正切S型tansig传递函数和对数S型logsig传递函数。一般隐藏层节点传递函数选用tansig函数或logsig函数，输出层节点传递函数选用tansig函数或purelin函数，如表2为4种组合的预测误差以及训练耗时。由表可知，当隐藏层选择tansig函数、输出层选择purelin函数时，其预测误差最小，且训练耗时最短，当隐藏层和输出层均选择tansig函数时，也可以取得较好的结果，但耗时较长。
-
3.2 模型结果与评估
-
通过上节确定了双隐藏层最优的节点个数和传递函数，得到了神经网络的结构为4-20-20-1，隐藏层传递函数为tansig函数、输出层传递函数为purelin函数，并通过训练PSO-BP神经网络模型得到了各层最优权值和阈值。选取四条剖面的训练结果与真实值对比结果如图6所示。为了直观显示预测的效果，图6展示了部分训练集中预测值与真实值的拟合曲线。训练集中误差主要集中1 km左右，说明网络很好地反映输入值与输出值的对应关系，具有较高的预测质量。
-
图5 不同隐藏层节点数下的训练误差
-
Fig.5 Training errors under different number of hidden layer nodes
-
图6 马尼拉俯冲带区域训练集真实值与预测值拟合曲线
-
Fig.6 Curve fitted to the true and predicted values of the training set in the Manila subduction zone region
-
（a）—2016-2剖面地壳厚度;（b）—NS8剖面地壳厚度;（c）—T2B剖面地壳厚度;（d）—T2933剖面地壳厚度
-
(a) —crustal thickness in profile 2016-2; (b) —crustal thickness in profile NS8; (c) —crustal thickness in profile T2B; (d) —crustal thickness in profile T2933
-
3.3 花东海盆及马尼拉俯冲带区域地壳厚度特征
-
基于BP神经网络的重力数据反演地壳厚度结果如图7所示。重力反演结果表明沿马尼拉俯冲带地壳厚度变化较大，分布并不均匀，较厚区域集中在海山和残余扩张脊附近，包括黄岩岛、中沙群岛东部的众多小海山。黄岩岛周边区域存在显著的局部地壳增厚，并向马尼拉海沟延伸。沿黄岩岛近EW向的残余扩张中心洋壳厚度在9 km以上，显著高于正常洋壳。残余扩张中心洋壳厚度呈现出南北不对称的特点，残余扩张中心北侧洋壳厚度约7 km略大于南侧约6 km。马尼拉俯冲带西北部可见地壳厚度自北向南经历先增大后减小，这反映了大陆边缘的废弃裂谷的分布和趋势，地壳厚度的强烈变化，反映了早期的强烈减薄运动。地壳厚度迅速减薄的洋陆过渡带包括南海东北部的笔架海山、东沙斜坡，南海东南部的礼乐斜坡。马尼拉海沟的几何形状与9 km洋壳厚度等值线基本一致。花东海盆洋壳厚度大部分为4~7 km，属于典型的洋壳特征，花东海盆东部的加瓜脊可见地壳厚度的显著增厚，可能代表了重要的构造边界。西菲律宾海盆地壳厚度主体为6 km，本哈姆隆起地壳厚度增加至约13 km，可能与早期岩浆供应充足有关（Okino et al.，2003）。
-
4 讨论
-
4.1 不同反演方法计算南海地壳厚度的有效性
-
为了评估网络的泛化能力，将我们的反演结果、地震剖面以及基于Parker-Oldenburg的重力反演结果进行对比（图8和图9a）。其中，Parker-Oldenburg的重力反演法与BP神经网络重力数据反演两者的输入重力异常一致。我们根据前人的研究工作分别为三条剖面设置了反演参数， Moho面深度参考Eakin et al.（2014） 的结果设置为19 km，壳幔密度差采用Huang Liang et al.（2023）的计算结果，地壳密度设置为2.65 g/cm3，地幔密度为3.3 g/cm3，壳幔密度差为0.65 g/cm3。本文随机选取分布在研究区域的3条剖面（图8），沿三条地震剖面，BP神经网络重力反演的结果与地震反演结果在数值上和总体趋势上更加吻合。如图9a所示，BP神经网络反演结果大部分（~52.3%）的误差在±1 km以内，大部分误差集中在±2 km以内。高斯分布拟合误差表明，其平均值为-0.36 km，标准偏差为 1.07 km，极值为7 km。而传统重力反演结果在过渡带处，只有小部分（~26.8%）的误差在±1 km以内，大部分（~55.2%）误差超过±2 km，这其中（~17.0%）超过了±5 km，其高斯分布平均值为 0.18 km，标准偏差为 3.22 km，极值为13 km。
-
图7 马尼拉俯冲带区域地壳厚度等值线图
-
Fig.7 Contour map of crustal thickness in the Manila subduction zone region
-
SCS—南海；WPB—西菲律宾海盆；MT—马尼拉海沟；GR—加瓜脊；HB—花东海盆；ESC—南海残余扩张脊
-
SCS—South China Sea; WPB—west Philippine basin; MT—Manila Trench; GR—Gagua Ridge; HB—Huatung basin; ESC—South China Sea extinct spreading center
-
Hung et al.（2021）在南海海盆残余扩张脊布设若干个海底地震仪站位来研究其增生和变形过程，获取了高精度的地壳厚度。部分海底地震仪站位与本文的研究区域重合，但BP神经网络训练过程并未纳入这些站位，图9b中反演结果的对比可以很好地说明BP神经网络重力反演对相邻研究区域的预测能力。结果表明，BP神经网络重力反演结果与地震反演结果比较一致，两者的差值大多数在1 km以内，只有个别站位如HY08反演结果差值超过2 km。由于训练集采用的卫星重力异常数据分辨率的影响，会在神经网络训练过程中造成部分细节的丢失，这也可能导致反演结果不准确的原因。
-
以上结果表明，在地壳厚度变化显著的过渡带区域，BP神经网络重力反演与OBS剖面更为吻合（图8和图9a），更符合实际地质情况，而传统重力反演通常会因为低通滤波的范围选取不当，会使反演结果过于光滑（Rao Weibo et al.，2023），导致对界面细节起伏特征的刻画能力较低；在Moho面深度起伏不大的区域如南海海盆残余洋中脊附近，BP神经网络与传统重力反演法具有 H比较好的反演结果（图9b）。因此，在该研究区域利用BP神经网络反演计算南海地壳厚度可能会得出更为准确的结果。但值得注意的一个现象是，在图8中OBS2012剖面的70~110 km的位置和图8c中NS8剖面130~170 km的位置，BP神经网络反演结果与OBS剖面出现较大差异。这是由于大陆坡以下地壳厚度存在明显的起伏变化，以及南海东北部陆缘的快速减薄等因素会增大反演的难度。而且Huang Liang et al. （2023） 在进行南海莫霍面精细反演时，也提到地震解释的不确定性。在 OBS 2011 和 OBS 2013-3的交点上，这两个剖面的莫霍面成像深度相差1.6 km，这可能是导致反演结果误差较大的重要因素。另外，本文采用的是卫星重力数据与地震勘探分辨率的不同可能也是造成结果不一致的原因之一（Wang Tingting et al.，2011）。
-
图8 马尼拉俯冲带区域BP神经网络预测结果与OBS剖面以及Parker-Oldenburg重力反演结果对比图
-
Fig.8 Comparison of deep learning inversion results with OBS profiles and Parker-Oldenburg gravity inversion results in the Manila subduction zone region
-
图9 马尼拉俯冲带区域反演结果误差对比的分布直方图（a，虚线表示高斯分布最佳拟合）及不同方法反演结果的比较（b）
-
Fig.9 Comparison of the error of the inversion results (a, the dashed line indicates the best fit of the Gaussian distribution) , and comparison of the inversion results of crustal thickness (b) for the Manila subduction zone region
-
4.2 马尼拉俯冲带、南海残余扩张脊和加瓜脊地壳厚度对称性的解释
-
为了揭示不同构造成因地形高的地壳结构特点，我们跨越马尼拉海沟、南海残余扩张脊和加瓜脊进行了剖面比较（图10），剖面展示了地表高度变化和通过PSO-BP神经网络方法反演得出的莫霍面深度变化（红色实线）。
-
横跨南海残余扩张脊的剖面（M1~M3），海底地形与Moho面深度表现出较好的镜像特征，即海底地形高对应着Moho面深度最大的区域。大洋钻探和岩石拖网定年结果表明（Li Chunfeng et al.，2014; He Enyuan et al.，2023），残留扩张脊附近的海山多数为扩张后岩浆作用。初步推测，南海残留扩张脊扩张作用不对称性弱，且裂后岩浆作用可能在洋脊岩石圈厚度仍较薄的情况下发生，因此表现出较好的镜像均衡特征。
-
沿着马尼拉俯冲带，从北向南（即剖面L1至L4），地壳厚度显示出明显向南增厚的趋势。此外，俯冲带附近地形高与Moho深度表现出非镜像不对称特征，Moho深度的最大值与地表地形最高点在垂向位置上不重合。特别是在海沟南段，随着俯冲角度变陡，剖面上（L3和L4）的地壳厚度变化更为剧烈。这种非镜像不对称的特征与很多俯冲带一致，推测与俯冲板片单向运动导致对上浮板片地幔造成了动力侵蚀（Stern，2002）。
-
加瓜脊的地壳上表面地形和Moho形态也表现出明显非镜像不对称性，这个现象在Eakin et al.（2015）的OBS地震探测中已有揭示。如剖面T2所示，其东侧的地壳厚度逐渐抬升，而西侧的地壳厚度则迅速减薄。通过对地壳厚度对称性的对比研究，加瓜脊的不对称性与马尼拉俯冲带的特征非常相似，这很可能是由于西菲律宾海盆向花东海盆俯冲的结果。此外，加瓜脊西侧地壳厚度的骤减可能是俯冲历史较短从而未造成地壳厚度的明显变大。Qian Shengping et al.（2021）通过地球化学分析指出，加瓜脊的熔岩显示出与俯冲作用相关的弧地球化学特征，进一步支持了过去存在俯冲事件的假说。
-
5 结论
-
本文提出了基于BP神经网络的重力数据反演的方法，构建了从水深、沉积物厚度、自由空间重力异常、延拓后的布格重力异常至地壳厚度的BP神经网络，开展了马尼拉俯冲带及邻区地壳厚度的计算，获得了以下结论：
-
（1）与传统重力反演法相比，在地壳厚度变化显著的区域，BP神经网络重力反演具有明显的优势，与OBS剖面较吻合，更符合实际地质情况。
-
（2）多影响指标控制的定量估计地壳厚度的模型重新估算了马尼拉俯冲带及其邻近区域的地壳厚度。东部次海盆残余扩张中心一直延伸至马尼拉海沟，地壳厚度局部增厚至10 km以上。花东海盆地壳厚度为4~7 km，为典型的洋壳。
-
图10 马尼拉俯冲带区域的Moho深度反演示意图（a）及横跨马尼拉海沟、南海残余扩张脊和加瓜脊的剖面图（b）
-
Fig.10 Schematic of Moho depth inversion in the Manila subduction zone region (a) and profiles across the Manila Trench, the South China Sea extinct spreading ridges and the Gagua Ridge (b)
-
图10a中白线为剖面的位置；T1、T2为Eakin et al.（2015）结果；SCS—南海；WPB—西菲律宾海盆；MT—马尼拉海沟；GR—加瓜脊；HB—花东海盆；ESC—南海残余扩张脊
-
The location of the profile is indicated by the white line in Fig.10a; the results of Eakin et al. (2015) are represented by T1 and T2; SCS—South China Sea; WPB—West Philippine basin; MT—Manila Trench; GR—Gagua Ridge; HB—Huatung basin; ESC—South China Sea extinct spreading center
-
（3）通过对比分析马尼拉俯冲带、南海残余扩张脊和加瓜脊地壳厚度的对称性，结果表明加瓜脊的不对称性与马尼拉俯冲带表现很相似，很可能是西菲律宾向花东俯冲的产物。
-
致谢：文章中使用了国家自然科学基金共享航次计划项目（编号NORC2022-08）采集的多道地震数据，该航次由“实验六”号科考船实施，在此致谢。
-
参考文献
-
Andrea V, Gianluca G, Maurizio F. 2023. Deep Learning to estimate the basement depth by gravity data using Feedforward neural network. Geophysics, 88(3): 1~36.
-
Bai Yongliang, Williams S E, Mueller R D, Liu Zhan, Hosseinpour M. 2014. Mapping crustal thickness using marine gravity data: Methods and uncertainties. Geophysics, 79(2): G27~G36.
-
Chen Zhihao, Li Jiabiao, Wu Ziyin, Wu Zhenli, Shang Jihong. 2009. Tectonic evolution implication of geometry shape characteristics for Manila Trench. Marine Geology & Quaternary Geology, 29(2): 59~65 (in Chinese with English abstract).
-
Christensen N I, Mooney W D. 1995. Seismic velocity structure and composition of the continental-crust—A global view. Journal of Geophysical Research-Solid Earth, 100: 9761~9788.
-
Deschamps A, Monié P, Lallemand S, Hsu S K, Yeh K Y. 2000. Evidence for Early Cretaceous oceanic crust trapped in the Philippine Sea Plate. Earth and Planetary Science Letters, 179(3~4): 503~516.
-
Eakin D H, McIntosh K D, Van Avendonk H J A, Lavier L, Lester R, Liu C S, Lee C S. 2014. Crustal-scale seismic profiles across the Manila subduction zone: The transition from intraoceanic subduction to incipient collision. Journal of Geophysical Research-Solid Earth, 119(1): 1~17.
-
Eakin D H, McIntosh K D, Van Avendonk H J A, Lavier L. 2015. New geophysical constraints on a failed subduction initiation: The structure and potential evolution of the Gagua Ridge and Huatung basin. Geochemistry Geophysics Geosystems, 16(2): 380~400.
-
Fan Chaoyan, Xia Shaohong, Cao Jinghe, Zhao Fang, Sun Jinlong, Wan Kuiyuan, Xu Huilong. 2019. Lateral crustal variation and post-rift magmatism in the northeastern South China Sea determined by wide-angle seismic data. Marine Geology, 410: 70~87.
-
Fan Jianke, Wu Shiguo, Spence G. 2015. Tomographic evidence for a slab tear induced by fossil ridge subduction at Manila Trench, South China Sea. International Geology Review, 57(5~8): 998~1013.
-
Gómez-Ortiz D, Agarwal B N P. 2005. 3DINVER. M: A MATLAB program to invert the gravity anomaly over a 3D horizontal density interface by Parker-Oldenburg's algorithm. Computers & Geosciences, 31(4): 513~520.
-
Hacker B R. 2008. H2O subduction beyond arcs. Geochemistry Geophysics Geosystems, 9.
-
Han Bo. 2008. Research on geophysical field and deep geological tectonics of the East China Sea. Graduate School of Chinese Academy of Sciences (Institute of Oceanography) (in Chinese with English abstract).
-
Haykin S. 1999. Neural Networks: A Comprehensive Foundation (2nd Edition).
-
He Enyuan, Zhao Minghui, Huang Haibo, Li Yuhan, Tan Pingchuan, Qiu Xuelin, Zhang Xubo, Liu Kang. 2023. Seismic structures of Zhenbei and Huangyan seamounts and their postspreading volcanic evolution. Tectonophysics, 869: 230102.
-
He Huiyou, Fang Jian, Chen Ming, Cui Ronghua. 2019. Moho depth of the East China Sea inversed using gravity data. Geomatics and Information Science of Wuhan University, 44(5): 682~689(in Chinese with English abstract).
-
He Mengxing. 2021. Genetic optimisation algorithm of BP neural network and its application to Moho interface inversion in the South China Sea. Master's thesis of China University of Geosciences (Beijing) (in Chinese with English abstract).
-
Hilde T W C, Lee C S. 1984. Origin and evolution of the west Philippine basin: A new Interpretation. Tectonophysics, 102: 85~104.
-
Huang Chiyue, Wang Pinxian, Yu Mengming, You Chenfeng, Liu C S, Zhao Xixi, Shao Lei, Zhong Guangfa, Yumul G P. 2019. Potential role of strike-slip faults in opening up the South China Sea. National Science Review, 6(5): 891~901.
-
Huang Liang, Wen Yonglin, Li Chunfeng, Peng Xi, Lu Zhezhe, Xu Liuna, Yao Yongjian. 2023. A refined Moho depth model from a joint analysis of gravity and seismic data of the South China Sea basin and its tectonic implications. Physics of the Earth and Planetary Interiors, 334(3): 106966.
-
Hung T D, Yang Ting, Le B M, Yu Youqiang, Xue Mei, Liu Baohua, Liu Chenguang, Wang Jian, Pan Mohan, Huong P T, Liu Fang, Morgan J P. 2021. Crustal structure across the extinct mid-ocean ridge in South China Sea from OBS receiver functions: Insights into the spreading rate and magma supply prior to the ridge cessation. Geophysical Research Letters, 48(3): e2020GL089755.
-
Ku Chiayen, Hsu Shukun. 2009. Crustal structure and deformation at the northern Manila trench between Taiwan and Luzon islands. Tectonophysics, 466(3): 229~240.
-
Li Jiabiao, Jin Xianglong, Ruan Aiguo, Wu Shimin, Wu Ziyin, Liu Jianhua. 2004. Extrusion tectonics of the middle part of the accretionary wedge in the Manila trench. Science Bulletin, (10): 1000~1008 (in Chinese with English abstract).
-
Lester R, McIntosh K, Van Avendonk H J A, Lavier L, Liu C S, Wang T K. 2013. Crustal accretion in the Manila Trench accretionary wedge at the transition from subduction to mountain-building in Taiwan. Earth and Planetary Science Letters, 375: 430~440.
-
Lester R, Van Avendonk H J A, McIntosh K, Lavier L, Liu C S, Wang T K, Wu F. 2014. Rifting and magmatism in the northeastern South China Sea from wide-angle tomography and seismic reflection imaging. Journal of Geophysical Research-Solid Earth, 119(3): 2305~2323.
-
Li Chunfeng, Xu Xing, Lin Jian, Sun Zhen, Zhu Jian, Yao Yongjian, Zhao Xixi, Liu Qingsong, Kulhanek D K, Wang Jian, Song Taoran, Zhao Junfeng, Qiu Ning, Guan Yongxian, Zhou Zhiyuan, Williams T, Bao Rui, Briais A, Brown E A, Chen Yifeng, Clift P D, Colwell F S, Dadd K A, Ding Weiwei, Almeida I H, Huang Xiaolong, Hyun S, Jiang Tao, Koppers A A P, Li Qianyu, Liu Chuanlian, Liu Zhifei, Nagai R H, Peleo-Alampay A, Su Xin, Tejada M L G, Hai Son T, Yeh Y C, Zhang Chuanlun, Zhang Fan, Zhang Guoliang. 2014. Ages and magnetic structures of the South China Sea constrained by deep tow magnetic surveys and IODP Expedition 349. Geochemistry Geophysics Geosystems, 15(12): 4958~4983.
-
Li Chunfeng, Li Jiabiao, Ding Weiwei, Franke D, Yao Yangjian, Shi Hesheng, Pang Xiong, Cao Ying, Lin Jian, Kulhanek D K, Williams T, Bao Rui, Briais A, Brown E A, Chen Yifeng, Clift P D, Colwell F S, Dadd K A, Hernández-Almeida I, Huang Xiaolong, Hyun S, Jiang Tao, Koppers A A P, Li Qianyu, Liu Chuanlian, Liu Qingsong, Liu Zhifei, Nagai R H, Peleo-Alampay A, Su Xin, Sun Zhen, Tejada M L G, Trinh H S, Yeh Y C, Zhang Chuanlun, Zhang Fan, Zhang Guoliang, Zhao Xixi. 2015. Seismic stratigraphy of the central South China Sea basin and implications for neotectonics. Journal of Geophysical Research-Solid Earth, 120(3): 1377~1399.
-
Liu Siqing, Zhao Minghui, Sibuet J C, Qiu Xuelin, Wu J, Zhang Jiazheng, Chen Chuanxu, Xu Ya, Sun Longtao. 2018. Geophysical constraints on the lithospheric structure in the northeastern South China Sea and its implications for the South China Sea geodynamics. Tectonophysics, 742: 101~119.
-
Lu Juanjuan, Chen Hong. 2006. Researching development on BP Neural networks. Control Engineering of China, 13(5): 449~451+456.
-
McIntosh K, van Avendonk H, Lavier L, Lester W R, Eakin D, Wu F, Liu C S, Lee C S. 2013. Inversion of a hyper-extended rifted margin in the southern Central Range of Taiwan. Geology, 41(8): 871~874.
-
McIntosh K, Lavier L, van Avendonk H, Lester R, Eakin D, Liu C S. 2014. Crustal structure and inferred rifting processes in the northeast South China Sea. Marine and Petroleum Geology, 58: 612~626.
-
Mutter C Z, Mutter J C. 1993. Varitions in thickness of layer-3 dominate oceanic crustal structure. Earth and Planetary Science Letters, 117: 295~317.
-
Okino K, Fujioka K. 2003. The central basin spreading center in the Philippine Sea: Structure of an extinct spreading center and implications for marginal basin formation. Journal of Geophysical Research-Solid Earth, 108(B1): 2040.
-
Oldenburg D W. 1974. Inversion and Interpretation of Gravity Anomalies. Geophysics 39(4): 526~536.
-
Parker R L. 1973. The rapid calculation of potential anomalies. Geophysical Journal International, 31(4): 447~455.
-
Qian Shenping, Zhang Xiaozhi, Wu J, Lallemand S, Nichols A R L, Huang Chiyue, Miggins D P, Zhou Huaiyang. 2021. First identification of a Cathaysian continental fragment beneath the Gagua Ridge, Philippine Sea, and its tectonic implications. Geology, 49(11): 1332~1336.
-
Rao Weibo, Yu Nan, Chen Gang, Xu Xinyu, Zhao Simin, Song Zhaoqi, Liu Jiazheng, Wan Longtu, Hu Changhong. 2023. Moho inversionby gravity anomalies in the South China Sea: Updates and improved iteration of the Parker-Oldenburg algorithm. Earth and Space Science, 10: e2022EA002794.
-
Roy A, Sharma R K, Jash D, Rao B P, Dev J A, Tomson J. 2024. Imaging of Moho topography with conditional generative adversarial network from observed gravity anomalies. Journal of Asian Earth Sciences, 265: 106093.
-
Sandwell D, Garcia E, Soofi K, Wessel P, Chandler M, Smith W H F. 2013. Toward 1-mGal accuracy in global marine gravity from CryoSat-2, Envisat, and Jason-1. Leading Edge, 32(8): 892~899.
-
Sibi P, Jones S A, Siddarth P. 2013. Analysis of different activation functions using backpropagation neural networks. Journal of Theoretical & Applied Information Technology, 17(3): 1264~1268.
-
Sibuet J C, Hsu S K. 1997. Geodynamics of the Taiwan arc-arc collision. Tectonophysics, 274(1): 221~251.
-
Stern R J. 2002. Subduction zones. Reviews of Geophysics, 40: RG1012.
-
Völker D, Stipp M. 2015. Water input and water release from the subducting Nazca Plate along southern Central Chile (33°S-46°S). Geochemistry Geophysics Geosystems, 16: 1825~1847.
-
Wan Kuiyuan, Xia Shaohong, Cao Jinghe, Sun Jinlong, Xu Huilong. 2017. Deep seismic structure of the northeastern South China Sea: Origin of a high-velocity layer in the lower crust. Journal of Geophysical Research-Solid Earth, 122: 2831~2858.
-
Wan Xiaoli, Li Chunfeng, Zhao Minghui, He Enyuan, Liu Siqing, Qiu Xuelin, Lu Yu, Chen Nan. 2019. Seismic velocity structure of the magnetic quiet zone and continent-ocean boundary in the northeastern south China Sea. Journal of Geophysical Research-Solid Earth, 124(11): 11866~11899.
-
Wang Tankin, Chen Mingkai, Lee Chaoshing, Xia Kanyuan. 2006. Seismic imaging of the transitional crust across the northeastern margin of the South China Sea. Tectonophysics, 412: 237~254.
-
Wang Tingting, Lin Jian, Tucholke B, Chen Y J. 2011. Crustal thickness anomalies in the North Atlantic Ocean basin from gravity analysis. Geochemistry Geophysics Geosystems 12: Q0AE02, 2.
-
Wessel P, Luis J F, Uieda L, Scharroo R, Wobbe F, Smith W H F, Tian D. 2019. The generic mapping tools Version 6. Geochemistry, Geophysics, Geosystems, 20(11): 5556~5564.
-
Wu Zhaocai, Gao Jinyao, Ding Weiwei, Shen Zhongyan, Zhang Tao, Yang Chunguo. 2017. Three-dimensional gravity-constrained inversion Moho surface depth and its characteristics in the South China Sea basin. Geophysical Journal, 60(7): 2599~2613(in Chinese with English abstract).
-
Yu Chuanhai, Zhao Junfeng, Shi Xiaobin, Yang Xiaoqiu, Ren Ziqiang, Chen Mei. 2017. Sediment density correction of gravity anomaly in the South China Sea and its significance to analyze regional tectonic characteristics. Chinese Journal of Geophysics, 60(8): 3151~3166(in Chinese with English abstract).
-
陈志豪, 李家彪, 吴自银, 吴振利, 尚继宏. 2009. 马尼拉海沟几何形态特征的构造演化意义. 海洋地质与第四纪地质, 29(2): 59~65.
-
韩波. 2008. 东海地球物理场及深部地质构造研究. 中国科学院研究生院（海洋研究所）.
-
何慧优, 方剑, 陈铭， 崔荣花. 2019. 利用重力数据反演中国东海海域莫霍面深度. 武汉大学学报(信息科学版), 44(5): 682~689.
-
何孟兴. 2021. 遗传优化算法的BP神经网络及其在南海莫霍界面反演中的应用. 中国地质大学(北京)硕士学位论文.
-
李家彪, 金翔龙, 阮爱国， 吴世敏， 吴自银， 刘建华. 2004. 马尼拉海沟增生楔中段的挤入构造. 科学通报, (10): 1000~1008.
-
鲁娟娟, 陈红. 2006. BP神经网络的研究进展. 控制工程, 13(5): 449~451+456.
-
吴招才, 高金耀, 丁巍伟, 沈中延, 张涛, 杨春国. 2017. 南海海盆三维重力约束反演莫霍面深度及其特征. 地球物理学报, 60(7): 2599~2613.
-
于传海, 赵俊峰, 施小斌, 杨小秋， 任自强， 陈梅. 2017. 南海重力异常的沉积层密度改正及其对区域构造特征分析的意义. 地球物理学报, 60(8): 3151~3166.
-
摘要
地壳厚度是地质体的重要参数，对认识区域地壳形变以及量化物质循环至关重要。由于变形破坏以及地壳厚度在较窄空间上的急剧变化，地震探测方法对俯冲增生带的莫霍面成像质量较差，很难获得准确的地壳厚度。学者们常利用重力数据反演俯冲带的地壳厚度，由于缺乏其他地球物理数据的约束，加上传统的 Parker-Oldenburg 算法在迭代过程中常需引入低通滤波以确保计算稳定。然而，如果低通滤波的范围选择不当，可能会使反演结果过于光滑。这些问题导致重力反演得到的地壳厚度结果与地震数据的一致性较差。本文基于卫星重力数据和大量海底地震仪（ocean bottom seismometer, OBS）剖面，运用优化的反向传播(back propagation, BP)神经网络的方法，建立了一个多影响指标控制的定量估计地壳厚度的模型，并利用该模型反演了马尼拉俯冲带与花东海盆区域地壳厚度。与传统重力反演（Parker-Oldenburg）法相比，在地壳厚度变化显著的区域，BP神经网络重力反演方法具有明显的优势，与OBS剖面较吻合，更符合实际地质情况。结果表明，南海海盆地壳厚度在4~9 km，东部海盆残余扩张中心呈SW-NE向一直延伸至马尼拉海沟，地壳厚度局部增厚至10 km以上。花东海盆地壳厚度为4~7 km，推测为洋壳属性。通过分析不同构造单元地壳厚度的对称性，本研究推测加瓜脊可能是西菲律宾海盆向花东海盆俯冲的产物。
Abstract
Determining crustal thickness is crucial for understanding regional crustal deformation and material cycling. However, in subduction accretionary zones, complex deformation and rapid changes in crustal thickness over narrow spaces often lead to poor imaging quality of Moho surfaces, thereby hindering the accurate determination of crustal thickness through seismic methods. Gravity data are frequently used to infer crustal thickness in these zones due to the dearth of alternative geophysical data and the fact that the conventional Parker-Oldenburg algorithm often necessitates low-pass filtering to ensure stability during the iterative calculation process. However, inappropriate filtering can lead to overly smooth inversion results, creating discrepancies between gravity inversion and seismic data regarding crustal thickness. In this paper, we utilize satellite gravity data and a substantial number of ocean bottom seismometer (OBS) profiles to develop a model for quantitatively estimating crustal thickness, taking into account a range of influencing factors. We apply this model to invert the crustal thicknesses of the Manila subduction zone and the Huatung basin. In comparison to the conventional Parker-Oldenburg method, the BP neural network gravity inversion method demonstrates notable advantages in regions characterized by significant crustal thickness variations. This approach aligns more closely with OBS profiles and is more consistent with the actual geological situation. Our results demonstrate that the crustal thickness of the South China Sea basin ranges from 4 to 9 km. A localized thickening to over 10 km occurs near the residual spreading center of the eastern basin, extending SW-NE towards the Manila trench. The crustal thickness of the eastern basin is 4 to 7 km, which is assumed to be oceanic crust. Furthermore, analysis of crustal thickness symmetry in different tectonic units suggests that the Gagua Ridge may be a product of subduction from the West Philippine basin to the Huatung basin.
Keywords
crustal thickness ； gravity inversion ； BP neural network ； Manila Trench
