引言

地球的圈层结构特征已从地震学、地质学和地球动力学等方面得到佐证，其中地球最外层长英质地壳与其下部镁铁质地幔之间的莫霍面是研究地球内部结构的一个重要间断面，其两侧介质的地震波速度、介质密度、各向异性、非均匀性尺度及化学组分均有明显变化，深入研究莫霍面的埋深或地壳厚度可为进一步研究地壳演化和现今构造运动的地球动力学过程提供重要的基础数据，尤其是现代地震学研究（危自根等，2019）显示浅地震震源深度分布与地壳厚度关系密切。全球岩石圈在各个地质时期的动力学演化过程不一致，造成这一全球性界面在各地不尽相同。东亚及西太平洋边缘海岩石圈是由欧亚板块、太平洋板块、印度板块及菲律宾板块拼合而成，这些板块在中新生代又经过强烈的演变，因此呈现出多种结构形式的复杂构造特征，其中中国大陆是晚古生代（距今约250—230 Ma）以来由多个微板块或地块（中朝、扬子、中南半岛、塔里木、哈萨克斯坦、土尔其、伊朗、伊拉克）快速汇聚而成的复合大陆。探测中国大陆复杂的壳幔结构特征，可为认识和理解该地区的地质演化和可能的地球动力学模式提供重要的深部约束（滕吉文等，2002；杨文采等，2017）。从莫霍面发现的100多年以来，中外科学家陆续开展了一系列研究，得到了不同区域不同尺度的地壳厚度模型，采用的方法主要集中在重力方法、主动源地震测深的地震波反射和折射波法、被动源天然地震体波接收函数和地震面波频散方法。典型工作如对区域重力异常数据进行密度界面约束反演，得到中国大陆（黄建平等，2006；郭良辉等，2012）、亚洲及邻区（Stolk et al，2013）的地壳厚度；早期如曾融生等（1995）采用深地震测深得到中国大陆的地壳厚度；近年来孙伟家等（2018）利用深地震反射、宽角地震折射和高温高压岩石物理实验结果确定了中国东部地区壳幔过渡带的厚度及速度。尤其是随着近些年宽频带地震仪的密集布设，利用接收函数研究中国大陆（Chen et al，2010；张攀等，2014；He et al，2014）或局部地区（张广成等，2013；刘启民等，2014；张攀等，2014；Liu et al，2014；Wei et al，2016；Wang et al，2017；Zhang et al，2019）的地壳厚度；或者综合深地震反射、宽角地震折射以及接收函数方法研究中国大陆的地壳厚度（熊小松，2010）。这些方法蓬勃发展并取得了丰硕成果，提供了中国大陆或局部区域的地壳厚度资料，但是不同方法所得地壳厚度模型的局部还不一致，尚有争议（Li et al，2014）。

面波由于其独特的传播特征，其频散性质非常适宜于研究地壳和上地幔的弹性结构。近年由于宽频带地震仪的广泛布设以及测量面波频散方法的改进，用面波对地壳结构和厚度进行成像的精度达到了新水平。大多数研究关注的焦点是基于不同周期的面波频散获得全球（Shapiro，Ritzwoller，2002；Zhou et al，2006）或区域（朱介寿等，2002；Huang et al，2003；易桂喜等，2008；Legendre et al，2015）不同精度的地下S波速度结构，其中代表性的有Huang等（2003）基于10—184 s的面波群速度反演得到了中国及邻区1°×1°高分辨率S波速度结构，易桂喜等（2008）反演得到中国大陆及邻区20—120 s基阶瑞雷面波的相速度分布，这些研究结果为探究中国大陆地下介质的结构特征及其动力学演化提供了重要资料。少数研究利用面波频散获得全球（Devilee et al，1999；Meier et al，2007；Lebedev et al，2013）或局部地区（Cheng et al，2019）的莫霍面深度分布。由于面波群速度和相速度采用的测量方法不同，且二者对深度的敏感性不同，探测的深度有所差别，因此同时采用相速度和群速度反演地壳厚度比单独使用其中之一的效果更为理想（Shapiro，Ritzwoller，2002）。反演地壳厚度的方法可以归纳为两类：一是基于模型驱动的方法，二是基于数据驱动的方法。基于模型驱动的方法考虑地壳厚度和地震面波频散之间的物理关系，其结果很大程度上依赖于初始模型的选择以及问题的线性化处理，而基于数据驱动的神经网络方法可用面波频散反演地壳厚度（Devilee et al，1999；Meier et al，2007）。神经网络适用于正演关系已知而反演问题未知或者很难用传统分析或数值方法得到解的情况（de Wit et al，2013），其目的在于通过学习发现反演问题的映射关系。神经网络在地球物理学领域应用非常广泛，van der Baan和Jutten （2000）总结了神经网络在测井、电阻率层析成像、地震信号处理等方面的应用。目前机器学习用于地学研究正如火如荼地展开，神经网络在解决重复性复杂问题方面具有绝对的优势，具体应用包括：用神经网络方法测定面波频散曲线从而建立地壳模型（Meier et al，2007）和过渡带模型（Meier et al，2009），用神经网络方法计算走时曲线和基阶分裂函数来获得一维地震模型（de Wit et al，2013，2014），基于深度网络来识别地震与爆破（隗永刚等，2019），利用神经网络进行震级分析（林彬华等，2019）。2018年2月于美国新墨西哥召开的固体地球科学中的机器学习（Machine Learning in Solid Earth Geosciences）会议展示了机器学习在地震学尤其是在震相提取、地震事件识别、地震预警、地震预测以及超高分辨率核幔边界图像等方面应用的研究成果（Marone，2018）。

