有限差分法地震正演模拟程序一．二阶公式推导 1．二维的弹性波动方程d2U _X + 22U 2U 九+ »Q2U炉 _ x + x + ZSt 2 p dx 2 p Qz 2 p dxdzS 2U 九 + 2 »S 2U »S 2U 九+ »S 2Ur _ l + l + xSt2 p Sz2 p Sx2 p SxSz2．对方程进行中间差分首先对时间进行二阶差分U (t + t)-2U (t)+ U (t - t) Un+1 - 2Un + Un-1_ —x x x _ —xi?k xi,k xi,k1)S 2UxSt22)3)4)S2U _zSt2U (t +勺)-2U (t)+ U (t-分) Un+1 - 2Un + Un-1z , _ z — _ zi _ k对x方向进行二阶差分S2U U (x + x)-2U (x)+ U (x - x) Un △- 2Un + Unl _ ——x x x _ ——xi+1,k xi^k xi-1,kSx2 ( x b ( x)2S2Uz _Sx2U (x+°x)-2U (x)+ U (x -抵)q / Q 、 QUn - 2Un + Un_ ——zi+1,k zi,k zi-1,k3对z方向进行二阶差分U (z + z)-2U (z)+ U (z - z) Un - 2Un + Un_ ——x x x _ ——xi,k+1 x^k xi ,k-1S2UxSz2S2Ur _Sz2U (z +△ z)-2U (z)+ U (z-々)Un - 2Un + U —zi, k+1 zi,k ,nzi，k-1S 2U _ S(SUSSxSzSxzI Sz丿Sx对应进行差分△A(Un—zi, k+12—Un ) zi, k -1 z丿Unzi +1,k+1 z—U n — U n + U n:i-1,k+1 zi+l,k-1 zi-1,k -14 z xU n - U n U n - U nzi+1,k+1 z，——1,k+1 — z，+1,k——1 z，——1,k 112 x 2 x2 zd 2U d(Qu ]_ q(Un — Un )xi k +1 xi k —1QxQz QxI Qz )Qx(2 z丿Un — Unxi 亠i,k—1 xi—1,k 112 xUn —Unxi + 1 k +1 xi—1 k +1Un 一 Un 一 Un + Unxi+l,k+1 xi—l,k+1 xi+l,k-1 xi—l,k 一 14 z x(5)把(1) (2)(3) (4)得到的结果带入波动方程3 •写出差分方程Un+1 — 2Un + Un-1 X + 2L! Un — 2Un + Un斗 弋 = xi+1k 吟 xi—1,k(t J2 P ( x J2! U n — 2U n + U n X + ! U n — U n — U n + U n+ xi, k+1 x^k xi ,k-1 + zi+l,k+1 zi—1,k+1 zi+l,k-1 zi-1,k-1P a ( z J2 Pa 4 z xUn+1 —2Un +Un—1 X +2X Un —2Un +Unzirk ¥ zir^ = zi+1,k 旳 zi—1,k A △(t J2 P ( x J2! U n — 2U n + U n+ zi ,k+1 zi^k zi, k-1 +PX + ! Un —Un —Un +Unxi+1,k+1 xi—1,k+1 xi+1,k—1 xi—1,k —14 z xU n +1 =2U n -U n —1 +xi,k xi,k xi,kt2△”X + 2LX Un — 2Un + Un ——7_x^k xi—1,k + (Q pX Un—2Un + Un )―?—— xi, k —1(z J2X+X Un —Un —Un +Un+ zi+1,k+1 zi—1,k+1 zi+1,k—1 zi—1,k—1P 4 z x即得到X + 2 X U n — 2U n + U n X U n — 2U n + U nzi+1k 《 zi^U + zi^Ui z^ zi^U!U n+1 =2U n -U n—1+ t2 zi ,k zi ,k zi ,k(Q p (z J2X+! Un —Un —Un +Un+ xi+1,k+1 xi—1,k+1 xi+1,k—1 xi—1,k—14 z x二 MATLAB 程序 clear;clc;Nx=200;Nz=300;fi=30;%%% 主频 t_step=300;%%%%时间米样点 dx=10.0;%空间米样间隔——x方向 dz=10.0;%空间采样间隔——z方向 dt=0.001;%时间采样间隔——1msIambda=66*le9;%砂岩拉梅常数 lamda mu=44*le9;%砂岩拉梅常数mu rho=2650;%砂岩密度%%%%%%模型扩展%%%%%vp=zeros(Nz+2,Nx+2);% 纵波速度 vs=zeros(Nz+2,Nx+2);% 横波速度c=zeros(Nz+2,Nx+2);% 交叉项系数for i=l:Nz+2 for j=l:Nx+2vp(i,j)=sqrt((lambda+2*mu)/rho); vs(i,j)=sqrt((mu)/rho);c(i,j)=sqrt((lambda+mu)/rho);endend%%%%%% 震源 %%%%%%%t_wavelet=[l:t_step]*dt-l.0/fi;source=((1-2*(pi*fi*t_wavelet).人2).*exp(-(pi*fi*t_wavelet).人2));% 雷克子波 amp=sqrt(2.0);% 振幅source_x=floor(Nx/2)+1;% 震源位置 x 坐标source_z=2;%震源位置 z坐标source_amp=zeros(Nz+2,Nx+2);% 震源振幅初始化(所有点处均为 0) source_amp(source_z,source_x)=amp;% 震源振幅，在位置 source_z,source_x 处振 幅为amp,其它位置为0%%%%%%%%%%%%%%%%%%%%%%%%%%%U=zeros(Nz,Nx);% 弹性波 x 分量V=zeros(Nz,Nx);% 弹性波 z 分量U0=zeros(Nz+2,Nx+2);% 前一时刻的 UU1=zeros(Nz+2,Nx+2);% 当前时刻的 UU2=zeros(Nz+2,Nx+2);% 下一时刻的 UV0=zeros(Nz+2,Nx+2);% 前一时刻的 VV1=zeros(Nz+2,Nx+2);% 当前时刻的 VV2=zeros(Nz+2,Nx+2);% 下一时刻的 V record_u=zeros(t_step,Nx);% x方向接收到的地震记录 U record_v=zeros(t_step,Nx);% x方向接收到的地震记录 V%%%%%% 有限差分计算 U V %%%%%%for it=1:t_stepfor x=2:Nx+1for z=2:Nz+1U2(z,x)=2*Ul(z,x)-U0(z,x)+(dt*dt)*(vp(z,x)人2*(1.0/(dx*dx))*(Ul(z,x+l)-2*Ul(z,x) +Ul(z,x-l))+vs(z,x)人2*(1.0/(dz*dz))*(Ul(z+l,x)-2*Ul(z,x)+Ul(z-l,x))+c(z,x)人2*(1. 0/(4*dx*dz))*(V1(z+1,x+1)-V1(z+1,x-1)-V1(z-1,x+1)+V1(z-1,x-1)))+source(it)*sour ce_amp(z,x)*dt*dt;V2(z,x)=2*Vl(z,x)-V0(z,x)+(dt*dt)*(vs(z,x)人2*(l.0/(dx*dx))*(Vl(z,x+l)-2*Vl(z,x) +Vl(z,x-l))+vp(z,x)人2*(l.0/(dz*dz))*(Vl(z+l,x)-2*Vl(z,x)+Vl(z-l,x))+c(z,x)人2*(l. 0/(4*dx*dz))*(Ul(z+l,x+l)-Ul(z+l,x-l)-Ul(z-l,x+l)+Ul(z-l,x-l)));endendU0=Ul;Ul=U2;V0=Vl;Vl=V2;record_u(it,:)=U2(2,[2:20l]);record_v(it,:)=V2(2,[2:20l]);U=U2(2:30l,2:20l);V=V2(2:30l,2:20l);end%%%%% 绘制 U 和 V 的波场图%%%%%%figureimagesc(U(l:300,l:200));title('U');xlabel('x');ylabel('z');figureimagesc(V(l:300,l:200));title('V');xlabel('x');ylabel('z');%%%%% 绘制 U 和 V 的地震记录%%%%%%figureimagesc(record_u(l:300,l:200));title('U');xlabel('x');ylabel('t_step');figureimagesc(record_v(l:300,l:200));title('V');xlabel('x');ylabel('t_step');U 分量波场图V 分量波场图U 分量地震记录V 分量地震记录。