《利用有限差分方法合成地震记录》由会员分享，可在线阅读，更多相关《利用有限差分方法合成地震记录（5页珍藏版）》请在金锄头文库上搜索。
1、利用有限差分方法合成地震记录program seismicimplicit none real:lamda1,mu1,r1,dt,dx,dz,vp1,vs1,freq,lamda2,mu2,r2,vp2,vs2 real:c11,c12,c21,c22,c31,c32,d1,d2,d3,pireal,dimension(101):x,z real,dimension(101,101):uz1,uz2,uz3,u3,ux1,ux2,ux3 integer:i,j,k,n,ccharacter*10:text,nfile1,nfile2,tpwrite(*,*)input ytpe(mid-ref
2、lect,len-direct)read(*,(a) tpif(tp.eq.mid) thenprint *, dx,dzread(*,*) dx,dzelse if(tp.eq.len) thenprint *, dx,dzread(*,*) dx,dzend ifvp1=1500.0;vs1=1100.0dt=0.00025freq=10.0r1=300.0lamda1=(vp1*2-2*vs1*2)*r1mu1=vs1*2*r1pi=atan(1.0)*4.00c11=0.5*(lamda1+2.0*mu1)/r1c21=0.5*(lamda1+mu1)/r1 c31=mu1/(2.0*
3、r1) d1=(dt/dx)*2 d2=dt*2/(4.0*dx*dz) d3=(dt/dz)*2 n=500ux1=0.0;uz1=0.0ux2=0.0;uz2=0.0ux3=0.0;uz3=0.0do i=2,100do j=2,100 ux2(i,j)=ux1(i,j)+c11*d1*(ux1(i+1,j)+ux1(i-1,j)-2*ux1(i,j) & +c21*d2*(uz1(i+1,j+1)-uz1(i-1,j+1)-(uz1(i+1,j-1)-uz1(i-1,j-1) & +c31*d3*(ux1(i,j+1)+ux1(i,j-1)-2*ux1(i,j)uz2(i,j)=uz1(
4、i,j)+c11*d1*(uz1(i,j+1)+uz1(i,j-1)-2*uz1(i,j) & +c21*d2*(ux1(i+1,j+1)-ux1(i-1,j+1)-(ux1(i+1,j-1)-ux1(i-1,j-1) & +c31*d3*(uz1(i+1,j)+uz1(i-1,j)-2*uz1(i,j)end doend doux2(51,51)=0.0 uz2(51,51)=0.0 ux2(50,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(50,50)=-3.0*exp(-2.0*pi*freq*1.*d
5、t)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(52,50)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(52,50)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(50,52)=-3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(50,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)
6、 ux2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(52,52)=3.0*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) ux2(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(51,50)
7、=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0) uz2(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*1.*dt)*sin(2*pi*freq*1.*dt+pi/2.0)do k=2,ndo i=2,100do j=2,100ux3(i,j)=-ux1(i,j)+2*ux2(i,j)+2*c11*d1*(ux2(i+1,j)+ux2(i-1,j)-2*ux2(i,j) & +2*c21*d2*(uz2(i+1,j+1)-uz2(i-1,j+1)-(uz2(i+1,j-1)-uz
8、2(i-1,j-1)& +2*c31*d3*(ux2(i,j+1)+ux2(i,j-1)-2*ux2(i,j)uz3(i,j)=-uz1(i,j)+2*uz2(i,j)+2*c11*d1*(uz2(i,j+1)+uz2(i,j-1)-2*uz2(i,j) & +2*c21*d2*(ux2(i+1,j+1)-ux2(i-1,j+1)-(ux2(i+1,j-1)-ux2(i-1,j-1)& +2*c31*d3*(uz2(i+1,j)+uz2(i-1,j)-2*uz2(i,j) u3(i,j)=sqrt(ux3(i,j)*2+uz3(i,j)*2)end doend doux3(51,51)=0.0
9、 uz3(51,51)=0.0 ux3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(50,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(52,50)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(52,50)=-3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(50,52)=-3.0*exp(-
10、2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(50,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(52,52)=3.0*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) ux3(50,51)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*fr
11、eq*k*dt+pi/2.0) ux3(52,51)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)uz3(51,50)=-3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0) uz3(51,52)=3.0*sqrt(2.0)*exp(-2.0*pi*freq*k*dt)*sin(2*pi*freq*k*dt+pi/2.0)write(text,(i3) k c=mod(k,10) if(c .eq. 0) then nfile1=text(
12、1:3)/ux.dat open(88,file=nfile1) do i=1,101 do j=1,101 x(i)=(i-1)*dx z(j)=(j-1)*dz write(88,(1x,i3,f10.4,3x,f10.4,3x,f16.4,3x,f16.4,3x,f16.4)k,x(i),z(j),ux3(i,j),& uz3(i,j),u3(i,j) end doend do close(88) end ifdo i=2,100do j=2,100 ux1(i,j)=ux2(i,j) ux2(i,j)=ux3(i,j) uz1(i,j)=uz2(i,j) uz2(i,j)=uz3(i,
13、j) end do end doend dostopend program seismic实例采用网格数为100*100,方向和z方向的间隔dx=dz=l,时间间隔dt=0.00025s, 模拟地质体的纵波速度vp=1500.0m/s,横波速度为vs=1100.0m/s,震源的频率为freq=10.0 Hz,地质体的密度为r=300.0震源放在模型区域的中心位置，其振动为爆炸震源。通过程序计算结果如下：当T=10*dt时间后，由于处于爆炸的初期，波刚刚向前传播了 10*dt*vp的距离。如下图 一所示。当 T=100*dt 时，波则传播到 100*dt*vp=100*0.00025*1500=37.5m，则上界达到100-37.5=12.5。下界达到50+37.5=87.5m。左边界到达12.5m，右边界为87.5m如下图 所示。I 口口日口-呂口-*1 口-31DBQ图二、当T=100*dt时的波场快照当 T=200*dt 时，波往外传播 200*dt*vp=500*0.00025*1500=75m，则波发生反射, 如下图三所示。图三、当T=200*dt时的波场快照
