cdp3d.m,v

具有特色的地震数据处理源码

head	3.4;
access;
symbols;
locks; strict;
comment	@// @;


3.4
date	2000.06.20.17.30.55;	author perron;	state Exp;
branches;
next	3.3;

3.3
date	2000.06.15.15.45.20;	author gilles;	state Exp;
branches;
next	3.2;

3.2
date	2000.06.14.18.52.44;	author gilles;	state Exp;
branches;
next	3.1;

3.1
date	2000.06.14.18.36.56;	author gilles;	state Exp;
branches;
next	3.0;

3.0
date	2000.06.13.19.18.02;	author gilles;	state Exp;
branches;
next	;


desc
@Release 3
@


3.4
log
@fix ntr = datain.fh{13}
@
text
@function [refpoints]=cdp3d(datain,v,strike,dip)

%[dataout]=cdp3d(datain,v,strike,dip)
%
%three dimensional CDP (common depth point) transform
%finds reflection points in 3 dimensions using rotations and ray 
%tracing method
%creates a 3 dimensional grid of bins, sums trace energy 
%corresponding in the bin  corresponding to its reflection point for 
%a given strike and dip of the geology, normalizes by the number of 
%traces contributing energy to each bin
%uses functions 'refpoints3d'
%
%use 'solve3d' for a picture of ray paths and reflection points
%use 'slicefmap' to extract 2D vertical or horizontal slices from the 
%3D volume these slices can then be displayed using 'dispseis' 
%(DSISoft) or else 'pcolor' or 'imagesc' (MATLAB).
%
%INPUT
%datain - dataset in DSI format
%     shot and receiver locations and travel time information are 
%     taken from the headers
%     trace headers 31, 29 and 33 are shot easting, northing and 
%     elevation respectively (x,y,z)
%     trace headers 37, 35, and 39 are receiver easting, northing, 
%     and elevation respectively
%     travel time startime:samplingrate:endtime (file headers 9, 8 
%     and 10 respectively)
%v - velocity of the waves (m/s)
%strike  - strike angle in degrees from North (or mine grid North)
%dip  - dip angle to the right, by convention, also in degrees
%
%OUTPUT
%dataout - the results of performing a 3D CDP transform on the input data
%    format is a modified DSI variable.  File and trace headers still 
%    exist but the .dat section has been replaced with a .fmap extension 
%    which contains a cell array (records) of 3 dimensional arrays.  
%    This is the results of the 3D CDP transform.  There are also .xsc, 
%    .ysc, and .zsc extensions, each containing vectors describing 
%    the coordinates of the centre of each of the bins.
%
%written by K.S. Beaty July, 1998
%last modified May, 1999
%by G. Perron
%Significantly modified May 2000
%by G. Bellefleur  
  
%$Id: cdp3d.m,v 3.3 2000/06/15 15:45:20 gilles Exp perron $
%$Log: cdp3d.m,v $
%Revision 3.3  2000/06/15 15:45:20  gilles
%Test to remove bad rx points
%
%
%Copyright (C) 1998 Seismology and Electromagnetic Section/
%Continental Geosciences Division/Geological Survey of Canada
%
%This library is free software; you can redistribute it and/or
%modify it under the terms of the GNU Library General Public
%License as published by the Free Software Foundation; either
%version 2 of the License, or (at your option) any later version.
%
%This library is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%Library General Public License for more details.
%
%You should have received a copy of the GNU Library General Public
%License along with this library; if not, write to the
%Free Software Foundation, Inc., 59 Temple Place - Suite 330,
%Boston, MA  02111-1307, USA.
%
%DSI Consortium
%Continental Geosciences Division
%Geological Survey of Canada
%615 Booth St.
%Ottawa, Ontario
%K1A 0E9
%
%email: dsi@@cg.nrcan.gc.ca

  
%disp('[dataout]=cdp3d(datain,v,strike,dip)')

