kirk_shot3dfz.m

用matlab语言写的地震偏移成像软件集,很有用的东东呵

function [shotmig,zmig,xmig,ymig,xgrid,ygrid]=kirk_shot3Dfz(shotrec,t,x,y,xshot,yshot,vrms,tv,params,zmig)
% KIRK_SHOT3Dfz: Kirchhoff prestack shot-record time migration in 3D, depth output   
% 
% [shotmig,zmig,xmig,ymig,xgrid,ygrid]=kirk_shot3Dfz(shotrec,t,x,y,xshot,yshot,velmod,tv,params,zmig)
%
% KIRK_SHOT3Dfz migrates a single 3D shot record using prestack time
% migration, and outputs in depth. This is especially useful if only select
% depth planes are desired. The 3D shot record is stored in a 2D 
% matrix and it is not necessary that the traces fall on a regular spatial 
% grid. The algorithm is simple traveltime path summation with a few
% options as described in the parameters vector. Traveltime is described by
% a single, laterally invariant, vrms function. Time shifts are
% accomplished by frequency domain phase shift where the frequency domain 
% operations are done on narrow fairways defined around the traveltime for
% each event. The choice of the fairway width is a key parameter and, for
% 2ms data, it has been found that 10ms fairway (params(9)) seems to work
% well. There is no topographic compensation at present. The
% source and all receivers are assumed to be on the same horizontal plane.
% The Kirchhoff integral is done as described in Bleistein, Cohen, and
% Stockwell p247 eqn 5.2.22. Although this derivation is for constant
% velocity, the adjustments for time migration are reasonably obvious.
%
% shotrec ... 2D matrix containing the shot record. One trace per column.
% t ... time coordinate vector for shotrec. 
%           Requirement: length(t)=size(shotrec,1);
% x ... row vector of the first space coordinate for shotrec. Let
%   [nt,ntr]=size(shotrec), then we must have ntr=length(x);
%           Requirement: length(x)=size(shotrec,2);
% y ... row vector of the first space coordinate for shotrec. Let
%   [nt,ntr]=size(shotrec), then we must have ntr=length(y);
%           Requirement: length(y)=size(shotrec,2)
% vrms ... vector of rms velocities as a function of time
% tv ... time coordinate vector for vrmsd
% xshot ... x coordinate of the shot. It must be in the same coordinate
%       system as vector x. 
% yshot ... y coordinate of the shot. It must be in the same coordinate
%       system as vector y. 
%
% params ... vector of migration parameters. An entry of nan gets the
% default.
%
%    params(1--3) : migration aperture and its taper
%    params(1) ... physical aperture in meters. This is the largest lateral
%               distance a trace may have from the output location and still be
%               allowed to contribute.
%	          default is the largest diagonal of the velocity model
%    params(2) ... width of the aperture taper
%                 default is 0.05*params(1)
%    params(3) ... = 0, linear taper
%                     = 1, cosine taper   
%                 default is 1 (cosine taper)
%    params(4-5) : scattering angle limit (degrees)
%    params(4) ... maximum scattering angle limit (degrees)
%                 default = 60
%    params(5) ... width of angle limit taper (linear in cosine)
%                 default = min([0.15*params(4), 90-params(4)])
% NOTE: The angle limit described by params(4-5) is applied to both the
% source angle (angle w.r.t vertical of the vector from imagepoint to source)
% and the receiver angle (similarly defined). If either angle exceeds this
% limit, then the response is either damped or attenuated altogether.
%    paqrams(6) ... fmin, minimum frequency to migrate
%           default = 10 Hz
%    params(7) ... fmax, maximum frequency to migrate
%            default = 70% of Nyquist
%    params(8) ... incf, frequency increment (an integer). Every incf'th
%           frequency will be imaged.
%            default = 1
% NOTE: frequencies below fmin and above fmax are ignored
%    params(9) : time slice width (seconds)
%           default = 5 dominant periods = 10/(fmin+fmax);
%                  (If already a power of 2, it will not be extended).
%    params(10--12) : defines migration target window
%	 params(10) ... zmin of migration target window
%	          default = 0 
%	 params(11) ... tmax of migration target window
%	          default = something appropriate
%    params(12) ... xmin of target window
%	          default = min(x)
%    params(13) ... xmax of migration target window
%	          default = max(x) 
%    params(14) ... ymin of target window
%	          default = min(y)
%    params(15) ... ymax of migration target window
%	          default = max(y) 
%    params(16) ... dx = spatial grid size of migrated volume. Same in x
%               and y.
%             default = min([mean(diff(x)) mean(diff(y))])
%    params(17) ... ievery: write a message after computing this many
%               migrated traces
%             default = 100
%    params(18) ... ttflag: 0 means calculated traveltimes for each point,
%               1 means build a traveltime table and do table lookup.
%             default = 0
%
% zmig ... if specified on input, then it overrides params(9,10). This is a
%       vector of explicit depths to image at. That is, if specified, then
%       only these depths will be imaged. (if params(6)=0, then these are
%       interpreted as times.
%
% OUTPUT arguments
%
%    shotmig ... 2D matrix of the output migrated shot record
%    tmig/zmig ... t/z coordinates of migrated data (depends on params(6)
%    xmig ... x coordinates of migrated data (length(xmig)=size(shotmig,2))
%    ymig ... y coordinates of migrated data (length(ymig)=size(shotmig,2))
%    xgrid ... the unique x coordinates in the output grid
%    ygrid ... the unique y coordinates in the output grid
%
% By G.F. Margrave
% CREWES Project, U of Calgary, 2007
%
% NOTE: It is illegal for you to use this software for a purpose other
% than non-profit education or research UNLESS you are employed by a CREWES
% Project sponsor. By using this software, you are agreeing to the terms
% detailed in this software's Matlab source file.
 
% BEGIN TERMS OF USE LICENSE
%
% This SOFTWARE is maintained by the CREWES Project at the Department
% of Geology and Geophysics of the University of Calgary, Calgary,
% Alberta, Canada.  The copyright and ownership is jointly held by 
% its author (identified above) and the CREWES Project.  The CREWES 
% project may be contacted via email at:  [e-mail]
% 
% The term 'SOFTWARE' refers to the Matlab source code, translations to
% any other computer language, or object code
%
% Terms of use of this SOFTWARE
%
% 1) Use of this SOFTWARE by any for-profit commercial organization is
%    expressly forbidden unless said organization is a CREWES Project
%    Sponsor.
%
% 2) A CREWES Project sponsor may use this SOFTWARE under the terms of the 
%    CREWES Project Sponsorship agreement.
%
% 3) A student or employee of a non-profit educational institution may 
%    use this SOFTWARE subject to the following terms and conditions:
%    - this SOFTWARE is for teaching or research purposes only.
%    - this SOFTWARE may be distributed to other students or researchers 
%      provided that these license terms are included.
%    - reselling the SOFTWARE, or including it or any portion of it, in any
%      software that will be resold is expressly forbidden.
%    - transfering the SOFTWARE in any form to a commercial firm or any 
%      other for-profit organization is expressly forbidden.
%
% END TERMS OF USE LICENSE