迄今为止，自编码网络主要应用于分类方面，少数用于回归，鲜有自编码网络用于地球物理反演，尤其是地壳厚度反演方面。本文拟考虑深度降噪自编码网络的优点以及根据频散反演地壳厚度的特点，对中国大陆的地壳厚度进行反演，以期通过神经网络表示地震面波频散与地壳厚度之间复杂的非线性关系。

1.   深度降噪自编码网络

地球物理由于其多尺度、多时段、多精度、多比例尺和多解性等的特点，使得联系地球物理观测数据与地球物理模型参数之间的数学关系具有非线性特征，被称为非线性问题。基于数据驱动的人工神经网络是一种高度非线性化的系统，可以模拟任意函数，适用于地球物理信息的提取和非线性处理，同时又是解决此类复杂非线性映射问题的有效手段，特别适用于参数变量与目标函数之间难以用明晰数学式表达的复杂反问题。第三代基于深度学习的机器学习算法在语音识别和机器视觉等领域已得到长足发展并取得了巨大成功。深度学习最近几年取得的进展如此之大主要在于以下几点：发展了基于卷积神经网络（convolutional neural networks，缩写为CNN）和递归神经网络（recursive neural network，缩写为RNN）的深层网络结构；采用了新的激活函数ReLU和新的参数初始化方法；改进了训练算法；实现了避免过拟合的新方法；采用GPU计算，极大地加速了网络训练；发展了一些新开源库（如Caffe，Keras，Theano以及Torch7），扩展了深度学习的应用范围。

对神经网络而言，深度（也称容量）指的是网络层次，是标识网络学习得到的函数与自变量非线性程度的指标。以前广泛采用的浅层神经网络一般指一个输入层、一个隐层和一个输出层的网络架构，只能表示非线性程度较低的函数关系。深度神经网络（也称深度学习）包含多个隐含层，目的在于建立、模拟人脑的多层次信息传递模式进行分析和学习，其主要思想是具有很多隐层的人工神经网络能很好地学习到数据之间的特征，这些特征是对数据更本质的刻画。因此深度网络能从复杂的数据中提取关键的抽象信息，进一步有效学习数据之间的复杂关系，进一步作出正确的预测或分类。

Rumelhart等（1986）提出自编码器神经网络，其主要思想是基于无监督学习实现数据压缩功能。自编码器由编码器网络和解码器网络两部分组成。数据输入编码器，编码器输出为该数据的编码。编码器的输出作为解码器的输入，而解码器的输出试图重建原始数据。由于编码器的维度通常远小于数据输入的维度，所以自编码器是一种数据压缩算法和数据降维方法的综合。早期的自编码器被用于深度神经网络初始化，通过逐层初始化和精调权重来训练一个深度神经网络（Hinton，Salakhutdinov，2006）。当前，自编码器作为一个基本神经网络结构及训练思想，已与其它深度神经网络结合，广泛应用于气体识别（Liu et al，2014，2015）、图像生成、语言翻译及语音检索等领域，并且取得了较好的效果。