tt=datain.fh{9}:datain.fh{8}:datain.fh{10}; %travel time vector
shot=[datain.th{1}(31,:)' datain.th{1}(29,:)' datain.th{1}(33,:)']; %shot x,y,z
rec=[datain.th{1}(37,:)' datain.th{1}(35,:)' datain.th{1}(39,:)'];  %receiver x,y,z

ntr=datain.fh{13};


for ii=1:ntr

   off(ii)=sqrt((shot(ii,1)-rec(ii,1))^2+(shot(ii,2)-rec(ii,2))^2+(shot(ii,3)-rec(ii,3))^2);
%  some trouble near the fb so we add +20 (to be fixed!)
   tt_fbsamp=ceil(off(ii)/(v*datain.fh{8}))+20;
   refpoints{ii}(1:tt_fbsamp,1)=0;
   refpoints{ii}(1:tt_fbsamp,2)=0;
   refpoints{ii}(1:tt_fbsamp,3)=0;
   refpoints{ii}(1:tt_fbsamp,4)=0; 
   refpoints{ii}(1:tt_fbsamp,5)=0;
   ii

   %find the reflection points in 3D 
   [cdps]=refpoints3d(shot(ii,:)',rec(ii,:)',tt(tt_fbsamp:length(tt)-1),dip,strike,v);
   
   refpoints{ii}(tt_fbsamp:length(tt)-1,1)=cdps(1,:);
   refpoints{ii}(tt_fbsamp:length(tt)-1,2)=cdps(2,:);
   refpoints{ii}(tt_fbsamp:length(tt)-1,3)=cdps(3,:);
   refpoints{ii}(tt_fbsamp:length(tt)-1,4)=tt(tt_fbsamp:length(tt)-1);
   refpoints{ii}(tt_fbsamp:length(tt)-1,5)=1.0;
   
% test to exclude cdp higher than sp

   test=refpoints{ii}(:,3)-shot(ii,3);
   ind=find(test<0);
   refpoints{ii}(ind,5)=0;
      
   ind1=find(imag(refpoints{ii}(:,1))~=0);
   ind2=find(imag(refpoints{ii}(:,2))~=0);
   ind3=find(imag(refpoints{ii}(:,3))~=0);
 
   refpoints{ii}(ind1,5)=0;
   refpoints{ii}(ind2,5)=0;
   refpoints{ii}(ind3,5)=0;
   
   clear ind ind1 ind2 ind3;
%   figure(1)
%   hold on
%   ind=find(refpoints{ii}(:,5)==1);
%   plot3(refpoints{ii}(ind,1),refpoints{ii}(ind,2),refpoints{ii}(ind,3),'.');
%   plot(refpoints{ii}(ind,1),refpoints{ii}(ind,2),'r.');
%   clear ind
end;

%   plot3(shot(1,1),shot(1,2),shot(1,3),'bv');
%   plot3(rec(:,1),rec(:,2),rec(:,3),'*k');
%   box on;
%   set(gca,'zdir','reverse');
%   rotate3d;

 
%end of function cdp3d
@


3.3
log
@Test to remove bad rx points
@
text
@d48 4
a51 2
%$Id:$
%$Log:$
d53 1
d88 1
a88 1
ntr=datain.fh{1};
@


3.2
log
@modified help text
@
text
@d87 1
d90 3
a92 1
   tt_fbsamp=floor(datain.th{1}(15,ii)/datain.fh{8});
d99 1
a99 1
      
d101 1
a101 1
   [cdps]=refpoints3d(shot(ii,:)',rec(ii,:)',tt(tt_fbsamp+1:length(tt)-1),dip,strike,v);
d103 5
a107 5
   refpoints{ii}(tt_fbsamp+1:length(tt)-1,1)=cdps(1,:);
   refpoints{ii}(tt_fbsamp+1:length(tt)-1,2)=cdps(2,:);
   refpoints{ii}(tt_fbsamp+1:length(tt)-1,3)=cdps(3,:);
   refpoints{ii}(tt_fbsamp+1:length(tt)-1,4)=tt(tt_fbsamp+1:length(tt)-1);
   refpoints{ii}(tt_fbsamp+1:length(tt)-1,5)=1;
a109 1
% Never found one so far!
d111 1
a111 1
   test=cdps(3,:)-shot(ii,3);
d114 8
d123 7
d132 7
@


3.1
log
@This version recovered from backup
@
text
@d1 1
a1 1
function [refpoints]=cdp3d_v2(datain,v,strike,dip)
d3 1
a3 1
%[dataout]=cdp3d_v2(datain,v,strike,dip)
d12 1
a12 1
%uses functions 'rotcoord' and 'findref3d'
d79 1
a79 1
disp('[dataout]=cdp3d_v3(datain,v,strike,dip)')
d115 1
a115 1
%end of function cdp3d_v3
@


3.0
log
@*** empty log message ***
@
text
@d1 1
a1 1
function [dataout]=cdp3d(datain,v,strike,dip,bins,records,limits,plt_flg)
d3 1
a3 1
%[dataout]=cdp3d(datain,v,strike,dip,bins,records,limits,plt_flg)
d6 1
a6 1
%finds reflection points in 3 dimensions using a rotation and ray 
d31 1
a31 6
%dip  - dip angle to the left, by convention, also in degrees
%bins - vector of bin widths in metres for fold map in x, y, and z 
%     directions.  should be set equal to sampling rate times velocity
%records - vector of records to work on eg. 1 or [1:3]
%limits - spatial limits of fold map [xmin xmax ymin ymax zmin zmax]
%plt_flg - use one or two sets of points (set = 1 or 2) 
a41 2
%DSISoft Version 1.0
%Customized VSP Processing Software
d43 34
a76 1
%last modified July 21, 1998
d78 2
a79 1
disp('[dataout]=cdp3d(datain,v,strike,dip,bins,records,limits,plt_flg)')
d82 2
a83 4
dataout=datain;
smp=datain.fh{8};
t1=datain.fh{9};
dataout=rmfield(dataout,'dat');
d85 1
a85 4
%vector of coordinates of centre of each fold map bin
xsc=limits(1):bins(1):limits(2);
ysc=limits(3):bins(2):limits(4);
zsc=limits(5):bins(3):limits(6);
d87 1
a87 6
for R=records %loop over records
 shot=[datain.th{R}(31,:)' datain.th{R}(29,:)' datain.th{R}(33,:)']; %shot x,y,z
 nshot=size(shot,1);
 rec=[datain.th{R}(37,:)' datain.th{R}(35,:)' datain.th{R}(39,:)']; %receiver x,y,z
 %find the reflection points in 3D using rotation/ray tracing method
 [refpoints]=findref3d(shot,rec,tt,v,strike,dip);
d89 7
a95 17
 fmap=zeros(length(xsc),length(ysc),length(zsc)); %initialize foldmap
 hits=fmap;

 for n=1:nshot %loop over shots
  if plt_flg==1
   fs=find(refpoints{n}(:,5)==1); %first set of reflection points
   rpts=refpoints{n}(fs,:);
  else
   rpts=refpoints{n}(:,:);
  end %if/else

  for i=1:size(rpts,1)
   %find bin associated with reflection points
   [z,f1]=min(abs(rpts(i,1)-xsc));
   [z,f2]=min(abs(rpts(i,2)-ysc));
   [z,f3]=min(abs(rpts(i,3)-zsc));
   if (round((rpts(i,4)-t1)/smp+1))>datain.fh{7}
d97 11
a107 8
   else
      amp=datain.dat{R}(round((rpts(i,4)-t1)/smp+1),n); 
      
      fmap(f1,f2,f3)=fmap(f1,f2,f3)+amp; %increment bin
      hits(f1,f2,f3)=hits(f1,f2,f3)+1; %increment bin
   end %(if)
  end %for i  
 end %loop over shots
d109 5
a113 16
 %set zeros on hits to NaN
 for a=1:size(hits,3)
  [i,j,v]=find(hits(:,:,a)==0);
  for b=1:length(i)
   hits(i(b),j(b),a)=NaN;
  end %for b
 end %for a
 %set dataout.fmap and normalize by number of hits to each bin
 dataout.fmap{R}=fmap./hits;
 for a=1:size(hits,3)
  [i,j,v]=find(isnan(hits(:,:,a)));
  for b=1:length(i)
   dataout.fmap{R}(i(b),j(b),a)=0;
  end %for b
 end %for a
end %loop over records
d115 1
a115 5
dataout.xsc=xsc;
dataout.ysc=ysc;
dataout.zsc=zsc;

%end of function cdp3d
@