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

1、本科生实验报告实验课程 油气勘探新方法 学院名称 地球物理学院 专业名称 勘查技术与工程（石油物探） 学生姓名 学生学号 指导教师 熊高君 实验地点 5417 实验成绩 2015年12月 成都理工大学油气勘探新方法实验报告实验时间2015年 12月开课单位地球物理学院指导教师熊高君实验题目：叠前地震记录的相移波动方程正演模拟实验姓名学号班级专业勘测技术与工程（石油物探）院（系）地球物理学院 地球探测与信息技术系单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总成绩实验报告一、 实验题目叠前地震记录的相移波动方程正演模拟实验二、实验目的掌握各向同性介质任意构造、水平层状速度结构地质模型的

2、叠前地震记录相移波动方程正演模拟基本理论、实现方法与程序编制，由正演记录初步分析地震信号的分辨率。三、原理公式1、地震波传播的波动方程设（x,z）为空间坐标，t为时间，地震波传播速度为v(x,z)，则二维介质中任意位置、任意时刻的地震波场为p(z,x,t)：压缩波纵波。则二维各向同性均匀介质中地震波传播遵循的声波方程： 2p(x,z,t)x2+2p(x,z,t)z2=1v2(x,z)2p(x,z,t)t2 (1)2、傅里叶变换的微分性质p(t)与其傅里叶变换的P()的关系： P=-pte-itdt 正傅里叶变换 pt=12-Peitdt 逆傅里叶变换 2则有时间微分性质： (i)P=-dptd

3、te-itdt 一阶微分 (i)2P=-d2ptdt2e-itdt 二阶微分 3为频率，=2T，T为周期。同理有空间微分性质： (ik)Pk=-dpxdxe-ikxdx 一阶微分 (ik)2Pk=-d2pxdx2e-ikxdx 二阶微分 (4)k为波数， k=2， 为波长。3、地震波传播的相移外推公式令速度v不随x变化，只随z变化，则利用傅里叶变换微分性质（3）和（4）式，把波动方程(1)式变换到频率波数域，得： ik2Pk,zi,+2P(k,z,)z2=i2vz2Pk,z, （5）或： 2P(k,z,)z2=-2vz2-k2Pk,z, (6)令：kz2=2vz2-k2则(6)式的解为： Pk

4、,z,=c1e-ikzz+c2eikzz (7)包括上行波和下行波两项，正演模拟取上行波： Pk,z,=c1e-ikzz (8)若Zj和Zj+1间隔为z，速度v(z)为在此间隔内不随Z变的常数，(8)式实现波场从Zj+1到Zj的延拓，即： Pk,zj,=c1e-ikzz （9）在深度Zj+1开始向上延拓到Zj，若延拓深度为零，即：Z=Zj+1-Zj，则 Pk,zj=zj+1,=c1e-ikz(zj+1-zj)=ce-ikz0=c (10)对于任意深度Zj+1到Zj的延拓，可得正演模拟中地震波的传播方程（延拓公式） Pk,zj,=Pk,zj+1,c1e-ikz(zj+1-zj) (11)4、初始

5、条件和边界条件按照爆炸界面理论，反射界面震源在 t=0 时刻同时起爆，此时刻的波场就是震源。根据不同情况，可直接使用反射系数脉冲或子波作震源。如果直接使用反射系数作震源脉冲，则初始条件可表示为： p0x,z,t=rx,z t=0 0 t=其他 (12)p0x,z,t对时间t和空间x做二维傅立叶变换，则得频率-波数域的初始波场p0k,z,。边界条件： p0x,z,t=rx,z t=0,xminxxmax,zminzzmax 0 t=其他,x=其他,z=其他 (13) 其他参数都是在xminxxmax,zminzzmax范围内定义的。5、边界处理（1）边界反射问题把实际无穷空间区域中求解波场的问题

