KH方法计算远震转换波.docx

《KH方法计算远震转换波.docx》由会员分享，可在线阅读，更多相关《KH方法计算远震转换波.docx（17页珍藏版）》请在三一办公上搜索。

1、用Kirchhoff-Helmholtz方法计算远震PS转换波1李红光，吴庆举中国地震局地球物理研究所，北京100081摘要本文将Kirchhoff-Helmholtz积分方法推广到远震转换波的合成地震图计算，其优点 是能够计算复杂界面的地震波。通过与反射率法及动力学射线追踪的对比，表明KH积分方 法能很好地模拟远震转换波震相，且精度较高。KH积分方法能够计算任意复杂界面的地震 波，是一种研究地壳上地幔结构的有效方法。关键词 KH方法，动力学射线追踪，反射率法，远震转换波1引言合成地震图是地震波形解释的一个重要工具。自Helmbergers提出了广义射线合成地 震图计算方法，并将其应用于实际地

2、震波形资料的对比和解释工作以后，水平层状介质合成 地震图的计算方法得到了迅猛发展，产生了广义射线、反射率法、WKBJ近似句、广义反 射一透射系数方法5,6、离散波数法7,8及全波理论9等合成地震图计算方法。随着地震台网 的密集，特别是现代大型地震科学台阵观测技术的问世，使地壳上地幔横向不均匀性的调查 和研究成为可能，因此，迫切需要一种能够快速计算横向不均匀介质的合成地震图方法。目 前横向不均匀介质中地震波响应的计算方法主要分为以下几类：直接求解弹性动力学波动方 程及边界条件的数值方法，例如有限差分方法，有限元方法13；将非均匀介质分解成参 考介质和微扰介质，把参考介质中的背景场与散射场相加得到

3、非均匀介质中总波场的扰动法 14；基于射线理论的高频渐近方法15；把地震波场表示为边界上的积分方程的边界积分类 方法16）等。有限差分方法虽然能很好地模拟复杂介质中的地震波传播过程，但由于有限差 基金项目国家自然科学基金项目（40574040）资助作者简介李红光.男.1982年12月出生.硕士研究生.专业:地球探测与信息技术Email:fengzhe122分方法非常耗时，限制了其在地壳上地幔尺度的地震波波形模拟和解释工作中的广泛应用。 扰动法只适用于横向非均匀性相对较弱条件下的波形模拟工作。当介质的非均匀性比地震波 波长大得多的时候，高频渐近方法效果最佳。最典型的高频渐近方法是传统的动力学射线

4、追 踪方法，但当存在奇点的时候，传统的射线理论方法不能正常工作。为了克服奇点带来的困 难，人们陆续发展了高斯光束方法m和Maslov方法mu”等。目前复杂介质中波形模拟应用 最广泛的是高频近似方法。边界积分方法是用边界上波场的积分方程来表示接收点处的地震波场，它很早就被应 用到物理学的其他领域，如固体力学、电磁学、流体动力学及热力动力学等。 Kirchhoff-Helmholtz积分方法（简称KH方法）近即是一种边界积分方法，它是从严格的 波动方程出发，将体积分转化为边界积分，当用于计算反射/透射波时，KH积分方法把反射 /透射界面上的每个点都看作一个点源，认为每个点源对地震波振幅都有一定的贡

5、献，把界 面上每个点的贡献相加就得到了地震波响应。KH方法不同于以往的Kirchhoff方法： Kirchhoff方法是在时间域内做积分，而KH方法是在频率域做积分然后通过傅氏变换得到 地震波响应。KH方法不仅计算效率高，而且还能计算介质中存在焦点和散射点的地震波响 应，因此，KH积分方法是一种非常好的计算横向非均匀介质合成地震图的方法。KH方法最早被Frazer：24应用在反射波的波场模拟，本文把KH方法推广到远震转换波的 计算。在推导出远震转换波的KH积分公式的基础上，编制了计算程序，给出了用KH积分 计算远震转换波的算例，并与反射率法及动力学射线追踪的结果进行了对比。2KH积分方法的原理

