声波方程有限差分法数值模拟

声波方程有限差分法数值模拟学院：海洋地球科学学院专业：地球信息科学与技术年级：2012级姓名：王昊指导教师：宋鹏声波方程有限差分法数值模拟2012级地球信息科学与技术12040032037王昊摘要：本实验应用声波方程作为正演模拟的波动方程，将所提供震源函数离散后绘图，并将给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图形。通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规律，在模拟水平层状速度模型中，体会地震波在两种介质分界面的传播规律，并能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。并通过模拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波形变化趋势。关键词：声波方程数值模拟地震波有限差分引言：地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的方法模拟该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野外数据采集过程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值模拟技术。地震波场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、震源和地下几何界面、物性参数（岩层密度、速度等）情况下，研究弹性波或声波的传播规律。地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。相对于上述几种方法，有限差分法是一种更为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。本实验给定大小两个模型，对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；对于大模型，定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。实验内容：对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222tSzuxuvtu（1）其中v(x,z)是介质在点(x,z)处的纵波速度，本实验中v为常数2500，u为描述速度位或者压力的波场，)(ts为震源函数，本实验选用雷克子波（零相位子波）fttftse2cos)/2()(22上式中，t为时间，fm为中心频率，一般取为20-40HZ，本实验取20，γ为控制频带宽度的参数，一般取3-5,这里取4，在实际计算过程中，需把此震源函数离散，参与波场计算。由图像可知子波能量主要集中在前部，大致上t35的时候S(t)都为零值，故在绘制S(t)图像时k取100，来比较清晰得呈现雷克子波。本实验不涉及吸收边界条件的使用。这样声波方程数值模拟所需的（1）震源函数（2）地层速度（波速）（3）边界条件均定义完全。为求（1）式的数值解，必须将此式离散化（包括时间离散和空间离散，这里dt取0.002，dh取4，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化，得到横向上x=i*dh,纵向上z=j*dh,t=n*dt，k时刻(i,j)点的波场值为kjiu,。实验中波场本应定义三维数组，但考虑到内存问题，将其定义为二维数组，采用循环交换赋值实现三维数组的功能，这里利用泰勒展式，展到二阶略去高阶小量，整理得(i,j)点k时刻的二阶时间微商，对于空间微分，采用四阶精度差分格式分别在(i,j)点k时刻展开到四阶小量，消除四阶小量并解出二阶微分，由时间和空间的二阶微商整理得到波场值。}25][34][121{2,,1,1,2,22221,,1,kjikjikjikjikjikjikjikjiuuuuuhtvuuu图1})(*)(*)(00jjiits（2）一般来讲，差分时差分格式阶数越高，得到的波场值精度越高，但稳定性会受影响，差分格式的稳定性不仅与差分格式本身有关，也随着差分精度的不同而不同，而且与网格步长之比的大小有关，值得指出的是差分格式的稳定性与微分方程无关，本实验所采用的（2）式为条件稳定，稳定条件：83/*maxhtv，（2）式由差分方程近似替代微分方程所得，所以如果空间和时间采样间隔不当，就会产生频散现象，导致波形畸变，甚至派生出多个同相轴，为减少频散，Δh需满足Dablain的经验公式：)/(minNGfvh𝑓𝑁为Nyquist频率，一般取为主频的两倍，G为每个波长所占的网格点数，时间、空间为两阶差分的情况G取8，而时间、空间为四阶差分的情况G取4。则由上面两式得v∈[1000，3061.8]由上面定义的数值，用Matlab将震源函数)(ts成图，下面是成图的源程序：clf;f=20;r=4;k=100;dt=0.002;fori=1:kw(i)=exp(-(2*pi*f/r)^2*i^2*dt^2)*cos(2*pi*f*dt*i);endplot(w,'k')axis([0,100,-1,1]);title('震源离散后图像')xlabel('t\rightarrow'),ylabel('s(t)\rightarrow');由上面的程序可得震源离散后的图像，如下图：kjikjikjikjikjiuuuuuhtv,1,1,2,2,22225][34][121{图2小模型波前快照选用一个小模型，即200*200的范围，视为完全弹性模型，为一层介质，在介质中取速度为定值，本实验中取为2500，将震源放在模型中央，分别记录两个时刻的波前快照，（即区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射。在小模型中，可将介质看做是单一介质，即在整个介质中，波的传播速度是相同的，在这里设地震波速度为2500m/s，以下是Matlab成图的源程序：clfforx=1:200f(x)=2500;endplot(f);axis([1,200,1,3000]);title('小模型的速度-深度模型');figure;x=[00200200];y=[02002000];patch(x,y,'g');title('小模型');axisij;用matlab编译后得图像，见下图：根据网格设计可知，网格间距为10m，则200*200的矩阵实际上大致是个2000*2000平方米的区域，这样震源在（100，100）即中心点，即在地质剖面上，宽度为1000m，深度为1000m，假设介质均匀各向同性，且边界为自由界面，所以在介质中速度v以一常数2500m/s的传播速度传播，应该在400ms左右时到达边界。通过以上的计算，可以将k=130和k=240时波场的数据进行导出，用Matlab绘出图像。其中，在Matlab里用到的成像程序如下（这里只列举第一张波前快照的程序，第二张的与之相同）：clearall;loadwavefront1.dat;图3m=wavefront1(1);n=wavefront1(2);len=length(wavefront1)-2;nn=len/m;ifn~=nn;disp('Thedatafilelengtherror!');return;enddx=1;dz=1;fori=1:mforj=1:nx(j)=(i-1)*dx+wavefront1((i-1)*n+j+2)*12*dx;y(j)=(j-1)*dz;endplot(x,y,'k');holdonendaxisij;title('小模型的波前快照（t=260）');通过调节增益（x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*12*dx，dx前面的12），对波场值放大一定倍数，可以使能量小的增大，使成像清晰。从下图中可清楚看到波前的图像，数层的地震波随时间的推移由震源开始向外扩散。k=130时，大约是时间为260ms时的图像，此时波传播了s=v*t=650m，地震波尚未到达边界，而在k=240时，大约是时间为480ms时，此时波传播了s=v*t=1200m1000m,所以形成的图像是到边界之后的反弹波形图像，即人工边界反射，相当于只能看到反射的情况。地震子波的能量有限，当它由震源向周围传播时，波前面愈来愈大，就使愈来愈远地离开震源前进着的地震波的振幅愈来愈小。图4图5以下为小模型波前快照的源程序：/***********************************声波方程正演模拟——小模型波前快照************************************/#includestdio.h#includemath.h#includestdlib.h#definePI3.14159265#defineFM20#defineR4#defineKMAX1000#defineDT0.002//时间间隔#defineI200//横轴点数#defineJ200//纵轴点数#defineSI100//震源横坐标#defineSJ100//震源纵坐标#defineDH10//网格大小double**dimension2(intx,inty)//建立二维动态数组函数{inti;double**m;m=(double**)malloc(x*sizeof(double*));for(i=0;ix;i++)m[i]=(double*)malloc(y*sizeof(double));returnm;}voidexchange(double**u1,double**u2,double**u3)//循环交换赋值{inti,j;for(i=0;iI;i++)for(j=0;jJ;j++)u1[i][j]=u2[i][j];for(i=0;iI;i++)for(j=0;jJ;j++)u2[i][j]=u3[i][j];}main(){FILE*fp;//定义文件inti,j,k;floatv[I][J],Dt[I][J];double**u1,**u2,**u3,**u4,**u5,w[KMAX];doublet1,t2;u1=dimension2(I,J);u2=dimension2(I,J);u3=dimension2(I,J);u4=dimension2(I,J);u5=dimension2(I,J);for(i=0;iI;i++)for(j=0;jJ;j++){v[i][j]=2500;u1[i][j]=0;u2[i][j]=0;u3[i][j]=0;//对速度，u1，u2，u3赋初值}for(k=0;kKMAX;k++){w[k]=exp(-pow(2*PI*FM/R,2)*pow(k*DT,2))*cos(2*PI*FM*k*DT);}//波场值计算for(i=0;iI;i++)for(j=0;jJ;j++)Dt[i][j]=0;Dt[SI][SJ]=1;k=1;while(k300){k++;for(i=2;iI-2;i++)for(j=2;jJ-2;j++){t1=-1.0/12*(u2[i-2][j]+u2[i+2][j])+4.0/3*(u2[i-1][j]+u2[i+1][j])-5.0/2*u2[i][j];t2=-1.0/12*(u2[i][j-2]+u2[i][j+2])+4.0/3*(u2[i][j-1]+u2[i][j+1])-5.0/2*u2[i][j];u3[i][j]=2*u2[i][j]-u1[i][j]+pow(v[i][j],2)*pow(DT,2)/pow(DH,2)*(t1+t2)+Dt[i][j]*w[k];}exchange(u1,u2,u3);if(k==130)