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

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

2、形式Bortfeld近似式（公式2）。 （2）其中Rpp表示纵波反射系数。vp1为入射层的纵波速度，vp2为透射层的纵波速度。1为入射角，2为透射角。1为入射层密度，2为折射层密度。Vs1为入射层横波速度，Vs2为折射层的横波速度。在本文中由于在临界角以外只有入射角，此时的反射系数可能会存在误差。在第六层与第七层的反射系数过小造成在最后的合成单炮记录上几乎肉眼不能看到。动校正的方法主要使用了常规的双曲动校正和无拉伸动校正，常规双曲动校正的原理是公式（3），在炮检距较大时动校正会有明显的动校拉伸，它的拉伸系数为公式（4）。 （3） 是自激自收的时间，x为炮检距，v为动校正的速度。在实际中动校正的

3、速度应该是变化的，但是为了计算的方便一般是设定成订制，这就导致了动校正在大偏移距时发生动校误差。 （4）t为反射点的时间与校正到点的时间的差，t为校正的点的时间。而无拉伸动校正的原理是在原始地震剖面上实现反射层位的部分实现整体搬离。本文中由于动校速度为定值，在大偏移距时存在比较明显的误差。在反射见面附近的动校速度为给定的动校速度，在该区域外运用的速度为插值得到的动校速度进行常规动校，该动校的好坏有直接关系的是双曲线反射界面的追踪的准确程度。生成雷克子波,频率40赫兹，采样时间为1ms利用时间进行时移相应的地震反射系数与子波的乘积，得到对应的地震合成单炮记录。计算得到反射系数利用Bortfeld

4、近似式编写反射系数的程序计算自激自收时间，和均方根速度记录各层传播的时间，并将得到的点的坐标连接得到射线路径。在炮点利用二分法进行试射，层间透射满足斯奈尔定理。定义观测系统利用厚度构造层位给定各层的初始数据包括纵横波速度，密度等。开始流程图图1合成地震记录的流程图循环进行上述过程，直到得到的点超出给定矩阵t0点的值赋值为0将计算得到的相应的点赋值到t0点是是否满足拉伸条件进行双曲线动校正的编写，并求出拉伸系数利用在单炮记录计算过程中的计算的叠加速度和时间t0插值得到相应的每层采样点的速度和时间建立用于常规动校正需要知道的参数否输出得到的新矩阵 图2 常规动校正的流程图，将上面动校过得的点赋值0

5、，对上面没有动校的区域进行常规动校，动校速度应重新插值。以中心线处的动校速度为整片区域的动校速度，进行该区域的无拉伸动校。在可能出现同相轴的中心线地方进行一定幅度的拓宽，在子波的一半多一点的长度为宜利用常规动校正的图在可能出现反射层位的地方 进行追踪。将校正后的数据输出图3 无拉伸动校的流程图该动校的好坏与中心线的确定有直接关系。图4 追踪的射线路径图5时间与炮检距的变化图4中在起始位置的由大到小为层位由小到大。图6 反射系数由于在临界角以后求得的反射系数为复数，在图中取的模长，反射系数第一层界面为深蓝线，第二层反射界面为红线，第三层反射系数为黄线，第四层的反射系数为紫线，第五层的反射系数为绿

6、线，第六层的反射系数为浅蓝线。图7反射系数在时间域图8 单炮记录图9 均方根速度的动校正图10均方根速度存在校正过量调整图11 切除拉伸系数30%图12基于图9的无拉伸动校正程序主程序部分clear all;%清除内存close all;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; %接

7、收道数目os=0; %最小炮检距；dj=40; %道间距 offset=zeros(1,nd);%偏移距的计算for i=1:nd; offset(i)=i*dj;endtnmo=zeros(1,nc);%各层垂直入射时间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 vnm

8、o(i)=sqrt(mmm/nnn);endvnmo(1)=2040;%由于均方根速度把存在明显地校正过量对速度调整vnmo(2)=2265;vnmo(3)=2400;vnmo(4)=2760;vnmo(5)=2976;vnmo(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);%

9、time用于存放追踪出的双程旅行时；jl=zeros(2400,nd);%用于存放形成的地震记录，记录长度4000m，nd道接收%产生雷克子波作为地震地波fm=40; %主频L=256; %L为时间域采样点数time_sample=0.001; %采样间隔为1mst=(-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_s

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

11、2,1,2);plot(fi,2*abs(fftw/wl);%绘制频谱图title(雷克子波的振幅谱);xlabel(频率(Hz);ylabel(幅值);axis(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;%画水平层状介质的各个界面endplot(0,0,r v);hold on;%画炮点title(射线路径图);xlabel(炮检距(m);ylabel(深度(m);recex=os:dj

12、:os+(nd-1)*dj;recey=zeros(1,nd);plot(recex,recey,r );hold on;%画检波点 %设置让每一层的入射角不会大于临界角d=zeros(nc,1);%用于存放让地震波入射到第i层的初始入射角（第一层的入射角）的最大值for i=2:nc d(i)=asin(v(1)/v(i);endC=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)

13、/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 on;%画第一层的射线路经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(

14、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+(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; s

15、x(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);%计算每一层射线射线路径的纵坐标 sy(2*i+2-m)=sy(m); end plot(sx,sy);hold on;%画射线路径 endendhold 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=a

16、bs(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,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);%反变换到时间域 sjseis

17、mic_wave=Re(i,j)* sjseismic_wave; nt=round(1000*time(i,j); %通过时移将震源子波放置到接收道上去 for k=nt-wl/2+1:nt+wl/2 jl(k,j)=jl(k,j)+sjseismic_wave(k-nt+wl/2); end endend%显示合成的单炮地震记录；也可以尝试用imagesc(jl)显示figure(3)wigb(jl);title(合成的单炮记录);xlabel(道间距(40m);ylabel(时间(ms);% % % figure(6)% % % imagesc(jl);MMM=nmo(jl,time_s

18、ample,offset,tnmo,vnmo,10);figure(4)wigb(MMM);title(常规动校正);% % % MMM1=zeros(4000,71);% % % for i=1:4000;% % % for j=1:71;% % % MMM1(i,j)=MMM(i,j); % % % end% % % endMMM1=nmo(jl,time_sample,offset,tnmo,vnmo,0.3);figure(6)wigb(MMM1);title(切除后的常规动校正);NNM=wula_nmo(jl,time_sample,offset,tnmo,vnmo);figure

19、(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; 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=reflac

20、t(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=(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,v

21、i = 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)*dt; vi = ones(1,nt)*vnmo;end; dout = zeros(size(d);for it = 1:nt; for ih = 1:nh; arg = ( ti(it)

22、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, vnmo(N); ti = (0:1:nt-1)*dt; vi = interp1(t1,

23、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