6、我们首先简要概括一下远震转换波的KH积分公式推导过程。如图1所示，V是一个2D 或3D弹性介质体，令3V为V的边界，n是V边界上任一点的外法线向量。f1是作用于空 间乂1处的体力，且在V和3V上为零，u和t分别是体力f产生的位移矢量和应力张量。f 11112是作用在空间X处的另一体力，它在V外和aV上为零，u和T分别是体力f产生位移矢 2222量和应力张量。设f为a方向的集中力，那么f = a 5（x- x）。我们有以下公式24 22222（1）a - u (x ) = j n - (t - u - T - u ) dA2121221图1 一个带有边界av的体积单元V，n是边界av的外法线方向

7、，我们希望计算出力f1在x 2点处引起的位移。Figure.1 A volume V in（E2 or E3） with boundary a V .The vector n is the outward-pointingunit normal to a V .We wish to calculate the displacement at x2 due to a force f1 at x根据公式（1），我们只要知道力f、f在av上所产生的位移矢量u、u与应力张量t、12121t，就能求出力f在空间V内任意一点x上的位移矢量u （x ）。21212下面，我们将进一步论述如何利用公式（1）计算

8、远震转换波场。图2用KH积分公式(12)来计算间断面Z上的透射波场Figure. 2 Application of the KH equations (12) to the calculation of transmission wavesfrom a material discontinuity .如图2所示，是一个复杂的透射界面，和3V组成一个闭合的空间V,震源x1位 于V外部，接收点x2位于V中。U1表示震源入射到界面上的透射波，U2表示界面上任 意一点向接收点辐射的波场，可视为接收点处虚拟震源在界面上所产生的位移场。在高频近似条件下，信号的波长远远小于积分面的变化尺度，因此，可以把附近

9、的u 1和u 2看作平面波。由平面波理论得到平面波的位移公式26:（2）u P = At P exp i （tP - r） / a ,其中A为平面波的振幅大小，为平面波的方向向量，3为平面波的频率，向量r为平（3）面波离开震源的距离向量，a为纵波速度。由公式（A1）我们得到:V uP = t PuP , V - uP = uP -1 P aa由于位移u与应力T存在如下关系:t = XV - uI + 闵V u + （V u）t （4）把公式（A2）代入（A3）就得到P波的应力张量t P: it P = X（uP - tP ）I + 2日uPtP .（5）a根据公式（A4），对于两个相互作用的P

10、波：n - （t p - u P - t P - u P ） = u P - X （t Pn 一 nt p ） + 2 日（nt P 一 t Pn） - u P（6）图3两个弹性半无限空间耦合在一起，边界是个不规则面Figure. 3 Two elastic half-spaces welded together along an irregular boundary.如图3所示，平面波U在界面Z上的透射波为U了、U尸。令U的位移为:uP = APtP11 1（7）那么透射波U了的位移可写为:uPP = Apptpp = uP -1 pPPtpp(8)111111其中：tPP =abG-疽1-

11、a 2b 2 n表示透射波U pp的方向向量，a p = tp - (I - nn) /以， 1111112a = | a p |，a = a p /a ； tP是入射波U p的方向向量；PP是平面波U p在界面上的透射系 11111数；n是z上的外法线向量。根据公式(6)得到upp和up的作用表达式：n - (t PP - uP -t P - uPP ) = uPP - 入(t PPn - nt P ) + 2日(nt PP - t Pn) - uP1221a 112122代入(A7)：=些 uP - tPPPtPP - 入(tPPn - nt P ) + 2 日(nt PP - tPn)