6、化为有穷区域求解时，左右两边使用零边界条件。物理上假设探区距Xmin与Xmax两个端点很远，在两个端点上收到的反射波很弱。但是，上述条件在实际中不能成立，造成零边界条件反而成为绝对阻止波通过的强反射面。在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。（2）边界强反射的处理镶边法、削波法、吸收边界都能有效消除边界强反射。削波法就是在波场延拓过程中，每延拓一次，在其两侧均匀衰减到零，从而消除边界强反射的影响。假设横向总长度为NX，以两边Lx道吸波为例，有以下吸波公式： AbsNx-Ix=AbsIx=sin(2LxLx-1) 0=IxLx AbsIx=1. Ix=其它 (14)6、定位原理某一

7、震源在某一固定检波点产生的记录， 就等于该固定检波点主动到所有可能反射点处接收的由这一震源激发的反射波的总和。形象地称为检波点下延记录原理。7、反射记录的形成 1）待激发的二次源的描述：介质的离散化把整个介质空间离散化，每个离散点都可以看成一个绕射点，某深度的反射界面与这一深度线的交点为这一反射界面的分布点。 每个离散点就是待激发的二次源 2）反射界面的分布某一深度层的反射界面分布在这一深度的离散点上，即反射界面与这一深度线的交点为这一反射界面的分布点。 3）反射界面二次源的形成 ：绕射波的产生震源从炸点开始向下传播，遍历介质的每一个离散点绕射点。每遇到一个绕射点，都会激发该点的绕射波。反射界

8、面与某深度离散点的 交叉点的反射系数不为零， 产生的绕射波也不为零， 而其他没有反射界面的地方，反 反射系数为零，是均匀介质， 产生的绕射波为零。 4）反射记录的形成把检波器的检波机制看成一个数学算子，或检波因子，可用脉冲或雷克子波。在延拓震源的同时，把检波器即检波因子也做完全相同的向下延拓，同样遍历每个离散点。震源和检波因子每延拓一个深度步长，震源、检波因子及反射系数三者做相关运算，则检波点就记录到了该深度每个绕射点的绕射波。所有二次源地震波都传播到测线某检波点，叠加而成该点的记录炮集中的某一记录道。测线排列多少道检波器，则都可记录道所在位置的地震模拟记录。8、数字化根据数字信号处理的采样定

9、理，把连续的信号变为计算机能处理的数字信号，使相移法正演模拟得以实现。频域抽样定理：一个频谱受限信号 f(t)，如果时间只占据 -tm-tm的范围，若在频域以不大于 1/2tm 频率间隔 f1/2tm 对信号f(t)的频谱F()采样，则抽样到的离散信号 F1 可以唯一表示原信号。时域抽样定理：一个时间受限信号 f(t)，如果频谱只占据-m-m的范围，则信号 f(t)可以用等间隔的抽样值唯一表示出来，而时间 t 抽样间隔必须不大于1/2fm，m=2fm，t=1/2fm。四、实验内容削波法相位移正演模拟（1）点绕射构造和水平层状速度模型（参数如图1所示）的正演数值模拟； 1）削波的正演；2）无削波

10、的正演；（2）计算中点和两个边界的信号位置，分析实验结果的正确性；（3）做同样模型的褶积模型数值模拟，对比分析两者的异同；（4）改变绕射点位置、速度，再做正演模拟；图1 点绕射构造与水平速度模型五、 方法路线1、参数初始化；2、二次源的形成：1）震源波场定义：Scr(x,z,t) Scrx,z,t=1 x=x0,z=z0,t=00 其它 （15）2）界面待激发爆炸源定义：Rsht(x,z,t) Rshtx,z,t=Rflctx,z t=00 其它 （16）3）形成削波数据：Absb( x )在所有深度都用同样的削波数据： AbsNx-Ix=AbsIx=sin(2LxLx-1) 0=IxLx A

11、bsIx=1. Ix=其它 （17）4）从测线深度开始向下做波场延拓 在所有深度都用同样的削波数据震源波场做削波处理 Scrx,zj,t=Scrx,zj,tAbsbx （18）震源做 x 和t 的傅里叶变换到频率-波数域： Scrx,zj,tFFTScrx,zj,t （19）计算相移延拓算子： eikzdzkx,zj,=e-i2vzj+12-kx2(zj+1-zj) （20）震源波场延拓计算 Scrkx,zj+1,= Scrkx,zj,eikzdzkx,zj, （21） 延拓后的震源波场变换到时间-空间域： Scrkx,zj+1,FFT-1 Scrx,zj,t （22）激发当前深度的二次源：震

