刘洋-研究生射线追踪及动校正编程报告中国石油大学.pdf

《刘洋-研究生射线追踪及动校正编程报告中国石油大学.pdf》由会员分享，可在线阅读，更多相关《刘洋-研究生射线追踪及动校正编程报告中国石油大学.pdf（18页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、 射线追踪及动校正编程报告 原理 报告中使用的地震子波为瑞克子波。射线追踪使用的二分法的方式，确定了极大值和极小值后，利用中间值进行试射逼近，然后以中间值为一个极值再次进行取中间值进行试射，循环多次，当射线得到的点与我们理想点之间的差距在给定的误差范围内时，记录此时射线经过各层的坐标和角度，对记录的坐标进行连接，得到射线的追踪路径。射线在各层之间的透射满足斯奈尔定理（公式1）。Sin(a)/v1=sin(b)/v2 （1）其中，a 是入射角，b 是透射角，v1 为入射层的层速度，v2 为透射层的层速度。在反射系数的计算中没有考虑每层的透射损失只考虑了最终层位的反射，反射系数使用的计算方法是做不

2、离子方程的简化形式 Bortfeld 近似式（公式 2）。222122211122211121112ln()cossin1ln()()()22cosln()ln()pPPssppsppppsvRvvvv vvvvv v （2）其中 Rpp 表示纵波反射系数。vp1 为入射层的纵波速度，vp2 为透射层的纵波速度。1 为入射角，2 为透射角。1 为入射层密度，2 为折射层密度。Vs1为入射层横波速度，Vs2 为折射层的横波速度。在本文中由于在临界角以外只有入射角，此时的反射系数可能会存在误差。在第六层与第七层的反射系数过小造成在最后的合成单炮记录上几乎肉眼不能看到。动校正的方法主要使用了常规的双

3、曲动校正和无拉伸动校正，常规双曲动校正的原理是公式（3），在炮检距较大时动校正会有明显的动校拉伸，它的拉伸系数为公式（4）。22220/ttxv （3）0t 是自激自收的时间，x 为炮检距，v 为动校正的速度。在实际中动校正的速度应该是变化的，但是为了计算的方便一般是设定成订制，这就导致了动校正在大偏移距时发生动校误差。/kt t （4）t 为反射点的时间与校正到点的时间的差，t 为校正的点的时间。而无拉伸动校正的原理是在原始地震剖面上实现反射层位的部分实现整体搬离。本文中由于动校速度为定值，在大偏移距时存在比较明显的误差。在反射见面附近的动校速度为给定的动校速度，在该区域外运用的速度为插值得

4、到的动校速度进行常规动校，该动校的好坏有直接关系的是双曲线反射界面的追踪的准确程度。流程图 利用厚度构造层位 开始 给定各层的初始数据包括纵横波速度，密度等。定义观测系统 在炮点利用二分法进行试射，层间透射满足斯奈尔定理。记录各层传播的时间，并将得到的点的坐标连接得到射线路径。利用Bortfeld 近似式编写反射系数的程序 计算自激自收时间，和均方根速度 计算得到反射系数 生成雷克子波,频率40 赫兹，采样时间为1ms 利用时间进行时移相应的地震反射系数与子波的乘积，得到对应的地震合成单炮记录。图 1 合成地震记录的流程图 图 2 常规动校正的流程图 建立用于常规动校正需要知道的参数 利用在单

5、炮记录计算过程中的计算的叠加速度和时间 t0 插值得到相应的每层采样点的速度和时间 进行双曲线动校正的编写，并求出拉伸系数 是否满足拉伸条件 是 否 将计算得到的相应的点赋值到 t0 点 t0 点的值赋值为 0 循环进行上述过程，直到得到的点超出给定矩阵 输出得到的新矩阵 图 3 无拉伸动校的流程图 该动校的好坏与中心线的确定有直接关系。图 4 追踪的射线路径 利用常规动校正的图在可能出现反射层位的地方 进行追踪。在可能出现同相轴的中心线地方进行一定幅度的拓宽，在子波的一半多一点的长度为宜 以中心线处的动校速度为整片区域的动校速度，进行该区域的无拉伸动校。，将上面动校过得的点赋值 0，对上面没

