地震波动方程

1、第三章地震波动方程现在，我们用前一章提出的应力和应变理论来建立和解在均匀全空间里弹性波传播的地震波动方程。这章涉及矢量运算和复数，附录2对一些数学问题进行了复习。3.1 运动方程（EquationofMotion）前一章考虑了在静力平衡和不随时间变化情况下的应力、应变和位移场。然而，因为地震波动是速度和加速度随时间变化的现象，因此，我们必须考虑动力学效应，为此，我们把牛顿定律（Fma）用于连续介质。3.1.1 一维空间之振动方程式质点面上由于应力差的存在而使质点产生振动。如图1-3所示，考虑一薄棒向x轴延伸，其位移量为u:Fig3-1则其作用力为“应力”X"其所在的质点面积”，所以其

2、两边的作用力差为dsxdxxdxdsx惯量(inertia)为2udxds-2t2所以得出2u，7.(3-1)ustress)=E°xtx其中p为密度（density）,<r为应力（3-1式表示，物体因介质中的应力梯度（stressgradient）而得到加速度。如果p与E为常数，则3-1式可写为22u1u-222XCt(3-2)运用分离变量法求解(3-2)式，设u=F(x)T(t),(3-2)式可以变为1XTXTc、.c2XT2设XT一i-则可得：Teit,Xec考虑欧拉公式：eitcos(t)isin(t),eitcos(t)isin(t)i_xcti_xcti_xctix

3、ctuAecBecCecDec(3-3)其中A,B,C,D为根据初始条件和边界条件确定的常数。考虑到可正可负，方程式的解具有ufxctgxct的形式，其中f及g为波的函数，以c的波行速度向+x与-x方向传递。我们可以采用如下程序模拟地震波的传播。剪切波的齐次微分方程可表达为：2_u下这里u是位移。对100公里的波长和假定差分法解这方程的计算机程序。用长度间距定在u(50公里)震源时间函数的形式为:u50tsin2t5用u(0公里)的应力自由边界条件和u差分图解来近似二次导数：2uui12x以4秒的间隔画出1-33秒的图。M=moviein(101);dx=1;dt=0.1;tlen=3;bet

