多道面波分析相移法频散曲线提取方法MATLAB程序

%%%相移法Park文章 function [EN,fre,v] = PhaseShift(rec,dt,offset,dx,vmin,dv,vmax,fmin,df,fmax) [Nt,Ntrace]=size(rec);%%读取单炮记录的时间采样数与道数 v=vmin:dv:vmax; % 速度向量 lv = length(v); % 相速度向量的长度 Fs = 1/dt; % 采样率 f = Fs*(0:(Nt/2))/Nt; % 频率向量 % 确定频散能量图中最小频率在频率向量中的索引 ind = find(f>fmin); fmin_ind = ind(1)-1; % 确定频散能量图中最大频率在频率向量中的索引 ind = find(f>fmax); fmax_ind = ind(1); freq=f(fmin_ind:fmax_ind); % 显示范围频率序列 for j=1:Ntrace rec_fft(:,j) = fft(rec(:,j),Nt,1); % 获得地震记录的频谱 end rec_fft_amp = abs(rec_fft); % 频谱振幅 rec_fft_n = rec_fft./rec_fft_amp; % 对地震记录的频谱作归一化 E = zeros(lv,fmax_ind); % 频散能量矩阵 %%%%%%%%%%%%%%%%%%%计算频散能谱 for j=1:fmax_ind for i=1:lv for k=1:Ntrace x=offset+(k-1)*dx;%%%%从第一道依次叠加相移后能谱 E(i,j) = E(i,j)+exp(1i*2*pi*f(j)*x/v(i))*rec_fft_n(j,k);%%%相移叠加过程 end end E(:,j) = abs(E(:,j)./max(abs(E(:,j)))); % 对每个频率的频散能量作归一化 end EF=E(:,fmin_ind:fmax_ind);%%截取频散能谱对应的频带范围内值 [V0,F0]=meshgrid(freq,v);%%%生成显示初始网格 fre=fmin:df:fmax;%%%根据给定频率间隔生成显示频率序列 [Vf,Ff]=meshgrid(fre,v);%%%生成最终初始网格 EN = interp2(V0,F0,EF,Vf,Ff,'cubic');%%%%插值生成最终的显示频散能谱 EN=EN.^10;%%%为突出能量峰值，进行幂函数计算。