实验要求：
1、应用声波方程作为正演模拟的波动方程；
3、给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图形来； 4、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻，体会其人工边界反射；
5、对于大模型，定义为水平层状速度模型（至少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分别对应哪些波。
实验目的：
1. 通过本次作业，加深对波动方程的理解，明白波动方程所代表的物理意义。
2. 通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规律。
3. 通过模拟水平层状速度模型，体会地震波在两种介质分界面的传播规律，并能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。
4. 通过模拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波形变化趋势。
需要的已知条件包括：
1）震源函数
2）地层速度（波速）
3）边界条件
2．弹性波方程：
声波方程的有限差分法数值模拟
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
（4-1）
是介质在点（x , z）处的纵波速度，为描述速度位或者压力的波场，为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。
设x、z方向的网格间隔长度为，为时间采样步长，则有：
（i为正整数）
(j为正整数) (n为正整数)
表示在(i,j)点，k时刻的波场值。
将在(i,j)点k时刻用Taylor展式展开：
(4-2)
将在(i,j)点k时刻用Taylor展式展开：
(4-3)
将上两式相加，略去高阶小量，整理得(i,j)点k时刻的二阶时间微商为：
(4-4)
对于空间微分，采用四阶精度差分格式，（以X方向为例）即将、、
分别在(i,j)点k时刻展开到四阶小量，消除四阶小量并解出二阶微分得：
(4-5)
同理可得：
(4-6)
这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得：
}
(4-7)铜绿微囊藻
式中为介质速度的空间离散值，是空间离散步长，为时间离散步长，为震源函数，关于一般使用一个理论的雷克型子波代替，即：
（4-8）
上式中，为时间， 为中心频率，一般取为20-40HZ，为控制频带宽度的参数，一般取3-5。在实际计算过程中，需把此震源函数离散，参与波场计算。确定震源位置。
稳定性条件： (4-8)
这里表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-8)式时，第推求解（4-7）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。
频散关系式： (4-9)
式中为最小速度，为Nyquist频率。一般取震源子波中的主频的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分取8，而对于空间为四阶差分的情况则取4方能有效减少频散。
二．实验步骤
1、应用声波方程作为正演模拟的波动方程，忽略转换波的产生、传播；
2、将所提供震源函数离散后绘图；
震源函数为雷克子波，离散绘图如下：
fm=30;
r=3;
t=0.002;
for n=1:200
w(n)=exp(-(2*pi*fm/r)^2*(t*n)^2)*cos(2*pi*fm*t*n);
end
figure(1),plot(w);
3、对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模型中间，分别记录两个时刻的波前快照（即区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻。
由稳定条件，设v为2000，可以令DT=0.001,DH=5,此时精度较高，且满足频散关系
程序，图形如下：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.141593
#define FM 30
#define R 3
#define KN 200
#define XN 101
#define ZN 101
#define DH 5
#define DT 0.001
void main()
{
FILE *fp;
int i,j,k,m,n;
float u1[XN][ZN],u2[XN][ZN],u3[XN][ZN],u4[XN][ZN],f[XN][ZN]; //不能直接初值为0
float u5[XN][ZN],v[XN][ZN],w[KN],uu0,uu1,uu2;
for(k=0;k<KN;k++)
w[k]=exp(-(2*PI*FM/R)*(2*PI*FM/R)*(k*DT)*(k*DT))*cos(2*PI*FM*k*DT);
for(i=0;i<XN;i++) //定义f函数，当且仅当i,j同时为50时，f为1，其余为0
for(j=0;j<ZN;j++)
{if(i==50&&j==50)
f[i][j]=1;
else f[i][j]=0;
}
for(i=0;i<XN;i++)
仿古建筑施工
for(j=0;j<ZN;j++)
{ u1[i][j]=0.0;
u2[i][j]=0.0;
透水混凝土做法 u3[i][j]=0.0;
u4[i][j]=0.0;
v[i][j]=2000; //速度相同表示同一介质
}
for(k=0;k<KN;k++)
{for(i=2;i<XN-2;i++)
for(j=2;j<ZN-2;j++)
{uu0=(v[i][j])*(v[i][j])*(DT/DH)*(DT/DH);
uu1=-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];
uu2=-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]+uu0*uu1+uu0*uu2+w[k]*f[i][j];
}
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
{u1[m][n]=u2[m][n];
u2[m][n]=u3[m][n];
}
if(k==100)
for(m=0;m<XN;m++)
for(n=0;n<ZN;n++)
u4[m][n]=u3[m][n];//记录波前快照,中间点阿比丹 艾山
}
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",u4[i][j]);
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议