www.pudn.com > logtools.rar > sonic2tz.m, change:2004-08-13,size:11548b


function [tz,zt,vins]=sonic2tz(sonic,z,nlegs,tnot,titlestring,flag)  % [tz,zt,vins]=sonic2tz(sonic,z,nlegs,tnot) % [tz,zt,vins]=sonic2tz(sonic,z,nlegs) % sonic2tz(sonic,z,nlegs,tnot,title,flag) ... for analysis display % sonic2tz(sonic,z,nlegs) ... for analysis display % % SONIC2TZ computes an approximate 2-way time depth curve from a sonic log for % use with depth conversion of logs. The approximation is a piecewise linear % curve whos points are taken directly from the exact t-z curve (which is % computed from directly integrating the sonic.) The sample points occur % at depths defined by the location of the largest peaks of the Hilbert  % envelope of the derivative of the sonic (the derivative of the sonic is % proportional to reflection coeficient at constant density.) At these sample % points the time depth approximation will give zero error compared to the % exact curve. The error between control points is reduced by computing % more legs on the approximate curve; however, many legged curves will lead % to long time-depth conversion run-times. To get help in choosing the number % of legs, run the function with no return arguments and it will produce a % display showing the fit of the approximation with estimated depth conversion % error. % % sonic = a vector containing the sonic log. Sonic should be in units of %   microseconds per length unit. % z = vector of depth's for the sonic log. % ******** length(z) must equal length(sonic) ********* % nlegs = number of legs on the piecewise linear approximation % NOTE: If sonic2tz is unable to find the requested number of legs, it %       will terminate with an error message. This can be avoided by %       entering nlegs as a negative number in which case execution will %       proced with a warning that fewer than nlegs were found. % tnot = starting time corresponding to z(1) % *********** default = 0.0 ********** % title = string with a title for the analysis displays % flag = 0 ... compute estimated error in time for analysis displays %      = 1 ... compute estimated error in depth for analysis displays % *********** default = 1 *********** % tz = vector of length nlegs+2 containing the 2-way times of the   %	approximate t-z curve % zt = vector of length nlegs+2 containing the depths of the approximate %	t-z curve % vins = vector of length nlegs+2 containing the instantaneous velocities  %	implied by the approximate t-z curve % % G.F. Margrave, March 1994 % % 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:  crewesinfo@crewes.org %  % 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 % % $Id: sonic2tz.m,v 1.3 2004/07/30 22:03:49 kwhall Exp $  if(nargin<4) 	tnot=0; end  if(nargin < 6) 	flag=1; end  iwarn=0; if( nlegs<0 ) 	iwarn=1; 	nlegs=abs(nlegs); end  %convert to column vectors [m,n]=size(sonic); sonic=sonic(:); z=z(:);   %compute detailed tz curve %compute delta z's nz=length(z); dz=diff(z); %compute the mean of each pair of sonic samples %s2=(sonic(1:nz-1)+sonic(2:nz))*.5; %integrate t=int_sonic(sonic,z)+tnot; if(length(t)<=abs(nlegs) & iwarn) 	disp(['Log too short for approximations, returning exact curve, length= '... 		int2str(length(t))]) 	tz=t; 	zt=z; 	vins= 1.e06 ./sonic; 	return; end   %compute approx rc's and then hilbert env and mask %(transpose is necessary to force gradient to give a real result) rcs=gradient(sonic',z'); [mask,henv]=hmask(rcs,.005); henv=abs(henv); spikes=mask.*henv;  %find the nlegs largest %kzone=nz/3; %nlegs=3*round(nlegs/3); %mlegs=nlegs/3; %ilegs=[]; %for k=1:3 % iz1=1+(k-1)*kzone; % iz2=min([iz1+kzone-1,nz]); %	spikezone=spikes(iz1:iz2); %	[spikezone,is]=sort(spikezone); %	spikezone=flipud(spikezone); %	is=flipud(is); %	ilegs=[ilegs (iz1+is(1:mlegs-1)-1)]; %end %[spikes,is]=sort(spikes); %spikes=flipud(spikes); %is=flipud(is); %ilegs=is(1:nlegs-1);  nlegsasked=nlegs; izone=nz*.06;      % Does anyone know what izone means? Please describe it here.                    % It -feels- like an overlap window or something to smooth                    % between legs. I don't really know though.                      ilegs=zeros(1,round(2*nlegs)); jlegs=1;  while(jlegs<nlegs)     spike2=spikes;       % take a copy of the spikes     ind=find(ilegs~=0);  % find all the bits of ilegs that are nonzero.                          % the first time through, there are none, so ind =                          % 0.                                for k=1:length(ind)          i1=max([1,ilegs(ind(k))-izone]);  % The minimum index of the ileg, including overlap         i2=min([ilegs(ind(k))+izone,nz]); % The max index of the ileg, including overlap         spike2(i1:i2)=zeros(size(i1:i2)); % set everything in this ileg range to zero on spike2.                                           % I think this part protects the                                           % ilegs that we have already                                           % determined. So it sort of masks                                           % them out on the spike2 set,                                           % which then constrains us to                                           % setting only those parts of                                           % ilegs which haven't been set                                           % before.     end          kmax=nlegs-length(ind);	              % we can only deal with the number of legs we have left.          for k=1:kmax         if( sum(spike2)==0.0) % If we have wiped out the entire spike2 in                               % that last operation, then we have filled                               % all of the slots in ilegs and so we're done. I                               % think so, anyhow.             break;         end         ii=find( spike2==max(spike2) ); % look for the biggest spike(s).                   ilegs(jlegs:jlegs+length(ii)-1)=ii; % set ilegs to be those highest spikes.                                              % If there are n spikes of                                             % the same height within the                                             % spikes, then n elements                                             % of ilegs are set.         for kk=1:length(ii)  % for each of the 'n' elements we just set            i1=max([1,ilegs(jlegs+kk-1)-izone]);  % find the minimum of the ileg including overlap            i2=min([ilegs(jlegs+kk-1)+izone,nz]); % find the max of the ileg including overlap            spike2(round(i1):round(i2))=zeros(size(round(i1):round(i2))); % now set that bit to zero                                                                           % so we ignore it next time around                     end         jlegs=jlegs+length(ii);     end     izone=floor(izone/2); % Now we shrink the smoothing window.      if (izone==0) break; end % If we have no more window, then stop. end nlegs=jlegs; % disp([ int2str(nlegs) ' legs found ']) if( nlegs<nlegsasked & ~iwarn) 	error(' Unable to find requested number of legs'); elseif( nlegs<nlegsasked & iwarn ) 	disp(['WARNING: SONIC2TZ found only ' int2str(nlegs) ' legs where '... 		int2str(nlegsasked) ' were requested']) end %make sure nlegs is exactly nlegsasked if( nlegs > nlegsasked ) 	nlegs=nlegsasked; end 	 nlegs=nlegs-1; ilegs=ilegs(1:round(nlegs));  	 %now the legged tz function i1=find(t(1)==t(ilegs)); i2=find(t(nz)==t(ilegs)); if(isempty(i1) & isempty(i2)) 	tz=[t(1);t(ilegs);t(nz)]; 	zt=[z(1);z(ilegs);z(nz)]; elseif(~isempty(i1)) 	tz=[t(ilegs);t(nz)]; 	zt=[z(ilegs);z(nz)]; elseif(~isempty(i2)) 	tz=[t(1);t(ilegs)]; 	zt=[z(1);z(ilegs)]; end  % sort into order [zt,is]=sort(zt); tz=tz(is);  % the instantaneous velocity % transpose for a real result, 2 for 2 way time vins=2*gradient(zt',tz'); vins=vins(:);  if(m==1) 	tz=tz'; 	zt=zt'; 	vins=vins'; end  if( nargout==0 ) % create analysis plots 	hfig=figure('visible','off','menubar','none'); 	% first the tz curve 	subplot(1,3,1); 	hax1=gca; 	plot(t,z); 	line(tz,zt,'color','r','marker','.','markersize',12); 	title('TZ curves'); 	set(gca,'ylabel',text(0,0,'Depth'),'xlabel',text(0,0,'Time'),'Ydir','reverse'); 	%text(0,-35,'MB1->draw box to zoom. MB1->single click to unzoom',... 	%		'units','pixels'); 	grid; 	%now the instantaneous velocity comparison 	subplot(1,3,2); 	hax2=gca; 	plot(1.e06 ./sonic,z) 	line(vins,zt,'color','r','linewidth',2); 	title('Vins'); 	set(gca,'xlabel',text(0,0,'Velocity'),'ylabel',text(0,0,'Depth'),'ydir','reverse'); 	grid; 	% now the error 	subplot(1,3,3); 	hax3=gca; 	if( flag ) 		tt=t(1):.008:t(length(t)); 		zz=pwlint(tz,zt,tt); 		t=linspace(min(t),max(t),length(t)); 		zzex=interp1(t,z,tt); 		plot(zz(:)-zzex(:),tt) 		grid 		title('Estimated error') 		set(gca,'xlabel',text(0,0,'Depth Error'),'ylabel',text(0,0,'Time'),'ydir','reverse'); 	else 		zz=linspace(z(1),z(length(z)),200); 		tt=pwlint(zt,tz,zz); 		ttex=interp1(z,t,zz); 		plot(tt(:)-ttex(:),zz) 		grid 		title('Estimated error') 		set(gca,'xlabel',text(0,0,'Time Error'),'ylabel',text(0,0,'Depth'),'ydir','reverse'); 	end 	simplezoom; 	tz=[]; 	zt=[]; 	vins=[];  	if(nargin>=5) 		%a figure title 		pos=get(hfig,'position'); 		ht=pos(4); 		Ht=pos(4)+41; 		htit=uicontrol('style','text','string',titlestring,'units','normalized', ... 			'position',[1/pos(3),(pos(4)+21)/Ht,(pos(3)-2)/pos(3),20/Ht]); 		htit2=uicontrol('style','text','string',... 			'Yellow=exact, Red=approx, MB1->draw box to zoom. MB1->single click to unzoom',... 			'units','normalized','position',[1/pos(3),pos(4)/Ht,(pos(3)-2)/pos(3),20/Ht]); 		%shift each axes down 		pa=get(hax1,'position'); 		set(hax1,'position',[pa(1) ht*pa(2)/Ht pa(3) ht*pa(4)/Ht]); 		pa=get(hax2,'position'); 		set(hax2,'position',[pa(1) ht*pa(2)/Ht pa(3) ht*pa(4)/Ht]); 		pa=get(hax3,'position'); 		set(hax3,'position',[pa(1) ht*pa(2)/Ht pa(3) ht*pa(4)/Ht]); 		%resize the figure 		set(hfig,'position',[pos(1:3) Ht]); 	end 	 	% a close button 	uicontrol('style','pushbutton','string','close','position',... 		[1,1,45,20],'callback','close(gcf)'); 		 	set(hfig,'visible','on'); end