kirk_shot3dfz.m

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

if( isnan(params(14)) ) 
		ymig1 = min(y);
else
		ymig1 = params(14);
end
if( isnan(params(15)) ) 
		ymig2 = max(y);
else
		ymig2 = params(15);
end
if ymig2 < ymig1
 	error(['the start y location of target trace range should be less than the end y location'...
        'i.e. params(13) < params(14)']);
end

if(isnan(params(16)))
    dx=max([median(diff(x)) median(diff(y))]);
else
    dx=params(16);
end
if(dx<=0)
    error('output grid size must be greater than zero')
end
if(isnan(params(17)))
    ievery=100;
else
    ievery=params(17);
end

if(isnan(params(18)))
    ttflag=0;
else
    ttflag=params(18);
end

%one way time
dt1=.5*dt;

vmin=min(vrms);
vmax=max(vrms);
%output samples targeted
if(nargin<10)
    %build a zmig vector
    dz=vmin/(4*fmax);
    %estimate a zmax if need be
    if(zmig2==-1)
        zmig2=vmax*max(t)/2; %should use average velocity here but it won't matter much
    end
    zmig=(zmig1:dz:zmig2)';
else
    zmig=zmig(:);%Force column vector
    if(zmig2==-1)
        zmig2=max(zmig);
    end
    zmig1=min(zmig);
end
nz=length(zmig);

