地震波的一维波动

数学文化课程报告地震波的一维波动方程正演问题求解学成绩________________________地震波的一维波动方程正演问题求解摘要：文中从地震勘探的实际问题出发，对实际勘探过程中运用到的人工激发地震波法进行简化，建立了一维波动方程，并由条件对波动方程的正问题进行求解，研究了非齐次边界条件转化成齐次边界条件对解定解问题的重要性。

关键词：地震波，一维波动方程，非齐次边界条件，正演一、引言：地震波是指从震源产生向四外辐射的弹性波。

地震发生时，震源区的介质发生急速的破裂和运动，这种扰动构成一个波源。

由于地球介质的连续性，这种波动就向地球内部及表层各处传播开去，形成了连续介质中的弹性波。

地震是照亮地下的明灯.为了用地震波进行远距离探测,地震勘探中通常使用大当量的炸药作为激发源,利用可控震源人工激发地震波是进行石油勘探的一种主要方式。

基本办法是用高爆炸力的TNT炸药在地面激起人工地震波，震波沿着与地面垂直的方向传播,在碰到质地相对致密的岩层以后,一部分波被反射回地面,预先,在地面上安置起许多呈现点阵的检波器,这些检波器能够把地面微弱的震动变成电子信号,通过连接线传输到接收机里,接收机的功能是分道记录不同位置的检波器的电信号，从而得到相关数据用以实际应用。

在此，根据查阅相关资料，建立地震波纵波的一维波动方程，并根据所得的边界条件和初始条件，求解波动方程，得到地震波振幅的解，预测出波场。