12、源波场与待激发波场相关 Rscndx,zj+1,t=-Scrx,zj+1,Rshtx,zj+1,t+d （23）将二次源变换到频率域存储： Rscndx,zj+1,tFFTRscndx,zj+1, 24从目前深度开始，重到步，震源波场继续向下延拓,直至介质最深处，生成介质所有点的绕射波- 整个介质的二次源形成。3、正演记录产生的数学过程 1）定义频率域延拓波场:wfld(x,z,)；2）调入削波数据：Absb(x)；3）在频率域延拓波场： 每给一个频率j，延拓波场初始化：wfldx,z,j=0；从最深处开始，波场向上延拓： a、取出当前频率和深度的二次源：Rscnd(x,zj,j)；b、形成新

13、延拓波场： wfldx,zj,j=wfldx,zj,j+Rscndx,zj,j （25）c、新延拓波场做削波处理 wfldx,zj,j=wfldx,zj,jAbsbx （26）d、新延拓波场变换到波数域： wfldx,zj,jFFTwfldkx,zj,j （27）e、计算相移延拓算子： eikzdzkx,zj,j=e-ij2vzj2-kx2(zj-zj-1) （28）f、波场延拓计算 wfldkx,zj-1,j=wfldkx,zj,jeikzdzkx,zj,j （29）g、波场变换到空间域： wlfdkx,zj-1,jFFT-1wlfdx,zj-1,j （30）h、重复a到g的计算，直至波场延

14、拓到测线深度，ZJ=Z1，得到测线上的频率域波场：wlfdx,z1,j，存储测线上的波场:wlfdx,z1,j;重复到的计算，直到完成所有频率的循环；4)把波场变换到时间域，得一炮的叠前正演模拟结果： wlfdx,z1,FFT-1wlfdx,z1,t （31）4、频率-空间域波场对频率做反傅里叶变换，得时间-空间波场；5、使用Fimage软件显示所得结果。六、实验结果1、结果显示及分析（点绕射结构）图2 削波（左）正演模拟结果与未削波正演模拟结果（右）由图2分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉

15、正常界面的反射。削波后，边界假反射的影响消除了，且曲线变得更光滑了。a、速度V1=5000、V2=5500，深度h=2000时，改变绕射点x的位置；图3 x=Nx/4-1时，反射系数（左）与削波正演模拟结果（右）图4 x=Nx/2-1时，反射系数（左）与削波正演模拟结果（右）由图3与图4对比分析可知：当速度，深度一定时，改变绕射点的位置，曲线随绕射点位置的改变而左右移动，且移动趋势与绕射点位置的移动一致。b、绕射点x=Nx/2-1，深度h=2000 ，速度V2=5500时，改变速度V1；V1=4000V2=5500图5 V1=4000时，速度模型（左）与削波正演模拟结果（右）V1=6000V2

16、=5500图6 V1=6000时，速度模型（左）与削波正演模拟结果（右）由图5与图6对比分析可知：当绕射点位置固定，深度一定时，随着介质速度增大，曲线变得越平缓。这是由于速度增大时，波的传播时间减小所致。c、速度V1=5000、V2=5500，绕射点x=Nx/2-1，改变深度h；图7 h=1000时，反射系数（左）与削波正演模拟结果（右）图8 h=3000时，反射系数（左）与削波正演模拟结果（右）由图7与图8对比分析可知：当介质速度一定，绕射点位置不变时，随着绕射点深度的增加，曲线变得越来越平滑。显然是由于深度变大，波的传播时间变小所致。2、结果分析（水平层速度模型）图9 未削波（左）正演模拟

17、结果与削波正演模拟结果（右）由图9分析可知：削波前，由于在两个端点上收到的反射波很弱，造成零边界条件反而成为绝对阻止波通过的强反射面，在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。削波后，边界假反射的影响消除了，且曲线变得更光滑了。a、 深度h=2000 ，速度V2=5500时，改变速度V1；V1=4000V2=5500图10 V1=4000时，速度模型（左）与削波正演模拟结果（右）V2=5500V1=6000图11 V1=6000时，速度模型（左）与削波正演模拟结果（右）由图10与图11对比分析可知：水平层深度不变时，随着速度的增加，结果显示的水平层向上移动。这是由于速度增大时，波的

