inv_taup.m

著名的seismiclab的代码 是地震学研究人员必备的工具

function [m] = inv_taup(d,dt,h,q,N,flow,fhigh,mu,type);
%INV_TAUP: Inverse Radon transform. 
%          Given a cmp gather this subroutine is used
%          to compute the  Radon panel. By Radon Transform
%          we mean the Linear or Parabolic Radon Transform.
%
%  [m] = inv_taup(d,dt,h,q,N,flow,fhigh,mu,type)
% 
%  IN   d:     seismic traces   
%       dt:    sampling in sec
%       h(nh): offset or position of traces in meters
%       q(nq): ray parameters  if N=1
%              residual moveout at far offset if N=2
%       N:     N=1 Linear tau-p  
%              N=2 Parabolic tau-p
%       flow:  freq.  where the inversion starts in HZ (> 0Hz)
%       fhigh: freq.  where the inversion ends in HZ (< Nyquist) 
%       mu:    regularization parameter 
%       type:  type='ls' least-squares solution
%              type='hr' non-iterative High Resolution solution
%
%  OUT  m:     the linear or parabolic tau-p panel
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi, last modified December 10,  1998.
%  [e-mail]
%
%  Copyright (C) 1998 Signal Analysis and Imaging Group
%                     Department of Physics
%                     The University of Alberta
%

[nt,nh]= size(d);
nq = max(size(q));

if N==2; h=h/max(abs(h));end;  

nfft = 2^nextpow2(nt);

D = fft(d,nfft,1);
M = zeros(nfft,nq);
i = sqrt(-1);

ilow  = floor(flow*dt*nfft)+1; if ilow<1; ilow=1;end;
ihigh = floor(fhigh*dt*nfft)+1;
if ihigh>floor(nfft/2)+1; ihigh=floor(nfft/2)+1;end

Q = eye(nq);
Qtr = nq;
Mtr = nq*nh;
beta = mu*Mtr/Qtr;

for ifreq=ilow:ihigh

 f = 2.*pi*(ifreq-1)/nfft/dt;
 L = exp(i*f*(h.^N)'*q); 
 y = D(ifreq,:)';
 x = L'*y;

 MATRIX = L'*L; 

 A = MATRIX + beta*Q;
 x = inv(A) *L'* y; 

 M(ifreq,:) = x';
 M(nfft+2-ifreq,:) = conj(x)';

if type=='hr'; 
   xx = x /max(abs(x)); Q = diag(1./(abs(xx).^2+0.0001)); 
end

end
 
M(nfft/2+1,:) = zeros(1,nq);
m = real(ifft(M,[],1));
m = m(1:nt,:);

return