劉洋-研究生射線追蹤及動校正編程報告中國石油大學(xué)

1、射線追蹤及動校正編程報告原理報告中使用的地震子波為瑞克子波。射線追蹤使用的二分法的方式，確定了極大值和極小值后，利用中間值進行試射逼近，然后以中間值為一個極值再次進行取中間值進行試射，循環(huán)多次，當(dāng)射線得到的點與我們理想點之間的差距在給定的誤差范圍內(nèi)時，記錄此時射線經(jīng)過各層的坐標(biāo)和角度，對記錄的坐標(biāo)進行連接，得到射線的追蹤路徑。射線在各層之間的透射滿足斯奈爾定理（公式1）。Sin(a)/v1=sin(b)/v2 （1）其中，a是入射角，b是透射角，v1為入射層的層速度，v2為透射層的層速度。在反射系數(shù)的計算中沒有考慮每層的透射損失只考慮了最終層位的反射，反射系數(shù)使用的計算方法是做不離子方程的簡化

2、形式Bortfeld近似式（公式2）。 （2）其中Rpp表示縱波反射系數(shù)。vp1為入射層的縱波速度，vp2為透射層的縱波速度。1為入射角，2為透射角。1為入射層密度，2為折射層密度。Vs1為入射層橫波速度，Vs2為折射層的橫波速度。在本文中由于在臨界角以外只有入射角，此時的反射系數(shù)可能會存在誤差。在第六層與第七層的反射系數(shù)過小造成在最后的合成單炮記錄上幾乎肉眼不能看到。動校正的方法主要使用了常規(guī)的雙曲動校正和無拉伸動校正，常規(guī)雙曲動校正的原理是公式（3），在炮檢距較大時動校正會有明顯的動校拉伸，它的拉伸系數(shù)為公式（4）。 （3） 是自激自收的時間，x為炮檢距，v為動校正的速度。在實際中動校正的

3、速度應(yīng)該是變化的，但是為了計算的方便一般是設(shè)定成訂制，這就導(dǎo)致了動校正在大偏移距時發(fā)生動校誤差。 （4）t為反射點的時間與校正到點的時間的差，t為校正的點的時間。而無拉伸動校正的原理是在原始地震剖面上實現(xiàn)反射層位的部分實現(xiàn)整體搬離。本文中由于動校速度為定值，在大偏移距時存在比較明顯的誤差。在反射見面附近的動校速度為給定的動校速度，在該區(qū)域外運用的速度為插值得到的動校速度進行常規(guī)動校，該動校的好壞有直接關(guān)系的是雙曲線反射界面的追蹤的準(zhǔn)確程度。生成雷克子波,頻率40赫茲，采樣時間為1ms利用時間進行時移相應(yīng)的地震反射系數(shù)與子波的乘積，得到對應(yīng)的地震合成單炮記錄。計算得到反射系數(shù)利用Bortfeld

4、近似式編寫反射系數(shù)的程序計算自激自收時間，和均方根速度記錄各層傳播的時間，并將得到的點的坐標(biāo)連接得到射線路徑。在炮點利用二分法進行試射，層間透射滿足斯奈爾定理。定義觀測系統(tǒng)利用厚度構(gòu)造層位給定各層的初始數(shù)據(jù)包括縱橫波速度，密度等。開始流程圖圖1合成地震記錄的流程圖循環(huán)進行上述過程，直到得到的點超出給定矩陣t0點的值賦值為0將計算得到的相應(yīng)的點賦值到t0點是是否滿足拉伸條件進行雙曲線動校正的編寫，并求出拉伸系數(shù)利用在單炮記錄計算過程中的計算的疊加速度和時間t0插值得到相應(yīng)的每層采樣點的速度和時間建立用于常規(guī)動校正需要知道的參數(shù)否輸出得到的新矩陣 圖2 常規(guī)動校正的流程圖，將上面動校過得的點賦值0

