弹性波场论基础

声波方程模拟 姓名_ 专业_ 年级_ 学号_ 任课教师_ 小模型 大模型 程序 分析 拓展 总分 小模型 整个区域的速度值设为常数，即只有一种介质，将震源点放在模型中间，分别记录两 个时刻的波前快照（即该时刻区域内所有网格点的波场值）。第一时刻为地震波还未 传播到边界上的某时刻，第二时刻为地震波已经传播到边界上的某时刻， (1)地震波未传播到边界的某时刻 时间采样间隔t 为 0.001 秒， 空间间隔长度h 为 5.0 中心频率 f=25 频带宽参数 r=4 波速 V=2000 满足稳定性条件即 Vt/ h=0.4 22 （时间二阶，空间二阶）； 满足减少频散经验公式hV/(G* )，其中 G=8(时间二阶，空间二阶)。 此即为时间 t=t*300=0.3 时刻的波前快照，震源在中心位置（150,150）激发， 传播 0.3 秒后的波前，由图可以看出，波前达到了（150,270）的位置，即地震波传播 了 112*h=560 的距离，即 V*t=600。 震源函数为 ，模拟图如图所示，为一衰减周期函数，模拟爆(2/) 22cos2 炸后震源附近的介质的 振动情况。 在震源处，当炸药 在介质内部爆炸时，产 生强烈冲击波，是周围 的岩石破碎。随着远离 爆炸中心，作用力急剧 衰减，引起岩石塑性应 变，当进一步远离震源， 作用力继续衰减时，岩 石表现出弹性性质，发 生弹性应变。从而形 成一瞬时膨胀点震源， 产生的波动向四周传播， 由于介质是各向同性的， 可以把爆炸中心周围的破碎带、塑性应变带与弹性应变带的分界面看成球面由于球对 称，质点只能发生径向位移，波动沿径向传播，波阵面位球面，所以可用纵波模拟地 震波。 区域大小为 300*300，网格点间隔为 5m，时间采样间隔为 t=0.001s ，即在地质 剖面上，宽度为 1500m，深度为 1500m，假设介质均匀各向同性，所以在介质中速度 v 为一常数，设为 2000m/s，震源在中心（150,150）激发，地震波是球对称的，波阵 面为球面，在空间二维上显示为深度，宽度的圆形波阵面。在时间 t=t*300=0.3s 内， 波传播了 s=v*t=600m,此时波前为一半径为 600m 的圆形波阵面，如图所示。 附：生成波场数据文件的程序： /波前快照 #include #include #include #define pi 3.1415926 #define Xn 300 #define Zn 300 #define Tn 400 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;/*频带宽参数*/ /*速度初赋值* for(i=0;i #include #include #define pi 3.1415926 #define Xn 600 /探测区域宽度 #define Zn 250 /探测区域深度 #define Tn 500 /记录时间长度 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); 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=100 else /第三层介质，深度为 150-250，波速为 v=4000 uijk+1=2*uijk-uijk- 1+(dt*dt/(dh*dh)*(2*vij)*(vij)*(ui+1jk-2*uijk+ui- 1jk+uij+1k-2*uijk+uij-1k)+wk*loc; /创建数据文件，导出 t=400 时刻的波场值* if(fp=fopen(“wavefront3.dat“,“w+“)!=NULL) fprintf(fp,“%dn“,Xn); fprintf(fp,“%dn“,Zn); for(i=0;iXn;i+) for(j=0;jZn;j+) fprintf(fp,“%fn“,uij400); fclose(fp); return 0; 其地质模型如图所示： 地震波的透射与反射与两层介质的性质有关，当地震波从波速较快的介质入射到 波速较慢的介质中时，反射波的能量要大于透射波的能量，反之，则反射波的能量要 小于透射波的能量。此模型为地震波从波速较慢的介质入射到波速较快的介质，透射 波的能量要大于反射波的能量。 地震波从波速较快的介质入射到波速较慢的介质 此模型为地震波从波速较快的介质入射到波速较快的介质，由图中可以看出，反 射波的能量要大于透射波的能量。 地震波 V 反射波 =2000m/s12 透射波 =3000m/s11 反射波 =3000m/s22 透射波 21 二、 合成地震记录 无边界吸收 此为一模拟合成地震记录，震源在深度为 70*h=350m 激发，检波器与震源在 同一深度，检波器接收到的不同深度的地震记录包括直达波，反射波，折射波等。 此地质模型为三层介质，第一层介质深度为 100*h=500m，第二层介质深度为 100-150，即 500m-750m，第三层介质深度为 150-200，即 750m-1000m。 地质模型为： 第一层介质 深度为 100，波速 V=2000m/s 第二层介质 深度为 100-150，波速 V=3000m/s 第三层介质 深度为 150200，波速 V=4000m/s = 其原理如图所示： 同相轴为 直线和双曲线，其中直线为直达波，双曲线为反射波，反射波 的同一相轴可带来一个地层信息。在此地震记录中： 呈直线为直达波，其时距方程为：x=vt,其中 v 为在介质中的波速， v=2000m/s， x 为检波器距震源的距离，t 为地震波传播的时间。 第一条双曲线为在第一层介质与第二层介质分界面上的反射波，其时距方程为： ,其中 v 为第一层介质中的波速， v=2000m/s，x 为检波器距震源的距离， 22424=2 h 为震源距第一层介质的深度，t 为地震波传播的时间。 第二条双曲线为地震波传播到自由界面时的反射波，由于没有透射波，即没有能 量损失，所以在地震记录中比较清晰。其时距方程为： ，其中 v 为第一 22424=2 层介质中的波速，x 为检波器距震源的距离，h 为震源的深度， t 为地震波传播的时间。 第三条双曲线为地震波透过第一层介质，透射波传播到第二层介质与第三层介质 分界面时的反射波， 其时距方程为： 21214214=21 t x 22224224=22 其中 为第一层介质中的波速， ， 为第二层介质中的波速，1 1=2000/2 ， 为震源距第一层介质的深度， 为第二层介质的深度， 为地震波2=3000/1 2 1 沿 x 轴方向传播的距离， 为地震波沿 x 轴方向传播的距离， 为地震波在第一层介质2 1 中传播的时间， 为地震波在第二层介质中传播的时间。2 第四条双曲线为二次反射波，即地震波在第一层介质与第二层介质分界面反射后 又在自由界面二次反射。 其时距方程为： 2214214=21 2224224=22 其中 v 为第一层介质中的波速，v=2000m/s ， 为震源据第一层介质的深度， 为震1 2 源深度，x 为检波器距震源的距离，t 为地震波传播的时间。 由于地质模型的有限性，当地震波传播到地质边界时会产生因边界内外速度的巨 大反差( 边界外的速度为零)，从而产生了不是由日标体界面产生的反射，这就是虚假 反射，即边界反射。 由此可见，在地震勘探中，边界反射的存在会形成干扰，影响成像效果，所以在 模拟时应去掉边界反射，即有边界吸收，生成无边界反射的地震记录。 边界吸收的效果： =5.0; =0.001; 此模型为 200 200，网格点间隔 5.0m，时间采样间隔 0.001s，震源在中心激发， 波速为 v=2000m/s，此为 t= *320=0.32s 时的波前快照。显然，在时间 t= 0.25s 时， 地震波已传播到边界，但是由于边界吸收，并没有反射波，边界吸收效果比较好，是 模拟更加接近真实地震记录。 边界吸收的方程为： P(0, j, k + 1)=(2 一 2A 一 )P(0, j，k)+2A(1+A)P(1, j, k)2 一 A2P(2, j, k) + (2A 一 1)P(0, j, k 一 1)一 2AP(1, j, k 一 1); P(M , j, k + 1)=(2 一 2A 一 A2)P(M, j, k)+2A(1+A)P(M 一 1, j, k) 一 A2P(M 一 2, j, k)+(2A 一 1)P(M, j, k 一 1)一 2AP(M 一 1, j, k 一 1); P(i, 0, k+1)=(2 一 2A 一 A2 )P(i, 0, k)+2A(1+A)P(i, 1, k) 一 A2P(i, 2, k)+(2A 一 1)P(i, 0, k 一 1)一 2AP(i, 1, k 一 1); P(i, N, k+1)=(2 一 2A 一 A2 )P( i, N, k)+2A(1+A)P(i, N 一 1, k) 一 A2P(i, N 一 2, k )+(2A 一 1)P(i, N, k 一 1)一 2AP(i, N 一 1, k 一 1); 其中 M,N 分别为区域的边界，A=V*dt/dh ，P( i , j , k )为 k 时刻的波场值。 绕射波（无人工边界反射）： 当地震波传播时，遇到地层中的伏潜山等地下突起物时，将在地下突起物为震源 处发生绕射，使之成为一个“震源”，反映在波前快照上，将有两个波阵面如图所示： 在各向同性均匀介质中的绕射波： =5.0; =0.001; 此区域为 300*300，震源在中心（150,150）处激发，波速为 v=2000m/s，地下 突起物在（150,100）处，大小为 5*5 一块区域。当地震波传播过程中遇到地下突起物 时，将形成第二个震源，只是能量要比中心震源小很多，此图为 t=t*350=0.35s 时 的波前快照，此时地震波传播到（150,10）点处，在 t=t*125=0.125s 时，地震波 传播到地下突起物处，在此形成震源，到 t=0.35s 时刻，地震波又传播了 0.225s，到 达（150,190）点处，从而显示在波前快照中为两个相切的圆，从中可以看出，绕射波 的能量要比地震波的能量小很多。 在分层介质中的绕射波： 当地震波传播时，遇到地层中的古潜山等地下突起物时，将在地下突起物为震源 处发生绕射，使之成为一个“震源”，反映在地震记录上，将有两个震源如图所示： 此区域为 300*250，震源在（5,60）处激发，介质分为三层，第一层深度为 120，波速为 v=2400m/s，第二层介质的深度为 120-160，波速为 1200m/s，第三层 介质的深度为 160-250，波速为 2400m/s，地下突起物在点（150,100）处，大小为 5*5 一块区域，在第一层介质中。当地震波传播过程中遇到地下突起物时，此处将成 为第二个震源，产生的绕射波被检波器接受，但绕射波有可能和发射波混迭，形成的 地震记录如图所示，从中可以看出，绕射波和反射波很接近，近乎重合。地震波在多 层介质中经过反射和透射、绕射后，能量已经衰减。 其地质模型为： 第一层介质 深度为 120，波速 V=2400m/s 第二层介质 深度为 120-160，波速 V=1200m/s 第三层介质 深度为 160250，波速 V=2400m/s = 地下突起物 同相轴为直线和双曲线，其中直线为直达波，双曲线为反射波，反射波的同一相 轴可带来一个地层信息。在此地震记录中： 呈直线为直达波，其时距方程为：x=vt,其中 v 为在介质中的波速， v=2400m/s， x 为检波器距震源的距离，t 为地震波传播的时间。 第一条双曲线为在第一层介质与第二层介质分界面上的反射波，其时距方程为： ,其中 v 为第一层介质中的波速， v=2000m/s，x 为检波器距震源的距离， 22424=2 h 为震源距第一层介质的深度，t 为地震波传播的时间。 第二条双曲线为地震波透过第一层介质，透射波传播到第二层介质与第三层介质 分界面时的反射波， 其时距方程为： 21214214=21 22224224=22 其中 为第一层介质中的波速， ， 为第二层介质中的波速，1 1=2400/2 ， 为震源距第一层介质的深度， 为第二层介质的深度， 为地震波2=1200/1 2 1 沿 x 轴方向传播的距离， 为地震波沿 x 轴方向传播的距离， 为地震波在第一层介质2 1 中传播的时间， 为地震波在第二层介质中传播的时间。2 在中间的为绕射波，由于地下突起物在（150,100 ）处，位于中心，所以绕射波将 成对称出现在地震记录中，其时距方程为： 2242=2 其中 v 为第一层介质中的波速，v=2400m/s ，x 为检波器距地下突起物的距离，h 为地 下突起物与检波器的垂直距离，h=40。 折射波： 当地震波以非零任意角度入射到界面时，情况比较复杂，这是有转换波产生，各 二次波的能量分配关系完全由做普利兹方程决定，也是 AVO 分析的理论基础。 当波由稀疏介质向密介质入射时，入射角 增大到一定角度，透射波的能量迅速 下降，表明能量发生再分配，会产生一种新的波动-折射波，该角度称为临界角， 在临界角附近反射能量很强，这叫广角反射。 其地震记录如图所示： 模拟原理 地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的方法模拟 该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野外数据采集过 程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值模拟技术。地震波 场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、震源和地下几何界面、 物性参数（岩层密度、速度等）情况下，研究弹性波或声波的传播规律。 在室内模拟时，常采用声波方程，弹性波方程模拟地震波。 声波方程： S(t)zpxvtp(222 能够描述且只能描述纵波的传播规律，包括直达波、反射波、透射波、折射波等， 但不能描述转换波传播规律。 应用声波方程需已知的条件为： 1）震源函数 2）地层速度 3）边界条件 弹性波方程： 能够描述纵、横波的传播规律，包括直达波、反射波、透射波、折射波以及转换 波等。其中 为纵波波速， 为横波波速。 为拉梅常数。= +2 = ， 进一步比较 和 ,求的纵波与横波速度比： =+2 =2(1)12 由于一般的材料的泊松比 1/4 1/3,故 。地下岩石的泊松比一般 32 为 1/4， ，由此可见，地层中纵波的速度要比横波的波速大，故在地震勘探=3 中纵波比横波先到达。地震勘探长期使用纵波，但目前横波日引起人们的注意，特别 是转换横波。横波具有较高的分辨率，能够对尖灭、小幅度构造、小断层、礁体、古 潜山等准确定位。对非构造油气藏勘探、真假亮点的识别、气筒内部成像、裂缝发育 分析、流体识别等十分有效，因为横波不受流体影响。 需要的已知条件包括： 1) 震源函数 2）地层速度或根据方程的类型需要提供的地层的其它弹性参数 3）边界条件 )()(222222zwxvtwtSutusp 在地震勘探中，地震波传播的实际介质是十分复杂的。在一定条件下，即震源作 用时间短，作用力微小，地球介质可以看作完全弹性模型，但随着地震勘探技术的发 展，勘探精度要求提高，面临复杂地质目标时，要求地震勘探采用更加符合实际的介 质模型进行研究。粘弹性介质模型更符合实际。 但是到目前为止，在地震资料反演处理中应用最多的还是声波方程，弹性波以及 粘弹性波方程的应用还只是停留在模拟层次上。 地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元 法和有限差分法等。 相对于上述几种方法，有限差分法是一种更为快速有效的方法。虽然其精度比不 上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的 重视与