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

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

2、的计算方法是做不离子方程的简化形式 Bortfeld 近似式（公式 2）。C 口ln(2)1Vp2工 cos 弓 siny222J=丿(=曲円1)他1讥2)2+v1v v2vp1HCOS2vp1ln(vp2)_ln(vp2vs1)vp1vp1vs2其中 Rpp 表示纵波反射系数。vp1 为入射层的纵波速度，vp2 为透射层的纵波速 度。91 为入射角，92 为透射角。p1 为入射层密度，p2 为折射层密度。Vs1 为入射层横波速度，Vs2为折射层的横波速度。在本文中由于在临界角以外只有 入射角，此时的反射系数可能会存在误差。 在第六层与第七层的反射系数过小造 成在最后的合成单炮记录上几乎肉眼不

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

4、为给定的动校速度，在该区域外运用的速度为插值得到的动校速度 进行常规动校，该动校的好坏有直接关系的是双曲线反射界面的追踪的准确程(3)度。图 2 2 常规动校正的流程图图 1 1 合成地震记录的流程图建立用于常规动校正需要知道的参数利用在单炮记录计算过加速度和时间to插值 样点的速度和时间【程中的计算的叠 宜得到相应的每层采进行双曲线动校正的编写，并 求出拉伸系数循环进行上述过程，直到得到 的点超出给定矩阵输出得到的新矩阵点赋值到:的相应的I扁点的值赋值厂图 4 4 追踪的射线路径该动校的好坏与中心线的确定有直接关系图 3 3 无拉伸动校的流程图H0imoism2KH3504ODZ2图 5 5

5、 时间与炮检距的变化图 4 中在起始位置的由大到小为层位由小到大。图 6 6 反射系数由于在临界角以后求得的反射系数为复数， 在图中取的模长，反射系数第一层界 面为深蓝线，第二层反射界面为红线，第三层反射系数为黄线，第四层的反射系 数为紫线，第五层的反射系数为绿线，第六层的反射系数为浅蓝线。图 9 9 均方根速度的动校正M40ST607DM图 7 7 反射系数在时间域arcjk?i40mi图 8 8 单炮记录0 100图 1212 基于图 9 9 的无拉伸动校正图 1010 均方根速度存在校正过量调整图 1111 切除拉伸系数 30%30%程序主程序部分clear all;% 清除内存clos

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

7、层垂直入射时间 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+t nm o(j)*(v(jF2); nnn=nnn+tnmo(j);endvnmo(i)=sqrt(mmm/nnn);end vnm o(1)=2040;%由于均方根速度把存在明显地校正过量对速度调整vnmo(2)=2265;vnmo(3)=2400;vnmo(4)=2760;vnm

8、o(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);%time 用于存放追踪出的双程旅行时； jl=zeros(2400,nd);% 用于存放形成的地震记录，记录长度 4000m，nd 道接收%产生雷克子波作为地震地波f

9、m=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).A2)*exp(-(pi*fm*t).A2); %雷克子波wl=length(seismic_wave);% 子波长度subplot(2,1,1);plot(1/time_sample*t,seismic_wave);title( 雷克子波 );xlabel( 时间(ms);ylabel( 幅值);a

10、xis(-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(2,1,2);plot(fi,2*abs(fftw/wl);%绘制频谱图title( 雷克子波

11、的振幅谱 );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+1ls=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:os+(nd-1)*dj;recey=zeros(1,nd);plot(recex,r

12、ecey,r A);hold on;% 画检波点%设置让每一层的入射角不会大于临界角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

13、,j-1)=a(1,j);figure(2) plot(sx,sy);hold on;%画第一层的射线路经end%追踪其余入射角和旅行时for j=2:nd+1for i=2:ncA=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-9c(1)=a(i,j); s=h(1)*tan(c(1);for k=2:i c(k)=asin(v(k)/v(k-1)*sin(

14、c(k-1); s=s+h(k)*tan(c(k);endif s-(os/2+(j-2)*dj/2)0B=a(i,j);else A=a(i,j);endC(i,j-1)=c(i);a(i,j)=(A+B)/2;endsx=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)=o

15、s+(j-2)*dj-sx(m);sy(m+1)=sy(m)-h(m);% 计算每一层射线射线路径的纵坐标 sy(2*i+2-m)=sy(m);endplot(sx,sy);hold on;% 画射线路径endendhold off;Re=zeros(nc,nd);% 反射系数的求取for i=1:nc;for j=1:ndRe(i,j)=reflact(vs(i),v(i),vs(i+1),v(i+1),C(i,j),mi(i),mi(i+1);endendRe=abs(Re);for i=1:nd;Re(4,i)=-Re(4,i);endDD=zeros(2400,101);for i=1

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

17、上去for k=nt-wl/2+1:nt+wl/2 jl(k,j)=jl(k,j)+sjseismic_wave(k-nt+wl/2);endendend%显示合成的单炮地震记录；也可以尝试用 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=

18、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(7)wigb(NNM);title( 无拉伸动校正 );for i=1:nc;w=0:(nd-1);% jj=10*w+1;% jj=

19、round(jj);FFF=time(i,w+1);显示figure(9);plot (40*w,FFF);title( 旅行时间 );hold on;endfor 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=reflact(vs1,vp1,vs2,vp2,thata,mi1,mi2)thata2=asin(vp2.*sin(thata)/vp1);s1=vp2.*mi

20、2.*cos(thata);s2=vp1.*mi1.*cos(thata2);a=1/2*log(s1/s2);s3=sin(thata)/vp1;s4=vs1.*s1-vs2.*s2;b=(s3)42.*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

21、, tnmo, (nt-1)*dt;v1 = vnmo(1), vnmo, vnmo(N);ti = (0:1:nt-1)*dt;vi = interp1(t1,v1,ti,linear);elseti = (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)A2 + (h(ih)/vi(it).A2 );time = sqrt(arg);stretch = (time-ti(it)/(ti(it)+1e-10);if stretchmax_s

22、tretch;%max_stretch/100; %M(it)= M(it) + 1;its = time/dt+1;it1 = floor(its);it2 = it1+1;a = its-it1; ifit2 1)t1 = 0, tnmo, (nt-1)*dt; v1 = vnmo(1), vnmo, vnmo(N); ti = (0:1:nt-1)*dt;vi = interp1(t1,v1,ti,linear); elseti = (0:1:nt-1)*dt;vi = ones(1,nt)*vnmo; end;dout1 = zeros(size(d);for j=1:N;for ih

23、 = 1:nh;c=m(j)-(ih)M.15/6;b=n (j)-(ih)A1.15/6; c=round(c); b=round(b); for it = c:b;% for ih = 1:nh;arg = ( ti(c).A2 + (h(ih)/vnmo(j).A2 );time = sqrt(arg);its = time/dt+1;it1 = floor(its);it2 = it1+1; a = its-it1; s=it-c;if it2=ntdout1(it,ih) = (1-a)*d(it1+s,ih)+a*d(it2+s,ih); k(it1+s,ih)=0;k(it2+s,ih)=0;end% endendend;end;for j=1:N-1;for ih = 1:nh;c=m(j+1)-(ih)A1.15/6;b=n(j)-(ih)A1.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