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

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

2、的計算方法是做不離子方程的簡化形式 Bortfeld 近似式（公式 2）。C 口ln(2)1Vp2工 cos 弓 siny222J=丿(=曲円1)他1譏2)2+v1v v2vp1HCOS2vp1ln(vp2)_ln(vp2vs1)vp1vp1vs2其中 Rpp 表示縱波反射系數(shù)。vp1 為入射層的縱波速度，vp2 為透射層的縱波速 度。91 為入射角，92 為透射角。p1 為入射層密度，p2 為折射層密度。Vs1 為入射層橫波速度，Vs2為折射層的橫波速度。在本文中由于在臨界角以外只有 入射角，此時的反射系數(shù)可能會存在誤差。 在第六層與第七層的反射系數(shù)過小造 成在最后的合成單炮記錄上幾乎肉眼不

3、能看到。動校正的方法主要使用了常規(guī)的雙曲動校正和無拉伸動校正， 常規(guī)雙曲動校正的 原理是公式（3），在炮檢距較大時動校正會有明顯的動校拉伸， 它的拉伸系數(shù)為 公式（4）。t2f2x2/v2to是自激自收的時間，x 為炮檢距，v 為動校正的速度。在實際中動校正的速度 應(yīng)該是變化的，但是為了計算的方便一般是設(shè)定成訂制， 這就導(dǎo)致了動校正在大 偏移距時發(fā)生動校誤差。k*：t/t（4） t 為反射點(diǎn)的時間與校正到點(diǎn)的時間的差，t 為校正的點(diǎn)的時間。而無拉伸動校正的原理是在原始地震剖面上實現(xiàn)反射層位的部分實現(xiàn)整體搬離。本文中由于動校速度為定值，在大偏移距時存在比較明顯的誤差。 在反射見面附 近的動校速度

4、為給定的動校速度，在該區(qū)域外運(yùn)用的速度為插值得到的動校速度 進(jìn)行常規(guī)動校，該動校的好壞有直接關(guān)系的是雙曲線反射界面的追蹤的準(zhǔn)確程(3)度。圖 2 2 常規(guī)動校正的流程圖圖 1 1 合成地震記錄的流程圖建立用于常規(guī)動校正需要知道的參數(shù)利用在單炮記錄計算過加速度和時間to插值 樣點(diǎn)的速度和時間【程中的計算的疊 宜得到相應(yīng)的每層采進(jìn)行雙曲線動校正的編寫，并 求出拉伸系數(shù)循環(huán)進(jìn)行上述過程，直到得到 的點(diǎn)超出給定矩陣輸出得到的新矩陣點(diǎn)賦值到:的相應(yīng)的I扁點(diǎn)的值賦值廠圖 4 4 追蹤的射線路徑該動校的好壞與中心線的確定有直接關(guān)系圖 3 3 無拉伸動校的流程圖H0imoism2KH3504ODZ2圖 5 5

5、 時間與炮檢距的變化圖 4 中在起始位置的由大到小為層位由小到大。圖 6 6 反射系數(shù)由于在臨界角以后求得的反射系數(shù)為復(fù)數(shù)， 在圖中取的模長，反射系數(shù)第一層界 面為深藍(lán)線，第二層反射界面為紅線，第三層反射系數(shù)為黃線，第四層的反射系 數(shù)為紫線，第五層的反射系數(shù)為綠線，第六層的反射系數(shù)為淺藍(lán)線。圖 9 9 均方根速度的動校正M40ST607DM圖 7 7 反射系數(shù)在時間域arcjk?i40mi圖 8 8 單炮記錄0 100圖 1212 基于圖 9 9 的無拉伸動校正圖 1010 均方根速度存在校正過量調(diào)整圖 1111 切除拉伸系數(shù) 30%30%程序主程序部分clear all;% 清除內(nèi)存clos

6、e 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; % 接收道數(shù)目 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;%由于均方根速度把存在明顯地校正過量對速度調(diào)整vnmo(2)=2265;vnmo(3)=2400;vnmo(4)=2760;vnm

8、o(5)=2976;vnmo(6)=3143;%程序所設(shè)中間變量及最終結(jié)果變量sx=zeros(1,105);sy=zeros(1,105);%sx 和 sy 是畫射線追蹤圖的各點(diǎn)的橫縱坐 標(biāo)a=zeros(nc,nd+1);c=zeros(nc,1);%a 的第 2 列到第 102 列存放入射角， 第 1 列 存放0 值，用于統(tǒng)一二分法射線追蹤的下限值； c 用于存放對地下每一界面的入 射角time=zeros(nc,nd);%time 用于存放追蹤出的雙程旅行時； jl=zeros(2400,nd);% 用于存放形成的地震記錄，記錄長度 4000m，nd 道接收%產(chǎn)生雷克子波作為地震地波f

9、m=40; % 主頻L=256; %L 為時間域采樣點(diǎn)數(shù) time_sample=0.001; % 采樣間隔為 1mst=(-L/2+1:L/2)*time_sample; % 子波延續(xù)時間從 -L/2ms 到 L/2ms， 共采集 L 個點(diǎn)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 離散，對應(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(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;%畫水平層狀介質(zhì)的各個界面endplot(0,0,r v);hold on;% 畫炮點(diǎn)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;% 畫檢波點(diǎn)%設(shè)置讓每一層的入射角不會大于臨界角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;%畫第一層的射線路經(jīng)end%追蹤其余入射角和旅行時for j=2:nd+1for i=2:ncA=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è)置炮點(diǎn)到入射點(diǎn)的橫向距離的初值； 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);% 初始化一下，防止由于前面道的計 算而遺留的坐標(biāo)值影響后面道的計算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);% 計算每一層射線射線路徑的橫 坐標(biāo)sx(2*i+2-m)=o

15、s+(j-2)*dj-sx(m);sy(m+1)=sy(m)-h(m);% 計算每一層射線射線路徑的縱坐標(biāo) sy(2*i+2-m)=sy(m);endplot(sx,sy);hold on;% 畫射線路徑endendhold off;Re=zeros(nc,nd);% 反射系數(shù)的求取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( 反射系數(shù) ); %根據(jù)射線追蹤的結(jié)果合成單炮記錄 for j=1:ndfor i=1:ncsjfftw=ampw.*exp(1i.*angw);%振幅譜結(jié)合相位譜， 重構(gòu)變換后的復(fù)信號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( 常規(guī)動校正 );% % % 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( 切除后的常規(guī)動校正 );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( 反射系數(shù) )hold on;end反射系數(shù)求取的程序： 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