This page contains the following errors:

Below is a rendering of the page up to the first error.

图  1  浅层降噪自编码网络结构图

Figure  1.  Structure of shallow denoising autoencoder

下载: 全尺寸图片 幻灯片

编码函数定义为

式中，f₁（x）表示编码激活函数，W ^（1）表示连接输入层与隐含层的权重；b^（1）表示输入层的偏置项。解码函数定义为

式中，f₂（x）表示解码激活函数，W ^（2）表示连接隐含层与输出层的权重；b^（2）表示隐含层的偏置项。

This page contains the following errors:

Below is a rendering of the page up to the first error.

This page contains the following errors:

Below is a rendering of the page up to the first error.

This page contains the following errors:

Below is a rendering of the page up to the first error.

深度去噪自编码（deep denoising autoencoder，缩写为DDAE）网络是由多个浅层自编码网络（DAE）或其变形网络堆叠而成，以上述DAE作为基本结构单元，在编码器和解码器中增加隐层数量而构成。输入层和输出层的神经元个数相同，均为数据的维数。用编码器的输出作为新特征向量代替原始特征向量输入到解码器（回归模型）中，可以大幅度提高模型预测的准确度，并增强鲁棒性（Rifai et al，2011）。与Cheng等（2019）的差别主要在于本文用深度自编码网络（DDAE）的同时训练正演和反演问题，所采用模型如图2所示。

图  2  深度降噪自编码网络（DDAE）结构

Figure  2.  Structure of DDAE

下载: 全尺寸图片 幻灯片

本文DDAE采用的目标函数为

This page contains the following errors:

Below is a rendering of the page up to the first error.

2.   基于深度降噪自编码网络反演中国大陆地壳厚度

神经网络主要用于解决从输入到输出的映射关系不是一一对应的反演问题。而针对影响面波频散的多因素原因，采用神经网络反演地球介质参数的思路框架：网络输入包括瑞雷波的群速度和相速度，网络输出为输入条件下的地壳厚度。根据瑞雷面波群速度和相速度的频散特性，建立多层次的深度神经网络模型，使得深度神经网络模型能准确地揭示面波频散特征与地球结构之间的关系，从而进一步优化网络，求解研究区地壳厚度。使用深度神经网络方法求解反演问题包括以下几步：

1） 对模型空间随机采样，求解正演问题。对采样的径向对称模型根据简正波理论计算理论群速度和相速度，地球模型和计算所得群速度和相速度值构成训练数据；

2） 训练神经网络。将群速度和相速度作为输入、地球模型作为输出来训练网络，即改变网络参数，使得网络输出与期望输出（地球模型）吻合，从而得到训练好的神经网络，这表示从群速度和相速度到地球模型结构的非线性映射关系建立起来了；

3） 反演。把实测的面波频散曲线输入到训练好的网络，得到接近于真实地球模型的输出。

采用基于自编码的深度学习方法来获得青藏高原东缘和扬子地台西缘的地壳厚度。首先，把群速度和相速度同时作为输入，把地壳厚度作为输出，对4万个理论模型进行训练，对1万个理论模型进行测试，从而构建深度神经网络；然后，采用这些神经网络对观测群速度和相速度进行反演，得到研究区的地壳厚度；最后，与其它研究结果的地壳模型进行测试误差和拟合误差的对比，从中选择最佳模型作为研究区的地壳厚度。

2.1   数据准备

根据基本地球参考模型（preliminary reference Earth model，缩写为PREM）把地球模型参数化为185个网格点，参数范围详见de Wit等（2014）一文中的表A.2−A.4；然后使用Minoes软件包进行正演计算（Masters et al，2007），得到5万个理论瑞雷波频散模型，将其作为训练神经网络的输入与输出。Shen等（2016）基于中国及邻区2 000多个地震台的背景噪声数据，使用背景噪声层析成像方法得到中国大陆8—50 s周期内的高分辨率瑞雷波群速度和相速度分布。采用丰富的地震资料和可靠的反演方法，该模型高精度地揭示了地壳和上地幔顶部丰富的地质和构造特征，因此本文采用Shen等（2016）的频散模型作为神经网络的预测输入。