二、模型建立：1、波动方程的建立：当有体力f作用于一个连续介质物体时，物体中每一点（无限小立方体）应力运动方程或动量方程是这是大多数地震学所依据的基本方程，其中σ为正应力，τ为剪应力，),,(zyxifi=为体力分量，),,(zyxiui=为位移分量，ρ为介质密度。

当没有体力是可得到齐次的运动方程，它反映了由震源区向外的地震波传播。

为了求解运动方程，需要利用位移u来表示应力τσ，：其中μλ和为拉梅参数，这两个拉梅参数描述了固体内部的应力-应变线性关系。

因为：所以方程（1-1）可以改写成矢量形式：()k Z j Y i X F++=ρρ，ρ为密度，X,Y,Z 为单位质量体力在坐标轴上的投影。

忽略体力后，式（1-4）为地震波波动方程的形式之一，适合于复杂结构的地质模型，但是该方程非常复杂而难以有效的求解。

因此主要有两种方法进行近似计算：①因为地下介质的速度与密度有一些经验关系，所以只考虑速度为密度的函数。

②忽略拉梅参数的梯度项，把介质近似为均匀介质；假设不考虑介质中的旋转形变，剪应力为0。

当地下介质为均匀介质时，式（1 -4）可以简化为根据弹性系数K与拉梅系数λ和μ的关系：μλ+=K,(2-4)式可以化为：假设地下介质是各向同性的水平层状介质，被震源激发的地震波是垂直传播的平面波，则描述地震勘探中地震波的一维波动方程为：由于K，ρ与深度x无关，所以可得地震波一维波动方程：其中：x为深度，t为时间，k为均匀介质的弹性参数，ρ为均匀介质的密度，u(x,t)为振幅,T为地面记录的最大时间，L为所要研究的最大深度。

2、边界条件和初始条件的求解：①初始条件：初始条件就时地震波在开始时刻的位移和速度，则：②边界条件：设g(t)为震源函数，则：三求解定解问题：）（4-3,)(0T t o tf u x≤≤==分析此定解问题，方程是齐次的方程，边界条件为非齐次的边界条件，所以要设法将边界条件化成齐次的边界条件。

具体做法是取一个适当的未知函数之间的代换，使其对新的未知函数，边界条件为齐次的边界条件：令),(),(),(t x W t x V t x u +=(3-4)选取),(t x W 使),(t x V 的边界条件化为齐次边界条件的方程，即：由(3-2)及(3-4)容易看出，要使(3-5)成立，只要为使得计算方便，取W 为x 的一次式，即：)()(),(t B x t A t x W +=，用方程(3-6)确定A ，B 得L t g t B t g t A )()(),()(-==，只要做代换：[]gL gx V u -+= (3-7)，就能使新的未知函数V 满足齐次的边界条件。

经过这个代换后，得到关于V 的定解问题为：其中：[][]⎩⎨⎧--=--=L g x g x L g x g x )0()0(0)()0()0(0)(ψϕ (3-9)接下来求解问题(3-8)：1注：3-8-3定解问题（3-8）的方程为式（3-8-1），边界条件为式（3-8-2），初始条件为式（3-8-3）.第一步，分离变量法：代入方程(3-8-1)得：)()()()(2t T xX a t T x X =，这个式子左端仅是x 的函数，右端仅是t 的函数，一般情况下二者不可能相等，，只有当它们均为常数时才相等。

令此常数为λ-，则有：这样得到两个常微分方程：)()(0)()(2=+=+x X x X t T a t T λλ(3-8-4) 0)()(0)()(2=+=+x X x X t T a t T λλ (3-8-5)利用边界条件(3-8-2)，由于)()(),(t T x X t x v =，故有0)()(,0)()0(==t T L X t T X但0)(≠t T ，所以0)()0(==L X X . (3-8-6)得常微分方程的边值问题：⎩⎨⎧0)()0(,0)()(===+l X X x X x X λ第二步，求解特征值：1°设λ=0，此时方程 (3-8-5)的通解为,)(B Ax x X +=因为 A x X =)(，由条件(3-8-6)得A=B=0，所以0≠λ。

2°设λ<0，此时方程 (3-8-5)的通解为由条件易得A=B=0，)(x X =0，所以舍去。

3°设λ<0，此时方程 (3-8-5)的通解为与这些特征值对应的方程(3-8-4)的通解为得的任意常数利用初始条件确定其中,,n n D C 因为)()(x x φϕ和是定义在[0,L]上的函数，所以只要选取n C 为()x ϕ的傅里叶正级数展开式的系数，将解代入(3-7)式，即可得原问题的解：[]=-+=gLgxVut xv),(Ltgxtg)()(-+四实际应用：录的回波），对地震参数),,(zyxρ与k(x,y,z)进行反演，实际中的地层为非均匀介质，所以ρ与k都不是常数，与空间位置有关。

五雷克子波频谱分析【fft傅里叶快速变换】（c语言）：#include "stdio.h"#include "math.h"#include "stdlib.h"#define PI 3.1415926void fft(float sr[],float sx[],int m0,int inv)；void main(){ void fft();float *xr,*xi;int i,np,nfft,nf,k;float t,dt,hf,df,f,z;FILE *fpar,*fp1,*fp2;char fil1[80],fil2[80];fpar=fopen("spectrum.par","r");fscanf(fpar,"%s",fil1);fscanf(fpar,"%s",fil2);printf("%s\n",fil1);printf("%s\n",fil2);fscanf(fpar,"%d",&np);printf("np=%d\n",np);fscanf(fpar,"%f",&dt);printf("dt=%8.3f ms\n",dt);fscanf(fpar,"%f",&hf);printf("hf=%8.3f Hz\n",hf);dt=dt/1000.0;// calculate fft pointk=log(np)/log(2);if(np>pow(2,k))k=k+1;nfft=pow(2,k);df=1.0/(nfft*dt);nf=hf/df;printf("nfft=%d k=%d\n",nfft,k);printf("dt=%8.4f df=%8.4f\n",dt,df);// allocate memoryxr=(float *)calloc(nfft,sizeof(double));xi=(float *)calloc(nfft,sizeof(double));fp1=fopen(fil1,"r");// input datafor(i=0;i<np;i++){ fscanf(fp1,"%f%f",&t,&z);xr[i]=z;}fclose(fp1);// fill zeroif(np<nfft){for(i=np;i<nfft;i++)xr[i]=0;}for(i=0;i<nfft;i++)xi[i]=0.0;// calculate fftfft(xr,xi,k,1);// output amplitude spectrum of the orignal signal fp2=fopen(fil2,"w");for(i=0;i<nf;i++){ f=i*df;fprintf(fp2,"%10.4f %12.4f\n",f,sqrt(xr[i]*xr[i]+xi[i]*xi[i])*dt);}fclose(fp2);free(xr);free(xi);}void fft(float sr[],float sx[],int m0,int inv) {int i,j,lm,li,k,lmx,np,lix,mm2;double t1,t2,c,s,cv,st,ct;if(m0<0)return;lmx=1;for(i=1;i<=m0;++i)lmx+=lmx; //form 2**m0cv=2.0*PI/(double)lmx;ct=cos(cv); st=-inv*sin(cv);np=lmx;mm2=m0-2;/* fft butterfly numeration */for(i=1;i<=mm2;++i){lix=lmx;lmx/=2;c=ct;s=st;for(li=0;li<np;li+=lix){j=li;k=j+lmx;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=t1;sx[k]=t2; ++j;++k;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;}for(lm=2;lm<lmx;++lm){cv=c;c=ct*c-st*s;s=st*cv+ct*s;for(li=0;li<np;li+=lix){j=li+lm;k=lmx+j;t1=sr[j]-sr[k];t2=sx[j]-sx[k];sr[j]+=sr[k];sx[j]+=sx[k];sr[k]=c*t1-s*t2;sx[k]=s*t1+c*t2;}}cv=ct;ct=2.0*ct*ct-1.0;st=2.0*st*cv;}参考文献：[1] 张丽琴,王家映,严德天利用同伦反演方法进行岩性油气藏研究[期刊论文]-天然气工业 20 05(08)[2] 范留明非均匀层状介质一维波动方程精确解的有限差分算法[期刊论文]-岩土力学 2013(0 9)[3]《数理物理方程与特殊函数》第三版王元明高等教育出版社[4]《信号分析与处理讲义》张军华中国石油大学（华东）未出版对于数学文化课程的认识和建议：老师讲课十分认真投入，内容纲举目分，条理性很强，而且特别善于举例，让同学联系理论与实际，学习起来十分轻松，而且印象深刻，收到良好的效果。