%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% SPECTRUM RATIO METHOD FOR Q ESTIMATION ON ONE CMP GATHER %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Qn=spec_ratio_method(u_nmo,Noffset,Nx,Dx,Nt,Dt,v_r,nf_low,nf_high,err)
%% 参数设置 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nf = Nt; %频率采样点数
Df = 1/(Dt*Nt); %频率采样间隔
Fs = Nf*Df; %采样频率
omg_h = pi*Fs; %Nyquist频率
f = (0:Nf-1)*Df;
omg = 2*pi*f;
%% Q值估算 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 根据波峰和零点在参考道寻找完整波形
n=0; %%%对找到的波形进行计数%%%
for j=2:Nt-1
if(u_nmo(j,1)>u_nmo(j-1,1) && u_nmo(j,1)>u_nmo(j+1,1) && abs(u_nmo(j,1))>(0.5*err))
n=n+1;
j_max(n) = j;
for j1 = 2:j_max(n)-1
if((u_nmo(j_max(n)-j1+1,1)>=u_nmo(j_max(n)-j1,1) && u_nmo(j_max(n)-j1+1,1)>u_nmo(j_max(n)-j1+2,1))||(u_nmo(j_max(n)-j1+1,1)>u_nmo(j_max(n)-j1+2,1) && abs(u_nmo(j_max(n)-j1+1,1))<=(0.005*err)));
j_left(n) = j_max(n)-j1+1;
break;
else
j_left(n) = 1;
end
end
for j2 = j_max(n)+1:Nt-1
if((u_nmo(j2,1)>u_nmo(j2-1,1) && u_nmo(j2,1)>=u_nmo(j2+1,1))||(u_nmo(j2,1)>u_nmo(j2-1,1) && abs(u_nmo(j2,1))<=(0.005*err)))
j_right(n) = j2;
break;
else
j_right(n) = Nt;
end
end
end
end
for j=2:Nt-1
if(u_nmo(j,1)>u_nmo(j-1,1) && u_nmo(j,1)>u_nmo(j+1,1) && abs(u_nmo(j,1))>(0.5*err))
n=n+1;
j_max(n) = j;
for j1 = 2:j_max(n)-1
if((u_nmo(j_max(n)-j1+1,1)>=u_nmo(j_max(n)-j1,1) && u_nmo(j_max(n)-j1+1,1)>u_nmo(j_max(n)-j1+2,1))||(u_nmo(j_max(n)-j1+1,1)>u_nmo(j_max(n)-j1+2,1) && abs(u_nmo(j_max(n)-j1+1,1))<=(0.005*err)));
j_left(n) = j_max(n)-j1+1;
break;
else
j_left(n) = 1;
end
end
for j2 = j_max(n)+1:Nt-1
if((u_nmo(j2,1)>u_nmo(j2-1,1) && u_nmo(j2,1)>=u_nmo(j2+1,1))||(u_nmo(j2,1)>u_nmo(j2-1,1) && abs(u_nmo(j2,1))<=(0.005*err)))
j_right(n) = j2;
break;
else
j_right(n) = Nt;
end
end
end
end
[st1,tt1,ff1] = st(u_nmo(:,1),0,Fs/2,Dt,Df); %%%对参考道进行S变换%%%
% 将参考道找到的完整波形对应采样点振幅谱相加
spec1=zeros(nf_high-nf_low+1,n);
for jn=1:n
for j=j_left(jn):j_right(jn)
spec1(:,jn)=spec1(:,jn)+st1(nf_low:nf_high,j);
end
end
for i=1:Nx-1
[st2,tt2,ff2] = st(u_nmo(:,i+1),0,Fs/2,Dt,Df); %%%对其它各道进行S变换%%%
% 在对应完整波形处计算Q值
spec2=zeros(nf_high-nf_low+1,n);
for jn=1:n
% 将参考道找到的完整波形对应采样点振幅谱相加
for j=j_left(jn):j_right(jn)
spec2(:,jn)=spec2(:,jn)+st2(nf_low:nf_high,j);
end
spec_r=log(abs(spec2(:,jn))./abs(spec1(:,jn)));
k1=abs(polyfit(f(nf_low:nf_high),spec_r',1)); %%%求斜率%%%
k(jn,i)=k1(1);
t0=j_max(jn)*Dt;
t1=sqrt(t0^2+4*((Noffset/2)*Dx/v_r(j_max(jn)))^2); %%%求对应道的走时%%%
t2=sqrt(t0^2+4*((Noffset/2+i)*Dx/v_r(j_max(jn)))^2);
Q(jn,i)=pi*(t2-t1)/k1(1); %%%计算对应位置的Q值%%%
end
end
% 计算零偏移距Q值
for jn=1:n
k2=abs(polyfit(((Noffset+(1:(Nx-1)))*Dx).^2,Q(jn,:),1));
Qe(jn)=k2(2);
end
% 将等效Q值转化为层Q值
Qm(1)=Qe(1);
for jn=2:n
Qa=j_max(1)*Dt/Qm(1);
for j=2:jn-1
Qa=Qa+(j_max(j)-j_max(j-1))*Dt/Qm(j);
end
Qm(jn)=(j_max(jn)-j_max(jn-1))*Dt/(j_max(jn)*Dt/Qe(jn)-Qa);
end
% 对Q值按时间采样间隔进行分层处理
for j=1:j_max(1)
Qn(j)=Qm(1);
end
for jn=2:n
for j=j_max(jn-1)+1:j_max(jn)
Qn(j)=Qm(jn);
end
end
for j=j_max(n)+1:Nt
Qn(j)=Qm(n);
end
end
- 1
- 2
前往页