2.2   训练神经网络

理论上，神经网络可以通过对输入数据和输出数据充分学习而挖掘出输入物理量与输出物理量之间的非线性关系。而如何实现充分学习，涉及到神经网络的搭建形式、参数的选择、超参数的调整等一系列复杂流程。本文设置训练集、验证集和测试集的比例为7 ∶ 2 ∶ 1，通过上百次测试和调整，经分析得出神经网络在解决面波频散反演地壳厚度过程中对拟合结果影响最大的参数为学习率、隐层数目、正则化因子、加噪率，优选这些参数的最佳组合得到较为理想的和稳健的神经网络结构。为了评价该网络的性能，本文提出除了考虑传统意义上的测试误差、训练误差之外，用已知物理原理的正演过程与网络预测的结果进行比较，即根据简正波理论正演得到的频散值和神经网络解码模型得到的频散值经主成分分析（principal component analysis，缩写为PCA）降维后进行相关分析（如图3）。在设计网络构架时，同时考虑该地学问题的正反演问题，即解码对应正演过程、编码对应反演过程。由此可见，当输入信号的加噪率不超过0.2时，二者相关度非常高，由此说明由神经网络基于数据信息探究地壳厚度与面波频散与根据正演原理得到的地壳厚度与面波频散关系的吻合度很高，神经网络能够学习到地壳厚度与面波频散之间非线性关系的函数。图4给出了DDAE预测的地壳厚度与实际地壳厚度的对比，可见对测试集的DDEA学习结果与测试集实际地壳厚度的相关性较好，表明本文所采用的网络结构能够很好地揭示地壳厚度与面波频散的关系。

图  3  深度降噪自编码的解码结果与基于频散、地壳厚度的物理关系的正演结果相关性

Figure  3.  Comparison of DDAE decoder result with forward result based on physical relationship between crustal thickness and dispersion

下载: 全尺寸图片 幻灯片

图  4  DDAE预测地壳厚度与实际地壳厚度的相关性

Figure  4.  Comparison of crustal thickness from DDAE with true crustal thickness

下载: 全尺寸图片 幻灯片

2.3   反演

根据已经训练好的神经网络，将Shen等（2016）的频散数据作为预测输入，几秒之内反演得到中国大陆的地壳厚度，如图5所示。CRUST2.0模型（Bassin et al，2000）是根据接收函数、地震波的反射和折射得到的地壳厚度模型（图6a）；CUB2.0模型（Shapiro，Ritzwoller，2002）是用蒙特卡洛方法反演得到的地壳厚度模型（图7a）；Feng模型（Feng et al，2017）首先基于中国大陆内1 669个台站的面波频散反演一维S 波波速模型，然后基于该模型利用接收函数（或与面波联合）反演一维波速初始模型，再根据所得到的波速模型确定台站下方的地壳厚度（图8a）。接下来绘制本文模型与CRUST2.0、CUB2.0以及Feng模型的交汇图（图6b，7b，8b），以定量地比较几个模型的吻合程度。可见：本文模型与CRUST2.0、CUB2.0以及Feng模型在地壳厚度小的地区吻合得较好，在地壳厚度大的区域吻合得较差。这可能与本研究在青藏高原腹地台站较少导致资料覆盖较差以及基于地下介质横向各向同性的假设有关。

图  5  基于深度降噪自编码神经网络反演的中国大陆地壳厚度

Figure  5.  Crustal thickness of Chinese mainland from inversion by DDAE method

下载: 全尺寸图片 幻灯片

图  6  CRUST2.0模型的地壳厚度图（Bassin et al，2000）（a）及其与DDAE反演地壳厚度的相关性（b）

Figure  6.  Crustal thickness from the model CRUST2.0 （Bassin et al，2000） （a） and its correlativity with the crustal thickness from DDAE inversion （b）

下载: 全尺寸图片 幻灯片

图  7  CUB2.0模型的地壳厚度图（Shapiro，Ritzwoller，2002）（a）及其与DDAE反演地壳厚度的相关性（b）

Figure  7.  Crustal thickness from the model CUB2.0 （Shapiro，Ritzwoller，2002） （a） and its correlativity with the crustal thickness from DDAE inversion （b）