5、，對上面沒有動校的區(qū)域進行常規(guī)動校，動校速度應(yīng)重新插值。以中心線處的動校速度為整片區(qū)域的動校速度，進行該區(qū)域的無拉伸動校。在可能出現(xiàn)同相軸的中心線地方進行一定幅度的拓寬，在子波的一半多一點的長度為宜利用常規(guī)動校正的圖在可能出現(xiàn)反射層位的地方 進行追蹤。將校正后的數(shù)據(jù)輸出圖3 無拉伸動校的流程圖該動校的好壞與中心線的確定有直接關(guān)系。圖4 追蹤的射線路徑圖5時間與炮檢距的變化圖4中在起始位置的由大到小為層位由小到大。圖6 反射系數(shù)由于在臨界角以后求得的反射系數(shù)為復(fù)數(shù)，在圖中取的模長，反射系數(shù)第一層界面為深藍線，第二層反射界面為紅線，第三層反射系數(shù)為黃線，第四層的反射系數(shù)為紫線，第五層的反射系數(shù)為綠

6、線，第六層的反射系數(shù)為淺藍線。圖7反射系數(shù)在時間域圖8 單炮記錄圖9 均方根速度的動校正圖10均方根速度存在校正過量調(diào)整圖11 切除拉伸系數(shù)30%圖12基于圖9的無拉伸動校正程序主程序部分clear all;%清除內(nèi)存close all;clc;%，關(guān)閉所有窗口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; %層數(shù)目nd=101; %接

7、收道數(shù)目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;%由于均方根速度把存在明顯地校正過量對速度調(diào)整vnmo(2)=2265;vnmo(3)=2400;vnmo(4)=2760;vnmo(5)=2976;vnmo(6)=3143;%程序所設(shè)中間變量及最終結(jié)果變量sx=zeros(1,105);sy=zeros(1,105);%sx和sy是畫射線追蹤圖的各點的橫縱坐標(biāo)a=zeros(nc,nd+1);c=zeros(nc,1);%a的第2列到第102列存放入射角，第1列存放0值，用于統(tǒng)一二分法射線追蹤的下限值；c用于存放對地下每一界面的入射角time=zeros(nc,nd);%

9、time用于存放追蹤出的雙程旅行時；jl=zeros(2400,nd);%用于存放形成的地震記錄，記錄長度4000m，nd道接收%產(chǎn)生雷克子波作為地震地波fm=40; %主頻L=256; %L為時間域采樣點數(shù)time_sample=0.001; %采樣間隔為1mst=(-L/2+1:L/2)*time_sample; %子波延續(xù)時間從-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離散，對應(yīng)頻譜周期性，周期為pfdf=1/(wl*time_sample); %譜線間隔df:1/NT,dw=2*pi*dfn=0:wl-1; %因為譜線除開零位置的只有N-1根，對稱性質(zhì)，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;%畫水平層狀介質(zhì)的各個界面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;%畫檢波點 %設(shè)置讓每一層的入射角不會大于臨界角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;%畫第一層的射線路經(jīng)end %追蹤其余入射角和旅行時for j=2:nd+1 for i=2:nc A=a(i,j-1);%設(shè)置二分法的下限 %設(shè)置二分法的上限 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;%設(shè)置炮點到入射點的橫向距離的初值； 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);%初始化一下，防止由于前面道的計算而遺留的坐標(biāo)值影響后面道的計算 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);%計算每一層射線射線路徑的橫坐標(biāo) sx(2*i+2-m)=os+(j-2)*dj-sx(m); sy(m+1)=sy(m)-h(m);%計算每一層射線射線路徑的縱坐標(biāo) sy(2*i+2-m)=sy(m); end plot(sx,sy);hold on;%畫射線路徑 endendhold off;Re=zeros(nc,nd);%反射系數(shù)的求取 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(反射系數(shù));%根據(jù)射線追蹤的結(jié)果合成單炮記錄for j=1:nd for i=1:nc sjfftw=ampw.*exp(1i.*angw);%振幅譜結(jié)合相位譜，重構(gòu)變換后的復(fù)信號 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(常規(guī)動校正);% % % 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(切除后的常規(guī)動校正);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(反射系數(shù)) hold on; end反射系數(shù)求取的程序：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