%tstart=clock;
[nsamp,ntr]=size(shotrec);
nvsamp=length(vrms);

% check the validity input arguments

if(length(t)~=nsamp)
    error('Incorrect time specification')
end
t=t(:);
dt=t(2)-t(1);
nt=length(t);

if(length(x)~=ntr)
    error('Incorrect x specification for input data')
end

if(length(y)~=ntr)
    error('Incorrect y specification for input data')
end



%  ---- test velocity info ----

if(length(tv)~=nvsamp)
    error('Time vector for velocity model is incorrect')
end

if(min(tv)>min(t) || max(tv)<max(t))
    error('Time vector for velocity model must span that for data')
end


%  ---- examine parameters ----
nparams=18; 				% number of defined parameters
					
if(nargin<9) 				% no parameters inputted
	params= nan*ones(1,nparams); 
end 	 	

if(length(params)<nparams) 
	params = [params nan*ones(1,nparams-length(params))];
end

%assign parameter defaults

if( isnan(params(1)) ) 
	aper = sqrt((max(x)-min(x))^2+(max(y)-min(y))^2);
else
    aper = params(1);
end

if( isnan(params(2)) )
    width1 = aper/20;
else
    width1 = params(2);
end

if( isnan(params(3)) )
    itaper1 = 1;
else
    itaper1 = params(3);
end

if( isnan(params(4)) )
	ang_limit = pi/3;
else
	ang_limit = params(4)*pi/180;
end

if( isnan(params(5)) )
	width2 = 0.15*ang_limit;
    if(width2>pi/2-ang_limit)
        width2=pi/2-ang_limit;
    end
else
	width2 = params(5)*pi/180;
end
anglemax = ang_limit + width2;
cosanglemax=cos(anglemax);
cosanglimit=cos(ang_limit);
if( isnan(params(6)) )
	fmin=10;
else
	fmin = params(6);
end
if( isnan(params(7)) )
	fmax=.7/(2*dt);
else
	fmax = params(7);
end
if( isnan(params(8)) )
	incf=1;
else
	incf = params(8);
end
if(incf<1 || round(incf)~=incf)
    error('frequency increment must be a positive integer')
end
if( isnan(params(9)) )
    tslicewidth=10/(fmin+fmax);
else
	tslicewidth=params(9);
end

if( isnan(params(10)) ) 
		zmig1 = 0;
else
		zmig1 = params(10);
end

if( isnan(params(11)) ) 
		zmig2 = -1;%this is a flag to be adjusted later
else
		zmig2 = params(11);
end
if zmig2 < zmig1
    if(zmig2~=-1)
        error(['the target time window start time should be smaller than the end time !'...
                'i.e. paraams(8) < params(9)']);
    end
end

if( isnan(params(12)) ) 
		xmig1 = min(x);
else
		xmig1 = params(12);
end
if( isnan(params(13)) ) 
		xmig2 = max(x);
else
		xmig2 = params(13);
end
if xmig2 < xmig1
 	error(['the start x location of target trace range should be less than the end x location'...
        'i.e. params(11) < params(12)']);
end