kirk_shot3d.m

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

end

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

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

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

if ibcfilter 
	% get a cumulative array from shotrec
	arycum=cumsum(shotrec);
end

%one way time
dt1=.5*dt;
t1=t/2;
nt=length(t);

%compute maximum time needed. This is the traveltime for a scatterpoint
%when the source and receiver are colocated a distance aper away.
vmin=min(velmod(:));
tmax=sqrt(tmig2^2 + (2*aper/vmin)^2);

%pad input to tmaxin
npad=ceil(tmax/dt1)-nsamp+5;
if( npad > 0)
	shotrec= [shotrec; zeros(npad,ntr)];
	t1 = [t1',(nsamp+1:nsamp+npad)*dt1]';
	if ibcfilter
        arycum=[arycum; ones(npad,1)*arycum(nsamp,:)];
	end
end
t2= t1.^2;

%output samples targeted
samptarget=near(t,tmig1,tmig2);
tmig=t(samptarget);

%build output grid
xgrid=min(xv):dx:max(xv);%row vector
ygrid=(min(yv):dx:max(yv))';%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);

%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);

%initialize output array
shotmig=zeros(length(samptarget),ntrmig);

%loop over migrated traces
%

disp(' ');
disp([' --- Total number of migrated traces : ' int2str(ntrmig) ' ---']);
disp(' ');

clock1=cputime;
steptimes=nan*ones(size(ntrmig));
ntimes=0;
ievery=20;%print a progress message every this many traces
%save1=zeros(size(ntrmig));
% loop over traces in aperture

for kmig=1:ntrmig 			% 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);
    
    %need average velocity at output point
    ivx=near(xv,xtr);
    ivy=near(yv,ytr);
    vrms_mig=velmod(:,ivx(1),ivy(1));
    vave_mig=vrms2vave(vrms_mig,tv);
    zmig=t1(samptarget).*vave_mig(samptarget);%depths at output point
    
	%shot offset and velocity 
    offsetshot2=(xtr-xshot)^2+(ytr-yshot)^2;
    xvshotside=(xshot+xtr)/2;
    yvshotside=(yshot+ytr)/2;
    ivshotx=near(xv,xvshotside);
    ivshoty=near(yv,yvshotside);
	vshot2 = velmod(:,ivshotx(1),ivshoty(1)).^2;
    
    %source vector: points from image point to source
    rs_vec=[(xshot-xtr)*ones(size(zmig)) (yshot-ytr)*ones(size(zmig)) -zmig];
    rs=sqrt(sum(rs_vec.^2,2));
    %cosine of the source angle
    cossource=zmig./(rs+100*eps);
        
    %gather=zeros(length(tmig),length(inaper));
	for kaper=1:length(inaper)
        xnow=x(inaper(kaper));
        ynow=y(inaper(kaper));
        %receiver offset and velocity
        offsetrec2=(xtr-xnow)^2+(ytr-ynow)^2;
        xvrecside=(xnow+xtr)/2;
        yvrecside=(ynow+ytr)/2;
        ivrecx=near(xv,xvrecside);
        ivrecy=near(yv,yvrecside);
        vrec2 = velmod(:,ivrecx(1),ivrecy(1)).^2;
        
        %receiver vector: points from image point to receiver
        rg_vec=[(xnow-xtr)*ones(size(zmig)) (ynow-ytr)*ones(size(zmig)) -zmig];
        rg=sqrt(sum(rg_vec.^2,2));
        %cosine of the receiver angle
        cosrec=zmig./(rg+100*eps);
        
		% source-receiver travel time via double square root equation
        % "t2" is the squared one-way vertical time at the output point
        tsmig=sqrt(offsetshot2./vshot2(samptarget) + t2(samptarget)); %time from shot to mig point
        trmig=sqrt(offsetrec2./vrec2(samptarget) + t2(samptarget));%time from receiver to mig point
        tsr= tsmig+trmig+100*eps;
        
        %cosine of opening angle
        costheta=sum(rg_vec.*rs_vec,2)./(rg.*rs+100*eps); %just the dot product formula

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

        ind1 = find( cosangle < cos(ang_limit));
        ind2 = find( cosangle(ind1) > cos(anglemax) );
        
        angleweight=ones(size(cosangle));
        angleweight(ind1)=0;
        angleweight(ind1(ind2))=(cosangle(ind1(ind2))-cos(anglemax))/(cos(ang_limit)-cos(anglemax));
        
