2011, Vol. 54
地震层析成像是研究地球内部结构和深部孕震环境的重要手段，自20世纪70年代发展以来，地震层析成象方法，包括体波走时成像和面波层析成像，主要都是基于波传播的几何射线理论.对于体波来说，认为波沿连接震源和台站的几何射线路经传播，对于面波来说，认为波沿连接震源和台站的大圆路经传播，这是目前绝大多数层析成象方法的理论假设.然而，射线理论仅仅是波动理论的一个高频近似，在实际中，地震记录通常含有有限的频率成分.其带宽限制表明，波的传播不是沿我们称为射线的无限窄的区域传播，而是扩展到沿射线周围的一个体积有限区域[1, 2]，该区域称为菲涅耳区(Fresnel zone)[3].
以体波走时成像为例，传统基于射线理论的方法只有射线穿过的介质参数扰动影响波的到时，这种扰动造成走时的提前或延迟，沿射线路径将记忆走时的这种改变，直到信号被接收.对于高频的波来说，正是利用这种走时的改变，通过对穿过非均匀体的多条射线进行拟合达到成像目的.对于有限频率的波场，或者非均匀体的尺度和波场可比，波将发生散射和衍射现象，这种现象将造成波前的复原(wavefronthealing)[4, 5]，即随着传播距离的增加，由于散射和衍射，来自非均匀体周围不同方向的波相消干涉，会自动补偿波在非均匀体内部传播引起的走时的提前或延迟，使得走时的改变不会一直记忆到信号被接收.如果此时仍然采用射线理论，将会造成走时测量的错误，从而影响成像结果.相反，偏离射线的小尺度非均匀体，则可能引起走时的改变[6~9]，这种改变在射线理论中并未考虑，这意味着，基于射线理论的方法，不能很好地分辨小尺度的三维非均匀结构.为此，人们提出了有限频率地震波层析成像理论[10~12]，对于速度大致呈梯度变化的地球介质来说，反演中用到的灵敏度核，即表征介质参数扰动对走时扰动敏感程度的物理量，类似香蕉-甜甜圈形状，因此这种理论又被称为香蕉-甜甜圈(banana-doughnut)理论[13, 14].很多研究者给出了有限频率体波走时和波形灵敏度核，并将其用于成像中[15~19]，和射线理论的结果进行了对比.本文主要关注有限频率面波成像问题.
由于面波的多模性和模式之间的耦合，及面波层析成像的二维特征，和体波走时成像不同，有限频率面波的三维灵敏度核，并不是香蕉-甜甜圈形状.针对球状地球模型，在远场和波恩近似下，Zhou(2004)研究了有限频率面波层析成像中的灵敏度核[20] ，并讨论了模式耦合的影响，而后又将这种研究扩展到黏弹性[21]和径向各向异性情形[22].YoshizawaandKennett(2005)则采用另一种方法[23]，基于波恩和Rytov 近似，利用面波势的表示，给出了有限频率面波灵敏度核.对于多次散射情形，Friderich(1999)[24]和Maupin(2001)[25]曾经考虑过各向同性和任意各向异性介质中面波的散射理论，而且包含了近场的影响，目前这种方法还没有用于有限频率面波成像中.另一种处理面波散射的方法是直接求解积分方程[26]，实现面波的全散射.对于具有21 个弹性常数的弱各向异性介质，Chenand Tromp(2007)的研究发现[27]，面波传播主要受13个参数的影响，这13 个参数是21 个独立弹性常数的不同组合.Sieminski等(2007)采用伴随谱元方法[28]，计算了各向异性介质中这13 个参数的面波灵敏度核，并分析了不同参数对不同模式面波影响的权重大小[29].
和体波相比，面波占据了地表附近能量的绝大部分，因此在浅层工程勘探中，面波方法也获得了广泛的应用[30~33].在面波浅层勘探中，通常假定探测对象是均匀水平分层介质，在地表布置线性排列的检波器接收面波信号，提取面波的频散曲线，然后对频散曲线反演获得二维的剪切波速度剖面.通过对多条测线进行同样的处理，可以获得工区简单的三维速度结构.和大尺度的天然地震层析成像类似，这个过程和均匀水平分层的前提假设相矛盾.要获得真正的三维速度结构，就需要考虑散射面波的三维反演，在面波浅层勘探领域，这个工作目前还主要集中在正演方面[24~26, 34, 35]，国外有学者开展过散射面波的反演工作[36~38]，但主要集中在一些简单的数值模型和超声实验研究，并且在反演中对模型做了相当简单的假设，比如仅考虑密度的扰动[38].散射面波反演的一个重要工作即是计算三维体扰动对散射场的灵敏度，本文针对面波在浅层勘探工程中的应用，基于单次散射和远场假设，计算各向同性介质情形介质密度和波速扰动的灵敏度核，通过一个数值例子，分析有限频率面波三维灵敏度核的特征，考察模式耦合对灵敏度核的影响.第二和第三部分对基于射线理论和有限频率理论面波层析成像，及模拟面波散射的模式耦合方法进行了较为详细的评述，主要侧重物理概念的理解.
2 面波层析成像:从射线理论到有限频率如图 1所示，在传统的面波层析成像中，假设面波沿连接两个台站的大圆路径传播[39]，通过处理研究区域内不同路径传播的面波，获得经过每一个大圆路径的面波相速度后，对研究区域进行分块，然后根据射线穿过不同块体的路径，对所有块体的面波相速度进行拟合.对于单个模式，整个路径的平均相速度的分数扰动和每个块体的局部相速度的相对扰动之间的关系，可以通过沿射线的路径积分用下式表示[40]
|
(1) |
c是面波沿整个路径传播的平均相速度，c是面波在每个分块体内的局部相速度，δc表示相速度c的扰动量(下同，变量前面加δ表示该量的扰动)，它们均是频率的函数，L是大圆路径的长度.对不同的周期，重复这个过程，就可以得到不同周期的面波相速度分布图.这个过程对整个介质而言，其分块是二维的，即仅在水平方向对介质进行剖分，沿深度方向，面波相速度的变化是由面波的穿透深度来反映的，实际由面波模式沿深度的特征函数分布决定.在图(1)中，沿大圆路径下方，给出了单个模式的特征函数示意图.
|
图 1 从震源S点，沿大圆路径传播到接收位置R 的面波射线，在大圆路径下方给出了面波模式的特征函数示意图，修改自(Maupin, 2007)[41] Fig. 1 The ray path along the great circle connected the source S and the receiver R.The mode eigenfunction is also plotted below the ray path.Adapted from (Maupin.2007)[41]. |
前述过程是面波层析成像的第一步，这一步得到了不同周期(反应深度方向的维度)的面波相速度在水平面内的(水平方向的维度)分布.从方程(1)可以看出，在这一步中，沿整个大圆路径，局部相速度的扰动对整个路径平均相速度扰动的灵敏度核是常数1/L，即灵敏度核K1Dc 是一维的，且K1Dc = 1/L，这意味着，忽略了偏离大圆路径的介质参数扰动.
第二步是根据第一步得到的不同块体的相速度频散曲线反演该块体的介质参数.这个过程通常假定该块体是横向均匀的，只反演介质参数随深度的变化.对于各向同性情形，介质参数的扰动对相速度扰动的灵敏度核是一维的，即
|
(2) |
这里α，β和ρ 分别是介质的纵波速度、横波速度和密度，K表示对应介质参数扰动的灵敏度核，下标1D表示该灵敏度核是一维的.方程(1)和(2)构成了传统两步法面波层析成像的基础.将方程(2)带入方程(1)，可以得到[40]
|
(3) |
利用方程(3)可以直接由大圆路径的相速度反演介质参数.
从前述过程可以看出，在大圆路径假设下，忽略了偏离大圆路径介质扰动对面波传播的影响，这仅在高频射线近似下成立.通常地震记录的频带总是有限频宽的，如果介质的非均匀尺度和波长可比，就必须考虑偏离大圆路径的介质扰动对面波传播的影响，方程(1)中的灵敏度核和所考虑的整个区域内(水平方向的维度)的介质扰动都有关系，此时方程(1)可以写为
|
(4) |
由于波场的干涉作用，对灵敏度核起主要作用的范围集中在沿射线附近分布的第一菲涅带.将方程(2)代入方程(4)，可以得到有限频率下关于介质参数扰动的三维灵敏度核:
|
(5) |
对于更一般的线性地震层析成像情形，观测数据和模型参数之间的关系可以广义地写成如下形式:
|
(6) |
δd/d表示观测数据的相对扰动，可以是波速，走时，或者波形等；δm/m表示模型参数的相扰动，Km称为Fréchet核，它表示模型参数的扰动对数据扰动的灵敏程度，这也是将其称为灵敏度核(Sensitivity kernel)的原因，在反演算法中，该量是控制模型扰动的重要参数.
对于有限频率面波成像，方程(5)构成了其理论基础.从基于射线近似的面波成像到有限频率的面波成像，只需用方程(5)中的三维灵敏度核代替方程(3)中的一维灵敏度核即可.求取三维灵敏度核的基础，是求解三维非均匀体对于面波的散射，下一节我们讨论模拟面波散射的模式耦合法.
3 面波的散射:模式与模式耦合模拟地震波的传播与散射，通常有三种典型的方法，即纯数值方法，比如有限元、有限差分、谱元等；射线追踪理论；以及简正模理论.纯数值方法几乎可以处理任意复杂结构中的波动理论，但是这种方法很难描述问题的物理意义；射线理论是高频近似下的一种结果，主要用于体波的射线追踪；简正模理论用于模拟地球的自由震荡，这种方法将地球视为有限的球体.类似于拨动的弦和振动的单摆，球状地球的任意震动可以用一系列简正模式的叠加来描述，这些简正模式在数学上构成了一组完备正交集，求得这些模式的特征值和特征向量之后，对于不同的激发方式，求解不同模式的激发幅度，通过叠加求和即可求得任意一点的振动.对于球状地球来说，这些模式称为简正模(NormalModes)，它描述了地球在一定震源的激发下引起的地球整体的震荡，即相当于地震波从震源沿球状地球反向传播叠加形成的驻波解，因此它由一系列离散的谱线构成.相对于纯数值方法和射线理论，这种方法具有明确的物理意义和完美的数学描述，作为一种半解析的方法，通常用于检测数值方法的可靠性.由于对于高频情形，需要计算的模式个数较多，比较适于模拟低频波的传播.
在浅层勘探领域，由于尺度较小，通常将地球视为水平分层介质.与球状地球的自由震荡理论类似，可以用地震面波模式(Surfacewavemodes)的求和来描述质点的振动，不同的是，此时的面波模式是沿介质表面传播的行波解，其谱特征是连续的.与简正模相比，面波模式不是完备集，即通过面波模式的求和，不能包含介质内传播的所有震相.这可以通过模式和射线之间的关系来说明:由震源激发的波，在介质不同分层之间来回反射和折射，大部分能量被限制于近地表分层介质中，面波模式就是用于描述这部分能量的，因此，频率域中面波的多模和空间域中来回反射的射线的多径，实际是波场的不同表述形式.和视为有限大小的球状地球不同，水平层状介质是无限的，除了限制于层间来回反射的能量外，还有少部分能量是辐射到最底部无限半空间中的，面波模式不能描述这部分能量，因此面波模式不是完备集.有人考虑锁模近似(Lock-modeapproximation)技术[42]，在最底层引入一个高速反射介面，将能量限制在层状波导中来解决这一问题，但这种方法引入的是一个非物理的边界条件.Maupin(1996)则通过引入辐射模式的概念[43]，用来描述辐射到半空间中这部分能量.
本文考虑均匀各向同性弹性分层半空间中的面波情形，只考虑传统实数特征值的面波模式，这对于面波在浅层勘探中的应用是足够的.对某个模式m，在接收点X的位移，可以表示为
|
(7) |
Sm(ω)为源的激发幅度，通常是求解面波激发问题的待求量，Em(z，ω)为模式的特征函数，它是深度z的函数，只和介质的几何和物理参数有关，和震源以及波前形状和传播方向无关.eikmr为远场水平方向的波函数，考虑近场情形时，可以用Hankel函数代替.公式(7)中深度方向和水平方向的变量是分离的，这是面波模式的一个重要性质，也是传统两步法面波层析成像的基础.另一个重要特征是面波模式特征函数的正交，这意味着在均匀分层介质中面波模式是独立传播的.
当存在横向非均匀体时，独立传播的模式将会被非均匀体散射，如果散射体尺度不大，或者相对背景介质，参数的扰动较弱，散射场仍然可以用背景介质中的面波模式来表示，只是这些模式将不再独立，相互耦合(coupling)在一起[41]，通过求解不同模式之间的耦合系数，可以求解面波的散射场，进而求解有限频率面波三维灵敏度核.从理论上来说，考虑横向非均匀的介质作为背景模型也是可能的，比如 Chen(2007)[44]曾研究了非规则层状介质中广义的拉夫波模式，相对横向均匀情形，其理论处理十分复杂，而且对于横向非均匀的介质，模式通常不能写成变量分离的形式.由于实际反演中，一般都假定参考模型是横向均匀的，然后在此基础上对模型进行扰动，因此我们选取横向水平分层介质作为参考模型考虑面波的散射问题.
模式耦合法(Modecoupling)模拟面波散射的经典工作是由Snieder(1986)奠定的[36]，他给出了远场和波恩近似下面波模式之间的耦合系数，随后又将均匀分层介质的情形推广到球状地球简正模的情形[45].Maupin(2001)将该方法推广到任意各向异性介质情形[25]，且未作远场假设和波恩近似，实现了面波的多次散射.这里我们采用Snider的方法，如图 2所示的水平方向的坐标关系，对于水平均匀分层介质，远场近似下，面波的格林函数在球坐标中可以表示成并矢形式[46]
|
图 2 波恩近似下(单次散射)的模式与模式耦合示意图及几何坐标关系，左图修改自(Snieder, 2002)[46] Fig. 2 The geometric variables for the surface wave modes under Born approximation.Left panel adapted from(Snieder, 2002)[46] |
|
(8) |
S表示源点，X表示场点，†表示共轭转置.为简洁，这里忽略了对模式的求和，仅给出了单个模式的形式，用km表示模式m的水平方向波数，对实际存在的多个模式的情形，需要对所考虑的模式进行求和.p为面波的极化矢量，对于瑞利波和拉夫波分别为[36, 46]
|
(9) |
V和U分别表示瑞利波径向和垂直方向的特征函数，W表示Love波横向特征函数.每个方向的单位矢量用向量e加相应方向的下标表示.这里我们采用归一化引子[36, 47]
|
(10) |
其中U为群速度，I1 为能量积分.
根据格林函数的物理意义，如果在源点由单向力F激发波场，则沿R传播的一次场，可以表示为[36, 46]
|
(11) |
如果在X0 处存在散射体，在波恩近似和远场假设下，散射体引起的位移的扰动可以表示为
|
(12) |
其中
|
(13) |
Ωmn称为模式耦合系数，cijkl为介质弹性常数.对于入射的面波模式m，由于介质参数的扰动，一部分能量会转换为模式n，耦合系数则表示了能量转换的大小.如图 2所示，转换后的模式n沿路径R2 传播至接收位置，和从大圆路径R传播的一次场构成了地震面波记录的总场.式(12)是水平方向的面积分，沿深度方向上的积分隐含在耦合系数(13)式中，这是(7)式变量分离的结果.Snieder(1986)给出了式(12)的完美数学描述[36]:从右至左，面波模式m由源以一定的辐射方式被激发，然后以平面波形式几何扩散传播至散射体附近，在散射体附近面波模式互相耦合，耦合之后以模式n的特征，几何扩散传播至接收位置，在接收位置以一定的极化方式被记录.
若令(12)式中的m=n，即仅考虑模式自身的耦合，称为JWKB 近似，不同模式之间的耦合称为交叉模式耦合.可以通过考虑处在大圆路径的非均匀体来考虑两种模式耦合的不同作用.假定非均匀体恰好处在R穿过的大圆路径上.如图 2 所示，即图中X0 处在大圆路径上，共有3个模式传播，其中一个模式n1 =m，另两个模式分别为n2 和n3 .若只考虑模式自身的耦合，模式m经非均匀体散射后，其波前仍然以相同的速度传播，波阵面保持相同，那么在接收位置处，因介质参数的扰动引起的面波相速度的扰动主要由面波在非均匀体内部的走时扰动引起，只要非均匀体处在大圆路径上，非均匀体的位置对面波相速度的扰动就没有贡献.但对于交叉模式耦合，模式m经散射体转换为模式n之后，通常具有不同的相速度，这样模式n的波前会发生变化，在接收位置处，模式m和模式n2 或者n3 的叠加会出现相位的相消或相长干涉，引起接收面波的相位发生改变，相位的改变程度和非均匀体处在大圆路径上的位置有关，此时，除了面波穿过非均匀体内部的走时对面波相速度的扰动有贡献外，交叉模式之间的耦合也会引入面波相速度的改变.
4 有限频率面波3D 灵敏度核求得(13)式中的耦合系数之后，即可给出有限频率面波的三维灵敏度核表达式.(13)式中耦合系数是作用于介质密度和弹性常数的扰动量，其具体表达式可以参见文献(Maupin, 2001)[25]和(钱美平等，2008)[35]，在实际的反演问题中，我们直接关注的是介质的波速和密度.对于均匀各向同性情形，根据波速与介质密度和弹性常数之间的关系，可以得到波速扰动和弹性常数扰动之间的关系[20]:
|
(14) |
根据(14)式，可以将密度和弹性常数扰动的耦合系数转换为关于介质纵波和横波速度的耦合系数.表 1给出了不同面波模式之间的耦合系数，其中特征函数上方的黑点表示对深度z的微分.Zhou(2004)采用类似的方法[20]，利用简正模的耦合，给出了球状地球介质参数扰动引起的耦合系数，她的方法主要用于全球面波层析成像.这里我们主要关注面波在工程勘察中的应用，考虑水平分层介质的面波模式耦合.
|
表 1 面波模式之间的耦合系数，m表示入射的面波模式，n表示散射的面波模式，θ=φ2-φ1 为散射角 Table 1 Coupling coefficients between surface wave modes, incident and scattered modes are denoted repectively by m and n，θ=φ2-φ1 the scattered angle |
一般的面波层析成像，比如基于(5)式的相速度成像，仅仅利用了面波的相位信息，对于散射波成像，面波的相位和幅度信息都是有用的.下面分析面波的相位和幅度扰动的灵敏度核.根据公式(7)的形式，可以将面波位移表示成如下形式[48]，
|
(15) |
A和ψ 分别为位移的幅度和相位.介质参数相对背景模型的扰动会引起幅度和相位的扰动，地震层析成像的目的就是利用地震记录的这种扰动，推断介质参数的扰动，参考模型通常是作为初始模型或先验信息输入反演算法.实际遇到的反问题，通常假定介质的扰动很弱，引起的位移扰动相对一次场也很弱，即δu$\ll $u，此时，波恩近似是适用的.对方程(15)进行扰动，可以得到
|
(16) |
根据(16)式，可知δu/u的实部是幅度的相对扰动，虚部是相位的扰动，即
|
(17) |
其中位移的相对扰动可根据方程(11)和(12)得到
|
(18) |
KA和Kψ 分别表示相对幅度扰动和相位扰动的灵敏度核，上标α，β，ρ 分别对应纵波速度，横波速度和密度的扰动.将表(1)中对应的耦合系数带入(18)式，即可得到有限频率面波3D 灵敏度核.比如，将表(1)中Rayleigh-Rayleigh模式耦合中横波速度扰动的耦合系数Ωmnβ带入(18)式，可以得到横波速度扰动引起的相对幅度和相位扰动的3D 灵敏度核:
|
(19) |
|
(20) |
本节利用前述的理论和方法，针对一个两层的水平分层背景介质模型，计算有限频率面波灵敏度核，模型的几何和物理参数如表 2所示.特征函数的计算采用Saito的程序[49, 50].
|
|
表 2 参考模型的几何结构和介质参数 Table 2 The geometry and physical parameters of the reference model |
图 3a给出了该模型瑞利波和拉夫波的相速度频散曲线，目前在面波的工程应用中，获得测量的相速度频散曲线后，在介质横向均匀的假设下，通过不断调整预测模型，使预测模型的相速度和测量的相速度满足一定的拟合误差，认为该预测模型就是要反演的结果.对于单个模式的反演来说，深度方向的分辨率主要取决于信号的频率，较低的频率具有较深的勘探深度，其灵敏度核，即相速度关于介质参数的偏微分沿更大的深度范围分布.对于多模式的联合反演，通常同一频率下，高阶模式具有更深的勘探深度.图 3(b, c)给出了336 Hz时，3个瑞利波模式和拉夫波模式相速度关于介质横波速度的偏微分，即∂c/∂β，和式(2)中一维灵敏度核的关系式K1Dc =(β/c)·(∂c/∂β).可以看出，基阶模式对近地表信息更为敏感，高阶模式则反映较深部的信息.正如前面讨论的，在大圆路径假设下(对于工程应用来说，震源和检波器通常沿直线布设，假定面波的传播沿该直线路径传播，测量的频散信息只反应该路径下的介质信息)，反演中用到的灵敏度核是如图 3(b, c)所示的一维灵敏度核，没有考虑偏离直线路径的介质扰动的影响.
|
图 3 两层水平分层介质模型中，面波模式的相速度频散曲线(a)和关于横波速度的偏微分沿深度的分布(b, c) Fig. 3 The dispersion curves of the surface wave for a two - layered model shown in Table 1(a)，and the normalized derivatives of the phase velocity with respect to the shear wave velocity (b, c) |
对于区域地震面波相速度成像，其第一步是根据公式(1)反演纯路径的平均相速度，第二步和工程应用中的面波反演方法一致，根据公式(2)调整预测模型对相速度进行拟合.利用数字处理技术获得面波的相速度之后，整个反演过程均不涉及面波的激发问题.而有限频率面波成像，考虑实际波的传播与散射，其灵敏度核是三维的，和震源的辐射模式有关，对于面波的工程应用，一般可以用爆炸点源或者垂直方向的单向力源模拟面波的激发.本文采用垂直方向的单位力激发，计算有限频率面波三维灵敏度核，即令公式(18)中的F=ez，除非特殊声明，以下图中给出的结果均是位移的垂直分量.
图 4是基阶模式瑞利波有限频率面波三维灵敏度核KAβ和Kψβ，频率为336Hz, 图中只考虑了基阶模式自身的耦合，即JWKB 近似，忽略了交叉模式之间的耦合.图 4(a, b)是10m 深度处，灵敏度核在水平面的分布，分别是相对位移扰动的实部和虚部，即相对幅度和相位的灵敏度核.可以明显看出，围绕连接震源S和接收位置R 的射线路径周围，灵敏度核呈现菲涅耳带状分布，不单受射线路经的影响，其主要影响区域是前几个菲涅耳带，在有些面波成像中，通常只考虑对面波传播影响最大的第一菲涅带.根据这一物理现象，Rizwolleret.al.(2002)[51]和 YoshizawaandKennett(2002)[52]在传统面波层析成像方法的基础上，利用动态射线追踪技术，考虑第一菲涅带的影响.另外，在传统的成像方法中，引入类似高斯窗的平滑函数的做法[53]，其物理作用和考虑实际的菲涅带类似，只是高斯窗是一个人为引入的函数，图 4(a, b)的结果表明，实际的灵敏度核并不是类似高斯窗那样将射线路径简单拓宽为一个影响区域，而是一个振荡的菲涅耳带.和速度梯度变化的地球介质中呈香蕉-甜甜圈形状的体波走时灵敏度核不同，对面波来说沿射线路径的灵敏度核通常不是0.图 4(c, d, e, f)分别是沿图 4(a, b)中所示的剖面AB以及震源和接收位置的连线SR、灵敏度核沿深度的分布.
|
图 4 不考虑交叉模式耦合时，频率为336 Hz的基阶模式瑞利波有限频率三维灵敏度核KAβ和Kψβ (a, b)在深度为10m 时，灵敏度核在沿水平方向的分布，(c, d)和(e, f)分别是沿AB和SR 的剖面. Fig. 4 Finite frequency surfacewave 3Dsensitivity kernel (KAβ and Kψβ) of the fundamental Rayleigh wave with 336 Hz without taking the cross-mode coupling into account (a) and (b) are the projection on horizontal plane at the depth 10 m.(c, d) and (e, f) are the profiles along AB and SR，respectively. |
JWKB 近似相当于表 1中的m=n，如果我们进一步假定前向散射占据主导地位，即令表 1 的散射角为0，瑞利波模式之间的耦合系数沿深度的积分可以表示为
|
(21) |
注意到公式(14)中波速扰动和弹性常数扰动之间的关系，以及本文所取的归一化引子(10)式，将(21)式与文献[47]中的(7.78)式比较[注:该式分母漏掉了k2 因子]，可以看出(21)式为-0.5k2δc/c，这表明我们可以将(17)式中关于密度和波速扰动的三维灵敏度核表示成关于面波相速度扰动的二维灵敏度核，即在水平方向的二维灵敏度核和深度方向的一维灵敏度核是解耦的，这种解耦关系可以通过面波相速度的扰动获得.二维灵敏度核的求取，只需将公式(18)中Ωmn关于深度的积分去掉，同时乘以-k2/2即可.
水平方向的二维灵敏度核和深度方向的一维灵敏度核的解耦关系表明，沿深度方向变化主要由面波特征函数和频率决定，这种解耦关系也是两步法面波层析成像的基础.在两步法面波层析成像中，基于射线理论，第一步反演得到的结果称为纯路径平均的相速度，对有限频率面波成像，只需利用图 4(a, b)所示的二维灵敏度核，代替公式(1)中的1/L，所得到的二维面波相速度的分布不再是一个纯路径的平均，可以看成是区域平均的结果，然后再利用公式(2)中所示的一维灵敏度核，对区域平均的结果进行反演，得到介质的横波速度，即在JWKB 近似下，有限频率面波成像仍可以和传统成像方法类似，分两步进行.如果考虑交叉模式耦合的影响，类似图 4(c, d, e, f)中，沿深度方向的灵敏度核不单受某一个模式的特征函数的影响，随深度的变化将受多个模式的影响.
图 5是考虑交叉模式耦合对基阶模式瑞利波影响的结果，频率为336 Hz.由于对基阶模式瑞利波来说，交叉模式耦合的影响较小，其三维灵敏度核和 JWKB近似下图 4 中所示的灵敏度核很难分辨，这也是一般基阶模式反演只考虑JWKB近似的原因.为此，我们在图 5中给出了两种情况下三维灵敏度核的差，KA和Kψ 表示JWKB 近似下的结果，KA′ 和Kψ′ 表示考虑交叉模式耦合结果.如图 3 所示，在336Hz时，共有3个瑞利波模式和3个拉夫波模式，由于这里考虑垂直方向的点力源激发，在图 5中，只考虑了3个瑞利波模式之间的交叉耦合，从定量的估计来看，模式耦合对灵敏度核的影响约在5%左右，这一点可以根据图中所示的色标数值进行估计.与图 4(c, d)进行对比，可以发现，图 5(c, d)中，沿深度方向出现了明显的波动变化特点，而且在20m 深度以下，灵敏度核和20 m 深度以上处在可比的数量级，这均不是图 4(c, d)中所示的基阶瑞利波模式的变化特点，更接近于图 3b中高阶瑞利波模式的变化特征.
|
图 5 考虑交叉模式耦合(KA′ ，Kψ′)，和未考虑交叉模式耦合时(KA，Kψ)，基阶模式瑞利波有限频率三维灵敏度核之差，频率为336 Hz.(a, b)是10m 深度处灵敏度核在沿水平方向的分布，(c, d)是沿AB的剖面 Fig. 5 Difference of the sensitivity kernels computed with (KA，Kψ) and without (KA′ ，Kψ′) taking the cross-mode coupling into account.(a) and (b) are the projection on horizontal plane at the depth 10 m； (c, d) are the protiles along AB.Both of them are for the fundamental Rayleigh mode with 336 Hz |
图 6(a, b)分别是JWKB 近似和考虑交叉模式耦合时第二阶模式瑞利波的三维灵敏度核，和图 5相比，可以发现，对于本文的模型，交叉模式的耦合对高阶模式的影响非常重要，这个结果并不奇怪.本文的模型是一个速度递增的介质，对于这样的模型，基阶模式的能量占据主导地位，高阶模式能量很小，从基阶模式耦合到高阶模式的能量也很小，如果考虑高阶模式的瑞利波入射，在接收点耦合到基阶模式的能量将很大，这主要有介质模型决定.
|
图 6 不考虑交叉模式耦合(a)和考虑交叉模式耦合(b)时，二阶瑞利波的灵敏度核，频率为336 Hz, 深度为10m Fig. 6 Sensitivity kernels computed without (a) and with (b) taking the cross-mode coupling into account for the tirst overtone Rayleigh mode with 336 Hz.The depth is 10 m |
模式耦合的影响说明，如果对于某个模式交叉模式耦合的作用不能忽略，那么有限频率面波成像将不能按照传统的两步法，先获得纯路径的平均相速度的二维分布，再反演横波速度结构，因为此时水平方向的二维灵敏度核和深度方向的一维灵敏度核是耦合在一起的，公式(5)所示的关于介质参数扰动的三维灵敏度核，不能简单演绎为公式(4)所示的关于面波相速度扰动的二维灵敏度核，而是传统两步法面波层析成像的理论基础.
需要指出，这里给出的灵敏度核均是单色波的计算结果，因此在对灵敏度核影响较大的振荡菲涅带周围仍有小幅的振荡存在.实际应用中处理的信号是具有一定频率范围的时域信号，其灵敏度核可以通过单色波灵敏度核的叠加求取[54]，或采用 multitaper技术[20]，此时，在单色波菲涅带周围的震荡现象，会因频率临近的信号相消干涉而消失，起主要作用的是沿射线路径周围分布的菲涅带.单色波三维灵敏度核的信息已包含了有限频率灵敏度核的所有信息，这里不再对时域三维灵敏度核进行讨论.另外，本文主要针对面波的工程应用，只分析了垂直点力源激发的情况，本文理论同样适用区域天然地震面波层析成像，只需要将公式(18)中的点力源用矩张量源代替，用(ikmeR+ez∂zs)pm(zs，φ1):M代替pm†(zs，φ1)·F即可，这里M为矩张量源，:为双点积[36, 47].
6 结论在工程应用领域中考虑面波有限频率影响的研究工作较少，Campman和Riyanti(2007)[38]等采用积分方程法给出了面波密度扰动的灵敏度核，利用散射面波对密度进行了初步的反演.本文在评述大尺度射线理论和有限频率面波层析成像的基础上，针对面波在工程应用中经常遇到的水平分层介质和垂直点源激发的情形，采用面波模式耦合的方法，分析了密度、纵波速度和横波速度扰动，引起的有限频率面波三维灵敏度核的特征.基于一个两层参考模型的数值例子，考察了模式耦合对三维灵敏核的影响.结果表明，在前向散射和JWKB 近似下，介质参数扰动引起的三维灵敏度核可以蜕化为面波相速度扰动的二维灵敏度核，水平方向的二维灵敏度核和深度方向的一维灵敏度核是解耦的，有限频率面波成像和基于射线理论的方法类似，可以分两步进行.如果考虑交叉模式耦合的影响，这种解耦关系将不存在，反演时需要引入三维面波灵敏度核.
| [1] | Dahlen F A, Hung S H, Nolet G. Frechet kernels for finite-frequency traveltimes -I. Theory. Geophys. J. Int. , 2000, 141: 157-174. DOI:10.1046/j.1365-246X.2000.00070.x |
| [2] | Hung S H, Dahlen F A, Nolet G. Frechet kernels for finite-frequency traveltimes -II. Example. Geophys. J. Int. , 2000, 141: 175-203. DOI:10.1046/j.1365-246X.2000.00072.x |
| [3] | Spetzler J, Snieder R. The Fresnel volume and transmitted waves. Geophysics , 2004, 69(3): 653-663. DOI:10.1190/1.1759451 |
| [4] | Hung S H, Dahlen F A, Nolet G. Wavefront healing: a banana-doughnut perspective. Geophys. J. Int. , 2001, 146: 289-312. DOI:10.1046/j.1365-246x.2001.01466.x |
| [5] | Nolet G, Dahlen F A. Wave front healing and the evolution of seicmic delay times. J.Geo.Res. , 2000, 105(B8): 19043-19054. DOI:10.1029/2000JB900161 |
| [6] | Spetzler J, Snieder R. The effect of small-scale heterogeneity on the arrival time of waves. Geophys. J. Int. , 2001, 145: 786-796. DOI:10.1046/j.1365-246x.2001.01438.x |
| [7] | Spetzler J, Snieder R. The formation of caustics in two- and three-dimensional media. Geophys. J. Int. , 2001, 144: 175-182. DOI:10.1046/j.1365-246x.2001.00308.x |
| [8] | Spetzler J, Trampert J, Snieder R. Are we exceeding the limits of the great circle approximation in global surface wave tomography?. Geophysical Research Letters , 2001, 28(12): 2341-2344. DOI:10.1029/2000GL012691 |
| [9] | Spetzler J, Trampert J, Snieder R. The effect of scattering in surface wave tomography. Geophys. J. Int. , 2002, 149: 755-767. DOI:10.1046/j.1365-246X.2002.01683.x |
| [10] | 杨峰, 黄金莉. 有限频率地震层析成像方法及研究进展. 地震 , 2009, 29(4): 52–62. Yang F, Huang J L. Finite-frequency seismic tomography: methods and progresses. Earthquake (in Chinese) , 2009, 29(4): 52-62. |
| [11] | 徐小明, 史大年, 李信富. 有限频层析成像方法研究进展. 地球物理学进展 , 2009, 24(2): 432–438. Xu X M, Shi D N, Li X F. Research progress of the finite frequency tomography method. Progress in Geophysics (in Chinese) , 2009, 24(2): 432-438. |
| [12] | 杨峰, 黄金莉, 杨挺. 应用远震有限频率层析成像反演首都圈上地幔速度结构. 地球物理学报 , 2010, 53(8): 1806–1816. Yang F, Huang J L, Yang T. Upper mantle structure beneath the Chinese acpital region from teleseismic finite frequency tomography. Chinese J. Geophys. (in Chinese) , 2010, 53(8): 1806-1816. |
| [13] | Marquering H, Dahlen H A, Nolet G. Three-dimensional sensitivity kernels for fnite-frequency traveltimes: the banana-doughnut paradox. Geophys. J. Int. , 1999, 137: 805-815. DOI:10.1046/j.1365-246x.1999.00837.x |
| [14] | van der Hilst R, de Hoop M. Banana-doughnut kernels and mantle tomography. Geophys. J. Int. , 2005, 163: 956-961. DOI:10.1111/gji.2005.163.issue-3 |
| [15] | Montelli R, Nolet G, Dahlen F A, et al. Finite frequency tomography reveals a variety of plumes in the mantle. Science , 2004, 303: 338-343. DOI:10.1126/science.1092485 |
| [16] | Montelli R, Nolet G, Masters G, et al. Global P and PP traveltime tomography: rays versus waves. Geophys. J. Int. , 2004, 158: 637-654. DOI:10.1111/gji.2004.158.issue-2 |
| [17] | Montelli R, Nolet G, Dahlen F A. Comment on 'Banana-doughnut kernels and mantle tomography' by van der Hilst and de Hoop. Geophys. J. Int. , 2006, 167: 1204-1210. DOI:10.1111/gji.2006.167.issue-3 |
| [18] | Ren Y, Shen Y. Finite frequency tomography in southeastern Tibet: Evidence for the causal relationship between mantle lithosphere delamination and the north-south trending rifts. J.Geo.Res. , 2008, 113: B10316. DOI:10.1029/2008JB005615 |
| [19] | Zhao L, Jordan T H, Chapman C H. Three-dimensional Frechet differential kernels for seismic delay times. Geophys. J. Int. , 2000, 141: 558-576. DOI:10.1046/j.1365-246x.2000.00085.x |
| [20] | Zhou Y, Dahlen F, Nolet G. Three-dimensional sensitivity kernels for surface wave observables. Geophys. J. Int. , 2004, 158: 142-168. DOI:10.1111/gji.2004.158.issue-1 |
| [21] | Zhou Y. Surface-wave sensitivity to 3-D anelasticity. Geophys. J. Int. , 2009, 178: 1403-1410. DOI:10.1111/gji.2009.178.issue-3 |
| [22] | Zhou Y. Multimodel surface wave sensitivity kernels in radially anisotropic earth media. Geophys. J. Int. , 2009, 176: 865-888. DOI:10.1111/gji.2009.176.issue-3 |
| [23] | Yoshizawa K, Kennett B L N. Sensitivity kernels for finite-frequency surface waves. Geophys.J.Int. , 2005, 162: 910-926. DOI:10.1111/gji.2005.162.issue-3 |
| [24] | Friederich W, Wielandt E, Stange S. Multiple forward scattering of surface waves: comparison with an exact solution and Born single scattering methods. Geophys. J. Int. , 1993, 112: 264-275. DOI:10.1111/gji.1993.112.issue-2 |
| [25] | Maupin V. A multiple-scattering scheme for modelling surface wave propagation in isotropic and anisotropic three-dimensional structures. Geophys. J. Int. , 2001, 146: 332-348. DOI:10.1046/j.1365-246x.2001.01460.x |
| [26] | Lu L Y, Maupin V, Zeng R S, et al. Scattering of surface waves modelled by the integral equation method. Geophys.J.Int. , 2008, 174: 857-872. DOI:10.1111/gji.2008.174.issue-3 |
| [27] | Chen M, Tromp J. Theoretical and numerical investigations of global and regional seismic wave propagation in weakly anisotropic earth models. Geophys, J.Int. , 2007, 168: 1130-1152. DOI:10.1111/gji.2007.168.issue-3 |
| [28] | Sieminski A, Liu Q, Trampert J, et al. Finite-frequency sensitivity of surface waves to anisotropy based upon adjoint methods. Geophys. J. Int. , 2007, 168: 1153-1174. DOI:10.1111/gji.2007.168.issue-3 |
| [29] | Sieminski A, Trampert J, Tromp J. Principal component analysis of anisotropic finite-frequency sensitivity kernels. Geophys. J. Int. , 2009, 179: 1186-1198. DOI:10.1111/gji.2009.179.issue-2 |
| [30] | 刘云祯, 王振东. 瞬态面波法的数据采集处理系统及其应用实例. 物探与化探 , 1996, 20(1): 15–18. Liu Y Z, Wang Z D. Data collection and processing system of transient surface wave method and examples of its application. Geophysical & Geochemical Exploration (in Chinese) , 1996, 20(1): 15-18. |
| [31] | Park C P, Miller R D, Xia J H. Multichannel analysis of surface waves. Geophysics , 1999, 64: 800-808. DOI:10.1190/1.1444590 |
| [32] | 张碧星, 鲁来玉, 鲍光淑. 瑞利波勘探中"之"字形频散曲线研究. 地球物理学报 , 2002, 45(2): 263–274. Zhang B X, Lu L Y, Bao G S. A study on zigzag dispersion curves in Rayleigh wave exploration. Chinese J. Geophys. (in Chinese) , 2002, 45(2): 263-274. |
| [33] | Lu L Y, Wang C H, Zhang B X. Inversion of multimode Rayleigh waves in the presence of a low-velocity layer: Numerical and Laboratory study. Geophys. J. Int. , 2007, 168: 1235-1246. DOI:10.1111/gji.2007.168.issue-3 |
| [34] | Riyanti C D, Herman G C. Three-dimensional elastic scattering by near-surface heterogeneities. Geophys. J.Int. , 2005, 160: 609-620. DOI:10.1111/gji.2005.160.issue-2 |
| [35] | 钱美平, 鲁来玉, 刘斌. 三维非均匀结构体对面波的多次散射. 物探化探计算技术 , 2008, 1: 1–9. Qian M P, Lu L Y, Liu B. Multi-scattering of surface waves from three-dimensional structure. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese) , 2008, 1: 1-9. |
| [36] | Snieder R. 3-D linearized scattering of surface waves and a formalism for surface wave holography. Geophys. J. R. astr. Soc. , 1986, 84: 581-605. DOI:10.1111/j.1365-246X.1986.tb04372.x |
| [37] | Kaslilar A. Inverse scattering of surface waves: imaging of near-surface heterogeneities. Geophys. J.Int. , 2007, 171: 352-367. DOI:10.1111/gji.2007.171.issue-1 |
| [38] | Campman X, Riyanti C D. Non-linear inversion of scattered seismic surface waves. Geophys.J.Int. , 2007, 171: 1118-1125. DOI:10.1111/j.1365-246X.2007.03557.x |
| [39] | 何正勤, 叶太兰, 丁志峰. 华北东北部面波相速度层析成像. 地球物理学报 , 2009, 52(5): 1233–1242. He Z Q, Ye T L, Ding Z F. Surface wave tomography for the phase velocity in the northeastern part of North China. Chinese J. Geophys. (in Chinese) , 2009, 52(5): 1233-1242. |
| [40] | Trampert J, Spetzler J. Surface wave tomography: finite frequency effects lost in the null space. Geophys. J. Int. , 2006, 164: 394-400. DOI:10.1111/gji.2006.164.issue-2 |
| [41] | Maupin V. Introduction to mode coupling methods for surface waves.In: Maupin V, Wu R S, eds. Advances in Wave Propagation in Heterogeneous Media. 2007, 48:127~155 |
| [42] | Nolet G, Sleeman R, Nijhot V, et al. Synthetic reflection seismograms in three dimensions by a locked-mode approximation. Geophysics , 1989, 54: 350-358. DOI:10.1190/1.1442660 |
| [43] | Maupin V. The radiation modes of a vertically varying half-space: a new representation of the complete Green's function in terms of modes. Geophys. J. Int. , 1996, 126: 762-780. DOI:10.1111/gji.1996.126.issue-3 |
| [44] | Chen X F. Generation and propagation of seismic SH waves in multi-layeredd media with irregular interfaces.In:Maupin V, Wu R S, eds. Advances in Wave Propagation in Heterogeneous Media.2007, 48:191~264 |
| [45] | Snieder R, Nolet G. Linearized scattering of surface waves on a spherical Earth. J. Geophys. , 1987, 61: 55-63. |
| [46] | Snieder R. Scattering of surface waves. In: Pike R, Sabatier P, eds. Scattering and Inverse Scattering in Pure and Applied Science. Academic Press, San Diedo, 2002, 562~577 |
| [47] | Aki K, Richards P G. Quantitative Seismology. 1980 . |
| [48] | Nolet G. A breviary of seismic tomography:Imaging the Interior of the Earth and Sun. 2008 . |
| [49] | Saito M. Disper80: a subroutine package for the calculation of seismic normal mode solutions. In: Durk J ed. Seismological Algorithms:Computational Methods and Computer Prorams. Doornbos, Academic Press, New York , 1988: 293-319. |
| [50] | Takeuchi H, Saito M. Seismic surface waves in seismology: surface waves and earth oscillations.In:Bolt B A ed.Methods in Computational Physics. New York:Academic Press , 1972, 11: 217-295. |
| [51] | Ritzwoller M H, Shapiro N M, Barmin M P, Levshin A L. Global surface wave diffraction tomography. J.Geophy. Res. , 2002: 107-B12. |
| [52] | Yoshizawa K, Kennett B L N. Determination of the influence zone for surface wave paths. Geophys. J. Int. , 2002, 149: 440-453. DOI:10.1046/j.1365-246X.2002.01659.x |
| [53] | Debayle E, Sambridge M. Inversion of massive surface wave data sets: Model construction and resolution assessment. J.Geo.Res. , 2004, 109: B02316. |
| [54] | 刘玉柱, 董良国, 李培明, 等. 初至波菲涅耳体地震层析成像. 地球物理学报 , 2009, 52(9): 2310–2320. Liu Y Z, Dong L G, Li P M, et al. Fresnel volume tomography based on the first arrival of the seismic wave. Chinese J. Geophys. (in Chinese) , 2009, 52(9): 2310-2320. |