弹性波场论基础

double**er_wei(int x,int y)
{ double **m; m=(double **)malloc(x*sizeof(double*)); for(int i=0;i<x;i++) { m[i]=(double*)malloc(y*sizeof(double)); } return m; } 将 t=△t*300=0.3s 时刻的波场值导出，生成 wavefront.dat 文件，用 mat lab 画出此 时的波前。
//*************速度初赋值**** for(i=0;i<Xn;i++) for(j=0;j<Zn;j++) { if((i==0||i==Xn-1)||(j==0||j==Zn-1)) v[i][j]=0.0; //在边界波速为零 else v[i][j]=2800.0; //在均匀介质中波速为一常数 } //波场初赋值 ***** for(i=0;i<Xn;i++) for(j=0;j<Zn;j++) for(k=0;k<Tn;k++) { u[i][j][k]=0.0; } //波场计算********* for(k=1;k<Tn;k++) { w[k]=exp((-4*pi*pi*f*f*k*k*dt*dt)/(r*r))*cos(2*pi*f*k*dt);//扰动函数，为一 周期衰减函数，周期是 1/f for(i=1;i<Xn-1;i++) for(j=1;j<Zn-1;j++) { loc=0.0; if(i==150&&j==150) loc=1.0; //震源位于中心 u[i][j][k+1]=2*u[i][j][k]-u[i][j][k1]+(dt*dt/(dh*dh))*v[i][j]*v[i][j]*(u[i+1][j][k]-2*u[i][j][k]+u[i1][j][k]+u[i][j+1][k]-2*u[i][j][k]+u[i][j-1][k])+w[k]*loc; } }
double**array_2(int x,int y);//建立动态数组函数 double***array_3(int x,int y,int z);//建立动态数组函数 int main() { FILE* fp;//文件指针 int i,j,k; double loc,dt,dh,r,f; double ***u,**v,*w; u=array_3(Xn,Zn,Tn);
声波方程模拟
姓名___________________ 专业___________________ 年级___________________ 学号___________________ 任课教师________________
小模型
Fra Baidu bibliotek
大模型
程序
分析
拓展
总分
小模型 整个区域的速度值设为常数，即只有一种介质，将震源点放在模型中间，分别记录两 个时刻的波前快照（即该时刻区域内所有网格点的波场值）。第一时刻为地震波还未 传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，
(2)地震波传到边界的某时刻
时间采样间隔△t 为 0.001 秒， 空间间隔长度△h 为 5.0m 中心频率 f=25Hz 频带宽参数 r=4 波速 V=2000m/s 满足稳定性条件即 V△t/△h=0.4≤ （时间二阶，空间二阶）； 满足减少频散经验公式△h≤V/(G*fN )，其中 G=8(时间二阶，空间二阶)。
(1) 地震波未传播到边界的某时刻 时间采样间隔△t 为 0.001 秒， 空间间隔长度△h 为 5.0
中心频率 f=25 频带宽参数 r=4 波速 V=2000 满足稳定性条件即 V△t/△h=0.4≤ （时间二阶，空间二阶）； 满足减少频散经验公式△h≤V/(G*fN )，其中 G=8(时间二阶，空间二阶)。 此即为时间 t=△t*300=0.3 时刻的波前快照，震源在中心位置（150,150）激发， 传播 0.3 秒后的波前，由图可以看出，波前达到了（150,270）的位置，即地震波传播 了 112*△h=560 的距离，即 V*t=600。 震源函数为 e(−2π f/γ ) 后震源附近的介质的振 动情况。 在震源处，当炸药 在介质内部爆炸时，产 生强烈冲击波，是周围 的岩石破碎。随着远离 爆炸中心，作用力急剧 衰减，引起岩石塑性应 变，当进一步远离震源， 作用力继续衰减时，岩 石表现出弹性性质，发 生弹性应变。。从而形 成一瞬时膨胀点震源， 产生的波动向四周传播， 由于介质是各向同性的， 可以把爆炸中心周围的破碎带、塑性应变带与弹性应变带的分界面看成球面由于球对 称，质点只能发生径向位移，波动沿径向传播，波阵面位球面，所以可用纵波模拟地 震波。 区域大小为 300*300，网格点间隔为 5m，时间采样间隔为△t=0.001s，即在地质 剖面上，宽度为 1500m，深度为 1500m，假设介质均匀各向同性，所以在介质中速度 v 为一常数，设为 2000m/s，震源在中心（150,150）激发，地震波是球对称的，波阵
dx=1; dz=1; for i=1:m for j=1:n x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*15*dx; y(j)=(j-1)*dz; end plot(x,y,'k'); hold on end axis ij; title('波场值的波前快照');
通过调节增益（x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*15*dx，dx 前面的 15）， 对波场值放大一定倍数，可以使能量小的增大，使成像清晰。
//创建文件********** if((fp=fopen("wavefront.dat","w+"))!=NULL) { fprintf(fp,"%d\n",Xn); fprintf(fp,"%d\n",Zn); for(i=0;i<Xn;i++) for(j=0;j<Zn;j++) { fprintf(fp,"%f\n",u[i][j][200]); //地震波传播 0.2 秒是的波前快照 } fclose(fp); } return 0; }
clear all;
load wavefront.dat; m=wavefront(1); n=wavefront(2); len=length(wavefront)-2; nn=len/m; if n~=nn; disp('The data file length error!'); return; end
v=array_2(Xn,Zn); w = (double*)malloc(Tn*sizeof(double)); //边界条件 dt=0.001;dh=5.0;//时间离散化 ，时间采样间隔为 0.001s；空间网格化，网格 间隔为 5m loc=0; r=3;/*滤波主频*/ f=25.0;/*频带宽参数*/ //*************速度初赋值**** for(i=0;i<Xn;i++) for(j=0;j<Zn;j++) { if((i==0||i==Xn-1)||(j==0||j==Zn-1)) v[i][j]=0.0; else v[i][j]=2000.0; } //波场初赋值 ***** for(i=0;i<Xn;i++) for(j=0;j<Zn;j++) for(k=0;k<Tn;k++) { u[i][j][k]=0.0; } //波场值计算********* for(k=1;k<Tn;k++) { w[k]=exp((-2*pi*pi*f*f*k*k*dt*dt)/(r*r))*sin(2*pi*f*k*dt);//震源函数，为 一雷克子波函数，衰减因子为 2πf/r 周期为 1/f for(i=1;i<Xn-1;i++)
大模型 定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任 一点，记录下某些时刻的波前快照；二是合成一个地震记录，即记录下与震源同一深 度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。 一、 震源在区域表层
时间采样间隔△t 为 0.001 秒， 空间间隔长度△h 为 5.0m 中心频率 f=25Hz 频带宽参数 r=4 震源置于区域表层中心，区域纵向深度为 250*△h=1250m，横向宽度为 600*△ h=3000m.其中纵向分为三层介质，第一层介质深度为 100*△h=500m，波速为 V=2000m/s，第二层介质深度为 100---150,即 500---750m，波速为 V=3000m/s，第 三层介质深度为 150----250，即 750m----1250m，波速为 V=4000m/s，图中为 t=△ t*450=0.45s 时刻的波前快照。地震波传播到第二层介质时，发生反射和透射，但透 射波波速为V11 =3000m/s ，反射波波速为V12 =2000m/s 。透射波传播到第三层介质时， 继续发生二次反射和透射,透射波波速为V21 =4000m/s，反射波波速为V22 =3000m/s。 生成波场值的程序为： #include<stdio.h> #include<stdlib.h> #include<math.h> #define pi 3.1415926 #define Xn #define Zn #define Tn 600 //探测区域宽度 250 //探测区域深度 500 //记录时间长度
double**er_wei(int x,int y);//建立动态数组 double***san_wei(int x,int y,int z); int main() { FILE* fp; int i,j,k; double loc,dt,dh,r,f; double ***u,**v,*w; //动态数组 u=san_wei(Xn,Zn,Tn); v=er_wei(Xn,Zn); w = (double*)malloc(Tn*sizeof(double)); //边界条件 dt=0.001;dh=5.0; loc=0; r=2.50;/*滤波主频*/ f=25.0;/*频带宽参数*/
2 2
此即为时间 t=△t*490=0.49 时刻的波前快照，震源在中心位置（150,150）激发，传 播 0.49 秒后的波前，由图可以看出，波前达到了边界后又发生反射。 区域大小为 300*300，网格点间隔为 5m，时间采样间隔为△t=0.001s，即在地质剖面 上，宽度为 1500m，深度为 1500m，假设介质均匀各向同性，且边界为自由界面，所 以在介质中速度 v 为一常数，设为 2000m/s，震源在中心（150,150）激发，地震波是 球对称的，波阵面为球面，在空间二维上显示为深度，宽度的圆形波阵面。波传播到 边界时将发生反射，在时间 t=△t*375=0.375s 内，波传播了 s=v*t=750m,此时波传 播到边界，开始反射，到时间 t=△t*490=0.49s 时间内，反射波又传播了 s=v*t=230m,形成如图所示的波阵面。
2 t2
2 2
cos 2πft，模拟图如图所示，为一衰减周期函数，模拟爆炸
面为球面，在空间二维上显示为深度，宽度的圆形波阵面。在时间 t=△t*300=0.3s 内， 波传播了 s=v*t=600m,此时波前为一半径为 600m 的圆形波阵面，如图所示。
附：生成波场数据文件的程序： //波前快照 #include<stdio.h> #include<stdlib.h> #include<math.h> #define pi 3.1415926 #define Xn #define Zn #define Tn 300 300 400
double***san_wei(int x,int y,int z) { double ***m; m=(double ***)malloc(x*sizeof(double**)); for(int i=0;i<x;i++) { m[i]=(double**)malloc(y*sizeof(double*)); for(int j=0;j<y;j++) m[i][j] = (double*)malloc(z*sizeof(double)); } return m; }