地震波动方程[教学知识]

第三章 地震波动方程现在，我们用前一章提出的应力和应变理论来建立和解在均匀全空间里弹性波传播的地震波动方程。这章涉及矢量运算和复数，附录2对一些数学问题进行了复习。3.1 运动方程（Equation of Motion）前一章考虑了在静力平衡和不随时间变化情况下的应力、应变和位移场。然而，因为地震波动是速度和加速度随时间变化的现象，因此，我们必须考虑动力学效应，为此，我们把牛顿定律（）用于连续介质。3.1.1一维空间之振动方程式质点面上由于应力差的存在而使质点产生振动。如图1-3所示，考虑一薄棒向x轴延伸，其位移量为u：Fig3-1则其作用力为“应力”X“其所在的质点面积”，所以其两边的作用力差为惯量inertia为所以得出. （3-1）其中为密度density，为应力stress=。3-1式表示，物体因介质中的应力梯度stress gradient而得到加速度。如果与E为常数，则3-1式可写为 （3-2）其中运用分离变量法求解（3-2）式，设u=F(x)T(t),(3-2)式可以变为设则可得：考虑欧拉公式： （3-3）其中A,B,C,D为根据初始条件和边界条件确定的常数。考虑到可正可负，方程式的解具有的形式，其中f及g为波的函数，以c的波行速度向+x与-x方向传递。我们可以采用如下程序模拟地震波的传播。平面波在均匀介质里沿方向传播，剪切波的齐次微分方程可表达为：这里是位移。对100公里的波长和假定的情况，我们写出用有限差分法解这方程的计算机程序。用长度间距，时间间距秒。假定在（50公里）震源时间函数的形式为： 05秒用（0公里）的应力自由边界条件和（100公里）的固定边界条件。用有限差分图解来近似二次导数：以4秒的间隔画出1-33秒的图。M = moviein(101);dx=1;dt=0.1;tlen=3;beta=4; %初始化变量,tlen为震源持续时间，beta为波传播的速度u1=zeros(101,1);u2=u1;u3=u1;%u1为前一个时刻的各点的位移，u2为当前时刻的位移，u3为下一个时刻的位移值，开始均假定为零t=0;jj=0;while (t=33) %模拟的最长时间为33秒 for ii=2:100 rhs=beta2*(u2(ii+1)-2*u2(ii)+u2(ii-1)/dx2; %方程的解 u3(ii)=dt2*rhs+2*u2(ii)-u1(ii); %对时间求导数 end %左边为自由边界条件，右边为固定边界条件 u3(1)=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).2; %地震震源时间函数 end for ii=1:101 u1(ii)=u2(ii);u2(ii)=u3(ii); %时刻的更新 end plot(u2); %绘制目前的波形图 ylim(-1.2 1.2); M(:,jj+1) = getframe; %获得当前的图像 t=t+dt; %时间延长endmovie(M) %演示波形传播3.1.2三维空间之振动方程式推导三维空间之振动方程式的过程，与上节中所采用的一维空间讨论方式类似，如图3-2所表示，先探讨在x方向之位移量u：Fig3-2在y-z面上的作用力差为：在x-z面上的作用力差为：在x-y面上的作用力差为：惯量为：得出. 3-4其中xx、yx及zx分別为stress tensor在xxx面方向、x力方向，yxy面方向、x力方向及zxz面方向、x力方向方向的分量。注意，在本讲义中有关stress tensor的两个下标indexes之定义，依序为面的方向与力的方向。将xx、yx及zx与其对应的应变之关系代入3-4式可推导得出三维空间之振动方程式如下：. 3-5a其中及为常数，而为Laplacian operator，代表。以相同的方法，可以得出在y及z方向的振动方程式，若其位移量分別为v与w，则其相对应之振动方程式可分別表示如下：. 3-5b. 3-5c若以向量形式来统一表示3-5a、b、c式，可改写如下：. 3-6其中为位移向量，在x、y与z方向的位移分量分別为u、v与w。其中为体力，只有在研究震源时，才考虑该体力。这是构成许多地震学理论基础的基本方程，称之为连续介质方程或运动方程。体力通常包括重力项和震源项。在正常模型地震学中，重力项是频率很低时的一个重要因子，但对所观测到的典型波长范围，即在体波和面波的计算中，通常可被忽略。在这本书后面我们将考虑震源项。在没有体力的情况下，有齐次运动方程： （3-7）在场论中考虑到： （3.8）将其变为更常用的形式，即： （3.8）将这个式子代入（3.12）得到：上式决定了在震源区以外，地震波的传播。解真实地球模型的上述方程是地震学的重要部分，这样的解给出了离震源某一距离的特定地点预期的地面运动，通常称为合成地震图。3.1.3体波纵波与橫波之振动方程式首先，我们考虑由介质伸缩所衍生的质点体积应变之振动方程式。从上节所描述的单一方向x、y、z上之位移量u、v、w所导出的振动方程式，可以进一步地推求体积应变所引发的振动方程式，由的基本定义可以很自然的联想到分別将3-4a、3-4b以及3-4c三式分別对x、y与z微分之后再相加，忽略体力，即可得到下式： 3-7另外，考虑剪切应变可能产生的振动方程式。若将3-5c式对y微分、3-5b式对z微分，然后相減，忽略体力可得到下式：. 3-7其中括弧內的项就是质点运动绕x轴的扭转角度。yFig3-3参考图3-3，一个质点Py、z向逆时针方向扭转到Py、z，扭转角度为x，若其扭转半径为r，根据几何关系可得到：，其位移形变为将其分別对y及z微分且相加，得出同理得到和，所以质点扭转的运动方程式可写为： 3-83-6式与3-8式可用通式描述如下：. 3-9其为典型之波动方程式。根据对3-8式而言，可得出. 3-10对3-6式而言，可得出. 3-113-10可视为纵波亦称为P波，因其质点运动方向与波的传播方向相同如图3-4。 质点运动方向 波传方向Fig3-41-24视为橫波亦称为S波，因为扭转应变，其质点的转动方向与波的传播方向成正交。S波依其质点振动方向的不同可分为SV及SH，如图3-5所示。 波传方向Fig3-5综合以上所得，在完全弹性介质中，当其受外力作用时，产生两种波相：纵波与橫波。由前节所述之各弹性系数的关系，我们可将3-10式以及3-11式写为： ； 其他弹性系数与速度的关系如下：. 3-12 3-13 3-14 3-15其中3-13式可化为. 3-16在地函內部，大部分的泊松比接近于1/4。若=1/4，则 ， 而且若=1/2，即介质为纯液体，则、及皆为零地震所产生之弹性波，穿过地球內部，藉由弹性波传播所产生的速度变化，参考弹性理论以及弹性系数关系，我们可以探索地球內部的情況。 3.1.4 地震波的势位移往往可以根据P波的标量势和S波的矢量： （3.25）那么有： （3.26）将其代入，得到： （3.27）将（3.27）代入可得： （3.29）P波的解由的标量波动方程给出，S波的解由的矢量波动方程给出。3.3 平面波 式（3.28）和（3.29）具有相同的形式，它们在直角坐标系可以表示为：我们用分离变量法来寻找形式的解。每个因子是仅仅一个变量的函数。由上式可得：这意味着是常数，令其为可得：同理，对于某常数,有应注意，因此解可由三个量，而不是四个量来表示。类似于一维形式的推导。该方程可以有如下形式的通解：其中，令下面我们看看的物理意义。令当t=t1时，当t=t2时，由平面解析几何知识可知第一式为离原点距离为的平面，第二式为离原点距离为的平面，并且两平面的法线方向都为。因此两平面之间的距离为，为波从t1时刻传播到t2时刻所传播的距离，传播的速度恰为c，这也是为什么我们在波动方程中将其称之为速度的原因。类似地，表示以速度c向-n方向传播的平面波。任意函数都可以写成简谐平面波叠加的形式根据Fourier叠加原理，可以把屋里上实际存在的平面波动，以数学形式分解成抽象的、覆盖整个频率范围的平面波的积分来表示：实际问题不考虑。因此通常取为方程的基本解。而为波传播的方向，由于c为波的传播速度，通常称为慢度矢量。对不同的做Fourier叠加即可得到任意函数形式的平面波。引进平面波的概念很有帮助。平面波是一个位移只在波的传播方向上变化，在与波传播方向相互垂直的方向上，位移为常数的波动方程的解。例如，沿轴传播的波，位移可表达为： （3.30） 这里是波的速度，是任意函数（矢量函数需表达出波的偏振），这波沿方向传播。位移不随变化。在方向上，波无限扩展。如果是离散的脉冲，那么假定有以平面波阵面传播的位移脉冲形式。更普遍地说，在位置矢量处，平面波在单位矢量方向传播的位移可表达为： （3.31） （3.32）这里是慢度矢量，它的值是速度的倒数。由于地震能量通常由局部的震源辐射出来，地震波阵面总有某种程度的弯曲。然而，在离震源足够大的距离，波阵面平坦到足以使平面波的近似在局部上是正确的。因此，许多解地震波动方程的方法总是把整个解表达为不同传播角度的平面波的和。往往通过变换到频率域，从方程中去掉与时间的依赖关系。在这种情况下，可以把特定角频率的位移表达为： （3.33） （3.34）这里叫做波数矢量。在这本书中，我们将用复数来表示谐波。其详细情况在附录2中作了复习。把谐波称为单色的平面波，有时也把它叫做调和的或稳态平面波解。用来描述这样的波的其他参数是波数，频率，周期和波长。波数为单位长度内波的震动次数。 在波的传播过程中，某一振动状态(周相)在单位时间内传播的距离为波速c，因此波速又叫做相速。应注意介质中各质点的振动速度和波的传播速度c是两个完全不同的概念。振动速度由震源确定，它是周期性变化的，而波速的大小只与介质性质有关。将不同的谐波参数归纳于表3.1。表3.1谐波参数角 频 率频 率周 期速 度波 长波 数3.4 P波和S波的偏振考虑沿方向传播的P波，根据（3.28）式有： （3.35）可以把（3.35）式的解写成： （3.36）这里减号相应于沿方向传播，加号相应于沿方向传播。因为，故有： （3.37）注意对沿方向传播的平面波，在和方向没有变化，所以空间导数和为零。对P波仅在沿轴波的传播方向上有位移。这样的波叫做纵波。而且因为，运动是不旋转的，或“无旋”的。由于P波使介质体积发生变化，所以P波也叫“压缩”波或“膨胀”波。然而，要注意的是P波包括剪切和压缩， 这是为什么P波速度对体积模量和剪切模量反应都灵敏的原因。实际的P波谐振运动可以用图3.2来说明。图3.2 沿页面水平传播的谐振平面P波（上面）和S波（下面）的位移。S波纯剪切，没有体积变化。而P波包括材料体积的变化和剪切（形状变化）。相对于地球实际的应变，这里应变被放大。现在考虑沿正方向传播的S波，矢量势为： （3.38）位移为： （3.39）这里我们再用，即给出： （3.40）运动在和方向，垂直于传播方向。S波的实际运动往往可以分成两个分量：在含传播矢量的垂直面里的运动（波）和取向与这个面垂直的水平运动（SH波）。因为，运动是纯剪切的，没有任何的体积变化（因此叫做剪切波）。在垂直方向偏振的剪切谐波（波）的质点运动如图3.2所示。3.5 球面波如果我们假定球对称，P波势中的标量波动方程（3.28）就可能有另外的解。在球坐标系里，拉普拉斯方程为： （3.41）因为球对称，这里去掉角的偏导数，由表达式（3.28），即得到： （3.42）在点以外，方程的解可表达为： （3.43）注意到除了因子外，这与平面波方程（3.30）是相同的。分别用+和-号表示向内和向外传播的波。因为这个表达式通常用来模拟从点源辐射的波，所以在正常情况下，项表示波的振幅随距离衰减的几何扩散因子，在第6章将进一步的探讨。在时，（3.43）不是方程（3.42）的正确的解。然而，这表明（例如Aki和Richards,4.1节），（3.43）可能是以下非齐次方程的解： （3.44）这里函数在以外的任何地方都为零，它的体积积分为1。因子表示在震源时间函数。在第9章讨论震源理论时，我们将回到这个方程上来。平面波的反射和折射 地壳及地球内部是成层结构，内部有不少分界面。地表也可看作一个界面，震源在各向同性的均匀介质中产生的地震波波阵面是成球形的一层一层向外传播，称为球面波。因此，严格来讲，我们应该讨论球面波遇到分界面时的情况。但当距离震源足够远时，也就是说震源到接收点的距离比波长大得多时，作为一种近似，可讨论平面波在分界面上的行为。同时当（为分界面的曲率半径），也可以将分界面看作平面，这样可使讨论大大简化而不影响对许多现象本质的揭示。同时，球面波在理论上可以看作是许多不同方向的均匀或不均匀的平面波的叠加，因而先弄清了平面波在分界面上的行为，也比较容易讨论球面波在分界面的行为。P波、SV波设平面波（指均匀的平面波）的传播方向在xz平面内，传播方向就是波阵面的法线方向，波的位移场可以表示为： （1）其中满足压缩波的波动方程，满足剪切波的波动方程.由于均匀平面波波阵面上的为常数，而这里平面波传播方向在xz平面内，因此垂直于xz平面的直线上的各点必在同一波阵面内，也就是：。P波产生的位移为：P波产生的应力为SV波的位移SV波产生的应力为：将上面两式代入（1）式得：分析界面条件，界面应力为：界面条件为界面两边应力相等，位移连续，即：分析位移场在y方向的分量,也就是v全部为横波场的分量。再由界面应力条件看，v只出现的表达式中，而u,w只出现在的表达式中。因此，SH波和P-SV波产生的波场是分离的。地球表面是一个特殊的分界面，它将无限介质划分为两个半空间，地面以上空气介质，其密度与地面以下的岩石或海平面以下的海水层相比可以忽略。地球表面可以看成是一个弹性半空间表面，表面下面视为理想弹性介质，表面上面为空气，这种界面称为自由界面，自由界面上的应力作用为零。本节中将介绍弹性波在自由表面上的反射。P波在自由界面的反射如图所示，取xoy平面为自由表面，设有一P波自下部介质入射到自由表面上，由于自由表面以上不存在介质，所以当波遇到自由表面时，只可能折回到原来的介质，而不会透过它，即只存在反射被而不存在透射波。当P波入射到自由表面上时，为满足自由表面处的边界条件，反射波中会同时产生P波和SV波两种成分，此时，SV波称为转换波。但是，由于SH波的振动方向与P被和SV波的振动方向是相互独立的，所以反射波中不会产生SH波。设入射P波为平面简谐波，入射面为xOz平面，法线为z轴，入射P波的入射角为,反射SV波的反射角为，由图中各波的传播方向与坐标轴方向的关系，它们的波函数可以写为：这里只考虑分量，这是由于只有产生xoz平面的振动。式中，由边界条件可知，在z=0处，方程为的线性组合（其中由于z=0，指数因子中的z因子全为零）。所以必有，因此必有：这就是Snell定律，回忆一下几何光学，可见上式与几何光学中的折射定律和反射定律完全一致，这是由于它们在本质上（波动性）有相同之处。而折射反射定律正是反映了物质的波动相关的一种规律。在光学中是从光学实验或惠更斯原理得到了折射反射定律，而这里我们从波动方程和边界条件出发也得到了它。我们在以后的推导中令上式为常数p。则波函数可以写为：则：P波产生的位移为：P波产生的应力为SV波的位移SV波产生的应力为：根据边界条件，可得：对于正应力：对于剪应力：将入射波反射波的势的表达式代入可得：由第二个式子可得：，代入第一式得到折射系数：反射系数为：由于、，代入上面的式子可得到：位移位的振幅并不表示质点的振幅，不具有实际物理意义，下面讨论作为位移振幅比的反射系数。对于稳态传播的P波，位移振幅为：；对于稳态传播的S波振幅为：。我们可以举例说明上面的式子成立，如对于上面所表示的入射波：其合成振幅为：对于上面提到的SV波： 其合成振幅为：由此可知，入射P波在做自由界面上的反射P波位移反射系数与势反射系数相同，而反射SV波的反射系数为势反射系数的倍，即假定SV波入射到自由表面上，其势振幅为A,入射角为，反射SV波的势振幅为B，由反射定律可知其反射角为，反射P波的势振幅为C，反射角为，则根据前面P波和SV波产生的势的定义式和表面应力条件可得：从而得到：由第一个式子可得：，代入第二式得到折射系数：将其带入上面的式子得由于、，代入上面的式子可得到：考虑势振幅和位移振幅之间的关系，可得SH波在自由界面上的反射设入射SH波的位移为：反射SH波位移分别表示为：边条件为：将其简化为：，即在自由表面SH波的反射系数为1.从前面的讨论可以看出，当一列P波入射到自由表面时，会产生一列反射P波和一列反射SV波；同样，如果一列SV波向自由表面入射，会产生一列反射SV波和一列反射P波。或者说，在一般反射问题中半空间内至少存在三列简谐平面波（纯SH波仅反射SH波）。如果我们令 式中的分子为零，则转换波的振幅为零，半空间中只存在一列反射波。即P波入射只反射SV波，SV波只反射P波，这种现象称为偏振交换。自由界面上的位移，视出射角地面测量得到的是地面的实际位移，也就是自由表面的位移。入射波射到自由表面后由于产生了反射波，因而自由表面上的位移并不等于入射波的位移，这是十分重要的。对于P波我们称自由表面位移向量与界面法线的夹角为视入射角。称自由表面上的位移向量与地面之间的夹角为视出射角。当P波入射时，有将P波入射反射为P波和SV波的势函数，并采用反射系数可得因此，由地震记录可得到P波入射到地面后地面位移的北南、东西与垂直分量，求北南、东西分量的平方和再开方得到地面的水平分量，而水平分量和垂直分量的比值就是 ，的一半即为SV波的反射角，根据折射定律即可求得，即真入射角：或：当SV波入射到自由表面时，其真入射角为，仿效这个式子，我们定义这种情况下的视入射角为，则在推导时应注意。当SH波入射到自由表面时，根据前面的推导，反射系数为1.我们同样可设入射SH波的位移为：总的位移为在z=0的面上即得到即自由表面的位移为入射SH波位移的2倍。全反射当入射波为SV波时，由折射定律有：由于，因此反射P波的反射角大于入射SV波的入射角，当入射角满足时，反射P波的反射角为90o。当时，根据有关复数的知识，这时必为复数，且有这里为双曲函数，并采用了欧拉公式的推论：。的数值由确定。则代入波的势函数可得反射P波的势函数为：易见这时反射P波为不均匀波，它的振幅随z的增大而按指数衰减，此不均匀波是沿x轴，也就是沿界面传播的，传播的相速度为。练习：3.1 当波长用以下参数表示时：波数，速度，频率，时间，不是以上任一个，给出与角频率的关系。3.2 考虑在均匀介质里沿方向传播的两类单色平面波：波，对波，在方向有位移的波，即。对每种情况，导出应力张量的非零分量的表达式。借助于（2.12）得到应变张量的分量，然后用（2.24）得到应力分量。3.3 假定谐振的波以在固体里传播。如果最大的应变是，那么周期为1秒，10秒，100秒的质点最大位移是多少？解：根据，则应变为因此，由此可以计算各个频率的最大位移A。3.4 离开点源传播的S波可能球对称吗？在什么条件下，爆炸会产生S波？3.5 说明（3.43）满足（3.42）。因此，（3.43）满足（3.42）3.6 （计算），平面波在均匀介质里沿方向传播，剪切波的齐次动量方程可表达为：这里是位移。对100公里的波长和假定的情况，写出用有限差分法解这方程的计算机程序。用长度间距，时间间距秒。假定在（50公里）震源时间函数的形式为： 05秒用（0公里）的应力自由边界条件和（100公里）的固定边界条件。用有限差分图解来近似二次导数：以4秒的间隔画出1-33秒的图。检验一下速度为4公里/秒的脉冲行程。在每个端点，反射脉冲出现什么情况？当脉冲交叉时，出现什么情况？3.3什么是偏振交换？3.4试给出P波入社到自由表面时的反射P波和反射SV波的相对振幅公式，并给出各个变量的意义。3.5试给出SV波入社到自由表面时的反射P波和反射SV波的相对振幅公式，并给出各个变量的意义。3.6试证明P波入射到自由表面时视出射角为反射SV波反射角的2倍。3.7试证明SH波入射到自由表面时，自由表面的位移为入射SH波位移的2倍。3.8 试证明SH波入射到自由表面时，反射系数为1.50练题