%compute maximum time needed. This is the traveltime for a scatterpoint
%when the source and receiver are colocated a distance aper away.
tmax=sqrt(zmig2^2 + 4*aper)/vmin;
%pad input to tmax
npad=ceil(tmax/dt1)-nsamp+5;
if( npad > 0)
	shotrec= [shotrec; zeros(npad,ntr)];
	t = [t;((nsamp+1:nsamp+npad)')*dt1];
end

%need average velocity
vave=vrms2vave(vrms,tv);
zv=tv.*vave/2;%depth axis for velocity function
%t1=pwlint(zv,tv,zmig)/2;%one way times to each output point
t1=interp1(zv,tv,zmig)/2;%one way times to each output point
t2=t1.^2;

%get rms and average velocities at each output depth
vrms_mig=pwlint(zv,vrms,zmig);
vrms_mig=vrms_mig(:);
vave_mig=pwlint(zv,vave,zmig);
vave_mig=vave_mig(:);
if(~ttflag) vrms_mig2=vrms_mig.^2; end
%tmig=2*zmig./vave_mig;


%build output grid
xgrid=min(x):dx:max(x);%row vector
ygrid=(min(y):dx:max(y))';%column vector
%find locations in target window
ixtarget= near(xgrid,xmig1,xmig2);
iytarget= near(ygrid,ymig1,ymig2)';
xmig=ones(size(iytarget))*xgrid(ixtarget);
ymig=ygrid(iytarget)*ones(size(ixtarget));
xmig=xmig(:)';
ymig=ymig(:)';
ntrmig=length(xmig);

if(ttflag)
    %compute traveltime table
    omax=sqrt((max(x)-min(x))^2+(max(y)-min(y))^2);%max offset
    do=dx/10;%sample increment in offset for table
    offs=(0:do:omax)';
    noffs=length(offs);
    offs=linspace(0,omax,noffs);
    do=offs(2);
    timetable=zeros(length(offs),length(zmig));
    for kz=1:length(zmig)
        timetable(:,kz)=sqrt(t2(kz)+offs.^2/vrms_mig(kz)^2);
    end
end

%Differentiate and phase rotate the shotrecord
nt2=2^nextpow2(nt);
[shotftmp,f]=fftrl(shotrec,t,0,nt2);
shotftmp=-2*pi*i*f(:,ones(1,length(x))).*shotftmp;
shotftmp(end,:)=0;%zero the nyquist
shotrec=ifftrl(shotftmp,f);

% ifuse=between(fmin,fmax,f);%indices of frequencies to use
% ifuse=ifuse(1:incf:end);
% omega=2*pi*f(ifuse);%column vector of frequencies
% shotf= omega(:,ones(1,length(x))).*shotftmp(ifuse,:).*exp(-i*pi/2); %multiply by omega
% eio=exp(i*omega);
% clear shotrec shotftmp f;
%initialize output array
shotmig=zeros(nz,ntrmig);

%
%loop over migrated traces
%

disp(' ');
disp([' --- Total number of input traces : ' int2str(ntr) ' ---']);
disp([' --- Total number of migrated traces : ' int2str(ntrmig) ' ---']);
if(ttflag)
    disp('Using traveltime table lookup')
else
    disp('Using traveltime calculations per sample')
end
disp(['Output grid spacing ' num2str(dx) ' length units'])
switch incf
    case {1}
        disp('Imaging every frequency')
    case{2}
        disp('Imaging every other frequency')
    case{3}
        disp('Imaging every third frequency')
    case{4}
        disp('Imaging every fourth frequency')
    case{5}
        disp('Imaging every fifth frequency')
    otherwise
        disp(['Trying something really ridiculous and imaging every '...
            int2str(incf) 'th frequency'])
end
clock1=cputime;
steptimes=nan*ones(size(ntrmig));
ntimes=0;
% loop over traces in aperture
% ntrmig=1;
% xmig=xshot;
% ymig=yshot;

x=x(:);%want coordinate vectors in columns
y=y(:);

for kmig=1:ntrmig 			% ktr--the index of the output trace
%for kmig=21 			% ktr--the index of the output trace
    xtr=xmig(kmig);%x coordinate of target
    ytr=ymig(kmig);%y coordinate of target

	%determine input traces in aperture
    inaperx=between(xtr-aper,xtr+aper,x,2);
    inapery=between(ytr-aper,ytr+aper,y(inaperx),2);
    inaper=inaperx(inapery);
    
    xnow=x(inaper);
    ynow=y(inaper);
    
	%shot offset
    if(ttflag)
        offsetshot=sqrt((xtr-xshot)^2+(ytr-yshot)^2);
    else
        offsetshot2=(xtr-xshot)^2+(ytr-yshot)^2;
    end
        
    %gather=zeros(length(zmig),length(inaper));
	for kz=1:length(zmig)
        %source vector: points from image point to source
        rs_vec=[(xshot-xtr)*ones(size(xnow)) (yshot-ytr)*ones(size(xnow)) -zmig(kz)*ones(size(xnow))];
        rs=sqrt(sum(rs_vec.^2,2));
        %cosine of the source angle
        cossource=zmig(kz)./rs;
        
        %receiver vector: points from image point to receiver
        rg_vec=[(xnow-xtr) (ynow-ytr) -zmig(kz)*ones(size(xnow))];
        rg=sqrt(sum(rg_vec.^2,2));%column vector of mag(receiver vector)
        %cosine of the receiver angle
        cosrec=zmig(kz)./rg;
        
        if(ttflag)
            %receiver offset
            offsetrec=sqrt((xtr-xnow)^2+(ytr-ynow)^2);
            % source-receiver travel time via double square root equation
            % look up times in table via offset
            ioffshot=round(offsetshot/do)+1;
            ioffrec=round(offsetrec/do)+1;
            tsr= timetable(ioffshot,kz)+timetable(ioffrec,kz)+100*eps;%column vector of traveltimes
        else
           %receiver offset
            offsetrec2=(xtr-xnow).^2+(ytr-ynow).^2; 
            tsr= sqrt(offsetshot2/vrms_mig2(kz) + t2(kz))+sqrt(offsetrec2/vrms_mig2(kz) + t2(kz))+100*eps;
        end
        
        %cosine of opening angle
        costheta=sum(rg_vec.*rs_vec,2)./rg.*rs; %just the dot product formula

        %angle limit and the taper
        cosangle=min([cossource cosrec],[],2);

        ind1 = find( cosangle < cosanglimit);
        ind2 = find( cosangle(ind1) > cosanglemax );
        
        angleweight=ones(size(cosangle));
        angleweight(ind1)=zeros(size(ind1));
        angleweight(ind1(ind2))=(cosangle(ind1(ind2))-cosanglemax)/(cosanglimit-cosanglemax);
        
		% aperture taper
		aperweight = 1.0/length(inaper);%divide by number of traces in aperture
        xtest=abs(aper-abs(xtr-xnow));%distance of trace from edge of aper
        if xtest < width1
            if(itaper1==1)
                aperweight=(.5+.5*cos(pi*(xtest-width1)/(180*width1)))/length(inaper);
            else
                aperweight=(xtest-width1)/(width1*length(inaper));
            end
        end
        
        %get a slice from the input data
        tsrrounded=round(tsr/dt)*dt;
        s=slicemat(shotrec(:,inaper),tsrrounded/dt+1,round(tslicewidth/(2*dt)),1);
        %transform the slice, windowed with a gaussian
        ntwid=size(s,1);%slice width in samples
        n2=2^nextpow2(ntwid);
        itrel=1:ntwid;%relative time index
        n0=(ntwid-1)/2+1; %relative index to center of slice
        trel=dt*(itrel-n0);
        gwin=exp(-((itrel-n0)/(n0/2)).^2)';
        [sftmp,f]=fftrl(s.*(gwin(:,ones(size(inaper)))),trel,0,n2);
        %get frequencies to use, etc....
        ifuse=between(fmin,fmax,f);%indices of frequencies to use
        ifuse=ifuse(1:incf:end);
        omega=2*pi*f(ifuse);%column vector of frequencies
        sf= sftmp(ifuse,:);
        %apply phase shift with times adjusted to the slice time.
        tmp=real(sum(exp(i*omega*(tsr-tsrrounded-trel(1))').*sf(:,inaper)))';    
        %see Bleistein, Cohen, and Stockwell p247 eqn 5.2.22
		shotmig(kz,kmig)= shotmig(kz,kmig)+sum(zmig(kz)*aperweight*tmp.*rs.*costheta...
            .*angleweight./(rg*vave_mig(kz)).^2);
        %gather(:,kaper)=tmp.*zmig.*rs.*costheta./(rg.*vave_mig).^2;
		
	end
	
    if(rem(kmig,ievery)==0)
	    disp([' Completed migrated trace no. ' ,int2str(kmig) ,' of ' int2str(ntrmig) ]);
        timenow=cputime-clock1;
        
        ntimes=ntimes+1;
        steptimes(ntimes)=timenow;
        if(ntimes>1)
            timeremaining = (length(x)/(length(inaper)+1))*(timenow-steptimes(ntimes-1))*(ntrmig-kmig)/ievery;
        else
            timeremaining = (length(x)/(length(inaper)+1))*timenow*(ntrmig-kmig)/ievery;
        end
        disp([' time so far ' int2str(timenow) ' S estimated remaining ' int2str(timeremaining) ' S' ]);
    end

end
ind=find(isnan(steptimes));
steptimes(ind)=[];

totaltime=cputime-clock1;
disp(['Total time required ' int2str(totaltime)])