example_taup.m

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

%
%  example_taup.m
%
%  Parabolic Radon Transform
%
%  SeismicLab
%  Version 1
%
%  written by M.D.Sacchi, last modified Sept., 2001
%  [e-mail]
%
%  Copyright (C) 1998 Seismic Processing and Imaging Group
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define parameters to generate data using 
% Forward Radon transform 

 np = 40;
 nt = 100;
 dt = 0.004;

 M = zeros(nt,np);
 M(nt/2+5,15) = 0.5;

 M(nt/2,1) = 0.5;

% Radon parameter axis (Curvature of the Parabola)

 pmax =0.000001;
 pmin =0;
 dp = (pmax - pmin)/(np-1);
 p = pmin:dp:pmax;
 
% Offset axis

 nh = 30;
 dh = 10;
 h = 0:dh:dh*(nh-1);

% Max and Min freq in data in Hz

flow = 5.;
fhigh = 86.;  

[D] = for_taup(M,dt,h,p,2,flow,fhigh);
Data=D;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Do the inversion (Find the RADON panel)
% mu is a regularization paramter

mu = 0.101;

type = 'hr'  % type ='ls'

[Minv] = inv_taup(D,dt,h,p,2,flow,fhigh,mu,type);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Normalize data and display

hm = max(h)^2;

M = M/max(max(M));
D = D/max(max(D));
Minv = Minv/max(max(Minv));

figure(1); clf;

t = 0:dt:(nt-1)*dt;

subplot(131); wigb(D,1,h,t);  
              title(' t - offset');
              xlabel('offset (m)'); 
              ylabel('t (s)');

subplot(132); wigb(Minv,1,p*hm,t); 
              title('Inverted \tau - q'); 
              xlabel('q (s)');
              ylabel('\tau (s)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Do filtering in the Radon panel and then came back
% to offset space



Minv(:,3:np)=0;

[D] = for_taup(Minv,dt,h,p,2,flow,fhigh);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Normalize data and display

figure(2); clf;

subplot(131); wigb(Minv,0.5,p*hm,t); 
              title('Filtered \tau  - p');
              xlabel('p (s/m)'); 
              ylabel('\tau (s)');

subplot(132); wigb(D,1,h,t); 
              title('Recovered signal in t - offset');
              xlabel('offset (m)');
              ylabel('t (s)');