6、有动校的区域进行常规动校，动校速度应重新插值。将校正后的数据输出 图 5 时间与炮检距的变化 图 4 中在起始位置的由大到小为层位由小到大。图 6 反射系数 由于在临界角以后求得的反射系数为复数，在图中取的模长，反射系数第一层界面为深蓝线，第二层反射界面为红线，第三层反射系数为黄线，第四层的反射系数为紫线，第五层的反射系数为绿线，第六层的反射系数为浅蓝线。图 7 反射系数在时间域 图 8 单炮记录 图 9 均方根速度的动校正 图 10 均方根速度存在校正过量调整 图 11 切除拉伸系数 30%图 12 基于图 9 的无拉伸动校正 程序 主程序部分 clear all;%清除内存 close a

7、ll;clc;%，关闭所有窗口 v=2000;2500;2900;4000;3600;4000;4200;%纵波速度 vs=700;1050;1350;2000;1900;2200;2400;%横波速度 mi=2.0;2.1;2.2;2.4;2.3;2.5;2.6;%密度 h=500;400;200;400;300;200;%各层厚度 nc=6;%层数目 nd=101;%接收道数目 os=0;%最小炮检距；dj=40;%道间距 offset=zeros(1,nd);%偏移距的计算 for i=1:nd;offset(i)=i*dj;end tnmo=zeros(1,nc);%各层垂直入射时间

8、tnmo(1)=0.5;for i=2:nc;tnmo(i)=tnmo(i-1)+2*h(i)/v(i);end%vnm=zeros(1,nc);%均方根速度 vnmo=zeros(1,nc);%动校正速度 mmm=0;nnn=0;for i=1:nc;for j=1:i;mmm=mmm+tnmo(j)*(v(j)2);nnn=nnn+tnmo(j);end vnmo(i)=sqrt(mmm/nnn);end vnmo(1)=2040;%由于均方根速度把存在明显地校正过量对速度调整 vnmo(2)=2265;vnmo(3)=2400;vnmo(4)=2760;vnmo(5)=2976;vnmo

9、(6)=3143;%程序所设中间变量及最终结果变量 sx=zeros(1,105);sy=zeros(1,105);%sx 和 sy 是画射线追踪图的各点的横纵坐标 a=zeros(nc,nd+1);c=zeros(nc,1);%a 的第 2 列到第 102 列存放入射角，第 1 列存放 0 值，用于统一二分法射线追踪的下限值；c 用于存放对地下每一界面的入射角 time=zeros(nc,nd);%time 用于存放追踪出的双程旅行时；jl=zeros(2400,nd);%用于存放形成的地震记录，记录长度 4000m，nd 道接收%产生雷克子波作为地震地波 fm=40;%主频 L=256;%