下载: 全尺寸图片 幻灯片

图  8  Feng模型的地壳厚度图（Feng et al，2017）（a）及其与DDAE反演地壳厚度的相关性（b）

Figure  8.  Crustal thickness from the model Feng （Feng et al，2017） （a） and its correlativity with the crustal thickness from DDAE inversion

下载: 全尺寸图片 幻灯片

3.   讨论与结论

基于深度自编码模型反演瑞雷波群速度和相速度得到中国大陆的地壳厚度（或莫霍界面埋深）（图5），从地壳厚度分布特点可以把中国大陆从西到东分为地壳逐渐减薄的三个区域，每个区域内随着距离南北向碰撞造山带越远，地壳厚度逐渐减薄。DDEA揭示的中国大陆构造宏观格架（图5）与刘光鼎（2007）提出的“三横、两竖、两个三角”的特点，以及滕吉文等（2002）、杨宝俊等（2005）和胡卫剑等（2014）提出的东亚大陆地区地壳厚度分布的总体趋势为大陆地区周边薄，向内陆逐渐增厚，其中央地域（青藏高原）最厚的特点吻合。

中国大陆地壳总体呈现西厚东薄的特点，厚度在20—80 km范围内变化，地壳厚度从西到东由青藏高原→西扬子克拉通→中国东部→中国大陆边缘海逐级减薄。按照地壳厚度差异自西向东分为三个区域：西部青藏高原地壳汇聚增厚区，其中主要有喜马拉雅地块、拉萨地块、羌塘地块和柴达木地块的地壳，厚度达60—80 km，为增厚型地壳；中西部鄂尔多斯—四川正常地壳区，地壳厚度变化较小，约45 km，主要包括一系列元古代以来地质活动相对较弱、变形相对较小的克拉通地块，如上扬子地块、塔里木地块、准噶尔地块和鄂尔多斯地块，该区域东西两侧的地壳厚度急剧变化，处于45—55 km之间，汶川地震发生所在地的龙门山断裂带正好处于这个地壳厚度陡变的交界区域；东部地壳减薄区，包括大陆东部的华北、华南、黄海和东海大陆架以及南海等，该区域地壳厚度减至20—35 km，这与华东的中生代拆沉有关，其中海南岛的地壳厚度处于25—30 km之间，反映了南海打开时岩石圈拉张减薄的结果（Li et al，2014）。

中国大陆由北向南的地壳厚度特征表现为北部的天山—阿尔金山—祁连山带的莫霍面深度在50—60 km范围内变化，该带以南为青藏高原周边增厚地壳，以北为正常地壳厚度区，塔里木盆地与柴达木盆地分别处于该带南北两侧，形成了沉积盆地的地壳厚度变薄、造山带区地壳厚度增厚的分布特点。南部的喜马拉雅带分布较宽，地壳厚度在40—70 km之间变化。与Li等（2014）给出的华东地壳厚度分布特征类似，本文的地壳厚度模型（图5）揭示出华东有两条窄带，一条为秦岭—大别造山带，另一条为大兴安岭和燕山—华北造山带。秦岭—大别造山带为中朝地块和扬子克拉通向东碰撞的结果，大兴安岭和燕山—华北造山带的地壳厚达40 km，是古代微大陆和弧形地块汇聚的结果。

地壳厚度在一定程度上反映了板块的活动性特征。青藏高原腹地的拉萨地块和羌塘地块的地壳厚达60—80 km，向东到喜马拉雅地块地壳厚度减薄至50 km左右。研究（滕吉文等，2002；杨宝俊等，2005；刘光鼎，2007；胡卫剑等，2014；Li et al，2014；Liu et al，2014；Feng et al，2017）认为，青藏高原巨厚地壳的厚度分布是由于印度板块与欧亚板块强烈的陆陆碰撞，地壳被置换为增厚型地壳结构所致。西扬子地块和鄂尔多斯盆地属于稳定的前寒武地块，在长期的地质时间内相对稳定，地壳物质基本保持原状，且较少经历剧烈的构造运动，其地壳厚度在40—50 km范围内。在大兴安岭—太行山—武陵山以东地区，由于地幔深部热物质上涌引起大规模的火山活动及岩浆活动，造成地壳拉张减薄，其中主要构造带如华北裂谷带、华夏地块、松辽裂谷带的地壳厚度约为30 km。