12、- uP (9) a 11112122令：T PP = t pP P tPP - 入(tPPn - nt p ) + 2 日(nt PP - tPn)121a 11212则公式（6）简化为:n - (t pp - uP -t p - uPP ) = uP - T PP - uP12211122(10)我们称Tpp为P波的透射相互作用系数。需要注意的是：上面公式中的Tpp只有在高频近似，信号的波长远远小于界面Z的变化尺度条件下才是正确的。在各向同性弹性介质中，假设震源位于点，a 8 p是与空间位置有关的光滑函数，并且信号的波长远远小于这些参数的变化尺度，那么我们得到x点处的P波的振动方程26八/

13、.一、()tP F P exp( i T P) 1J pap a. Bp(11)公式（11）中，tp = tP （X,X ）是由震源x发出的P波射线在x点的切向量，TP = TP （X,X ） 1111111是P波从x1到x的走时。p 1和a 1分别是x1处的介质密度和纵波速度，p和a分别是x处的介质密度和纵波速度。FP是震源参数，Bp是扩展因子。(a)图4 (a)三维中的射线束；(b)二维中的射线束Figure. 4 (a) Ray tube in 3-D;(b) ray tube in 2-D3D介质中，扩展因子B可以表示为：BP = 4兀a -JdA /dQ( 12)112D介质中，扩展

