2010, Vol. 53
2. Department of Earth, Atmospheric and Planetary Science, MIT, Cambridge, MA 02139, USA
, LI Yu1, CHEN Jiu-Hui1, van er Hilst dR. D.2, GUO Biao1, WANG Jun1, QI Shao-Hua1, LI Shun-Cheng1
2. Department of Earth, Atmospheric and Planetary Science, MIT, Cambridge, MA 02139, USA
1979年Langston[1]证明，在等效震源假定下，从长周期远震P波可以获得台站下方地壳的脉冲响应(接收函数).在此基础上，接收函数方法得到了迅速发展，并在地壳上地幔速度结构研究中获得了广泛的应用.但是，对接收函数反演方法的评价，人们尚有不同的认识.Ammon等[2]根据数值检验的结果认为接收函数反演结果依赖初始模型(非惟一性)，因为接收函数仅对介质速度结构的相对信息(波阻抗)较为敏感.Wu等[3]证明采取恰当的反演方法并非不能正确地估计壳内间断面深度及其速度结构.
事实上，评价接收函数反演对实际数据的解释能力要比模型数值检验复杂得多.地壳结构的复杂性常常使依赖于一维模型的接收函数反演具有多解性.地球物理反问题的本质决定了其必然存在某种程度的不确定性，因为任何测量仅能构成一个有限维的数据空间，而反演寻求的模型参数却具有无限的自由度[4].减少反演非惟一性的出路在于尽可能减小模型参数的自由度.对于接收函数反演来说，与其他独立观测的联合反演无疑应成为接收函数反演研究发展的主要方向.
事实上，Özalaybey等[5]已率先发展了接收函数与面波相速度的线性联合反演方法.Julia等[6]进一步研究了接收函数与面波群速度频散的线性阻尼最小二乘反演.Sung等[7]及Lawrence和Wiens[8]则分别发展了基于遗传算法和基于生态位遗传算法(Niching Genetic Algorithm)的接收函数与面波的非线性联合反演.胡家富等[9]曾利用接收函数与面波群速度的阻尼最小二乘反演研究了滇西地区地壳上地幔速度结构.他们的结果表明:引入面波观测数据对于约束接收函数反演给出的地壳速度参数非常有效.但是，联合反演需要考虑较完整的周期范围.为此，有必要在联合反演中包括周期较短的面波频散数据[6].这意味着接收函数与面波联合反演在某种程度上会受到研究区地震活动性的限制.
近年来，地震环境噪声研究取得了突破性进展.研究表明，利用互相关方法，从散射波场中(包括环境噪声和散射尾波)，不仅可以提取台站间的格林函数，而且可以得到高精度的相速度或群速度频散数据[10~14].即使考虑到噪声源的不均匀性，其引起的相速度误差也能达到小于1%的水平[15].川西台阵环境噪声成像的结果表明[16]，从密集台阵观测的环境噪声数据甚至可以提取周期2s的相速度.这意味着环境噪声数据不仅可以提供地壳尺度的速度结构信息，而且可以对地表浅层结构给予更充分的约束.因此，接收函数与环境噪声的联合反演将为被动源密集台阵地震波形成像开辟新的途径.
接收函数-环境噪声的联合反演与接收函数-面波的联合反演理论上并没有本质的不同.但是，已有的接收函数-面波联合反演存在某些不足:(1)线性反演技术要求尽可能接近真实的地壳初始模型，以至于它并不适于缺乏前期研究的区域; (2)理论上全局搜索方法有助于避免反演陷于局部极小值，但实际上它并不能保证反演收敛到全局最优解; 特别是，它们的效率太低，以至于对密集流动台阵观测数据的解释来说，这类方法常常是难以忍受的.
针对上述问题，本文将在已有的接收函数非线性复谱比反演方法的基础上[17]，发展基于贝叶斯反演理论的接收函数与环境噪声的联合反演方法.它不但有助于形成相应的非线性反演方法，而且可以有效避免低效率的全局搜索.与此同时，我们还将引入对地壳泊松比的全局性搜索，以便获得对地壳泊松比的估计.这对优化接收函数与环境噪声面波频散数据的拟合具有重要价值.我们的讨论将从简要介绍本文的接收函数和瑞雷波相速度的正演算法开始.
2 正演算法 2.1 接收函数接收函数及其微分地震图的正演算法对反演的效率起着决定性的作用.类似文献[17]，我们采用Müller[18]给出的反射率法计算接收函数的复谱比.我们将假定，在接收台站下方，远震P波以平行平面波入射，以避免对慢度的积分运算.Müller[18]给出的反射率法包含了地震波上行和下行的迭代递归算法，这为构筑接收函数微分地震图的高效算法创造了有利条件[17~19].
目前，从远震P波波形数据中分离接收函数大多依据等效震源假定[1].数值实验表明，对于宽频带记录的情况，当接收区上部存在明显间断面时，接收函数的垂直分量并不能简单地看成σ函数[17, 20].克服这个困难最简易的办法是用径向与垂向接收函数的复谱比代替原来意义下的接收函数.除特别声明，本文的接收函数均指径向与垂向接收函数的复谱比.实际上，复谱域与时间域的接收函数是等效的.相对于时间域的接收函数反演，接收函数复谱域反演避免了Fourier变换的重复计算，同时减少了接收函数正演的计算量(相当于数据压缩).接收函数复谱域的正演算法有助于自然形成从低频向高频的扩张式分频段接收函数反演算法.
2.2 相速度频散的正演算法经典面波频散的正演大多为基于Thomson-Haskell矩阵的传播矩阵理论[21~24].其中，Herrmann等[24]的计算程序得到了广泛的认可和应用.计算面波频散的另外一类方法是Kennett[25]及Luco与Apsel[26]发展的广义反射透射系数(R/T)法.虽然广义R/T方法具有在所有频率上保持数值稳定的优点，但它的效率相对较低，以至于在面波成像中应用较少.针对上述不足，Chen[27]对广义R/T方法作了改进，使得理论公式更为简洁，并对高频及厚层情况具有更好的精度和稳定性.在此基础上，Pei等[28, 29]进一步针对面波频散计算提出了快速广义R/T方法.由于无需计算修正R/T矩阵，并直接采用层矩阵求逆的解析解，该方法较原有的广义R/T方法约节省了50%以上的计算时间[28].对于接收函数与面波联合反演来说，快速广义R/T方法为在统一理论体系下构筑高效的计算程序创造了条件.
根据文献[28]，面波相速度可通过求解本征方程(注:文献[27, 28]给出的特征方程有误):
|
(1) |
得到.这里，det表示行列式，E211和E221为分层地壳模型第1层层矩阵的子阵(第1层的上界面位于地表).若E1表示第1层的层矩阵，则有
|
其中，elk为层矩阵的元素(l=1，2，3，4;k=1，2，3，4)，这里，上标表示地壳模型的层序号，下标表示矩阵元素.式(1)中Λu1(0)为的地壳模型第1层深度为零时的相位延迟矩阵，
接收函数与面波的联合反演通常涉及分层更细的地壳模型(通常，模型层数N＞30).因此，层矩阵及其逆矩阵的计算对于面波频散正演算法的效率有很大影响.在Pei等[28]工作的基础上，潘佳铁等[30]曾给出了形式上与其不同的层矩阵逆的解析表达式.与上述研究不同，根据文献[31]，我们则导出了以相速度为独立变量的层矩阵及其逆矩阵的解析公式.对于瑞雷波，层矩阵的元素可写作
|
层矩阵的逆矩阵可写作
|
(2) |
并有
|
这里，α和β分别为介质层的P波和S波速度，ρ为密度，c为相速度，ω为角频率，并有
|
为简明起见，上述各式中均略去了层序号.
不难发现，以相速度为独立变量的层矩阵及其逆矩阵公式不但有助于进一步减少计算量，也便于我们理解相速度与层矩阵及其逆矩阵的关系.由(2)式的矩阵元素可见，当相速度与介质层内的S波速度相同或者非常接近时，层矩阵的逆将会失稳(这是为提高计算效率而付出的代价).为此，在求解式(1)给出的本征方程时，我们须避免这种情况的发生.对于非完全弹性介质，引入复速度参数是必要的.这同时可有效避免上述失稳情况的发生.
3 接收函数和相速度的联合反演基于贝叶斯定理的地球物理反演理论最早由法国Tarantola教授创立[31, 32].贝叶斯反演建立在概率论的基础上，形成了把数据信息与模型先验信息联系起来的理论框架，从而通过模型的先验信息可以约束模型的后验参数.这是贝叶斯反演理论得到广泛应用的一个重要原因[31~36].但是，如何具体确立数据及模型的先验概率密度(PDF)是一个引起争论的问题[33, 34].多年的争论已有效推进了相关问题的研究[32].
对于接收函数与环境噪声数据的联合反演来说，上述问题在某种程度上得到了缓解.这是因为远震及环境噪声观测提供了大量重复测量的机会，以至于对数据的先验信息，我们可以做出客观的参数估计.虽然模型先验概率分布的估计仍然是一个有待解决的问题，但是环境噪声数据的引入大大减少了先验模型参数的自由度.对于一维线性模型空间(不考虑地壳横向非均匀散射及各向异性)的情况，我们将假定观测数据和模型参数的先验PDF符合正态分布.这对我们讨论的问题应是一个足够好的估计.
根据贝叶斯反演理论[32, 33]，模型空间的后验概率密度可表示为
|
(3) |
其中，L(m)称为或然性函数，它反映了模型预测结果与观测数据的拟合程度，ρ(m)为独立于观测数据的先验模型PDF.在同时考虑接收函数与环境噪声数据的情况下，
|
(4) |
并有
|
(5) |
这里，g(m)为正演算子，d为观测数据矢量，m为模型参数矢量，*表示复共轭，T表示转置，而下标RF和SF分别表示接收函数和面波相速度频散，下标p表示先验信息.CRF和CSF分别为接收函数和面波相速度频散的数据协方差矩阵.CM为先验模型的协方差矩阵.
根据式(4)和(5)，我们可以得到相应的目标函数
|
(6) |
而求解反问题意味着要求模型空间后验概率密度σ(m)具有极大值，或者说需要寻求(6)式的极小值.共轭梯度法是优化标函数S(m)的一个较为稳健高效的方法.其关键在于如何计算目标函数的梯度。根据式(6)，我们可得到目标函数的梯度:
|
(7) |
其中，
由式(7)可知，计算目标函数梯度的核心是计算接收函数及面波频散对模型参数的偏导数.对于接收函数微分地震图计算，我们采用了文献[17]的方法.关于面波频散对模型参数偏导数的计算已有多种(包括解析和数值的)方法可供选择.文献[31]和[37]曾分别对此作了较全面的综述.由于面波频散对模型参数的偏导数计算与相应的正演算法密不可分.针对快速广义R/T方法，我们将采用Novotiny [38, 39]的隐函数算法.根据隐函数微分定理，对于给定的周期，面波频散对模型速度参数的偏导数
|
(8) |
这里，
另外，根据贝叶斯反演理论[32, 33]，我们可以导出模型的后验协方差
|
(9) |
而模型空间解的分辨率可用后验分辨率算子
|
(10) |
估计.这里，m∞表示模型空间中的最大或然性点，I为单位矩阵.
理论上，利用预条件共轭梯度法可以加速反演的收敛[17, 32].但是，我们的实际计算表明，对于接收函数与环境噪声数据的联合反演，预条件共轭梯度法并非是必要的，而且当求解的问题具有很强的非线性时，由于问题已经远离正态分布的假定，预条件共轭梯度法很可能是不收敛的.
图 1给出了本文算法的流程框图.我们对图 1可作以下几点说明:
|
图 1 接收函数与环境噪声联合反演的流程框图 Fig. 1 Diagram for the joint inversion of receiver function and ambient noise |
(1)如同文献[17]，本文算法引入了接收函数垂直分量与径向分量的初至振幅比，并包含了逐步扩展的接收函数分频段反演，具体的频段分配是:0~0.5Hz，0~0.75 Hz，0~1.0 Hz，0~1.25 Hz，0~1.5Hz.理论和实际研究表明，接收函数反演包括1.0Hz以上的信息是必要的[17].缺少1.0Hz以上的短周期信息是导致接收函数多解性的重要原因之一.
(2)本文算法包含了对泊松比的全局性扫描，扫描范围为0.2~0.35，扫描步长可以自由调节.但是，对泊松比的调整仅限于接收区壳内的局部区域，而地壳的平均泊松比可通过加权平均得到.这意味着本文给出的地壳平均泊松比是在某些先验假定下给出的估计值.具体地，对P波与S波的速度比，我们规定:在地幔区域为1.8015;在地壳区域(当S波速度大于3.57km/s时)为1.74;当S波速度小于2.5km/s时，P波速度按经验公式α=2.1+1.83(β-0.7)计算.相应的密度模型根据公式ρ=0.32α+0.77估算.这里，α和β分别为P波和S波的速度.上述假定主要基于全球模型和前人研究的结果，并根据我们的实际经验作了调整.对于接收函数和地震环境噪声的联合反演来说，相对于对整个地壳采用单一的泊松比，本文给出的地壳平均泊松比的估计方法应是一个更为合理的办法.
(3)反演收敛判据主要依据接收函数及相速度频散拟合的均方根误差.显然，接收函数和相速度频散具有不同的量纲.它们的预测值和观测值的拟合误差不但在量级上是不同的，而且它们的收敛也是不协调的.为此，可由接收函数及相速度频散拟合的总体方差公式
|
(11) |
和
|
(12) |
估计，并把(11)式的极小值作为收敛判据.其中，式(12)将由迭代的初值给定.这里，vr和vp分别表示接收函数及相速度频散拟合的均方根误差，而avr和avp分别表示观测接收函数振幅谱的算术平均和各周期观测相速度的算术平均.
(4)与接收函数与面波线性联合反演方法不同[6]，本文方法避免使用接收函数与面波之间人为规定的平衡系数，这是贝叶斯反演理论带来的另外一个优点.因为数据协方差的引入使得不同量纲数据的同化变得十分容易.相对于其他反演方法，贝叶斯反演理论的另外一个优点则是数据的协方差协调了不同信度的数据拟合.这对于采用一维速度参数模型的接收函数反演具有特别的优势.
4 数值检验表 1给出了对本文方法数值检验所用的接收区模型.反演中它将被细化为2km层厚的层状模型.壳幔界面的深度假定为40km.接收区以外的地球介质采用PREM模型[40].我们假定震中距为40°，震源深度为20km.与模型相应的“真实”接收函数和瑞雷波相速度频散由本文的正演算法得到.接收函数的先验均方根误差为谱振幅的1%，各周期相速度的先验均方根误差为0.05km/s.
|
表 1 数值检验所用的接收区模型 Table 1 Model parameters in the receiver area used for the numerical tests |
图 2给出了对本文方法的数值检验结果.由图 2可见，在“真实”接收区模型中，地壳和地幔区内均包含有低速层和高速层，相频散的周期为3~30s.我们假定了两个垂向均匀的接收区初始模型:对于第1个初始模型，我们令α=5.239km/s，β=3.009km/s，ρ=2.5g/cm3; 对于第2个初始模型，我们令α=8.82km/s，β=4.90km/s，ρ=2.5g/cm3.这里，α和β分别为介质层的P波和S波速度，ρ为密度.虽然这两个初始模型远离“真实”模型，但图 2表明，本文的方法可以相当好地恢复深度60km以上的接收区速度模型.在更大的深度上，虽然预测模型的误差随着深度增加而变大，但对80km界面深度仍能给出正确的推断.
|
图 2 接收函数与环境噪声联合反演的数值检验 (a)接收区模型，黑色实线为“真实”模型，分别相应于P波和S波速度; 红色虚线分别相应于反演得到的P波和S波速度; 蓝色的点线分别相应于初始的P波和S波速度.(b)接收函数，黑色实线为“真实”的接收函数; 红色虚线为反演得到的接收函数; V和R分别表示垂直分量和径向分量，右边的数字为相关系数.(c)频散曲线，黑色实线为“真实”的瑞雷波的相速度频散曲线(3~30s)，短线表示误差; 红色的圆点表示反演给出的各周期的相速度.上图和下图分别为第1个和第2个初始模型及其检验结果. Fig. 2 Numerical Tests on joint inversion of receiver function and ambient noise (a) Receiver model. Black solid lines represent the "true" model, relevant to P and S velocity; Red dash lines are the inversion model, relevant to P and S velocity; Blue dot lines are the initial model of P and S velocity, respectively.(b) Receiver function. Black lines are "true" receiver functions; red dash lines are the results after inversion; V and R represent the vertical and radial component, respectively; The digits on the right are correlation coefficient.(c) Dispersion curve. Black solid lines with error bar represent the "true" phase dispersion curve of Rayleigh waves (3~30s); The red dots are the phase velocity at each period after inversion. In the upper part and lower part are shown the results by the numerical tests on the first and second initional model, respectively. |
为了进一步考察相频散对联合反演结果的影响，图 3给出了相频散周期为2~40s的数值检验结果.由图 3可见，初始模型为图 2中的“低速”模型.相对于图 2，反演得到的60~80km的速度结构更加接近“真实的”模型参数，表明相频散周期向长周期方向的扩展有助于进一步约束深层的速度结构.另外，与图 2不同的是，对于近地表2 km范围内的速度结构，我们采用了1km层厚的分层.图 3的结果表明，由于包含了周期2s的相频散数据，本文的反演方法可以很好地预测浅表速度结构.实际上，数值检验结果表明，对于相频散数据限于3s以上的情况，浅表结构的反演存在多解性.这是因为在这种情况下，无论接收函数，还是面波数据均不能严格限制浅表结构.这个结果对于指导实际观测数据的解释具有重要参考价值.
|
图 3 接收函数与环境噪声联合反演的数值检验，所用符号与图 2相同. Fig. 3 Numerical Tests on joint inversion of receiver function and ambient noise. The notation used is same with those in Fig.2. |
图 4给出了反演得到的接收函数及瑞雷波相速度拟合误差随地壳泊松比的变化.反演采用的初始模型为图 2中的“高速”模型.图 4表明:无论接收函数，还是相速度频散，预测与“观测”数据的拟合程度均具有单峰极值; 反演得到的地壳平均泊松比为0.2648，而“真实”模型(表 1)的地壳平均泊松比为0.2639.其相对误差约为0.3%.
|
图 4 接收函数(b)及瑞雷波相速度(a)拟合误差随地壳泊松比σ的变化 Fig. 4 Variations of fitting errors of the receiver function (b) and phase velocity (a) with averaged Poisson′s ratio over the crust |
2006~2009年中国地震局地质研究所地震动力学国家重点实验室在川西地区(26°N~32°N，100°E~105°E)布设了由297台站组成的密集流动宽频带观测台阵(以下简称川西台阵)[41].刘启元等[42]曾利用接收函数方法研究了川西台阵沿31°N线台站下方地壳上地幔速度结构.李昱等[16]给出了川西台阵环境噪声面波层析成像的结果.这为我们进一步检验本文方法创造了条件.我们将直接采用他们给出的相频散数据，进行接收函数与环境噪声的联合反演.有关环境噪声数据的处理可见文献[16].
作为一个例子，我们利用本文的反演方法重新研究了川西台阵台站KWC05(位于31°N线内)下方100km深度范围内的地壳上地幔S波速度结构.图 5给出了台站KWC05不同方位的时间域接收函数.与文献[42]不同，图 5给出的接收函数是采用迭代反褶积方法由55个远震事件的P波波形数据得到[43].其震中距范围为35°~60°，台站方位角(station azimuth)为46°~207°.需要说明的是，图 5所示的接收函数是按方位分组叠加的结果.每组接收函数的方位角变化限制在1°以内.接收函数的分组叠加有助于等权重地压制不同方位介质横向非均匀引起的散射.
|
图 5 台站KWC05不同方位的接收函数 (a)接收函数，Baz为台站方位角(单位:(°))，SUM表示叠加平均; (b)台站位置，LMS表示龙门山断裂带. Fig. 5 Receiver functions at Station KWC05 over different azimuths (a) Receiver function; Baz denotes the back azimuth (°); SUM denotes the summation over azimuths.(b) Station map; LMS denotes the Longmen Shan faults. |
图 6给出了利用本文方法得到的反演结果.图 6表明，无论接收函数，还是环境噪声的相频散数据都得到了很好的解释.台站KWC05下方100km深度范围内的地壳上地幔S波速度结构表明，其地壳厚度为44km，上地壳为高速结构，而深度24~42km的中下地壳形成了相对的低速区.我们的重新解释结果表明，台站KWC05下方地壳的平均泊松比为0.262，壳内低速区泊松比为0.270.
|
图 6 台站KWC05下方地壳上地幔速度结构反演 所用符号与图 2相同.环境噪声相速度频散的周期为2~35s.左图:黑色实线为接收函数与环境噪声联合反演的结果，蓝色虚线为接收函数单独反演的结果. Fig. 6 Inversion of the crust and upper mantle velocity structure beneath Station KWC05 The notations are same with those in Fig.2. The period of the phase dispersion of ambient noise is 2~35s. Left side:the black solid line shows the crust and upper mantle velocity structure by the joint inversion, and blue dash line shows that obtained by the receiver function inversion. |
图 6中同时给出了利用非线性接收函数复谱比反演得到的地壳上地幔S波速度结构[42].两者比较的结果表明它们之间存在明显的差别.据此可以推断，环境噪声的相频散数据对于限定地壳及上地幔顶部的S波的绝对速度发挥了重要作用.事实上，在单独反演接收函数的情况下，判断反演结果的惟一依据只能是接收函数的波形拟合程度，以至于在人们追求不断改善波形拟合时并不能判断介质横向变化对接收函数的影响.
6 讨论和结论在已有的接收函数非线性复谱比反演方法的基础上[17]，本文进一步发展了基于贝叶斯反演理论的接收函数与环境噪声的非线性联合反演方法.数值检验和实际数据解释的结果表明:
(1)本文方法给出的反演结果与初始模型无关.即使对台站下方地壳速度结构仅有非常粗略的估计，利用本文方法仍可对模型参数作出很好的预测，表明环境噪声的相频散数据可以对地壳的S波绝对速度参数提供有力的约束.因此，对于环境噪声与接收函数的联合反演来说，全局反演方法并非是绝对必要的.
(2)本文方法可同样应用于地震面波和接收函数联合反演.但是，由于环境噪声能够提供与接收函数频带范围大体相同的面波频散数据，环境噪声与接收函数的联合反演能够更好地约束台站下方地壳上地幔，特别是近地表的速度结构.当环境噪声相速度频散的周期范围为2~40s时，环境噪声与接收函数的联合反演能够可靠地预测台站下方0~80km深度范围的S波速度结构，浅表速度结构的空间分辨率可达到1km.
(3)与已有的地震面波-接收函数联合反演方法不同，本文方法同时引入了地壳泊松比的全局性搜索.这样，我们不但可以更合理地解释接收函数和环境噪声的面波频散数据，改善数据拟合的质量，而且可在某些先验假定下同时获得地壳泊松比的估计.
(4)对川西台阵KWC05台站观测数据的联合反演结果表明，该台站下方地壳厚度为44km，上地壳具有明显的高速结构，而在24~42km深度范围的中下地壳具有低速结构.其地壳的平均泊松比为0.262，壳内低速带的泊松比为0.270.
(5)对于地壳结构包含低速层的情况，介质的衰减吸收参数对环境噪声与接收函数联合反演具有重要影响.在实际数据解释中，本文对介质衰减吸收参数的估计采用了简单的尝试法.进一步同时反演介质衰减吸收参数对于本文的方法来说并不存在原则的困难.
(6)环境噪声与接收函数的联合反演明显优于单独的接收函数反演，也优于地震面波与接收函数联合反演.但是，限于环境噪声的观测条件，本文方法主要适用于密集地震台网或台阵观测系统下方高分辨率地壳上地幔速度结构研究.
| [1] | Langston C A. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. , 1979, 84: 4749-4762. DOI:10.1029/JB084iB09p04749 |
| [2] | Ammon C J, Randall G, Zandt G. On the nonuniqueness of receiver function inversions. J. Geophys. Res. , 1990, 95(B10): 15303-15318. DOI:10.1029/JB095iB10p15303 |
| [3] | Wu Q, Li Y, Zhang R, et al. Wavelet modeling of broad-band receiver functions. Geophys. J. Int. , 2007, 170(2): 534-544. DOI:10.1111/gji.2007.170.issue-2 |
| [4] | Scales J A, Snieder R. The anatomy of inverse problems. Geophysics , 2000, 65(6): 1708-1710. DOI:10.1190/geo2000-0001.1 |
| [5] | Özalaybey S, Savage M K, Sheehan A F, et al. Shear-wave velocity structure in the northern Basin and Range province from the combined analysis of receiver functions and surface waves. Bull. Seism. Soc. Am. , 1997, 87(1): 183-199. |
| [6] | Julia J, Ammon C J, Herrmann R B, et al. Joint inversion of receiver function and surface wave dispersion observations. Geophys. J. Int. , 2000, 143(1): 1-19. DOI:10.1046/j.1365-246x.2000.00197.x |
| [7] | Sung J C, Chang-Eob B, Langston C A. Joint analysis of teleseismic receiver functions and surface wave dispersion using the Genetic algorithm. Bull. Seism. Soc. Am. , 2004, 94(2): 691-704. DOI:10.1785/0120030110 |
| [8] | Lawrence J F, Wiens D A. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm:application to Patagonia. Bull. Seism. Soc. Am. , 2004, 94(3): 977-987. DOI:10.1785/0120030172 |
| [9] | 胡家富, 朱雄关, 夏静瑜, 等. 利用面波和接收函数联合反演滇西地区壳幔速度结构. 地球物理学报 , 2005, 48(5): 1069–1076. Hu J F, Zhu X G, Xia J Y, et al. Using surface wave and receiver function to jointly inverse the crust-mantle velocity structure in the West Yunnan area. Chinese J. Geophys. (in Chinese) , 2005, 48(5): 1069-1076. DOI:10.1002/cjg2.750 |
| [10] | Lobkis O I, Weaver R L. On the emergence of the Green's function in the correlations of a diffusive field. J. Acoust. Soc. Am. , 2001, 110(6): 3011-3017. DOI:10.1121/1.1417528 |
| [11] | Weaver R L, Lobkis O I. Diffuse fields in open systems and the emergence of the Green's function. J. Acoust. Soc. Am. , 2004, 116(5): 2731-2734. DOI:10.1121/1.1810232 |
| [12] | Campillo M, Paul A. Long-Range correlations in the diffuse seismic coda. Science , 2003, 299(5606): 547-549. DOI:10.1126/science.1078551 |
| [13] | Shapiro N M, Campillo M. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise. Geophys. Res. Lett. , 2004, 31(7): 1615-1619. |
| [14] | Shapiro N M, Campillo M, Stehly L, et al. High-resolution surface-wave tomography from ambient seismic noise. Science , 2005, 307: 1615-1618. DOI:10.1126/science.1108339 |
| [15] | Yao H J, van der Hilst R D. Analysis of ambient noise energy distribution and phase velocity bias in ambient noise tomography, with application to SE Tibet. Geophys. J. Int. , 2009, 179(4): 1113-1132. |
| [16] | 李昱, 姚华建, 刘启元, 等. 川西地区台阵环境噪声瑞利波相速度层析成像. 地球物理学报 , 2010, 53(4): 842–852. Li Y, Yao H J, Liu Q Y, et al. Phase velocity array tomography of Rayleigh waves in western Sichuan from ambient noise. Chinese J. Geophys. (in Chinese) , 2010, 53(4): 842-852. |
| [17] | 刘启元, KindR, 李顺成. 接收函数复频谱比的最大或然性估计及非线性反演. 地球物理学报 , 1996, 39(4): 502–513. Liu Q Y, Kind R, Li S C. Maximal likelihood estimation and nonlinear inversion of the complex receiver function spectrum ratio. Chinese J. Geophys. (in Chinese) , 1996, 39(4): 502-513. |
| [18] | Müller G. The reflectivity method:a tutorial. J. Geophys. , 1985, 58(3): 153-174. |
| [19] | Randall G E. Efficient calculation of differential seismograms for lithospheric receiver functions. J. Geophys. Int. , 1989, 99(3): 469-481. DOI:10.1111/gji.1989.99.issue-3 |
| [20] | Scherbaum F. Seismic imaging of the site response using micro-earthquake recordings, part I:method. Bull. Seism. Soc. Am. , 1987, 77: 1905-1923. |
| [21] | Thomson W T. Transmission of elastic waves through a stratified solid medium. J. Appl. Physics , 1950, 21(2): 89-93. DOI:10.1063/1.1699629 |
| [22] | Haskell N A. The dispersion of surface waves on a multi-layered medium. Bull. Seism. Soc. Am. , 1953, 43(1): 17-34. |
| [23] | Dunkin J W. Computation of modal solution in layered, elastic media at high frequencies. Bull. Seism. Soc. Am. , 1965, 53(2): 335-358. |
| [24] | Herrmann R B, Ammon C J.Computer programs in seismology version 3.20:Surface waves, receiver functions, and crustal structure. St. Louis University, Missouri. 2002, Available at http://mnw.eas.slu.edu/People/RBHerrmann |
| [25] | Kennett B L N. Seismic Wave Propagation in Stratified Media. Cambridge University Press, 1983 |
| [26] | Luco J E, Apsel R J. On the Green's function for a layered half-space, part I. Bull. Seism. Soc. Am. , 1983, 73(4): 909-929. |
| [27] | Chen X. A systematic and efficient method of computing normal modes for multilayered half-space. Geophy. J Int. , 1993, 115(2): 391-409. DOI:10.1111/gji.1993.115.issue-2 |
| [28] | Pei D H, John N L, Pullammanappallil S K, et al. Improvements on computation of phase velocities of Rayleigh waves based on the generalized R/T coefficient method. Bull. Seism. Soc. Am. , 2008, 98(1): 280-287. DOI:10.1785/0120070057 |
| [29] | Pei D H, John N L, Pullammanappallil S K, et al. Erratum to improvements on computation of phase velocities of Rayleigh waves based on the generalized R/T coefficient method. Bull. Seism. Soc. Am. , 2009, 99(4): 2610-2611. DOI:10.1785/0120090003 |
| [30] | 潘佳铁, 吴庆举, 李永华, 等. 基于快速广义反射透射系数方法的面波群速度计算. 地球物理学进展 , 2009, 24(6): 2030–2035. |
| [31] | Aki K, Richards P G. Quantitative Seismology, Second Ed. University Science Books, Sausalito, Californian, 2002 |
| [32] | Tarantola A. Inverse Problem Theory, Methods for Data Fitting and Model Parameter Estimation. Elsevier, Amsterdam, 1987 |
| [33] | Tarantola A. Inverse Problem Theory, and Methods for Model Parameter Estimation. SIAM, PA 19104-2688 USA, 2005 |
| [34] | Scales J A, Snieder R. To Bayes or not to Bayes. Geophysics , 1997, 62(4): 1045-1046. DOI:10.1190/1.6241045.1 |
| [35] | Gouveia W P, Scales J A. Bayesian seismic waveform inversion:parameter estimation and uncertainty analysis. J. Geophys. Res. , 1998, 103(2): 2759-2779. |
| [36] | Gao S. A Bayesian nonlinear inversion of seismic body-wave attenuation factors. Bull. Seism. Soc. Am , 1997, 87(4): 961-970. |
| [37] | Cercato M. Computation of partial derivatives of Rayleigh-wave phase velocity using second-order subdeterminants. Geophys. J. Int. , 2007, 170(1): 217-238. DOI:10.1111/gji.2007.170.issue-1 |
| [38] | Novotiny O. Partial derivatives of dispersion curves of Love waves in layered medium. Studia Geoph. et Geod. , 1970, 14(1): 36-50. DOI:10.1007/BF02585549 |
| [39] | Novotiny O. Methods of computing the partial derivatives of dispersion curves. Pure. Appl. Geophys. , 1976, 114(5): 765-775. DOI:10.1007/BF00875786 |
| [40] | Dziewonski A M, Anderson D L. Preliminary reference Earth model. Phys. Earth Planet. Int. , 1981, 25(4): 297-356. DOI:10.1016/0031-9201(81)90046-7 |
| [41] | 刘启元, 陈九辉, 李顺成, 等. 汶川Ms8.0地震:川西流动地震台阵观测数据的初步分析. 地震地质 , 2008, 30(3): 584–596. Liu Q Y, Chen J H, Li S C, et al. The Ms8.0 Wenchuan earthquake:Preliminary results from the western Sichuan mobile seismic array observations. Seismology and Geology (in Chinese) , 2008, 30(3): 584-596. |
| [42] | 刘启元, 李昱, 陈九辉, 等. 汶川Ms8.0地震:地壳上地幔S波速度结构的初步研究. 地球物理学报 , 2009, 52(2): 309–319. Liu Q Y, Li Y, Chen J H. Wenchuan Ms8. 0 earthquake:S-wave velocity structure of the crust and upper mantle. Chinese J. Geophys. (in Chinese) , 2009, 52(2): 309-319. |
| [43] | Ligorria J P, Ammon C J. Iterative deconvolution of teleseismic seismograms and receiver function estimation. Bull. Seism. Soc. Am. , 1999, 89(5): 1395-1400. |