为分析反演结果的可靠性，本次研究计算了基于DDAE模型采用Minoes 软件包（Masters et al，2007）所得理论群速度与Shen等（2016）的观测群速度的相对误差，结果如图9所示。可见：在地壳厚度较厚的西部区域（如青藏高原），群速度拟合误差较大；在地壳厚度减薄的中西部地区（如华南地块），群速度拟合误差较小。其原因之一可能是在西部地区地震台站分布相对稀疏，射线覆盖率相对不高；之二是在测量群速度时假定基于横向各向同性介质射线沿大圆路径传播的前提，事实上西部地区构造活动强烈、结构复杂，地下介质具有各向异性特征，这使得地震波在传播过程中产生散射从而偏离大圆路径。本文根据DDAE模型计算的理论群速度与观测群速度之间的相对误差结果也显示104°E以东地区的相对误差小于西部地区，这与易桂喜等（2010）给出的中国大陆各向异性总体特征是大概在104°E以东地区地壳上地幔各向异性弱于东部地区的结果相一致。由于介质各向异性以及射线偏离大圆路径传播，射线沿偏离大圆路径传播通常比沿大圆路径传播的路径要长，因此所测得的群速度偏小，由这样的面波频散资料反演的地壳厚度尤其是在西部地区会偏厚（图6b，7b，8b）。

图  9  DDAE模型计算的理论群速度与Shen等（2016）所得观测群速度的相对误差

Figure  9.  The relative error between group velocity calculated by our model DDAE and observed by Shen et al （2016）

下载: 全尺寸图片 幻灯片

基于深度降噪自编码神经网络能有效地揭示出面波频散与地壳厚度的关系，因此采用深度降噪自编码神经网络反演地壳厚度是一种可行的和可信的方法。对神经网络架构体系的评价，本文提出除了考虑传统意义上的测试误差、训练误差之外，用已知物理原理的正演过程与网络预测的结果进行比较，通过二者的相关程度来评价。在设计网络构架时，需要该网络能实现由模型求取瑞雷波相速度和群速度的正演问题，以及由瑞雷波相速度和群速度求模型（即地壳厚度）的反演过程，即解码对应正演过程、编码对应反演过程。同时，针对观测频散数据包含噪声的特点，把训练样本用噪声污染，使解码器解码出无噪声输入，以达到对观测数据降噪的目的。对网络各种参数多次调试、分析再优化组合，最终获得稳健的神经网络，并据此反演出了中国大陆地壳厚度，并经与其它方法得到的地壳厚度模型进行相关分析对比，结果显示本研究结果与其吻合度较高。

最近一些研究通过面波和接收函数联合反演结果（Liu et al，2014）以及结合深地震反射、宽角地震折射和高温高压实验的结果（孙伟家等，2018）显示，莫霍界面在青藏高原东缘以及中国东部的壳幔地震波速度逐渐增大，接收函数反演结果（孙伟家等，2018）也显示出其下不具有尖锐界面的特征，由此推断壳幔分界的不是一个间断面，而是反映地壳与地幔物质交换的“过渡带” 。如何利用神经网络方法深入研究壳幔过渡带的厚度而不仅仅是莫霍界面埋深，是下一步工作需要思考和探讨的方向。

近年来人工智能技术发展迅猛，并成功应用于各行各业。另一方面，地震学经过百余年的发展，对地球精细结构的研究也取得了重大进展。本文目的在于尝试将人工智能的深度学习方法用于地震学领域，试图让机器学习到地震学参量之间的函数关系。本文结果表明基于神经网络用面波频散反演地壳厚度切实可行，且反演结果与其它方法得到的结果吻合度较高。但是，通过神经网络学习到的面波频散与地壳厚度之间的非线性关系，其精度尚需进一步提高，尤其在反演细结构的过程中还需要进一步研究网络结构、优化参数的选择和组合，以及在网络训练过程中的过拟合现象。另外，长期以来，神经网络由于其黑盒子性质而饱受诟病，而可解释神经网络具有内部工作原理透明的特点，并且能对输出结果进行解释，因此可解释神经网络的提出与深入研究应用，将为解决地球物理反演问题另辟蹊径。