2020, Vol. 63
2. 中国科学技术大学地球和空间科学学院地震与地球内部物理实验室, 合肥 230026;
3. 中国科学技术大学蒙城地球物理国家野外科学观测研究站, 安徽蒙城 233527;
4. Nanyang Technological University, 50 Nanyang Ave 639798, Singapore;
5. University of Toronto, Toronto, Ontario, M5S 1A1, Canada
2. Laboratory of Seismology and Physics of Earth's Interior, School of Earth and Space Sciences, University of Science and Technology of China, Hefei 230026, China;
3. Mengcheng National Geophysical Observatory, University of Science and Technology of China, Mengcheng Anhui 233527, China;
4. Division of Mathematical Sciences, School of Physical and Mathematical Sciences, Nanyang Technological University, Singapore;
5. Department of Physics, University of Toronto, Toronto, Ontario, M5S 1A1, Canada
地震体波和面波成像方法是目前地球岩石圈内部结构成像的重要手段.由于它们的频率特性以及传播路径存在差异，各自的成像深度范围及精度有所不同.其中，面波数据具有较强的频散特性，主要对地壳速度结构敏感，因此常常用于岩石圈地壳上地幔结构的研究.但地震面波周期通常大于20 s，难以准确约束浅部的地壳速度结构.前人研究发现利用背景噪声数据互相关可以提取出中短周期的面波成分(Shapiro and Campillo, 2004; Shapiro et al., 2005)，从而可以约束浅层的S波速度结构，有效地推动了面波成像技术的发展和应用(Yang et al., 2008; Yao et al., 2006, 2008; Lin et al., 2008).但面波的频率成分偏低，只能获得平滑的S波速度结构，对于精细结构的分辨能力有限.同时面波成像难以对Moho面形态进行精确约束.远震体波由于经过地幔向上传播可以穿透岩石圈被台站接收，因此可用于岩石圈上地幔和下地壳的结构研究，例如：基于远震体波数据的接收函数成像方法可以利用速度间断面(如Moho面、410 km间断面等)所产生的转换波(或多次波)信号获取台站下方的S波速度结构、地壳厚度以及各向异性分布等(Zhu and Kanamori, 2000; Ammon et al., 1990; Galetti and Curtis, 2012; Owens et al., 1984).但接收函数反演中存在强烈的非唯一性，对初始模型依赖程度较高.只利用面波或远震体波数据单独反演均存在缺陷，因此前人提出了面波和体波联合反演的思路，利用两种数据的高度互补性，获得岩石圈的精细结构.其中，基于射线理论的面波频散和体波接收函数联合成像方法已经发展成熟(Julià et al., 2000; Liu et al., 2010; Li et al., 2017a)，并应用于中国大陆岩石圈(Chen and Niu, 2016)、美国中西部(Shen et al., 2013)及区域尺度的3-D横波速度结构的构建(Liu et al., 2018; Bao et al., 2015; Li et al., 2017b; 郑晨等, 2018).尽管如此，目前联合反演方法大多基于单点下方的1-D模型假设，其分辨率会受到一定的限制；另外最终获得的3-D模型通常由单点反演的1-D模型通过插值产生，准确性有待提高，因此迫切需要发展更为高分辨的地壳上地幔成像方法.
得益于数值建模技术的发展以及计算机存储和计算能力的提升，数值求解全波方程来模拟地震波在大区域复杂介质中的传播已经变成了现实(Komatitsch et al., 2004).直接利用波动方程理论，求解弹性波方程或者标量波动方程来建立层析成像反演理论成为一种自然的选择.因此，基于全波方程数值模拟的波形层析成像得到发展，被认为是最有潜力的新一代成像方法(Liu and Gu, 2012).目前通常有两种实现形式:一种是散射积分法(Scattering-integral Method)，通过计算和存储每个数据泛函的Fréchet导数来求解反演问题(Zhao et al., 2005; Chen et al., 2007)；另一种是伴随方法(Adjoint-state Method)，通过从接收点处反传波场的形式对地下结构成像，称为伴随层析成像(Tarantola, 1982; Tromp et al., 2005).目前，在天然地震结构成像中应用最为广泛的是基于波形走时的伴随方法.它介于走时层析成像和全波形反演之间，既克服了射线理论的高频近似假设，又降低了全波形反演过程中非线性化程度.这种方法在区域结构和全球结构成像中得到成功应用(Tape et al., 2009; Zhu et al., 2017; Chen et al., 2017; Bozdaǧ et al., 2016; Fichtner et al., 2018; Tao et al., 2018).而在小尺度区域结构(2°~5°)成像中，由于地震事件的缺乏，主要依靠远震体波数据进行成像，但目前谱元法模拟全球尺度周期8 s以下的远震波形仍比较困难.为了降低计算成本，地震学家提出可以采用一个较小尺寸的局部三维区域用于三维波场建模，即采用一个混合方法来计算高频远震地震波场.目前主流的混合方法有两类：一类为直接解-谱元方法(简称DSM-SEM, e.g., Monteiller et al., 2013)；另一类为频率波数域-谱元方法(简称FK-SEM, e.g., Tong et al. 2014; Liu et al., 2017).在混合方法大幅提升计算效率之后，远震体波伴随层析成像方法也得到发展和应用(Wang et al., 2016; Beller et al., 2017).除了天然地震数据，伴随层析成像也应用到背景噪声数据.Chen等(2014)首次提出背景噪声伴随层析成像方法，并成功应用到青藏高原南部与四川盆地边缘地区，获得了地壳及上地幔顶部横波速度结构，反映出该地区地壳存在3-D低速异常体的布局.Zhu(2018)和Wang等(2018)先后将背景噪声伴随层析成像应用到美国北德克萨斯和俄克拉荷马州以及南加州区域地壳结构成像.
在获得高精度3-D速度结构的同时，不能忽略的是3-D伴随层析成像方法庞大的计算消耗和用时(Bozdaǧ et al., 2016).相对而言，2-D成像方法对计算资源需求更低，并且便于程序的测试和维护.在过去的二三十年内，2-D全波形反演方法在工业界的发展较为成熟(Virieux and Operto, 2009; Wu et al., 2014)，但在天然地震学中的应用还处于探索阶段.随着全球数字地震台网的不断建设，越来越多的密集线性台阵完成布设，尤其是在中国华北克拉通地区.例如中国地震局地球物理研究所布设的“华北科学地震台阵”(North China Seismic Array, 缩写为NCSA)(陈顒等, 2005)，以及中国科学院地质与地球物理研究所布设的“华北地区地球内部结构探测计划”(North China Interior Structure Project, 缩写为NCISP)地震台阵(朱日祥等, 2011; Zheng et al., 2017).基于它们记录的波形数据，线性台阵伴随方法开始得到发展.Tong等(2014)提出了基于线性台阵远震体波的2-D伴随层析成像方法，首先发展了频率波数域-谱元(FK-SEM)混合方法：在区域场外使用FK算子计算1-D介质的平面波响应，在区域场内再使用谱元法模拟远震体波在2-D复杂介质中的传播，因此在保持高效计算效率的同时还能精确地捕捉到线性台站下方非均匀结构体对远震体波的影响，同时混合方法能产生和Moho面相关的各种震相(及尾波)，对速度间断面也有很好的约束.Zhang等(2018)提出了基于线性台阵背景噪声伴随层析成像方法，并成功应用到中国地震局地球物理所在华北地区布设的一条线性台阵(48个地震台)实际成像工作.和传统的背景噪声成像结果相比，该方法获得的模型分辨率得到明显提升，同时仅利用196个CPU核线程，耗时4.5个小时便完成14次迭代，显示出计算的高效性.
如前文所述，由于背景噪声和远震体波数据的频率成分不同，对速度结构敏感性存在差异，仅仅利用其中单一数据的伴随方法同样存在不足.具体表现为：背景噪声伴随层析成像虽然相比传统方法分辨率有所提高，但总体而言，频率成分较低的面波只能获得平滑的S波速度结构，无法有效地约束速度间断面形态.远震体波伴随层析成像，能恢复更多精细的结构，同时能约束精细的速度间断面，但其对初始模型高度依赖，反演存在强烈的非唯一性(Tong et al., 2014).基于此，在伴随方法理论下发展背景噪声和远震体波数据联合成像方法十分必要.首先，两类数据能够照明地下不同深度的区域，具有高度互补性：背景噪声对浅部S波速度结构敏感，主要约束地壳波速结构；远震体波从岩石圈底部近垂直入射，可以约束下地壳及上地幔波速结构.其次，面波频率较低，能提供较平滑的背景速度，减小伴随层析成像对初始模型的高度依赖性；而远震体波的频率成分较高，能恢复精细的结构，提高成像的分辨率，同时在界面产生的转换波对速度间断面也起很好的约束.
基于此，利用两类数据的高度互补性，本文提出基于线性台阵背景噪声和远震体波数据的联合伴随层析成像方法，同时约束台阵下方S波速度结构及速度间断面形态.并建立符合华北克拉通岩石圈典型结构特征的三个理论模型开展成像测试，验证了方法的可行性和稳定性.
1 线性台阵数据伴随层析成像原理 1.1 背景噪声面波伴随层析成像根据Zhang等(2018)提出的线性台阵背景噪声面波伴随层析成像方法，面波深度敏感核随周期变化，其目标函数χS可以表示为多个周期频段面波数据的加权总和(张超，2018)，如下所示：
|
(1) |
其中m表示S波速度模型，Nc表示周期频段个数，Mi表示某频段下波形走时差测量总数，ΔTj(ω)表示经验格林函数EGFs与模拟的格林函数SGFs波形走时差，σj(ω)表示测量误差.在伴随状态法中，目标函数的扰动量和纵波、横波波速扰动量δlnα、δlnβ以及密度扰动量δlnρ线性相关，表达式如下：
|
(2) |
面波对横波速度最为敏感，但短周期数据对浅部的纵波速度与密度也有一定的敏感性.本研究采用Brocher(2005)提出的经验关系式获得纵波速度和密度，可推导出：
|
(3) |
|
(4) |
χn[α]和χn[ρ]表示拟合多项式的系数，将公式(3)和(4)代入公式(2)，可以得到
|
(5) |
KρS, KαS和KβS分别表示面波频散对密度(ρ)、纵波速度(α)和横波速度(β)的敏感核函数.可以看出最终构建的联合敏感核，既包含了横波速度的敏感核，同时也充分考虑了面波频散数据对纵波速度和密度的敏感性.
1.2 远震体波伴随层析成像Tong等(2014)发展了2-D的FK-SEM混合方法用于模拟远震波场，该方法能产生和速度间断面相关的各种转换波震相及尾波，因此能有效地约束速度间断面.在远震体波伴随层析成像中，目标函数B表达式如下：
|
(6) |
其中i表示反演使用的频率，Mi表示该频率下的测量总数，wi(t)表示波形走时差的测量窗口，u和d分别表示模拟和观测的远震体波波形.除了对密度(ρ)、纵波速度(α)和横波速度(β)敏感外，远震体波还对速度间断面(d)有一定的敏感性，因此目标函数的扰动δχB可以表示为
|
(7) |
KρB, KαB和KβB分别表示远震体波对密度(ρ)、纵波速度(α)和横波速度(β)的敏感核函数，KdB表示对速度间断面深度敏感核函数.
2 联合反演思路以上两种数据单独反演的目标函数与敏感核表达式，均建立在伴随方法理论之上.因此可以直接将体波和面波数据整合到同一个波形反演框架中.同时注意到，它们成像过程中使用了不同的正演算子，其中：面波基于谱元法(SEM)，远震体波基于FK-SEM混合方法，因此联合反演首先需要通过各自的正演算子分别产生两种数据的敏感核，然后进行加权建立联合的敏感核函数，最后通过Limited-memory BFGS(简称L-BFGS)(Nocedal, 1980)等优化反演算法对横波速度结构及Moho界面进行迭代更新，具体反演框架如图 1所示.联合反演的目标函数可以表示为
|
(8) |
|
图 1 基于伴随方法的背景噪声面波经验格林函数和远震体波联合成像流程图 Fig. 1 The workflow for joint inversion of ambient noise surface-wave empirical Green′s functions and teleseismic body-wave data based on the adjoint-state method |
其中χS(m)和χB(m)分别对应背景噪声面波和远震体波单独反演的目标函数(公式(1)和公式(6)).相应的联合横波速度及Moho面深度敏感核表达如下：
|
(9) |
其中ω为数据的权重系数，用于调节两种数据的贡献占比，防止最终反演结果被某一数据主导.目前在联合反演方法中，利用“L-curve”方法确定权重系数的做法最为普遍.但在伴随层析成像中，这种需要进行反复大量的权重测试的方法并不是很有效.本文借鉴Fang等(Fang et al., 2016)利用数据误差量化不同数据集权重的思路，定义权重参数：
|
(10) |
其中σS和σB分别表示反演前面波和体波两种数据的波形走时残差测量的标准差，NS和NB分别表示各自数据的测量总数.利用数据和质量控制权重系数，从而避免人为的选择.另外，由于面波数据对速度间断面约束非常有限，联合反演中界面信息仅由远震体波数据进行约束.Moho面深度敏感核表达如下：
|
(11) |
为了测试联合反演方法的有效性，本节设计了三个模型进行成像测试.
3.1 地壳低速体模型首先设计一地壳低速体模型作为真实模型来验证方法的可靠性，即在典型的地壳-地幔双层模型中置入一个长方形低速异常体(图 2a).其中双层模型尺寸为400 km×100 km，上层为地壳，下层为地幔，黑色虚线表示Moho面，具体的模型参数设置见表 1.低速异常体的尺寸为100 km×30 km, 相对于背景地壳S波速度的扰动为6%.沿着地表设置21个线性地震台站(图 2a中的三角形)，台间距为10 km.为了体现联合反演的优势，我们分别进行了三种类型的伴随层析成像测试用于对比：(1)远震体波数据单独成像(后文统称为体波反演)；(2)背景噪声数据单独成像(后文统称为面波反演)；(3)远震体波和背景噪声数据联合成像(后文统称为联合反演).为了减小反演对初始模型的依赖，我们将真实模型移除低速异常体后的地壳-地幔双层模型作为初始模型.表 2分别给出了体波和面波反演中震源参数的设置，具体的在体波反演中，设置8个远震事件，P波入射角分别为16°，20°，24°，28°，332°，336°，340°，344°，这些平面波入射角严格控制在[-28° 28°]区间范围内从而满足远震体波的属性，正负值分别表示远震波场从模型左下角和右下角入射进模型.震源子波类型设置为Ricker子波(中心频率：0.4 Hz).在面波反演中，我们利用了8个台站作为虚拟震源(从左边第一个台站起，每间隔2个台站分布，具体坐标见表 2).子波设置为高斯型(中心频率：0.2 Hz).面波反演使用四个周期频段的数据，分别滤波到：6~15 s，10~20 s，15~30 s，和20~35 s；体波反演使用单频率的数据，滤波到0.4 Hz.对以上两种数据均不添加噪声，根据前文联合反演中权重因子的计算方法(公式(10))，得到ω为0.4，三种类型的伴随层析成像均采用L-BFGS反演方法(Nocedal, 1980).
|
图 2 三种类型伴随层析成像第1次迭代的S波速度敏感核函数及反演模型对比. (a)真实模型：地壳低速体模型，黑色虚线表示Moho面，上下层分别为地壳和地幔，模型参数参见表 1；(b)-(d)分别表示三种类型反演第1次迭代对应的S波速度敏感核函数；(e)-(g)分别表示三种类型反演第1次迭代产生的S波速度模型 Fig. 2 The S-wave velocity sensitivity kernel and inverted model of the 1st iteration for 3 types of adjoint inversion. (a) True model: crust-over-mantle model with low-velocity anomaly in the crust; the black dashed line gives the Moho interface; the model parameters are given in Table 1; (b)-(d) The S wave velocity sensitivity kernel of the 1st iteration for 3 types of inversion; (e)-(g) The inverted S wave velocity model of the 1st iteration for 3 types of inversion. |
|
表 1 地壳-地幔双层模型参数 Table 1 Parameters in the two layer crust-mantle model |
|
|
表 2 体波和面波反演中震源参数设置 Table 2 Source parameter settings for body wave and surface wave tomography |
为了直观展示联合反演的优势，图 2给出了三种类型反演第1次迭代的S波速度敏感核函数(图 2b-d)及反演模型(图 2e-g).从图中可以看出，黑色方框内敏感核数值均为正数，表示模型需要向更低的速度进行更新从而逼近真实模型，这与真实模型中此处存在低速异常体的情况相符.另外，它们的敏感核各具特点：(1)在体波反演中，混合方法产生的远震波场近垂直入射到模型中，因此垂向分辨率比较高，敏感核函数的分布很好地对应了低速异常体的位置(图 2b).但远震体波的频率较高，遇到界面和异常体时容易发生散射等情况，在反演过程中会带来一些高频的假象，例如在Moho面位置及台站附近(图 2e).这种现象在Tong等(2014)研究中也能观察到，可能使反演陷入局部最小值；(2)在面波反演中，由于震源位于地表加之面波存在较强的频散特性，它的敏感核函数在横向上的分辨率较高(图 2c).和体波相比，它的敏感核函数比较平滑且主要分布在地壳浅层；同时面波的频率较低，很少发生散射等效应，因此更新的模型中不会产生高频的假象(图 2f)，这样的特性刚好和体波互补；(3)联合反演兼具了两种数据的优势，其敏感核函数(图 2d)既保留了体波垂向方向的刻画能力，又保留了面波在横向上的分辨率.和体波反演的结果(图 2e)对比，整合了面波数据的联合反演，能提供较平滑的背景速度，因此Moho面位置及台站附近的高频假象得到了很好的压制(图 2g).和面波反演的结果(图 2f)对比，整合了体波数据的联合反演，分辨率得到了极大的提高，反演获得的模型(图 2g)更接近真实模型.
随着迭代次数的增加，模型逐步得到更新.图 3a-c分别显示了联合反演第1，4，8次迭代后获得的S波速度模型.可以看出联合反演方法仅通过8次迭代已经获得和真实模型相近的结果.作为对比，分别使用体波和面波数据单独反演获得的模型如图 3d和图 3e所示.和之前的结论类似，相比于单独反演，联合反演得到的模型具有更高的分辨率，能更好地恢复异常体结构，同时在界面处及浅层处的假象压制更好.
|
图 3 三种类型伴随层析成像最终反演模型对比. (a)-(c)不含噪声合成数据联合反演第1，4，8次迭代更新获得的S波速度模型；(d)不含噪声的体波反演的最终速度模型；(e)不含噪声的面波反演的最终S波速度模型; (f)含噪声合成数据联合反演成像结果：其中面波、体波噪声水平分别为2%和5%; (g) (a)-(f)速度模型中X=200 km处的S波速度垂直剖面 Fig. 3 The final S wave velocity models for 3 types of inversion. (a)-(c) The inverted model after 1、4、5 and 8 iterations for joint inversion using synthetic data without noise; (d) The final inverted model for body wave tomography without noise; (e) The final inverted model for surface wave tomography without noise; (f) The inverted model for joint inversion with 2% random noise in surface waves and 5% random noise in body waves; (g) The vertical profile of S wave velocity located at X=200 km in the inverted model in (a)-(f) |
为了显示反演前后波形拟合的情况，我们以台站STA(图 2a中红色三角形)为例，展示了不含噪声合成数据的联合反演第1次，第4次，第8次迭代时，面波和体波的波形拟合(图 4a-f).可以看出随着迭代的增加，面波和体波都向真实观测数据逐渐逼近，在第8次迭代后基本完全重合.同时我们给出了体波、面波及联合反演的数据误差下降曲线(图 4g).可以看出，相比于单独反演需要的12次迭代，联合反演仅仅8次迭代后就能达到稳定，并且数据残差更小.这是因为联合反演利用了两种数据的高度互补性，面波成分提供了平滑的背景速度模型，降低了体波反演中强烈的非唯一性，从而加快反演的进程.综上分析，我们可以认为，联合反演不仅能更好地恢复S波速度异常体形态，而且能加快反演的收敛速度，具有更好的计算效率.
|
图 4 联合反演(不含噪声数据)第1，4，8次迭代波形拟合及数据误差下降曲线. (a)(c)(e)面波波形拟合; (b)(d)(f)体波波形拟合, 蓝、红线分别表示观测和理论波形；(g)数据误差下降曲线, 黑色、蓝色和红色五角星分别代表体波、面波和联合反演 Fig. 4 Waveform fitting and reduction of data misfit for joint inversion without noise in the data. (a)(c)(e) Waveform fitting for surface waves; (b)(d)(f) Waveform fitting for teleseismic body waves; The blue and red waveforms are the observed and synthetic waveforms, respectively; (g) Reduction of data misfit values as a function of iteration number; The black, blue, and red stars represent results for body wave, surface wave, and joint inversion, respectively |
为了测试联合反演方法的稳定性，我们对图 4中的合成数据加入了高斯随机噪声进行了联合反演成像的测试.根据实际研究中的数据质量，面波数据和远震体波数据的噪声水平分别设置为2%和5%.联合反演获得的S波速度结构如图 3f所示，可以看出反演结果和不含噪声数据的反演结果类似(图 3c), 仍然恢复出了真实模型(图 2a).我们提取出各自模型中x=200 km处的S波速度垂直剖面曲线进行对比(图 3g)，均能较好地拟合.因此测试结果表明，该方法即使在加入随机噪声的真实情况中，也可以获得较为合理和稳定的反演结果.
3.2 含起伏Moho面的地壳-地幔模型为了验证联合反演方法能有效反演S波速度的同时还能刻画速度间断面的能力，我们设计了含起伏界面的地壳-地幔模型.如图 5a所示，黑色实线表示Moho面的位置，Moho面在X=130 km和X=270 km处发生突降，深度由35 km逐渐递增到60 km，形成梯形状形态.同时将Moho面变化处的S波速度扰动设置为6%.初始模型设置为Moho面深度固定在35 km(图中白色虚线)的地壳-地幔模型双层模型.震源和台站参数与3.1节测试中的设置保持一致.
|
图 5 含起伏Moho界面的地壳-地幔模型联合反演. (a)真实模型，黑色实线表示Moho面；(b)-(d)分别为第1，4，7次迭代获得的S波速度模型; (e)-(g)分别为第9，12，15次迭代获得的S波速度模型；白色虚线表示更新的Moho界面 Fig. 5 Joint inversion for the crust-mantle model with an undulating Moho interface. (a) True model, with the black solid line showing the location of Moho interface; (b)-(d) The inverted S wave velocity model after 1、4 and 7 iterations; (e)-(g) The inverted S wave velocity model after 9、12 and 15 iterations; The white dashed line denotes the updating Moho interface |
联合反演过程中对体波和面波反演的策略有所差别.对于远震体波，本文采用Tong等(2014)的多尺度反演方法，即按照低频数据到高频数据分步反演的策略，从而减小反演中的非线性化程度.第一阶段(Step-1)先采用低频数据(0.4 Hz)进行反演；当数据拟合较好时，再加入高频的数据(0.8 Hz)，并利用第一阶段反演的结果作为初始模型，此为第二阶段(Step-2)；对于面波，在以上两个阶段中均采用Zhang等(2018)中使用4个周期频段数据(6~15 s, 10~20 s, 15~30 s, 20~35 s)进行加权反演的策略，从而有效约束0~60 km深度范围的波速结构.基于以上思路，我们同时对S波速度以及速度间断面进行了反演.图 5b-d给出了第一阶段对应的联合反演第1，4，7次迭代更新的模型，图 5e-g给出了第二阶段对应的联合反演第9，12，15次迭代更新的模型; 其中白色虚线对应每次更新的Moho面.可以看出，第一阶段从第1迭代起，Moho面就开始得到更新，通过7次迭代就接近真实界面位置(图 5d)，体现了远震体波对速度间断面较强的约束能力.同时，低速异常体开始进行恢复，由于现阶段使用的面波和低频体波，反演获得的模型比较平滑，尤其在Moho面起伏变化的边界处速度扰动值离真实值6%仍有一定的偏差.尽管如此，第一阶段获得的平滑模型(图 5d)刚好为第二阶段的反演提供了可靠的初始模型.第二阶段加入高频体波数据后，界面处的速度扰动很快接近了6%，边界的轮廓也越来越清晰(图 5e)，同时Moho面进一步向真实形态收敛.在第15次迭代后，低速异常体和Moho面形态的恢复和真实模型基本保持一致(图 5g).
图 6展示了联合反演第一阶段(第1，7次迭代)和第二阶段(第9，15次迭代)，台站STA的面波和体波的波形拟合情况.可以看出：第一阶段面波很快就得到了拟合(图 6c)，这是因为在反演的初期面波占主导地位.同时低频的远震体波也实现了较好的拟合(图 6d).第二阶段由于加入高频成分，显示体波波形仍然有很大的拟合空间(图 6f).在第15次迭代后，高频的体波也得到了很好的拟合(图 6h).同理，图 7给出了联合反演的误差下降曲线，可看出两个阶段的误差随着迭代的进行都有显著的下降.注意到第9次迭代时，误差有很明显的突升，这是因为反演中开始加入高频体波成分.
|
图 6 联合反演第一阶段和第二阶段反演波形拟合.蓝色实线表示观测数据，红色实线表示理论正演数据 Fig. 6 Waveform fitting of the 1st and 2nd step for joint inversion. The blue solid line denotes the observed data, while the red solid line denotes the synthetic data |
|
图 7 联合反演数据误差下降曲线.蓝色和红色五角星分别代表第一阶段和第二阶段反演 Fig. 7 Reduction of data misfit values as a function of iteration number for the 1st step (blue star) and the 2nd step (red star) |
通过以上分析，说明联合反演方法对于含起伏界面的模型，无论对速度结构的恢复还是界面矫正，都能起到很好的约束.
3.3 中国华北克拉通地壳-地幔模型华北克拉通(North China Craton，缩写为NCC)作为稳定古老大陆遭受破坏的典型实例，是目前国际地球科学研究的热点(朱日祥等, 2011, 2012; 陈凌等, 2010).为更深入认识该区域地质构造和演化过程，获得精细且全面的华北克拉通岩石圈内部结构成为地震学家研究的重点.为展示本文中联合反演方法应用到华北岩石圈结构成像上的潜力，根据Zhang等(2018)的S波成像结果，我们建立了具有典型华北地区岩石圈特征的地壳模型.如图 8a所示，该模型在中下地壳设置一椭圆状的低速异常体，相对于背景模型的扰动量为6%.低速体下方设计为速度递增的异常体，从而造成Moho面在X=200 km和X=350 km处出现突降.呈现出Moho面东高西低的分布变化，最大的深度差异达30 km，较好的符合华北克拉通东部岩石圈被破坏，西部保持稳定的地质特征.
|
图 8 华北克拉通地壳-地幔模型联合反演. (a)真实模型，中下地壳中存在一个低速体，相对速度扰动6%，低速体下方在X=200 km和X=350 km之间出现Moho面突降；为了分析联合反演方法对初始模型的依赖性程度，测试了两种不同的初始模型; (b)含起伏Moho面的地壳-地幔双层模型(简称model-1)；(e)含水平Moho面的地壳-地幔双层模型(简称model-2); (c) (f)分别为基于model-1和model-2的联合反演结果；(d) (g)分别为基于model-1和model-2的面波单独反演结果 Fig. 8 Joint inversion for the typical crust-mantle model of North China Craton. (a) True model, where a low velocity anomaly locates in the mid-lower crust with 6% velocity perturbation and there is a Moho depression between X=200 km and X=350km; In order to test the impact of initial model selection on the final model, we set up two different initial models: (b) the crust-mantle model with undulating Moho (model-1); (e) the crust-mantle model with flatting Moho (model-2); (c) (f) The final model for joint inversion based on model-1 and model-2; (d) (g) The final model for surface wave tomography based on model-1 and model-2 |
由于该模型区域较大，也相对复杂，我们增加了远震事件数和台站数量，在入射角[-28° 28°]区间内按照2°的间隔产生16个事件.台站数量增加至41个，从50 km到450 km沿地表分布，间距为10 km.前人研究表明，无论背景噪声还是远震体波伴随层析成像对初始速度模型都有很强的依赖性，为了分析联合反演方法对初始模型的依赖程度，我们设计了两种不同的初始模型：(1)含起伏Moho面的地壳-地幔双层模型(model-1，见图 8b)；(2)含水平Moho面的地壳-地幔双层模型(model-2，见图 8e).和上一节的反演思路相同，仍然采用多尺度反演方法，反演结果见图 8c和图 8f；作为对比，我们还显示了使用面波单独反演的结果，见图 8d和图 8g.我们定义真实和反演模型的相对误差为
|
其中x, z表示格点的空间坐标，mtrue(x, z)和msyn(x, z)分别表示真实模型和反演模型，根据相对误差再计算出模型的RMS误差, 从而定量评估反演结果的可靠程度.首先分析model-1作为初始模型的结果，因为Moho面形态已经和真实模型相近，降低了反演的难度.因此即使使用面波反演，低速异常以及下方的速度递增异常体均有一定程度的恢复.但由于面波频率较低，无法完整地勾勒出低速体的形态，对速度间断面也无法进行很好的约束(图 8d).而联合反演的成像结果(图 8c)有了显著地提高，反演模型的RMS误差由面波反演的12.63%下降到5.41%(见表 3).首先低速异常体的形态得到了恢复，其次低速异常体下方的速度递增的异常体也和真实模型接近，更为显著地是，不同于面波的模糊界面，联合反演精细地刻画出了Moho面的真实形态.model-2的Moho面设置为水平，增加了反演的难度.由于面波对速度间断面的约束能力本身有限，因此面波反演结果离真实值偏差更大(模型的RMS误差为15.64%)，无论是低速异常体还是Moho面形态都无法得到很好的约束(图 8g).而联合反演的结果(图 8f)仍然能恢复和真实模型极度相似的速度异常特征及速度间断面(模型的RMS误差下降至6.68%)，这得益于远震体波对速度界面较强的约束能力，也充分说明联合反演方法对初始模型的依赖性降低.
|
|
表 3 数值模拟参数及反演计算时间统计 Table 3 Parameters of numerical modeling and computational time for inversion |
在波形拟合方面，我们展示了以model-1作为初始模型反演前后的波形对比.图 9给出了联合反演第一阶段(第1，10次迭代)和第二阶段(第11，18次迭代)面波和体波波形拟合的情况，图中蓝色实线表示观测数据，红色实线表示理论正演数据.可以看出，随着反演次数的增加，面波和体波的波形逐渐达到拟合，直至反演终止.结合此次反演的误差下降曲线(图 10)，说明反演效果较好.
|
图 9 联合反演第一阶段和第二阶段反演波形拟合.蓝色实线表示观测数据，红色实线表示理论正演数据 Fig. 9 Waveform fitting of the 1st and 2nd step for joint inversion. The blue solid line denotes the observed data, while the red solid line denotes the synthetic data |
|
图 10 数据误差下降曲线.蓝色和红色五角星分别代表第一阶段和第二阶段反演 Fig. 10 Reduction of data misfit values as a function of iteration number for the 1st step (blue star) and the 2nd step (red star) |
本文提出的联合反演方法作为线性台阵背景噪声伴随层析成像方法(Zhang et al., 2018)的拓展，同样保留了计算的高效性.以3.3节的模型1测试为例，该模型尺寸为500 km×100 km (谱元法中网格剖分为：NX=200, NZ=40)，计算步长和步数分别为0.02 s和12000步.面波单独反演共有16个正演和敏感核函数计算.每个虚拟震源使用4个进程(Processors)，一共消耗64个CPU计算资源，完成14次迭代仅需3.65 h.联合反演由于增加了体波数据的敏感核函数计算，CPU耗费增加了一倍(128个)，但计算耗时上并没增加太多，完成18次迭代也仅需4.8 h(见表 3).
理论测试验证了联合反演方法的可行性，对于符合华北地区典型地壳结构的三个模型都进行了较好的恢复.和单独反演结果对比，体现出联合反演的优势：相比于面波反演，联合反演的成像分辨率有显著提升，能约束更多精细的速度结构；相比于体波反演，联合反演对初始模型的依赖程度更低，同时由于整合了低频的面波成分，反演过程中的非线性化程度得到极大的减小.另外，反演过程较为稳定，即使在加入随机噪声的真实情况中，也可以获得合理的反演结果.
联合反演方法除了对速度结构进行更新外，还能约束Moho面精细形态.地震成像领域对Moho面的约束主要依赖于接收函数方法，目前常规的接收函数偏移成像方法目前主要采用一维背景速度模型进行偏移成像，忽略了地下结构的强烈横向非均匀性影响，导致偏移结果对地震学界面的定位不够准确性，同时它对初始速度模型的依赖程度较高.本文提出的联合反演方法由于整合了远震体波数据，产生的和Moho面相关转换震相(及尾波)能有效的约束速度间断面.并且降低了对初始模型的依赖，即使初始模型偏离真实界面较远，联合反演框架也能较好的恢复到真实界面.目前华北地区和世界很多其他地区布设了大量的密集线性地震台阵，为该方法的应用提供了充足的数据.虽然本文给出的是二维背景噪声和远震体波的联合反演，如果能解决计算资源和效率的问题，这一联合反演框架也可直接推广至三维情形.
除了以上的优点，该方法同样存在一些缺陷，例如发生在同一测线方位上的远震事件数量可能有限，当地震事件偏离线性台站时，需要通过一定的矫正方法对事件进行分析，如何评估偏离程度，找到适合的校正方法需要更深入的研究.本文成像测试中面波和远震体波在反演中的权重主要根据数据的拟合误差来评估，但在实际反演时真实的数据误差可能才是影响权重的重要因素，如何选择合适的权重因子需要进一步探讨.另外，当数据缺少低频信号时，联合反演过程中非线性问题仍然可能存在，因此需要找到更合适的反演策略.
5 结论本文发展了线性台阵背景噪声和远震体波联合伴随层析成像方法，利用两类数据的互补性，高效、精准地约束台阵下方不同深度的S波速度结构及Moho面形态.该方法将Tong等(2014)的2-D远震体波波形反演方法和Zhang等(2018)的线性台阵面波数据波形反演方法相结合，建立了体波和面波联合反演的框架.我们将其应用到具有典型华北克拉通岩石圈结构特征的理论模型，测试结果表明联合反演在成像精度和效率方面都具有优势.该方法的提出有望提供一种更为高效精准的线性台阵成像方法，提升岩石圈成像分辨率，并为后续其它类型波形数据的引入提供思路和方法.
致谢 感谢三位评审专家提出的宝贵意见，感谢中山大学的王易博士在文章编写过程中的建议和指导，研究中使用了的2-D谱元法为开源程序包，网址见https://geodynamics.org/cig/software/specfem2d.
Ammon C J, Randall G E, Zandt G. 1990. On the nonuniqueness of receiver function inversions. Journal of Geophysical Research:Solid Earth, 95(B10): 15303-15318. DOI:10.1029/JB095iB10p15303 |
Bao X W, Sun X X, Xu M J, et al. 2015. Two crustal low-velocity channels beneath SE Tibet revealed by joint inversion of Rayleigh wave dispersion and receiver functions. Earth and Planetary Science Letters, 415: 16-24. DOI:10.1016/j.epsl.2015.01.020 |
Beller S, Monteiller V, Operto S, et al. 2017. Lithospheric architecture of the South-Western Alps revealed by multiparameter teleseismic full-waveform inversion. Geophysical Journal International, 212(2): 1369-1388. |
Bozdaǧ E, Peter D, Lefebvre M, et al. 2016. Global adjoint tomography:first-generation model. Geophysical Journal International, 207(3): 1739-1766. DOI:10.1093/gji/ggw356 |
Brocher T M. 2005. Empirical relations between elastic wavespeeds and density in the Earth's crust. Bulletin of the Seismological Society of America, 95(6): 2081-2092. DOI:10.1785/0120050077 |
Chen L, Wei Z G, Cheng C. 2010. Significant structural variations in the Central and Western North China craton and its implications for the craton destruction. Earth Science Frontiers (in Chinese), 17(1): 212-228. |
Chen M, Huang H, Yao H, et al. 2014. Low wave speed zones in the crust beneath SE Tibet revealed by ambient noise adjoint tomography. Geophysical Research Letters, 41: 334-340. DOI:10.1002/2013GL058476 |
Chen M, Niu F L, Tromp J, et al. 2017. Lithospheric foundering and underthrusting imaged beneath Tibet. Nature Communications, 8: 15659. DOI:10.1038/ncomms15659 |
Chen P, Zhao L, Jordan T H. 2007. Full 3D tomography for the crustal structure of the Los Angeles region. Bulletin of the Seismological Society of America, 97(4): 1094-1120. DOI:10.1785/0120060222 |
Chen Y, Zhou H W, Ge H K. 2005. Seismic array in North China. Journal of Geodesy and Geodynamics (in Chinese), 25(4): 1-5. |
Chen Y L, Niu F L. 2016. Joint inversion of receiver functions and surface waves with enhanced preconditioning on densely distributed CNDSN stations:Crustal and upper mantle structure beneath China. Journal of Geophysical Research:Solid Earth, 121(2): 743-766. DOI:10.1002/2015JB012450 |
Fang H J, Zhang H J, Yao H J, et al. 2016. A new algorithm for three-dimensional joint inversion of body wave and surface wave data and its application to the Southern California plate boundary region. Journal of Geophysical Research:Solid Earth, 121(5): 3557-3569. DOI:10.1002/2015JB012702 |
Fichtner A, Van Herwaarden D P, Afanasiev M, et al. 2018. The collaborative seismic earth model:generation 1. Geophysical Research Letters, 45(9): 4007-4016. DOI:10.1029/2018GL077338 |
Galetti E, Curtis A. 2012. Generalised receiver functions and seismic interferometry. Tectonophysics, 532-535: 1-26. DOI:10.1016/j.tecto.2011.12.004 |
Julià J, Ammon C J, Herrmann R B, et al. 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophysical Journal International, 143(1): 99-112. DOI:10.1046/j.1365-246x.2000.00217.x |
Komatitsch D, Liu Q Y, Tromp J, et al. 2004. Simulations of ground motion in the Los Angeles Basin based upon the Spectral-Element Method. Bulletin of the Seismological Society of America, 94(1): 187-206. DOI:10.1785/0120030077 |
Li J T, Song X D, Zhu L P, et al. 2017a. Joint inversion of surface wave dispersions and receiver functions with P velocity constraints:Application to Southeastern Tibet. Journal Geophysical Research:Solid Earth, 122(9): 7291-7310. DOI:10.1002/2017JB014135 |
Li Y H, Wang X C, Zhang R Q, et al. 2017b. Crustal structure across the NE Tibetan Plateau and Ordos Block from the joint inversion of receiver functions and Rayleigh-wave dispersions. Tectonophysics, 705: 33-41. DOI:10.1016/j.tecto.2017.03.020 |
Lin F C, Moschetti M P, Ritzwoller M H. 2008. Surface wave tomography of the western United States from ambient seismic noise:Rayleigh and Love wave phase velocity maps. Geophysical Journal International, 173(1): 281-298. DOI:10.1111/j.1365-246X.2008.03720.x |
Liu C L, Chen H P, Xie J. 2018. Progress in the studies of the joint inversion of surface wave dispersion and receiver functions. Progress in Geophysics (in Chinese), 33(2): 479-488. DOI:10.6038/pg2018BB0189 |
Liu Q, Gu Y J. 2012. Seismic imaging:From classical to adjoint tomography. Tectonophysics, 566-567: 31-66. DOI:10.1016/j.tecto.2012.07.006 |
Liu Q Y, Li Y, Chen J H, et al. 2010. Joint inversion of receiver function and ambient noise based on Bayesian theory. Chinese Journal of Geophysics (in Chinese), 53(11): 2603-2612. |
Liu S L, Yang D H, Dong X P, et al. 2017. Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling. Solid Earth, 8(5): 969-986. DOI:10.5194/se-8-969-2017 |
Monteiller V, Chevrot S, Komatitsch D, et al. 2013. A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model. Geophysical Journal International, 192(1): 230-247. DOI:10.1093/gji/ggs006 |
Nocedal J. 1980. Updating quasi-Newton matrices with limited storage. Mathematics of Computation, 35(151): 773-782. DOI:10.1090/S0025-5718-1980-0572855-7 |
Owens T J, Zandt G, Taylor S R. 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee:A detailed analysis of broadband teleseismic P waveforms. Journal of Geophysical Research:Solid Earth, 89(B9): 7783-7795. DOI:10.1029/JB089iB09p07783 |
Shapiro N M, Campillo M. 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophysical Research Letters, 31(7): L07614. DOI:10.1029/2004GL019491 |
Shapiro N M, Campillo M, Stehly L, et al. 2005. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715): 1615-1618. DOI:10.1126/science.1108339 |
Shen W S, Ritzwoller M H, Schulte-Pelkum V. 2013. A 3-D model of the crust and uppermost mantle beneath the Central and Western US by joint inversion of receiver functions and surface wave dispersion. Journal of Geophysical Research:Solid Earth, 118(1): 262-276. DOI:10.1029/2012JB009602 |
Tao K, Grand S P, Niu F L. 2018. Seismic structure of the upper mantle beneath eastern asia from full waveform seismic tomography. Geochemistry, Geophysics, Geosystems, 19(8): 2732-2763. DOI:10.1029/2018GC007460 |
Tape C, Liu Q Y, Maggi A, et al. 2009. Adjoint tomography of the southern California crust. Science, 325(5943): 988-992. DOI:10.1126/science.1175298 |
Tarantola A, Valette B. 1982. Generalized nonlinear inverse problems solved using the least squares criterion. Reviews of Geophysics, 20(2): 219-232. DOI:10.1029/RG020i002p00219 |
Tong P, Chen C W, Komatitsch D, et al. 2014. High-resolution seismic array imaging based on an SEM-FK hybrid method. Geophysical Journal International, 197: 369-395. DOI:10.1093/gji/ggt508 |
Tromp J, Tape C, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. Geophysical Journal International, 160(1): 195-216. |
Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367 |
Wang K, Yang Y J, Basini P, et al. 2018. Refined crustal and uppermost mantle structure of southern California by ambient noise adjoint tomography. Geophysical Journal International, 215(2): 844-863. DOI:10.1093/gji/ggy312 |
Wang Y, Chevrot S, Monteiller V, et al. 2016. The deep roots of the western pyrenees revealed by full waveform inversion of teleseismic P waves. Geology, 44(6): 475-478. DOI:10.1130/G37812.1 |
Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model. Geophysics, 79(3): WA13-WA24. DOI:10.1190/geo2013-0294.1 |
Yang Y J, Ritzwoller M H, Lin F C, et al. 2008. Structure of the crust and uppermost mantle beneath the western United States revealed by ambient noise and earthquake tomography. Journal of Geophysical Research:Solid Earth, 113(B12): B12310. DOI:10.1029/2008JB005833 |
Yao H J, Van Der Hilst R D, De Hoop M V. 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-I. Phase velocity maps. Geophysical Journal International, 166(2): 732-744. DOI:10.1111/j.1365-246X.2006.03028.x |
Yao H J, Beghein C, Van Der Hils t R D. 2008. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-Ⅱ. Crustal and upper-mantle structure.. Geophysical Journal International, 173(1): 205-219. DOI:10.1111/j.1365-246X.2007.03696.x |
Zhang C, Yao H J, Liu Q Y, et al. 2018. Linear array ambient noise Adjoint tomography reveals intense crust-mantle interactions in North China Craton. Journal of Geophysical Research:Solid Earth, 123(1): 368-383. DOI:10.1002/2017JB015019 |
Zhang C. 2018. Linear array ambient noise and body wave waveform adjoint tomography: methodology and applications[Ph. D. thesis] (in Chinese). Hefei: University of Science and Technology of China.
|
Zhao L, Jordan T H, Olsen K B, et al. 2005. Fréchet kernels for imaging regional earth structure based on three-dimensional reference models. Bulletin of the Seismological Society of America, 95(6): 2066-2080. DOI:10.1785/0120050081 |
Zheng C, Ding Z F, Song X D. 2018. Joint inversion of surface wave dispersion and receiver functions for crustal and uppermost mantle structure beneath the northern north-south seismic zone. Chinese Journal of Geophysics (in Chinese), 61(4): 1211-1224. DOI:10.6038/cjg2018L0443 |
Zheng T Y, Duan Y H, Xu W W, et al. 2017. A seismic model for crustal structure in North China Craton. Earth and Planetary Physics, 1(1): 26-34. DOI:10.26464/epp2017004 |
Zhu H J, Komatitsch D, Tromp J. 2017. Radial anisotropy of the North American upper mantle based on adjoint tomography with USArray. Geophysical Journal International, 211(1): 349-377. DOI:10.1093/gji/ggx305 |
Zhu H J. 2018. Crustal wave speed structure of North Texas and Oklahoma based on ambient noise cross-correlation functions and adjoint tomography. Geophysical Journal International, 214(1): 716-730. DOI:10.1093/gji/ggy169 |
Zhu L P, Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 105(B2): 2969-2980. DOI:10.1029/1999JB900322 |
Zhu R X, Chen L, Wu F Y, et al. 2011. Timing, scale and mechanism of the destruction of the North China Craton. Science China Earth Sciences, 54(6): 789-797. DOI:10.1007/s11430-011-4203-4 |
Zhu R X, Xu Y G, Zhu G, et al. 2012. Destruction of the North China craton. Science China Earth Sciences, 55(10): 1565-1587. DOI:10.1007/s11430-012-4516-y |
陈凌, 危自根, 程骋. 2010. 从华北克拉通中、西部结构的区域差异性探讨克拉通破坏. 地学前缘, 17(1): 212-228. |
陈顒, 周华伟, 葛洪魁. 2005. 华北地震台阵探测计划. 大地测量与地球动力学, 25(4): 1-5. |
刘成林, 陈浩朋, 谢军. 2018. 面波频散与体波接收函数联合反演研究回顾及展望. 地球物理学进展, 33(2): 479-488. DOI:10.6038/pg2018BB0189 |
刘启元, 李昱, 陈九辉, 等. 2010. 基于贝叶斯理论的接收函数与环境噪声联合反演. 地球物理学报, 53(11): 2603-2612. DOI:10.3969/j.issn.0001-5733.2010.11.008 |
张超. 2018.基于线性台阵背景噪声与体波波形伴随成像方法研究及应用[博士论文].合肥: 中国科学技术大学.
|
郑晨, 丁志峰, 宋晓东. 2018. 面波频散与接收函数联合反演南北地震带北段壳幔速度结构. 地球物理学报, 61(4): 1211-1224. DOI:10.6038/cjg2018L0443 |
朱日祥, 陈凌, 吴福元, 等. 2011. 华北克拉通破坏的时间、范围与机制. 中国科学:地球科学, 41(5): 583-592. |
朱日祥, 徐义刚, 朱光, 等. 2012. 华北克拉通破坏. 中国科学:地球科学, 42(8): 1135-1159. |