14、因子表示为：Bp =、：8氏淼a dl /d0 exp( is /4)(13)各参数如图7所示。在各向同性的介质中，有;dA / dQ =|xxJI, dl / d0 =|xx。方程(12)和(13)是假设射线没有出现焦散的情况。如果实际中射线存在一个或多个 焦散点，那么扩展因子B就必须乘上exp i sgn()。(x, x尸/2，其中a (x,七)是KMAN 指数27在2D问题中一条射线的KMAN指数初始值为零，每遇到一个焦散点KMAN指数 就自增一。对于集中力f、= a()8 (x - x1)，震源参数Fp为：FP = a -1P(x )(14)111其中，P(号是从x1到x的P波射线在震

15、源x1的切向量。对于一个双力偶 f = M V5 (x - x )，其中M=M(3 )是二阶对称张量矩阵，震源参数Fp为：F P = -1P (x ) M tP (x )(15)1 a 1111对点爆炸源f、= P()Vb(x - x1)，有:isF P = P ()i前面我们已经得到:(16)a - u (x ) = j n - (t - u - t - u ) dA212评1221(17)把(10)、(11)代入公式(17)中就得到了远震转换波的波场up :a u pp (x ) = j u p (x) - T pp - u p (x) d Z (x)2121122Z (x )=J M=U

16、P .玖上业心 d (x) Jp a z(x)Jpa 1122 B p(18)其中Fp = a tp (x )，tp (x )是从x到x的射线在点x上的单位切向量。由于公式(18)22222222两边都与a 2成线性关系，两边消掉 a，得到pp透射波的位移公式:u pp (x ) = ：= J t p (x ) up T pp t p- :d Z (x)12 Jp a(x) 221122 B ppa(19)其中:t pp = t pPP迪tpp (入(tpp121 a 11n 一 nt p) + 2日(n pp 一 tpn)212tpp = a qq 11q p = t p (i 一 nn)

17、/ a，q = ii q p ii1121tp , 3分别为源x和x的入射方向。p、a表示界面Z上反射点处介质的密度与纵波速度，o、12121a1表示震源x1处介质的密度与纵波速度，p2、a2表示接收点x2处介质的密度与纵波速度。Bp 是扩展因子。PS透射波的积分公式:ups (x ) = J u p T ps b b (x ) + c c (x ) P( d Z (x)12 Jp2P2 z(x) 1122 2 2222 B s pp2(20)对于远震情况，我们把u1看作平面波入射，因此u1的方向是不变的。3数值检验3.1计算参数的选取由公式（9）和（20）我们可以看出：转换波场是通过边界积分

18、得到的，所以积分区间 和积分步长的选取是制约KH积分方法计算效率的两个主要因素。KH积分是个收敛的积分， 而我们不可能取无穷大的积分区间，所以必须找一个合适的积分区间，使得结果即能满足我 们的要求，又有较高的计算效率。积分步长取得越小，得到的结果越精确，但计算效率会很 低；如果积分步长增大，计算效率会得到提高，但会损失计算精度。因此，需要选择合适的 积分区间和尽可能大的积分步长，使得在保证计算精度的前提下提高计算效率。我们首先来考虑积分步长的选取。研究表明：积分步长需要满足AL. WAT，其中 L是积分步长，Vmax是介质中最大波速，AT是采样间隔。根据采样定理：T =2f max所以L 2

19、fmax厂，其中fmax是地震波最大频率，人mk是指地震波最大频率maxfmaxfmax在最大波速介质中的波长。也就是说，我们选取积分步长不能大于最大波长4ax的一半。图5给出了不同的积分步长对波形的影响。图5不同的积分步长对波形的影响Figure 5 Traces showing the effect of using different steps of the integral图5是用不同的积分步长计算的结果，图中可以看出：当积分步长大于xq2时，波形 会出现震荡，积分步长越大，震荡越明显。为了提高计算效率，我们在以后的计算中取积分 步长为入ma/2，这样即保证了计算精度，也兼顾了计算效

20、率。)下面我们来考虑积分区间大小的选取。首先我们采用频率域高斯函数“0 ) = e-以20作 为入射平面波的频谱，其中心频率(优势频率)为，在其两侧，高斯函数的振幅谱视a 的大小而有不同的衰减特性：a越小，衰减越快；a越大，衰减越慢。图6是不同中心频率的高斯函数振幅谱，其中a取2.5。图6不同优势频率的入射频谱图Figure 6 The frequency spectrum of incidence waves with the different main frequency根据菲涅尔原理，我们要计算的地震波能量主要集中在菲涅尔第一半径内。那么我们现在就尝试找出积分区间的选取与菲涅尔第一半径

21、的关系。菲涅尔第一半径可以表示为：.TvhRf = vf 2 cos 0其中，Rf是菲涅尔第一半径，T是地震波周期，v是地震波波速，h是界面的深度，e 是入射角度。图7是一个计算菲涅尔半径的示意图。其中射线2是根据射线理论得到的射线 路径，射线1和3是菲涅尔第一半径区的边界射线。图7菲涅尔第一半径区示意图Figure 7 Schematic illustration of the first Fresnel zone下面我们用图6中的10种入射波分别计算远震地震波的振幅与积分区间的关系。如图8, 图中横坐标是积分区间大小（我们用对数坐标表示，能更清楚的看出振幅的变化趋势），纵 坐标是地震波振幅

22、大小。折射波积茶区间与振幅的关系.0.5T1030潮务.区.间成图8折射波积分区间大小与菲涅尔第一半径的关系Figure 8 the relation the range of integral interface of refraction with Fresnel first radii图8我们可以看出，随着积分区间的变大，振幅将趋于一个定值，这说明了我们的积 分区间取到一定大小就可以（比如例中取10Km就足够了）。图8中，菲涅尔第一半径（黑点） 与地震波最大振幅出现的位置基本重合，这说明我们选取KH积分区间可以以菲涅尔第一半 径为参考，比它略大最好。我们实际中用的入射波信号并不是单一频率

23、的，而菲涅尔半径是与频率的平方根成反比 的，所以我们需要通过试验来确定积分区间的大小。如图9,横坐标代表积分区间的大小（用 最小的波长数目来表示），纵坐标是地震波折射振幅大小。我们分别用积分步长dn=M2, dn=M4, dn/8得到三条曲线（人是信号最小波长）。图中很明显看出，三条曲线基本重合 在一起，而且在积分区间达到20个人的时候，地震波的振幅就趋于一个定值。Figure 9 the relation of maximum amplitudes of refraction with the range of integral interface3.2数值计算我们首先计算单层水平界面的远震

24、转换波，也是为了检验程序。模型参数见表1。我们 分别用KH方法和反射率法5计算了不同慢度的远震转换波响应。图7给出了两种方法在水 平界面上的透射转换波PP和PS响应。横坐标表示地震波的到时，纵坐标是射线的慢度。表1单层水平介质的结构模型参数Table 1 Parameters used in a single flat horizontal interface界面深度(km)P波速度(km/s)S 波速度(km/s)密度(g/cm3)306.243.62.777.794.53.26从图10的结果中我们发现：两种方法在震相、到时以及相对振幅的变化上都吻合的很 好，而且对于不同慢度的远震波都能得到

25、相近地结果。表明KH积分方法能很好地模拟水平 界面上的远震转换波。0 060.05horizoncal component140.vertical component(与语落 S1DUMO-S3 O O.(a)(b)图10 KH方法与反射率法在水平界面上的比较结果a：水平分量；b：垂直分量Figure. 10 Comparison of results of KH synthetics and Reflectivity method in a flat horizontal interface. (a) the horizontal component of wave respond;(b)

26、the vertical component of wave respond.进一步地，我们用KH方法模拟了倾斜界面(倾角为10)上的远震转换波(介质参 数同上)，并与动力学射线追踪法28的结果进行了对比，结果见图11，横坐标表示地震波的 到时(相对到时)，纵坐标是射线的慢度。不难发现：KH积分方法所得到的远震转换波与 动力学射线追踪法的结果非常相似，表明KH积分方法能很好地模拟横向变速介质中的远震 转换波。8O0=horizoncal component8OO.020.123456Time/s.070=6 6 4 o o O 0.0.0. (LUXafssCDUMO-s3 O 0.7 6 5

27、 4 o o o O 0.0.0.0.(专学50)如。03 O 0.02vertical componentTime/s(a)(b)图11 KH方法与动力学射线追踪法在倾斜界面上的比较结果：a,水平分量;b,垂向分量Figure 11 Comparison of results of KH synthetics with Dynamic ray tracing method on an incline.(a) the horizontal component of wave respond;(b) the vertical component of wave respond 4结果与讨论KH积

28、分方法通过边界积分来计算任意点处的波场，避免了焦散点处位移直接计算的困 难。本文将KH积分方法应用于远震转换波的计算，数值计算分析表明：积分区间与第一菲 涅尔半径相近时，波场振幅达到最大，为了保证计算的稳定，我们一般取积分区间大于第一 菲涅尔半径。积分步长满足公式 L/VmaxAT，其中AL时积分步长，Vmax是介质中的最 大速度，AT是采样间隔。模型数值检验表明用KH积分方法所得到的远震转换波与反射率 及动力学射线追踪的结果非常相似。本文仅是用KH积分计算远震转换波的初步尝试。我们的工作表明KH积分方法是一种 比较好的、能模拟复杂界面透射波的合成地震图方法。通过对其作进一步的发展，可望应用

29、于地壳上地幔尺度的波形模拟工作。参考文献1 HelmbergerD V 1968. Crust-mantle transitioninthe Bering seaJ.Bull.Seism.Soc.Am , 58:1792142 HelmbergerDV1974. Generalized ray theoryforshear dislocationsJ.Bull.Seism.Soc.Am , 64:45543 Fuchs K, Muller G 1971. Computation of synthetic seismograms with the reflectivity method and

30、comparison with observationJ. Geophy.J.R.astr.Soc, 23(2):4174334 Chapman C H. 1978. A new method for computing synthetic seismogramsJ. Geophy.J.R.astr.Soc, 54:4815185 Kennett BLN. 1983. Seismic waves propagation in a stratified mediaM. Cambridge:Cambridge University, 3426 Chen X F. 1990. Seismograms

31、 synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method.Part I.Theory of 2-D SH caseJ. Bull.Seism.Soc.Am,80:169617247 Bouchon M, Aki K. 1977. Discrete wave-number representation of seismic-source wave fieldsJ. Bull.Seismol.Soc.Am, 67

32、:2592778 Bouchon M. 2003. A review of the discrete wavenumber methodJ. Pure.Appl.Geophy, 160:4454659 Cormier V F, Richards PG. 1977. Full wave theory applied to a discontinuous velocity increase:the inner core boundaryJ. J.Geophys, 43:32210 Levander R. 1988. Fourth-order finite-difference P-SV seism

33、ogramsJ.Geophysics, 53:1425143611 Igel H,Mora RRiollet B. 1995. Anisotropic wave propagation through finite-differencegridsJ.12 Zhang W, Chen X F. 2006. Traction image method for irregular free surface boundaries in finite difference seismic wave simulationJ. Geophys.J.Int,167:33735313 Smith W D. 19

34、75. The application of finite element analysis to body-wave propagation problemsJ. Geophys.J.R.Astr.Soc,42:74776814 Wu R S. 1989. The perturbation method in elastic wave scatteringJ. Pure.Appl.Geophy, 131:60563715 Cerveny V Molotkov I A. 1977. ray method in seismologyM. Karlova University: Praha,798

35、616 Dmowska R. Advances in wave propagation in heterogeneous earth M. 2007. London: Academic Press, 160617 Cerveny V 1983. Synthetic body wave seismograms for laterally varying layered structures by Gaussian beam methodJ. Geophy.J.R.astr.Soc, 73:38942618 Chapman C H, Drumnood R. 1982. Body wave seim

36、ograms in inhomogeneous media using Maslov asymptotic theoryJ . Bull.Seism.Soc.Am , 72(6):27731719 朱良保.1993. Maslov渐近理论地震图D.北京：中国地震局地球物理研究所,178Zhu L B. Maslovs asymptotic theoretical seismogramPh.D.thesis(in Chinese).Beijing:Institute of Geophysics, China Earthquake Administration,199320 Scott P, He

37、lmberger D. 1983. Applications of the Kirchhoff-Helmholtz integral to problems in seismologyJ. Geophys.J.R.astr.Soc, 72:23725421 Hilterman F J. 1970. Three dimensional seismic modelingJ. Geophysics, 35:10201037 Geophysics ,60:1203121622 Berryhill J R. 1979. Wave-equation datumingJ. Geophysics, 44:13

38、29134423 Carter J A, Frazer L N. 1983. A method for modeling reflection data from media with lateral velocity changesJ. J.geophys.Res, 88:646964724 Frazer L N, Sen M K. 1985. Kirchhoff-Helmholtz reflection seismograms in a laterally inhomogeneous multi-layered elastic medium - I .TheoryJ. Geophys.J.

39、R.astr.Soc, 80:12114725 Sen M K, Frazer L N. 1987. Synthetic seismograms using multifold path integrals-II .TheoryJ. Geophys. J. R. astr. Soc,88:64767126 Aki K, Richards P G. 1980. Quantitative Seismology Theory and MethodsM. 1 &2, San Francisco: W H Freeman, 1 932.27 Ziolkowski R W, Deschamps G A.

40、1980. The Maslov method and the asymptotic Fourier transform:caustic analysis, Electro,Lab.Sci.Rep.No.80-9,University of Illinois at Urbana-Champain28 Langston C A. 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameterJ. Bull.Seism.Soc,Am, 67:102910

41、50Computing Teleseismic Conversion PS-Waves withKirchhoff-Helmholtz TheoryLI Hong-guang, WU Qing-juAbstractInstitute of Geophysics,China Earthquake Administration,Beijing 100081,ChinaWe synthesize the teleseismic conversion waves seismograms in a laterally inhomogeneous elastic medium with an asympt

42、otic method namedKirchhoff-Helmholtz integral which do not break down when there is a caustic on the reflector . The accuracy of this approximation is justified by comparing the numberical results with the method of reflectivity and Dynamic ray tracing. KH method can work quickly on curved or irregular interface, so it will be an efficacious method for researching the structure of the crust and upper mantle.Keywords KH method, Dynamic ray tracing, Reflectivity method, Teleseismic conversion waves