%         i1=1;
%         if(~isempty(ind))
%             i1 = ind(end);
%         end
%         ind = find( cosangle < cos(ang_limit) );
%         i2=1;
%         if(~isempty(ind))
%             i2 = ind(end);
%         end
% 
%         if i1 < i2
%             if itaper2  ==  0
%                 coef2 = lin_taper(i2,i1);
%             else
%                 coef2 = cos_taper(i2,i1);
%             end
%             cosangle(1:i1) = zeros(i1,1);
%             cosangle(i1+1:i2) = coef2(i2-i1:-1:1)'.*cosangle(i1+1:i2);
%         end
	
		% boxcar anti-aliasing filter
		if ibcfilter
			lt0=round((dx*tanalpha./velmod(samptarget,ktr)/dt));
			indt = round((tsr/dt))+1;
			lentr = nsamp+npad;
			lt = ones(lentr,1)*max(lt0);
			lt(indt)=lt0;
			lt(max(indt)+1:lentr) = ones(lentr-max(indt),1)*min(lt0);
			it = (1:lentr)';
			l1=it-lt-1;
			l2=it+lt;
			ind = find(l1 < 1);
			l1(ind) = ones(length(ind),1);
			ind = find(l2> lentr);
			l2(ind)=ones(length(ind),1)*lentr;
			tmp0=t;
			tmp0(1) = arycum(1,inaper(kaper));
			ind = 2:lentr;
			tmp0(ind) = (arycum(l2(ind),inaper(kaper))-arycum(l1(ind),inaper(kaper)))...
				    ./(l2(ind)-l1(ind));
		else
			tmp0 = shotrec(:,inaper(kaper));
		end	

		%interpolation
		% Linear
		if interp_type == 1
			tnumber = tsr/dt;
			it0 = floor( tnumber ) + 1;
			it1 = it0+1; 
			xt0 = tnumber - it0 + 1;
			xt1 = it0-tnumber;
			tmp = xt1.*tmp0(it0)+xt0.*tmp0(it1);
		end
		% Spline
		if interp_type == 2
			tmp = interp1(t,tmp0,tsr,'spline');
		end
		% Cubic
		if interp_type == 3
			tmp = interp1(t,tmp0,tsr,'cubic');
		end
		% Sinc
		if interp_type == 4
			tmp = sinci(tmp0,t,tsr);
        end


		% aperture taper
		aperweight = 1.0/length(inaper);
        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
		tmp = tmp .* aperweight.*angleweight;
		
        %see Bleistein, Cohen, and Stockwell p247 eqn 5.2.22
		tmp = tmp.*zmig.*rs.*costheta./(rg.*vave_mig+100*eps).^2;
        %save1(kmig)=save1(kmig)+sum(costheta);
        %gather(:,kaper)=tmp;
        ind=find(isnan(tmp),1);
		shotmig(:,kmig)= shotmig(:,kmig)+tmp;
		
	end
	
	% scaling and 90 degree phase shift
% 	scalemig = sqrt(vrec2(samptarget)).*sqrt(pi.*(tmig+0.0001)) ;
% 	shotmig(:,kmig) = shotmig(:,kmig)./scalemig ;
    %$$$$$ Don't I need to divide by number of traces summed?
    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 ' num2str(timenow) ' estimated remaining ' num2str(timeremaining) ]);
    end

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

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