18、传播时间减小的。b、 速度V1=5000、V2=5500，改变深度h；图12 h=1000时，反射系数（左）与削波正演模拟结果（右）图13 h=3000时，反射系数（左）与削波正演模拟结果（右）由图12与图13对比分析可知：当波的传播速度不改变时，随着水平层深度的增加，结果显示的水平层向下移动。这是由于深度增加，波的传播时间增大。七、讨论建议1、实验收获通过本次实验，对叠后相移波动方程正演模拟的过程有了更深入的理解，以及对其公式和意义都有了进一步的认识。编程能力又有了一定提高。2、存在问题图14 x=3*Nx/4-1时，反射系数（左）与削波正演模拟结果（右）如图14，当速度，深度一定时，绕射点

19、位置x=3*Nx/4-1，其反射系数不能用Fimage正确显示，由显示的削波正演模拟结果可知，使用削波函数并未起到削波的作用。3、心得体会本次实验，并不容易。在实验中，必须要理解理论知识，才能够正确的编写程序。在实验过程中，遇到许多问题，通过和同学讨论，向老师求解，最终才得出实验结果。附程序代码：/Include #include #include #include #include #include /Parameters set #define Nx 128 /Trace Number #define Nt 256 /Record Number #define Nz 100 /Depth

20、Number #define Labs 15 /Length Of Boundary Absorbing #define Dx 20 /Trace Interval #define Dt 0.004 /Record Interval #define Dz 20 /Depth Interval #define Pai 3.1415926 /function:Judge the power of 2 int Po2Judge(int N) int k=0; long Ln=0; for(k=0;N-Ln0;k+) Ln=(long)pow(2,k); Ln=(long)pow(2,k-1); if