10、L 为时间域采样点数 time_sample=0.001;%采样间隔为 1ms t=(-L/2+1:L/2)*time_sample;%子波延续时间从-L/2ms 到 L/2ms，共采集 L 个点 seismic_wave=(1-2*(pi*fm*t).2).*exp(-(pi*fm*t).2);%雷克子波 wl=length(seismic_wave);%子波长度 subplot(2,1,1);plot(1/time_sample*t,seismic_wave);title(雷克子波);xlabel(时间(ms);ylabel(幅值);axis(-80 80-2 10);%求雷克子波的频谱

11、fftw=fft(seismic_wave);ampw=abs(fftw);angw=angle(fftw);%求雷克子波的相位谱 pf=1/time_sample;%时域 dt 离散，对应频谱周期性，周期为 pf df=1/(wl*time_sample);%谱线间隔 df:1/NT,dw=2*pi*df n=0:wl-1;%因为谱线除开零位置的只有 N-1 根，对称性质，n=0 这个谱线是不考虑的 fi=df*n;subplot(2,1,2);plot(fi,2*abs(fftw/wl);%绘制频谱图 title(雷克子波的振幅谱);xlabel(频率(Hz);ylabel(幅值);axi

12、s(0 150 0 0.5);figure;xx=0;os+(nd-1)*dj;ls=0;hh=0;h(1:nc);for i=1:nc+1 ls=ls+hh(i);yy=-1*ls*ones(1,2);plot(xx,yy,k);hold on;%画水平层状介质的各个界面 end plot(0,0,r v);hold on;%画炮点 title(射线路径图);xlabel(炮检距(m);ylabel(深度(m);recex=os:dj:os+(nd-1)*dj;recey=zeros(1,nd);plot(recex,recey,r);hold on;%画检波点%设置让每一层的入射角不会大于

13、临界角 d=zeros(nc,1);%用于存放让地震波入射到第 i 层的初始入射角（第一层的入射角）的最大值 for i=2:nc d(i)=asin(v(1)/v(i);end C=zeros(nc,nd);%追踪第一层的入射角和旅行时 for j=2:nd+1 a(1,j)=atan(os/2+(j-2)*dj/2)/h(1);time(1,j-1)=2*(h(1)/cos(a(1,j)/v(1);sx(2)=(os+(j-2)*dj)/2;sx(3)=os+(j-2)*dj;sy(2)=-1*h(1);C(1,j-1)=a(1,j);figure(2)plot(sx,sy);hold o

14、n;%画第一层的射线路经 end%追踪其余入射角和旅行时 for j=2:nd+1 for i=2:nc A=a(i,j-1);%设置二分法的下限%设置二分法的上限 if a(i-1,j)d(i)B=d(i);else B=a(i-1,j);end a(i,j)=(A+B)/2;s=0;%设置炮点到入射点的横向距离的初值；while abs(s-(os/2+(j-2)*dj/2)1e-9 c(1)=a(i,j);s=h(1)*tan(c(1);for k=2:i c(k)=asin(v(k)/v(k-1)*sin(c(k-1);s=s+h(k)*tan(c(k);end if s-(os/2+

15、(j-2)*dj/2)0 B=a(i,j);else A=a(i,j);end C(i,j-1)=c(i);a(i,j)=(A+B)/2;end sx=zeros(1,11);sy=zeros(1,11);%初始化一下，防止由于前面道的计算而遗留的坐标值影响后面道的计算 for m=1:i temp=2*(h(m)/cos(c(m)/v(m);time(i,j-1)=time(i,j-1)+temp;sx(m+1)=sx(m)+h(m)*tan(c(m);%计算每一层射线射线路径的横坐标 sx(2*i+2-m)=os+(j-2)*dj-sx(m);sy(m+1)=sy(m)-h(m);%计算每

16、一层射线射线路径的纵坐标 sy(2*i+2-m)=sy(m);end plot(sx,sy);hold on;%画射线路径 end end hold off;Re=zeros(nc,nd);%反射系数的求取 for i=1:nc;for j=1:nd Re(i,j)=reflact(vs(i),v(i),vs(i+1),v(i+1),C(i,j),mi(i),mi(i+1);end end Re=abs(Re);for i=1:nd;Re(4,i)=-Re(4,i);end DD=zeros(2400,101);for i=1:nc;for j=1:nd;D=round(1000*time(i

17、,j);DD(D,j)=Re(i,j);end end figure(5)wigb(DD);title(反射系数);%根据射线追踪的结果合成单炮记录 for j=1:nd for i=1:nc sjfftw=ampw.*exp(1i.*angw);%振幅谱结合相位谱，重构变换后的复信号 sjseismic_wave=real(ifft(sjfftw);%反变换到时间域 sjseismic_wave=Re(i,j)*sjseismic_wave;nt=round(1000*time(i,j);%通过时移将震源子波放置到接收道上去 for k=nt-wl/2+1:nt+wl/2 jl(k,j)=j

18、l(k,j)+sjseismic_wave(k-nt+wl/2);end end end%显示合成的单炮地震记录；也可以尝试用 imagesc(jl)显示 figure(3)wigb(jl);title(合成的单炮记录);xlabel(道间距(40m);ylabel(时间(ms);%figure(6)%imagesc(jl);MMM=nmo(jl,time_sample,offset,tnmo,vnmo,10);figure(4)wigb(MMM);title(常规动校正);%MMM1=zeros(4000,71);%for i=1:4000;%for j=1:71;%MMM1(i,j)=MM

19、M(i,j);%end%end MMM1=nmo(jl,time_sample,offset,tnmo,vnmo,0.3);figure(6)wigb(MMM1);title(切除后的常规动校正);NNM=wula_nmo(jl,time_sample,offset,tnmo,vnmo);figure(7)wigb(NNM);title(无拉伸动校正);for i=1:nc;w=0:(nd-1);%jj=10*w+1;%jj=round(jj);FFF=time(i,w+1);figure(9);plot(40*w,FFF);title(旅行时间);hold on;end for i=1:nc

20、;w=0:(nd-1);%jj=10*w+1;%jj=round(jj);FF=Re(i,w+1);figure(10);plot(40*w,Re(i,w+1);title(反射系数)hold on;end 反射系数求取的程序：functionRF=reflact(vs1,vp1,vs2,vp2,thata,mi1,mi2)thata2=asin(vp2.*sin(thata)/vp1);s1=vp2.*mi2.*cos(thata);s2=vp1.*mi1.*cos(thata2);a=1/2*log(s1/s2);s3=sin(thata)/vp1;s4=vs1.*s1-vs2.*s2;b

21、=(s3).2.*s4;m=log(mi2/mi1);nn=log(vp2/vp1);yy=log(vp2.*vs1/(vp1.*vs2);RF=a+b.*(2+m/(nn-yy);return 动校正的程序：function dout,ti,vi=nmo(d,dt,h,tnmo,vnmo,max_stretch)nt,nh=size(d);N=length(vnmo);if(N1)t1=0,tnmo,(nt-1)*dt;v1=vnmo(1),vnmo,vnmo(N);ti=(0:1:nt-1)*dt;vi=interp1(t1,v1,ti,linear);else ti=(0:1:nt-1)

22、*dt;vi=ones(1,nt)*vnmo;end;dout=zeros(size(d);for it=1:nt;for ih=1:nh;arg=(ti(it)2+(h(ih)/vi(it).2);time=sqrt(arg);stretch=(time-ti(it)/(ti(it)+1e-10);if stretchmax_stretch;%max_stretch/100;%M(it)=M(it)+1;its=time/dt+1;it1=floor(its);it2=it1+1;a=its-it1;if it2 1)t1=0,tnmo,(nt-1)*dt;v1=vnmo(1),vnmo,v

23、nmo(N);ti=(0:1:nt-1)*dt;vi=interp1(t1,v1,ti,linear);else ti=(0:1:nt-1)*dt;vi=ones(1,nt)*vnmo;end;dout1=zeros(size(d);for j=1:N;for ih=1:nh;c=m(j)-(ih)1.15/6;b=n(j)-(ih)1.15/6;c=round(c);b=round(b);for it=c:b;%for ih=1:nh;arg=(ti(c).2+(h(ih)/vnmo(j).2);time=sqrt(arg);its=time/dt+1;it1=floor(its);it2=

24、it1+1;a=its-it1;s=it-c;if it2=nt dout1(it,ih)=(1-a)*d(it1+s,ih)+a*d(it2+s,ih);k(it1+s,ih)=0;k(it2+s,ih)=0;end%end end end;end;for j=1:N-1;for ih=1:nh;c=m(j+1)-(ih)1.15/6;b=n(j)-(ih)1.15/6;c=round(c);b=round(b);G=zeros(1,nt);G(b)=vi(j);G(c)=vi(j+1);G=interp1(G,ti,linear);for it=b:c;arg=(ti(it)2+(h(ih

25、)/G(it).2);time=sqrt(arg);its=time/dt+1;it1=floor(its);it2=it1+1;a=its-it1;if it2=nt;dout1(it,ih)=(1-a)*k(it1,ih)+a*k(it2,ih);end;end end end for ih=1:nh;c=m(1)-(ih)1.15/6;c=round(c);G=zeros(1,nt);%G(b)=vi(j);G(c)=vi(1);G=interp1(G,ti,linear);for it=1:c;arg=(ti(it)2+(h(ih)/G(it).2);time=sqrt(arg);it

26、s=time/dt+1;it1=floor(its);it2=it1+1;a=its-it1;if it2=nt;dout1(it,ih)=(1-a)*k(it1,ih)+a*k(it2,ih);end;end end for ih=1:nh;b=n(6)-(ih)1.15/6;b=round(b);G=zeros(1,nt);G(b)=vi(6);%G(c)=vi(j+1);G=interp1(G,ti,linear);for it=b:nt;arg=(ti(it)2+(h(ih)/G(it).2);time=sqrt(arg);its=time/dt+1;it1=floor(its);it2=it1+1;a=its-it1;if it2=nt;dout1(it,ih)=(1-a)*k(it1,ih)+a*k(it2,ih);end;end end return;