4、a=4;%初始化变量,tlenu1=zeros(101,1);u2=u1;u3=u1;平面波在均匀介质里沿x方向传播，22u2x4公里/秒的情况，我们写出用有限dx1公里，时间问距dt0.1秒。假0Vt<5秒(100公里)的固定边界条件。用有限25ui1dx2为震源持续时间，beta为波传播的速度%u伪前一个时刻的各点的位移，u2为当前时刻的位移，u3为下一个时刻的位移值，开始均假定为零t=0;jj=0；while(t<=33)%模拟的最长时间为33秒forii=2:100rhs=betaA2*(u2(ii+1)-2*u2(ii)+u2(ii-1)/dxA2;%方程的解u3(ii)

5、=dtA2*rhs+2*u2(ii)-u1(ii);%对时间求导数end%左边为自由边界条件，右边为固定边界条件u36=u3(2);%左边为自由边界条件u3(101)=0.0;%右边为固定边界条件%左右两边为自由边界条件%u3(1)=u3(2);%左边为自由边界条件%u3(101)=u3(100);%右边为自由边界条件%左右两边为固定边界条件%u3(1)=0.0;%左边为固定边界条件%u3(101)=0.0;%右边为固定边界条件if(t<=tlen)u3(51)=(sin(pi*t/tlen).A2;%地震震源时间函数endforii=1:101u1(ii)=u2(ii);u2(ii)=

6、u3(ii);%时刻的更新endplot(u2);%绘制目前的波形图ylim(-1.21.2);M(:,jj+1)=getframe;%获得当前的图像t=t+dt;%时间延长endmovie(M)%演示波形传播Fig3-2在y-z面上的作用力差为:xxdxxxXdydz在x-z面上的作用力差为:yxdyyxydxdz在x-y面上的作用力差为:zxdzzxzdxdy惯量为:dxdydz得出2ut2xxyxzxfx.(3-4)其中xxx、（Tyx及（Tzx分别为stresstensor在xx（x面方向、x力方向），yx（y面方向、x力方向）及zx（z面方向、x力方向）方向的分量。注意，在本讲义中有

7、关stresstensor的两个下标（indexes）之定义，依序为面的方向与力的方向。将（Txx、byx及bzx与其对应的应变之关系代入3-4式可推导得出三维空间之振动方程式如下:2ut22ufx.(3-5a)其中入及科为常数，而22、为Laplacianoperator,代表2xy2以相同的方法，可以得出在y及z方向的振动方程式，若其位移量分别为v与w,则其相对应之振动方程式可分别表示如下:1.3 三维空间之振动方程式推导三维空间之振动方程式的过程，与上节中所采用的一维空间讨论方式类似，如图3-2所表示，先探讨在x方向之位移量u:dx2vt2fy(3-5b)若以向量形式来统一表示3-5a、

8、b、graddivufz.(3-5c)c式，可改写如下:2uf.(3-6)其中u为位移向量，在x、y与z方向的位移分量分别为u、v与w。其中fx,fy,fz为体力，只有在研究震源时，才考虑该体力。这是构成许多地震学理论基础的基本方程，称之为连续介质方程或运动方程。体力f通常包括重力项fg和震源项fs。在正常模型地震学中，重力项是频率很低时的一个重要因子，但对所观测到的典型波长围,即在体波和面波的计算中，通常可被忽略。在这本书后面我们将考虑震源项fs在没有体力的情况下,有齐次运动方程:在场论中考虑到:graddivu(3-7)2u(3.8)将其变为更常用的形式，即:(3.8)将这个式子代入（3.

9、12)得到:2ut22graddivurot(rotu)上式决定了在震源区以外，地震波的传播。解真实地球模型的上述方程是地震学的重要部分，这样的解给出了离震源某一距离的特定地点预期的地面运动，通常称为合成地震图。1.3 体波（纵波与横波）之振动方程式首先，我们考虑由介质伸缩所衍生的质点体积应变之振动方程式。从上节所描述的单一方向（x、v、z）上之位移量（u、v、w）所导出的振动方程式，可以进一步地推求3-4a、3-4b体积应变所引发的振动方程式，由的基本定义可以很自然的联想到分别将以及3-4c三式分别对x、y与z微分之后再相加，忽略体力，即可得到下式:t2(3-7)另外，考虑剪切应变可能产生的

10、振动方程式。若将3-5c式对y微分、3-5b式对z微分，然后相减，忽略体力可得到下式:t2y其中括弧内的yv项就是质点运动绕x轴的扭转角度。z参考图3-3,一个质点P（y、z）向逆时针方向扭转到P'（y、z）,扭转角度为x,若其扭转半径为r,根据几何关系可得到:rsin其位移形变为vrxsinwrxcos将其分别对y及z微分且相加,得出同理得到z1,所以质点扭转的运动方程式可写2xy为:t22yF2(3-8)zF3-6式与3-8式可用通式描述如下:(3-9)2t7其为典型之波动方程式。根据对3-8式而言,可得出cV1(3-10)(3-11)对3-6式而言,3-10可视为纵波（亦称为P波

11、），因其质点运动方向与波的传播方向相同（如图3-4）。Particlemotion质点运动方向-*Ray波传方向Fig3-41-24视为横波（亦称为S波），因为扭转应变，其质点的转动方向与波的传播方向成正交。S波依其质点振动方向的不同可分为SV及SH,如图3-5所示。SHwave®qSVwave波传方向*RayFig3-5波。由前节所述之各弹性系数的关系，我们可将3-10式以及3-11式写为:综合以上所得，在完全弹性介质中，当其受外力作用时，产生两种波相：纵波与横ViV2其他弹性系数与速度的关系如下:_223Vl24V222.(3-12V1V2122(3-13VV2-V3(3-14其

12、中3-13式可化为2V1V2V22(3-15.(3-16在地函内部，大部分的泊松比b接近于1/4。若er=1/4,则2V1V2V11.732V2而且5V；65V2225V125V22若（7=1/2,即介质为纯液体，则、E及V2皆为零V12地震所产生之弹性波，穿过地球内部，藉由弹性波传播所产生的速度变化，参考弹性理论以及弹性系数关系，我们可以探索地球内部的情况。3.1.4地震波的势位移U往往可以根据P波的标量势和S波的矢量也：(3.25)那么有:(3.26)将其代入t2t2(-y(-z(根据3.13)(因为0)VUQ次()?(3.27)将(3.27)代入x2xt212P波的解由3.3平面波式（3

13、.28）和（3.29）具有相同的形式，它们在直角坐标系可以表不为:2f2'f2f2f2c22-2txyz我们用分离变量法来寻找X(x)Y(y)Z(z)T(t)形式的解。每个因子是仅仅一个变量的函数。由上式可得:2,222222cdXcdYcdZ1dT2"222Xdx2Ydy2Zdz2Tdt2这意味着19是常数，令其为2可得:Tdt2d2Tdt22T0,Te同理，对于某常数kx，ky,有0,ikxtd2Xdx2d2Ydy2kyY0,ikyyed2Zdz20,ik7ze应汪息,k22-2ck；,因此解可由三个量，kx,ky,而不是四个量来表示。类似于一维形式的推导该方程可以有如下

14、形式的通解:fx,y,z,tfitnxxnyy“zf2nxxnvynzzxyzc其中kx%,kycny.nz222一,kz一，nxnynz1,令nxxnyyn%znrcc下面我们看看fit%xnyy”z的物理xXnyy型tconstac当t=t1时，nXxnyynzzc(t1a)当t=t2时，nXxnyynzzc(t2a)由平面解析几何知识可知第一式为离原点距离为c(t1a)的平面，第二式为离原点距离为c(t2a)的平面，并且两平面的法线方向都为nx,ny,nz。因此两平面之间的距离为c(t2ti),为波从ti时刻传播到t2时刻所传播的距离，传播的速度恰为c,这也是为什么我们在波动方程中将其称

15、之为速度的原因。类似地，f2t”n"一生表示以速度c向-n方向传播的平面波。c任意函数都可以写成简谐平面波叠加的形式根据Fourier叠加原理，可以把屋里上实际存在的平面波动，以数学形式分解成抽象的、覆盖整个频率围的平面波的积分来表示：rijXjnjXjdt实际问题不考虑。it”因此通常取f(r,t,)Aec?由于c为波的传播速度，通常称n为慢度矢量。c到任意函数形式的平面波。为方程的基本解。而？为波传播的方向，对不同的A,做Fourier叠加即可得引进平面波的概念很有帮助。平面波是一个位移只在波的传播方向上变化，在与波传播方向相互垂直的方向上，位移为常数的波动方程的解。例如，沿X轴

16、传播的波，位移可表达为：uX,tft%(3.30)这里c是波的速度，f是任意函数(矢量函数需表达出波的偏振)，这波沿X或X方向传播。位移不随y或z变化。在丫、Z方向上，波无限扩展。如果f(t)是离散的脉冲，那么假定u有以平面波阵面传播的位移脉冲形式。更普遍地(3.31)(3.32)说，在位置矢量x处，平面波在单位矢量§方向传播的位移可表达为:u(x,t)f(t?x/c)ftsx这里s8/c是慢度矢量，它的值是速度的倒数。由于地震能量通常由局部的震源辐射出来，地震波阵面总有某种程度的弯曲。然而，在离震源足够大的距离，波阵面平坦到足以使平面波的近似在局部上是正确的。因此，许多解地震波动方

17、程的方法总是把整个解表达为不同传播角度的平面波的和。往往通过变换到频率域，从方程中去掉与时间的依赖关系。在这种情况下，可以把特定角频率的位移表达为：ux,tAeitsx(3.33)Aeitkx(3.34)这里ks/c?叫做波数矢量。在这本书中，我们将用复数来表示谐波。其详细情况在附录2中作了复习。把谐波称为单色的平面波，有时也把它叫做调和的或稳态平面波解。用来描述这样的波的其他参数是波数kk/c,频率f/,周期T%和波长cTo波数为单位长度波的震动次数。在波的传播过程中，某一振动状态(周相)在单位时间传播的距离为波速c,因此波速又叫做相速。应注意介质中各质点的振动速度和波的传播速度c是两个完全

18、不同的概念。振动速度由震源确定，它是周期性变化的，而波速的大小只与介质性质有关。将不同的谐波参数归纳于表3.1表3.1谐波参数2T1TckcT,22fkcc3.4P波和S波的偏振考虑沿x方向传播的P波，根据（3.28）式有:(3.35)2xxtt0tx可以把（3.35）式的解写成:(3.36)这里减号相应于沿x方向传播，加号相应于沿x方向传播。因为U,故有:(3.37)UxxUy0Uz0注意对沿x方向传播的平面波，在y和z方向没有变化，所以空间导数y和z为零。对P波仅在沿x轴波的传播方向上有位移。这样的波叫做纵波。而且因为0,运动是不旋转的，或“无旋”的。由于P波使介质体积发生变化，所以P波也

19、叫“压缩”波或“膨胀”波。然而，要注意的是P波包括剪切和压缩，这是为什么P波速度对体积模量和剪切模量反应都灵敏的原因。实际的P波谐振运动可以用图3.2来说明。波傅方向二图3.2沿页面水平传播的谐振平面P波（上面）和S波（下面）的位移。S波纯剪切，没有体积变化。而P波包括材料体积的变化和剪切（形状变化）。相对于地球实际的应变，这里应变被放大。(3.38)现在考虑沿正x方向传播的S波，矢量势为：巾Xtx?ytx?位移为：(3.39)UxWxzyyUyWyxzzUzWzyxxy(3.40)这里我们再用yz0,即给出：Uxz?xy今运动在y和z方向，垂直于传播方向。S波的实际在含传播矢量的垂直面里的运

20、动（SV波）和取向与这个面垂直的水平运动（SH波）。因为uW0,运动是纯剪切的，没有任何的体积变化（因此叫做剪切波）。在垂直方向偏振的剪切谐波（SV波）的质点运动如图3.2所示。3.5球面波(3.41)如果我们假定球对称，P波势中的标量波动方程（3.28）就可能有另外的解。在球坐标系里，拉普拉斯方程为：2r因为球对称，这里去掉角的偏导数,由表达式（3.28),即得到:2在0(3.42)在点r0以外，方程的解可表达为:r,t(3.43)、1注意到除了1因子外，这与平面波方程（3.30）是相同的。分别用+和-号表小向r和向外传播的波。因为这个表达式通常用来模拟从点源辐射的波，所以在正常情况下，1项

21、表示波的振幅随距离衰减的几何扩散因子，在第6章将进一步的探r讨。在r0时，（3.43）不是方程（3.42）的正确的解。然而，这表明（例如Aki和Richards,4.1节），（3.43）可能是以下非齐次方程的解：212一/、r=-24rft（3.44）2t2这里r函数在r0以外的任何地方都为零，它的体积积分为1。因子4rft表示在震源时间函数。在第9章讨论震源理论时，我们将回到这个方程上来。平面波的反射和折射地壳及地球部是成层结构，部有不少分界面。地表也可看作一个界面，震源在各向同性的均匀介质中产生的地震波波阵面是成球形的一层一层向外传播，称为球面波。因此，严格来讲，我们应该讨论球面波遇到分界

22、面时的情况。但当距离震源足够远时，也就是说震源到接收点的距离比波长大得多时，作为一种近似，可讨论平面波在分界面上的行为。同时当（为分界面的曲率半径），也可以将分界面看作平面，这样可使讨论大大简化而不影响对许多现象本质的揭示。同时，球面波在理论上可以看作是许多不同方向的均匀或不均匀的平面波的叠加，因而先弄清了平面波在分界面上的行为，也比较容易讨论球面波在分界面的行为。P波、SV波设平面波（指均匀的平面波）的传播方向在xz平面，传播方向就是波阵面的法线方向，波的位移场可以表示为：uUpUs6（1）其中满足压缩波的波动方程，小满足剪切波的波动方程.由于均匀平面波波阵面上的，x,y,z为常数，而这里平

23、面波传播方向在XZ平面，因此垂直于xz平面的直线上的各点必在同一波阵面，也就是:0。xyzyyyyp波产生的位移为:upvpwpzP波产生的应力为zzzxzyezySV波的位移usvusxusyuszsv波产生的应力为:zzsvzxsvzysvezy将上面两式代入zxz2uw:2zxzvwv0zyz1)式得：2wyezz22yyy-2xupus一xvpxvs一zwp分析界面条件，界面应力为:zz2ezzw2zzx2ezx222yy2-Txzxzzy2ezy界面条件为界面两边应力相等，位移连续，即:zzzz',zxzx',zyzyuu',vv',ww'分析

24、位移场在y方向的分量v,vvpVsvs,也就是v全部为横波场的分量。再由界面应力条件看，v只出现xy的表达式中，而u,w只出现在zz,zx的表达式中。因此，SH波和P-SV波产生的波场是分离的。地球表面是一个特殊的分界面，它将无限介质划分为两个半空间，地面以上空气介质，其密度与地面以下的岩石或海平面以下的海水层相比可以忽略。地球表面可以看成是一个弹性半空间表面，表面下面视为理想弹性介质，表面上面为空气，这种界面称为自由界面，自由界面上的应力作用为零。本节中将介绍弹性波在自由表面上的反射。P波在自由界面的反射SV波（出,）如图所示，取xoy平面为自由表面，设有一P波自下部介质入射到自由表面上，由

25、于自由表面以上不存在介质，所以当波遇到自由表面时，只可能折回到原来的介质，而不会透过它，即只存在反射被而不存在透射波。当P波入射到自由表面上时，为满足自由表面处的边界条件，反射波中会同时产生P波和SV波两种成分，此时，SV波称为转换波。但是，由于SH波的振动方向与P被和SV波的振动方向是相互独立的，所以反射波中不会产生SH波。设入射P波为平面简谐波，入射面为xOz平面，法线为z轴，入射P波的入射角为id，反射SV波的反射角为is,由图中各波的传播方向与坐标轴方向的关itkYxk7zBe系，它们的波函数可以写为:1Adtkx这里只考虑y分量，这是由于只有y?产生xoz平面的振动。kx-sinid

26、,kx一sinid,kxsinis式中，kzcosid,kzcosid,kzcosis由边界条件可知，在z=0处，方程为etkxxkzz,eitkxxkzz,eitkxxkzz的线性组合（其中由于z=0,指数因子中的z因子全为零）。所以必有kxkxkx,因此必有：sinidsinidsinis这就是Snell定律，回忆一下几何光学，可见上式与几何光学中的折射定律和反射定律完全一致，这是由于它们在本质上（波动性）有相同之处。而折射反射定律正是反映了物质的波动相关的一种规律。在光学中是从光学实验或惠更斯原理得到了折射反射定律，而这里我们从波动方程和边界条件出发也得到了它。我们在以后的推导中令上式为

27、常数p则波函数可以写为：itpx空zAeA2ecos'tpxSzBecoshpxSz则：P波产生的位移为:upupxWpUpzP波产生的应力为zzp22xz2.2cosid2-2.2cosid2-222cosid-222一2.Sinid2-2cosid2-2cosid2.2cosid2sin2id2sin2id一.一2.sinidzxp2ezx2.ip-z22sinidSV波的位移usvusxSV波产生的应力为:zzsvezz2zxsv2sinis2xz2yTz2y-2xp2.ip2.cosis2-2py根据边界条件，可得:对于正应力：zz|z0zzpzzsv对于剪应力：zx|z0zx

28、p2*2*4ipzz22p22222八2.yc12p122ip0z将入射波反射波的势的表达式代入可得:222COSie124p12p反射系数为:p2AA222psB0c2COSid22.2pA1A212pB022pcosid由第二个式子可得:A；E-21，代入第一式得到折射系数:,42COSisCOSid/c222A2A14p12p42cosiscosidp222122p2sinidsinis122p2sin2i.一2cos2is,代入上面的式子可得到:A2.八.22八.Asin2idsin2iscos2is2二.二22A1sin2idsin2iscos2is2B2sin2idcos2is7-

29、222Asin2idsin2iscos2is位移位的振幅并不表示质点的振幅,移振幅比的反射系数。不具有实际物理意义，下面讨论作为位对于稳态传播的P波，位移振幅为:势振幅一；对于稳态传播的S波振幅为：势振幅一。我们可以举例说明上面的式子成立，如对于上面所表示的入射波:upxvp1.ix0Wpupzcosis其合成振幅为:upwp对于上面提到的SV波:usvuscosisyWsvUs其合成振幅为：Usv由此可知，入射P波在做自由界面上的反射P波位移反射系数与势反射系数相同，而反射SV波的反射系数为势反射系数的一倍，即a2a1b2.八.八.sin2idsin2is2sin2idsin21s22-co

30、s2is22cos2is2sin2idcos2is2.22a1sin2idsin2iscos2is假定SV波入射到自由表面上，其势振幅为A,入射角为is,反射SV波的势振幅为B,由反射定律可知其反射角为is,反射P波的势振幅为C,反射角为ip,则根据前面P波和SV波产生的势的定义式和表面应力条件可得:zx2-ip2222py1y20zz2y1ip-z从而得到:cosis-八p1c0cosissA由第一个式子可得:22p222、1-,代入第二式得到折射系数:cosipA2cosiscosiD4p2spA1c22,2cosiscosip12P4p-将其带入上面的式子得cosis1_24p-2p21

31、22cosiscosiD-22p24p2sp2sin2L一2Acos2is,代入上面的式子可得到:BACA2八八sin2ipsin2issin21Psin2is2cos22is2cos22is22sin2idcos2is222sin2ipsin2iscos2is考虑势振幅和位移振幅之间的关系,可得2-c22八sin2ipsin2iscos2is222sin2ipsin2iscos2isfvp2sin2iscos2is222sin2ipsin2iscos2isSH波在自由界面上的反射设入射SH波的位移为：ixsiniszcosisitSHe反射SH波位移分别表示为：ixsini'szco

32、si'sitS'H'e边条件为：SS'yz|z00z将其简化为：HH,即在自由表面SH波的反射系数为1.从前面的讨论可以看出，当一列P波入射到自由表面时，会产生一列反射P波和一列反射SV波；同样，如果一列SV波向自由表面入射，会产生一列反射SV波和一列反射P波。或者说，在一般反射问题中半空间至少存在三列简谐平面波（纯SH波仅反射SH波）。如果我们令式中的分子为零，则转换波的振幅为零，半空间中只存在一列反射波。即P波入射只反射SV波,SV波只反射P波，这种现象称为偏振交换。自由界面上的位移，视出射角地面测量得到的是地面的实际位移，也就是自由表面的位移。入射波射到自

33、由表面后由于产生了反射波，因而自由表面上的位移并不等于入射波的位移，这是十分重要的。对于p波我们称自由表面位移向量与界面法线的夹角ip为视入射角。称自由表面上的位移向量与地面之间的夹角有e为视出射角。当P波入射时,12yxxz12yzzxtan2is或：cosep90°cosepuupusvtanip|z0wWpWsv将P波入射反射为P波和SV波的势函数，并采用反射系数可得4P2COSidCOSis22P2122P2tan：pxxz_s1n21s12,2122p2cos1122p222p28s21szzx因此Tp2is,由地震记录可得到P波入射到地面后地面位移的北南、东西与垂直分量，

34、求北南、东西分量的平方和再开方得到地面的水平分量，而水平分量和垂直分量的比值就是tanip,ip的一半即为SV波的反射角，根据折射定律即可求得ip,即真入射角：sinipsinis-sinip当SV波入射到自由表面时，其真入射角为is,tanisws,仿效这个式子,Us我们定义这种情况下的视入射角为，则ip(1cosipCAcosiip1BCAp-AwWsiWs2Wptan-UUsiUs2Up222.2cotipp-72cotrs1-sin2issinipcosipsinis""2"7：cos2iscosis在推导时应注意psinippsiniss1nlp,。si

35、nis当SH波入射到自由表面时，根据前面的推导，反射系数为入射SH波的位移为：.我们同样可设xsiniszcosisitSHe总的位移为vSS2ixsiniszcosisitHeixsiniszcosisite在z=0的面上即得到vS1i-xsinisiti-xsinisit$Heei-xsinisit2He即自由表面的位移为入射SH波位移的2倍。全反射当入射波为SV波时，由折射定律有:sinipsinisps由于，因此反射P波的反射角大于入射SV波的入射角is，当入射角is满足issin1一时，反射P波的反射角为90°。issin1一时，sinip1,根据有关复数的知识，这时ip必为复数，且有ip-p20,sinipsin一2sincosi2cos-sinicosicosh2cosipcos2coscosi