叠加地震记录的相移波动方程正演模拟

本科生实验报告实验课程数值模型模拟学院名称地球物理学院专业名称勘测技术与工程（石油物探）学生姓名学生学号指导教师熊高君实验地点5417实验成绩2015年5月成都理工大学《地震数值模拟》实验报告实验时间2015年5月开课单位地球物理学院指导教师熊高君实验题目：叠加地震记录的相移波动方程正演模拟姓名学号班级专业勘测技术与工程（石油物探）院（系）地球物理学院地球探测与信息技术系单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总成绩实验报告一、实验题目叠加地震记录的相移波动方程正演模拟二、实验目的掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论、实现方法与程序编制，由正演记录初步分析地震信号的分辨率。三、原理公式1、地震波传播的波动方程设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二维介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波——纵波。则二维各向同性均匀介质中地震波传播遵循的声波方程：2、傅里叶变换的微分性质p(t)与其傅里叶变换的P()的关系：{∫正傅里叶变换∫逆傅里叶变换则有时间微分性质：{∫一阶微分∫二阶微分ω为频率，⁄，T为周期。同理有空间微分性质：{∫一阶微分∫二阶微分k为波数，k=⁄，为波长3、地震波传播的相移外推公式令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质（3）和（4）式，把波动方程(1)式变换到频率‐波数域，得：或：()令：则(5)式的解为：包括上行波和下行波两项，正演模拟取上行波：若和间隔为，速度v(z)为在此间隔内不随Z变的常数，(7)式实现波场从到的延拓，即：()在深度开始向上延拓到，若延拓深度为零，即：，则()对于任意深度到的延拓，可得正演模拟中地震波的传播方程（延拓公式()()4、初始条件和边界条件按照爆炸界面理论，反射界面震源在t=0时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为：{其他对时间t和空间x做二维傅立叶变换，则得频率-波数域的初始波场。边界条件：{其他其他其他其他参数都是在范围内定义的。5、边界处理（1）边界反射问题把实际无穷空间区域中求解波场的问题化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距与两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。削波法就是在波场延拓过程中，每延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式：{√其它6、数字化根据数字信号处理的采样定理，把连续的信号变为计算机能处理的数字信号，使相移法正演模拟得以实现。频域抽样定理：一个频谱受限信号，如果时间只占据的范围，若在频域以不大于频率间隔对信号的频谱采样，则抽样到的离散信号可以唯一表示原信号。时域抽样定理：一个时间受限信号，如果频谱只占据的范围，则信号可以用等间隔的抽样值唯一表示出来，而时间抽样间隔必须不大于，，。四、实验内容削波法相位移正演模拟（1）点绕射构造和水平层状速度模型（参数如图1所示）的正演数值模拟；1）削波的正演；2）无削波的正演；（2）计算中点和两个边界的信号位置，分析实验结果的正确性；（3）做同样模型的褶积模型数值模拟，对比分析两者的异同；（4）改变绕射点位置、速度，再做正演模拟；五、方法路线1、参数初始化；2、形成边界削波数据；3、波场初始化；3、Zmax层波场延拓到深度Zmax-1；5、Zi+1层波场延拓到深度Zi；6、重复5，从Iz=Nz-1开始，直到Iz=1，得测线上的频率—空间域波场；7、频率-空间域波场对频率做反傅里叶变换，得时间-空间波场；8、使用Fimage软件显示所得结果。六、实验结果1、结果显示及分析（点绕射结构）图1点绕射构造与水平速度模型图2削波（左）正演模拟结果与未削波正演模拟结果（右）由图2分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。削波后，边界假反射的影响消除了，且曲线变得更光滑了。a、速度V1=5000、V2=5500，深度h=2000时，改变绕射点x的位置；图3x=Nx/4-1时，反射系数（左）与削波正演模拟结果（右）图4x=Nx/2-1时，反射系数（左）与削波正演模拟结果（右）由图3与图4对比分析可知：当速度，深度一定时，改变绕射点的位置，曲线随绕射点位置的改变而左右移动，且移动趋势与绕射点位置的移动一致。b、绕射点x=Nx/2-1，深度h=2000，速度V2=5500时，改变速度V1；图5V1=4000时，速度模型（左）与削波正演模拟结果（右）V1=4000V2=5500图6V1=6000时，速度模型（左）与削波正演模拟结果（右）由图5与图6对比分析可知：当绕射点位置固定，深度一定时，随着介质速度增大，曲线变得越平缓。这是由于速度增大时，波的传播时间减小所致。c、速度V1=5000、V2=5500，绕射点x=Nx/2-1，改变深度h；图7h=1000时，反射系数（左）与削波正演模拟结果（右）V2=5500V1=6000图8h=3000时，反射系数（左）与削波正演模拟结果（右）由图7与图8对比分析可知：当介质速度一定，绕射点位置不变时，随着绕射点深度的增加，曲线变得越来越平滑。显然是由于深度变大，波的传播时间变小所致。2、结果分析（水平层速度模型）图9未削波（左）正演模拟结果与削波正演模拟结果（右）由图9分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。削波后，边界假反射的影响消除了，且曲线变得更光滑了。a、深度h=2000，速度V2=5500时，改变速度V1；图10V1=4000时，速度模型（左）与削波正演模拟结果（右）图11V1=6000时，速度模型（左）与削波正演模拟结果（右）由图10与图11对比分析可知：水平层深度不变时，随着速度的增加，结果显示的水平层向上移动。这是由于速度增大时，波的传播时间减小的。V1=4000V2=5500V2=5500V1=6000b、速度V1=5000、V2=5500，改变深度h；图12h=1000时，反射系数（左）与削波正演模拟结果（右）图13h=3000时，反射系数（左）与削波正演模拟结果（右）由图12与图13对比分析可知：当波的传播速度不改变时，随着水平层深度的增加，结果显示的水平层向下移动。这是由于深度增加，波的传播时间增大。七、讨论建议1、实验收获通过本次实验，对叠后相移波动方程正演模拟的过程有了更深入的理解，以及对其公式和意义都有了进一步的认识。编程能力又有了一定提高。2、存在问题图14x=3*Nx/4-1时，反射系数（左）与削波正演模拟结果（右）如图14，当速度，深度一定时，绕射点位置x=3*Nx/4-1，其反射系数不能用Fimage正确显示，由显示的削波正演模拟结果可知，使用削波函数并未起到削波的作用。3、心得体会本次实验，并不容易。在实验中，必须要理解理论知识，才能够正确的编写程序。在实验过程中，遇到许多问题，通过和同学讨论，向老师求解，最终才得出实验结果。在以后的学习工作中，当一个人的力量不够时，要积极和其他人求教，共同学习，共同进步。附程序代码：#includestdlib.h#includeconio.h#includemath.h#includestdio.h#includestring.h#defineNx128//TraceNumber#defineNt256//RecordNumber#defineNz100//DepthNumber#defineLabs10//LengthOfBoundaryAbsorbing#defineDx20//TraceInterval#defineDt0.004//RecordInterval#defineDz20//DepthInterval#definePai3.1415926//function01:Judgethepowerof2intPo2Judge(intN){intk=0;longLn=0;for(k=0;N-Ln0;k++)Ln=(long)pow(2,k);Ln=(long)pow(2,k-1);if(fabs(Ln-N)=1)return(0);return(1);}//function02:AbsorbBounderyintAbsorb(){FILE*fp_Abs;intIx;floatAbs[Nx];if((fp_Abs=fopen(Absb.dat,wb))==NULL)printf(ConnotopenfileAbsb);for(Ix=0;IxNx;Ix++)Abs[Ix]=1.0;//1.AbsorbBoundaryInitializing？for(Ix=0;IxLabs;Ix++){Abs[Ix]=sqrt(sin((Pai/2)*Ix/(Labs-1)));Abs[Nx-Ix-1]=Abs[Ix];//2.AbsorbBoundaryCompute？}for(Ix=0;IxNx;Ix++)fwrite(&Abs[Ix],sizeof(Abs[Ix]),Nx,fp_Abs);//3.ByteNumberoffclose(fp_Abs);return(1);}//function03:FormReflectStructureModelintRflct(){FILE*fp_Rflct;intIx,Iz;floatRflct[Nz];if((fp_Rflct=fopen(Rflct.dat,wb))==NULL)printf(ConnotopenfileReflection);for(Ix=0;IxNx;Ix++){for(Iz=0;IzNz;Iz++){Rflct[Iz]=0.;if(Ix==Nx/4-1&&Iz==Nz/2-1)Rflct[Iz]=1;fwrite(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct);//4.}}fclose(fp_Rflct);return(1);}//function04:FormVelocityModelintVlcty(){FILE*fp_Vlcty;intIx,Iz;floatVlcty[Nz];if((fp_Vlcty=fopen(Vlcty.dat,wb))==NULL)printf(ConnotopenfileVlcty);for(Ix=0;IxNx;Ix++){for(Iz=0;Iz(int)(3*Nz/4);Iz++){Vlcty[Iz]=5000.;fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);//5.}for(Iz=(int)(3*Nz/4);IzNz;Iz++){Vlcty[Iz]=5500.;fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);//6.}}fclose(fp_Vlcty);return(1);}//FourierTransform:FFT(i=0)orIFFT(l=1)intkkfft(floatpr[],floatpi[],intn,intl){intit,m,is,i,j,nv,l0,il=0;floatp,q,s,vr,vi,poddr,poddi;floatfr[4096],fi[4096];intk=0;longLn=0;for(k=0;n-Ln0;k++){Ln=(long)pow(2,k);}k=k-1;for