21、(fabs(Ln-N)=1)return(0); return(1); /function:Form Absorb Bounderyint Absorb( ) FILE *fp_Abs; int Ix; double AbsNx; if(fp_Abs=fopen(Absb.dat,wb)=NULL)printf(Connot open file Absb); for(Ix=0;IxNx;Ix+) AbsIx=1.; for(Ix=0;Ix=Labs-1;Ix+) AbsIx=sqrt(sin(Pai/2)*(Ix/(Labs-1);AbsNx-Ix-1=AbsIx; for(Ix=0;IxNx

22、;Ix+) fwrite(&AbsIx,sizeof(AbsIx),1,fp_Abs); fclose(fp_Abs); return(1); /function:Form Reflect Structure Model int Rflct( ) FILE *fp_Rflct; int Ix,Iz; float RflctNz; if(fp_Rflct=fopen(Rflct.dat,wb)=NULL)printf(Connot open file Reflection); for(Ix=0;IxNx;Ix+) for(Iz=0;IzNz;Iz+) RflctIz=0.; if(Iz=Nz/5

23、&Ix=Nx/2) RflctIz=1.; /if(Iz=3*Nz/5&fmod(Ix+1,16)=0.) RflctIz=5.; /if(Iz=4*Nz/5&fmod(Ix+1,32)=0.) RflctIz=7.; fwrite(&RflctIz,sizeof(RflctIz),1,fp_Rflct); fclose(fp_Rflct); return(1); /function:Rsht(x,z,t)int Rsht()int Ix,Iz,It;float RshtNt,RflctNz;FILE *fp_Rflct=fopen(Rflct.dat,rb);FILE *fp_Rsht=fo

24、pen(Rsht.dat,wb);for(Ix=0;IxNx;Ix+)for(Iz=0;IzNz;Iz+)fread(&RflctIz,sizeof(RflctIz),1,fp_Rflct);for(It=0;ItNt;It+)RshtIt=0.;Rsht0=RflctIz;for(It=0;ItNt;It+)fwrite(&RshtIt,sizeof(RshtIt),1,fp_Rsht);fclose(fp_Rflct);fclose(fp_Rsht);return(1); /function:Form Velocity Modelint Vlcty( ) FILE *fp_Vlcty; i

25、nt Ix,Iz; float VlctyNz; if(fp_Vlcty=fopen(Vlcty.dat,wb)=NULL) printf(Connot open file Vlcty); for(Ix=0;IxNx;Ix+) for(Iz=0;Iz(int)(3*Nz/4);Iz+) VlctyIz=5000.; fwrite(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty); for(Iz=(int)(3*Nz/4);Iz0;k+)Ln=(long)pow(2,k); k=k-1; for (it=0; it=n-1; it+)m = it; is = 0; for

26、(i=0; i=k-1; i+)j = m/2; is = 2*is+(m-2*j); m = j; frit = pris; fiit = piis; pr0 = 1.0; pi0 = 0.0; p = 6.283185306/(1.0*n); pr1 = (float) cos(p); pi1 = -(float)sin(p); if (l!=0)pi1=-pi1; for (i=2; i=n-1; i+)p = pri-1*pr1; q = pii-1*pi1; s = (pri-1+pii-1)*(pr1+pi1); pri = p-q; pii = s-p-q;for (it=0;

27、it=0; l0-)m = m/2; nv = 2*nv; for(it=0; it=(m-1)*nv; it=it+nv)for (j=0; j=(nv/2)-1; j+)p = prm*j*frit+j+nv/2; q = pim*j*fiit+j+nv/2; s = prm*j+pim*j; s = s*(frit+j+nv/2+fiit+j+nv/2); poddr = p-q; poddi = s-p-q; frit+j+nv/2 = frit+j-poddr; fiit+j+nv/2 = fiit+j-poddi; frit+j = frit+j+poddr; fiit+j = f

28、iit+j+poddi; if(l!=0)for(i=0; i=n-1; i+)fri = fri/(1.0*n); fii = fii/(1.0*n); if(il!=0)for(i=0; i=n-1; i+)pri = sqrt(fri*fri+fii*fii); if(fabs(fri)0)pii = 90.0;elsepii = -90.0;elsepii = atan(fii/fri)*360.0/6.283185306; for(i=0;in;i+)pri=fri; pii=fii; return ( 1 );/function:Form Inition Wave Field:Re

29、flectCoefficiency int Scr0()/Function05.1:kkfft(Wfld0r,Wfld0i,Nt,0) FILE *fp_Scr0r,*fp_Scr0i; int Ix,Iz,It; float Scr0rNt,Scr0iNt; if(fp_Scr0r=fopen(Scr0r.dat,wb)=NULL)printf(Connot open Wfld0r.dat); if(fp_Scr0i=fopen(Scr0i.dat,wb)=NULL)printf(Connot open Wfld0i.dat); /if(fp_Rflct=fopen(Rflct.dat, r

30、b)=NULL) printf(Connot open Rflct.dat); for(Ix=0;IxNx;Ix+) printf(Scr0_FFT:Ix=%dn,Ix); for(Iz=0;IzNz;Iz+) for(It=0;ItNt;It+) Scr0rIt=0; Scr0iIt=0; if(Ix=64&Iz=0)Scr0r0=1.; if(kkfft(Scr0r,Scr0i,Nt,0)!=1) printf(FFT is error); exit(0); for(It=0;ItNt/2+1;It+) fwrite(&Scr0rIt,sizeof(Scr0iIt),1,fp_Scr0r)

31、; fwrite(&Scr0iIt,sizeof(Scr0iIt),1,fp_Scr0i);fclose(fp_Scr0r); fclose(fp_Scr0i); return(1); /Function06.1.1:Read In Velocity Data and Absorb Boundary Data int ReadVlctyAbsb(float Vlcty,float Absb) FILE *fp_Vlcty,*fp_Absb; int Iz,Ix; if(fp_Vlcty=fopen(Vlcty.dat,rb)=NULL)printf(Connot open Rflct.dat

32、);for(Iz=0;IzNz;Iz+)fread(&VlctyIz,sizeof(VlctyIz),1,fp_Vlcty);fclose(fp_Vlcty); if(fp_Absb=fopen(Absb.dat,rb)=NULL)printf(Connot open Absb.dat); for(Ix=0;IxNx;Ix+)fread(&AbsbIx,sizeof(AbsbIx),1,fp_Absb);fclose(fp_Absb);/remove(Absb.dat);return(1);/06.1.2.1:Read Data From (Ix,Iz,Iw)Oder to (Iw,Iz,Ix

33、)Oderint ReadIxIzIwToIwIzIx(FILE *fp_Scr0r,FILE *fp_Scr0i,float Scr0r,float Scr0i,int Iz,int Iw)int Ix,Nw=Nt; long AddfrmStrt; for(Ix=0;IxNx;Ix+)/Data Numbers From start Position To Current Position AddfrmStrt=(Ix*Nz+Iz)*(Nt/2+1)+Iw; /Byte Numbers From start Position To Current Position fseek(fp_Scr

34、0r, sizeof(Scr0rIx)*AddfrmStrt,0); fread(&Scr0rIx,sizeof(Scr0rIx),1,fp_Scr0r); fseek(fp_Scr0i, sizeof(Scr0iIx)*AddfrmStrt,0); fread(&Scr0iIx,sizeof(Scr0iIx),1,fp_Scr0i);return(1);/function06.1.2:Form New Wave Field int FrmNewWfld(FILE *fp_Scr0r,FILE *fp_Scr0i,float Scrr,float Scri,float Absb,int Iz,

35、int Iw) /function06.1.2.1 ReadIxIzIwToIwIzIx(.) int ReadIxIzIwToIwIzIx(FILE *fp_Scr0r,FILE *fp_Scr0i,float Scr0r,float Scr0i,int Iz,int Iw); int Ix,Nw=Nt; float Scr0rNx,Scr0iNx; /Take out Initial Wave Field Data With Individual Depth if(ReadIxIzIwToIwIzIx(fp_Scr0r,fp_Scr0i,Scr0r,Scr0i,Iz,Iw)!=1)prin

36、tf(exp_ikzDz is error);exit(0); /Form New Wave Field for(Ix=0;Ix0) kz=sqrt(kz); eikzdz0=(float)cos(-kz*Dz); eikzdz1=(float)sin(-kz*Dz); return(1); /function06.1.3: Extrapolate One Depth Stepint MoveOneDz(float Scrr,float Scri,float Vz,float Dkx,float Dw,int Iw)/06.1.3.1:exp_ikzDz(kz,Ix,Vz,Iw,Dw,Dkx)

37、 /06.1.3.2:kkfft(Wfldr, Wfldi,Nx,j) int Ikx,Nkx=Nx; float kz2,Scr_r,Scr_i; if(kkfft(Scrr,Scri,Nx,0)!=1) printf(FFT is error);exit(0); for(Ikx=0;IkxNx/2+1;Ikx+) /4.2.3.1 Computing Phaseshift Function if(exp_ikzDz(kz,Ikx,Vz,Iw,Dw,Dkx)!=1)printf(exp_ikzDz is error);exit(0); /4.2.3.2 WaveField multiply

38、Phaseshift Function /Compute WaveField Phaseshift Scr_r=ScrrIkx*kz0-ScriIkx*kz1; Scr_i=ScriIkx*kz0+ScrrIkx*kz1; ScrrIkx=Scr_r; ScriIkx=Scr_i; if(Ikx!=0&Ikx!=Nkx/2) Scr_r=kz0*ScrrNkx-Ikx-kz1*ScriNkx-Ikx; Scr_i=kz1*ScrrNkx-Ikx+kz0*ScriNkx-Ikx; ScrrNkx-Ikx=Scr_r; ScriNkx-Ikx=Scr_i; /4.2.4 IFFT of New W

39、ave field From WaveNumber to Space if(kkfft(Scrr,Scri,Nkx,1)!=1)/37.FFT or Inverse FFT printf(FFT is error);exit(0); return(1);/function06.1:PhaseShift calculate int PhaseShift( ) int ReadVlctyAbsb(float *,float *); int FrmNewWfld(FILE *,FILE *,float *,float *,float *,int,int); int MoveOneDz (float *,float *,float,float,float,int); /1.2 Define Varibles FILE *fp_Scrr,*fp_Scri,*fp_Scr0r,*fp_Scr