add traveltime functions

code/traveltime/#readme.txt

@@ -0,0 +1,21 @@
+dijkstra.m - Dijkstra's algorithm (Shortest path problem)
+getstartmodel.m - Function to create gradient starting model
+grid2mesh.m - Converts FD grid into FE mesh
+gridconst2d.m - Constraint matrix for 2d FD mesh
+meshsmoothness - Constraint matrix for unstructured mesh
+plotshot.m - Plot RA data (and responses)
+readgli.m - read GLI3D data file
+readgrm.m - read GREMIX data file
+readsensvectors.m - read sensitivity matrix by vectors
+showrays.m - show rays starting from one shot position
+tripatchmod.m - show mesh attributes (m/s)
+waymatrix.m - compute way matrix for traveltime calculation
+writepoly.m - write poly file as input for dctriangle
+dctriangle(.exe) - create (para) mesh from poly file
+
+ratomo.fig/.m - GUI for refraction tomography
+Options to be set:
+  Mesh: depth, quality of Mesh, smoothing
+  Start: vmin/vmax, Gradient or 2layer-model
+  Inversion: Parameter (log,v,s), lower/upper bound, smoothness
+             Flatness factor, robust, blocky, ...

code/traveltime/astern.m

@@ -0,0 +1,25 @@
+function [dist,prec]=astern(Mesh,start)
+
+U=1:Mesh.nnodes;
+dist=ones(Mesh.nnodes,1)*999;dist(start)=0;
+prec=zeros(Mesh.nnodes,1);prec(start)=start;
+while ~isempty(U),
+    [umin,iu]=min(dist(U));
+    u=U(iu);U(iu)=[];
+    fi1=find(Mesh.bound(:,1)==u);
+    fi2=find(Mesh.bound(:,2)==u);fi=[fi1;fi2];
+    allv=[Mesh.bound(fi1,2);Mesh.bound(fi2,1)];
+    for iv=1:length(allv),
+        v=allv(iv);
+        if dist(u)+Mesh.edgetime(fi(iv))<dist(v),
+            dist(v)=dist(u)+Mesh.edgetime(fi(iv));
+            prec(v)=u;
+        end
+    end
+end
+% way=[stop prec(stop)];
+% while way(end)~=start,
+%     way=[way prec(way(end))];
+% end
+% way=fliplr(way)
+% time=dist(stop);

code/traveltime/dijkstra.m

@@ -0,0 +1,31 @@
+function [dist,prec]=dijkstra(Mesh,start)
+
+% DIJKSTRA - perform dijskstra routing algorithm on mesh
+% [dist,prec] = dijkstra(Mesh,start)
+
+U=1:Mesh.nnodes;
+dist=ones(Mesh.nnodes,1)*999;dist(start)=0;
+prec=zeros(Mesh.nnodes,1);prec(start)=start;
+while ~isempty(U),
+    [umin,iu]=min(dist(U));
+    u=U(iu);U(iu)=[];
+    fi1=find(Mesh.bound(:,1)==u);
+    fi2=find(Mesh.bound(:,2)==u);fi=[fi1;fi2];
+    allv=[Mesh.bound(fi1,2);Mesh.bound(fi2,1)];
+%     [fi,fi2]=find(Mesh.bound==u);
+%     allv=Mesh.bound(fi,:);allv(allv==u)=[];
+    for iv=1:length(allv),
+        v=allv(iv);
+        newt=dist(u)+Mesh.edgetime(fi(iv));
+        if newt<dist(v),
+            dist(v)=newt;
+            prec(v)=u;
+        end
+    end
+end
+% way=[stop prec(stop)];
+% while way(end)~=start,
+%     way=[way prec(way(end))];
+% end
+% way=fliplr(way)
+% time=dist(stop);

code/traveltime/getshotrec.m

@@ -0,0 +1,27 @@
+function rec = getshotrec(Shot)
+
+% GETSHOTREC - Get shot reciprocity
+% rec = getshotrec(Shot)
+% Shot..shot structure with pos,ns,nx and t
+% rec..relative reciprocity
+
+l=0;
+for i=1:length(Shot.ns),
+    for j=1:length(Shot.nx{i}),
+        l=l+1;
+        nn=sort([Shot.ns{i} Shot.nx{i}(j)]);
+        po(l)=nn(1)*10000+nn(2);
+    end
+end
+rec=zeros(l,1);
+l=0;
+for i=1:length(Shot.ns),
+    for j=1:length(Shot.nx{i}),
+        l=l+1;
+        fi=find(po==po(l));
+        if any(fi), 
+            tt=Shot.t(fi);
+            rec(l)=std(tt)/mean(tt); 
+        end
+    end
+end

code/traveltime/getstartmodel.m

@@ -0,0 +1,36 @@
+function velocity=getstartmodel(Mesh,vmin,vmax,zz)
+
+% GETSTARTMODEL - Get Start model
+% velocity=getstartmodel(Mesh,v) - uniform
+% velocity=getstartmodel(Mesh,vmin,vmax) - gradient model
+% velocity=getstartmodel(Mesh,vmin,vmax,zlayer) - 2-layered
+% velocity=getstartmodel(Mesh,Shot) - determine from data
+
+if nargin<4, zz=0; end
+if nargin<3, vmax=0; end
+if nargin<2, error('Must specify mesh and Shot or min/max!'); end
+if (~isstruct(Mesh))||(~isfield(Mesh,'node'))||(~isfield(Mesh,'cell')),
+    error('1st argument must be a valid Mesh!'); end
+if isstruct(vmin), %automatic detection
+    error('not yet implemented!');
+    %vmin=
+    %vmax=
+end
+% vmin and vmax specified
+velocity=ones(Mesh.ncells,1)*vmin;
+if vmax>0, % not uniform
+    cellmids=zeros(Mesh.ncells,Mesh.dim);
+    for i=1:Mesh.ncells, cellmids(i,:)=mean(Mesh.node(Mesh.cell(i,:),:)); end
+    A=ones(Mesh.ncells,2);A(:,2)=cellmids(:,1);
+    ab=A\cellmids(:,Mesh.dim); % remove linear trend by regression
+    cellmids(:,Mesh.dim)=cellmids(:,Mesh.dim)-ab(1)-ab(2)*cellmids(:,1);
+    if zz>0, % 2-layered case
+        velocity(abs(cellmids(:,Mesh.dim))>zz)=vmax;
+    else, % gradient model
+%         deprel=abs(cellmids(:,Mesh.dim)-max(Mesh.node(:,Mesh.dim)));
+        deprel=abs(cellmids(:,Mesh.dim)-max(cellmids(:,Mesh.dim)));
+        deprel=deprel/max(deprel);
+        velocity=(vmax/vmin).^deprel*vmin;
+        % velocity=vmin+vmax*abs(cellmids(:,2))/max(abs(cellmids(:,2)));
+    end
+end

code/traveltime/getva.m

@@ -0,0 +1,20 @@
+function [va,di]=getva(Shot)
+
+% GETVA - Get apparent velocity of a Shot struct
+% va = getva(Shot)
+
+di=[];
+if isfield(Shot,'s')&&isfield(Shot,'g'),
+    if isfield(Shot,'pos'),
+        di=sqrt(sum((Shot.pos(Shot.s,:)-Shot.pos(Shot.g,:)).^2,2));
+    elseif isfield(Shot,'elec'),
+        di=sqrt(sum((Shot.elec(Shot.s,:)-Shot.elec(Shot.g,:)).^2,2));
+    else di=1;
+    end
+else
+    for i=1:length(Shot.ns),
+        rel=Shot.pos(Shot.nx{i},:)-Shot.pos(repmat(Shot.ns{i},length(Shot.nx{i}),1),:);
+        di=[di;sqrt(sum(rel.^2,2))];
+    end
+end
+va=di./Shot.t;

code/traveltime/gridconst2d.m

@@ -0,0 +1,35 @@
+function C = gridconst2d(M)
+
+% GRIDCONST2D - Grid constraint matrix
+% C = gridconst(M)
+
+% M=ones(3,2);
+alphaz=0.3;
+nx=size(M,1);nz=size(M,2);
+hor=(nx-1)*nz;ver=nx*(nz-1);
+% C=spalloc(hor+ver,nx*nz,(hor+ver)*2);
+l=0;
+mm=[];nn=[];vv=[];
+ox=ones(nx-1,1);oz=ones(nz-1,1)*alphaz;
+ii=(1:nx-1)';kk=(1:nz-1)';
+for k=1:nz,
+    ik=(k-1)*nx+ii;
+    mm=[mm;l+ii];
+    nn=[nn;ik];
+    vv=[vv;ox];
+    mm=[mm;l+ii];
+    nn=[nn;ik+1];
+    vv=[vv;-ox];
+    l=l+nx-1;
+end
+for i=1:nx,
+    ik=(kk-1)*nx+i;
+    mm=[mm;l+kk];
+    nn=[nn;ik];
+    vv=[vv;oz];
+    mm=[mm;l+kk];
+    nn=[nn;ik+nx];
+    vv=[vv;-oz];
+    l=l+nz-1;
+end
+C=sparse(mm,nn,vv);

code/traveltime/korrrad.m

@@ -0,0 +1,22 @@
+function Shot1=korrad(Shot)
+
+% KORRAD - corrigate Shot positions to ly on radius if feasible
+% Shot_new = korrad(Shot)
+
+mid=mean(Shot.pos);
+if max(abs(mid))<0.05, mid=zeros(1,size(Shot.pos,2)); end
+rel=Shot.pos-repmat(mid,size(Shot.pos,1),1);
+rad=sqrt(sum(rel(:,1).^2+rel(:,2).^2,2));
+ang=atan2(rel(:,2),rel(:,1))*180/pi;
+Shot1=Shot;
+if min(rad)/max(rad)>0.99, % on circle
+   ra=rndig(mean(rad),2);
+   da=2.5;
+   [ang,ii,jj]=unique(round(ang/da)*da);
+   Shot1.pos=[ra*cos(ang*pi/180) ra*sin(ang*pi/180)];
+   Shot1.pos=Shot1.pos+repmat(mid,size(Shot1.pos,1),1);
+   for i=1:length(Shot.ns),
+      Shot1.ns{i}=jj(Shot.ns{i});
+      Shot1.nx{i}=jj(Shot.nx{i});
+   end
+end

code/traveltime/linesearch2.m

@@ -0,0 +1,14 @@
+function [tauopt,appnew]=linesearch(soll,old,new)
+
+% LINESEARCH - Applies line-search procedure
+% tauopt = linesearch(old,new)
+
+for i=1:40,
+    tau=i*0.05;
+    inbet=old+tau*(new-old);
+    ch(i)=sqrt(mean((soll-inbet).^2));
+%     ch(i)=rms(soll,inbet);
+end
+[xx,nn]=min(ch);
+tauopt=nn*0.05;    
+appnew=old+tauopt*(new-old);

code/traveltime/mknewshot.m

@@ -0,0 +1,31 @@
+function Shot=mknewshot(npos,dshot,dx)
+
+% MKNEWSHOT - Make new shot
+% Shot = mknewshot(npos,dshot,dx)
+%        pos  - number of shot/geophone positions
+%               or positions itself
+%        dshot - shot every dshot positions
+%        dx    - geophone distance
+
+if nargin<1, npos=41; end % default values
+if nargin<2, dshot=5; end
+if nargin<3, dx=1; end
+if length(npos)>1, % positions given
+    Shot.pos=npos;
+    npos=size(Shot.pos,1);
+else % number of 
+    Shot.pos=(0:npos-1)';Shot.pos(1,2)=0;
+end
+Shot.t=[];
+ishot=0;nshot=1-dshot; % counters
+while(nshot<npos), % shot while not at end
+    ishot=ishot+1;nshot=nshot+dshot;
+    if nshot>npos, nshot=npos; end
+    Shot.ns{ishot}=nshot;
+    aa=(1:npos)';aa(nshot)=[];Shot.nx{ishot}=aa;
+    Shot.loc(ishot)=Shot.pos(nshot,1); %historical
+    Shot.x{ishot}=Shot.pos(aa); % only for plotting
+    Shot.tt{ishot}=abs(aa-nshot);
+    Shot.t=[Shot.t;Shot.tt{ishot}/1000];
+end
+if nargout<1, plotshot(Shot); end

code/traveltime/mksinkhole1.m

@@ -0,0 +1,80 @@
+xmax=300; % profile length
+zmax=100; % model depth
+dg=5; % geophone spreading
+dsdg=2; % 2 shots per geophone
+z1=10;
+Poly=[];
+Poly.node=[0 0;xmax 0;xmax z1;xmax zmax;0 zmax;0 z1];
+Poly.node(:,3)=0;
+Poly.edge=[1 2;2 3;3 4;4 5;5 6;6 1;6 3];
+Poly.edge(:,3)=0;
+Poly.region=[1 1 1;1 zmax-1 2];
+%%
+pos=50;len=20;dlen2=5;dep=10;ln=size(Poly.node,1);
+Poly.node(end+1:end+4,1:2)=[pos z1;pos+len z1;pos+len-dlen2 z1+dep;pos+dlen2 z1+dep];
+Poly.edge(end,2)=ln+1;
+Poly.edge(end+1:end+5,1:2)=ln+[1 2;2 3;3 4;4 1;2 0];
+Poly.edge(end,2)=3;
+Poly.region(end+1,:)=[pos+len/2 z1+1 size(Poly.region,1)+1];
+%%
+pos=120;len=50;dlen2=10;dep=20;ln=size(Poly.node,1);
+Poly.node(end+1:end+4,1:2)=[pos z1;pos+len z1;pos+len-dlen2 z1+dep;pos+dlen2 z1+dep];
+Poly.edge(end,2)=ln+1;
+Poly.edge(end+1:end+5,1:2)=ln+[1 2;2 3;3 4;4 1;2 0];
+Poly.edge(end,2)=3;
+Poly.region(end+1,:)=[pos+len/2 z1+1 size(Poly.region,1)+1];
+%%
+pos=220;len=10;dlen2=2;dep=5;ln=size(Poly.node,1);
+Poly.node(end+1:end+4,1:2)=[pos z1;pos+len z1;pos+len-dlen2 z1+dep;pos+dlen2 z1+dep];
+Poly.edge(end,2)=ln+1;
+Poly.edge(end+1:end+5,1:2)=ln+[1 2;2 3;3 4;4 1;2 0];
+Poly.edge(end,2)=3;
+Poly.region(end+1,:)=[pos+len/2 z1+1 size(Poly.region,1)+1];
+show2dpoly(Poly);
+%%
+Poly.node(1,3)=-99;
+dg=5;gpos=(Poly.node(1,1)+dg:dg:Poly.node(2,1)-dg)';
+% Poly.node(2,3)=-99;
+gpos(:,2)=0;gpos(:,3)=-99;Poly.node=[Poly.node;gpos];
+Poly.node(:,2)=-Poly.node(:,2);Poly.region(:,2)=-Poly.region(:,2);
+show2dpoly(Poly);
+%%
+epos=gpos;epos(:,2)=epos(:,2)-1;epos(:,3)=0;
+Poly1=Poly;Poly1.node=[Poly.node;epos];
+if 1, return; end
+%%
+writepoly2d('sinkhole1.poly',Poly1);
+system('dctriangle -v -q34 sinkhole1.poly');
+Mesh=loadmesh('sinkhole1.bms');
+tripatchmod(Mesh)
+%%
+Shot=[];Shot.pos=Mesh.node(Mesh.nodemarker==-99,1:2);
+Shot.ns={};Shot.nx={};Shot.tt={};Shot.t=[];
+ng=sum(Mesh.nodemarker==-99);
+all=(1:ng)';
+i=1;
+while i<ng,    
+    Shot.ns{end+1}=i;
+    Shot.nx{end+1}=all;
+    Shot.nx{end}(Shot.nx{end}==i)=[];    
+    Shot.tt{end+1}=zeros(size(Shot.nx{end}));
+    Shot.t=[Shot.t;Shot.tt{end}];
+    i=i+2;
+end
+savesgtfile('synth1.sgt',Shot);
+%%
+veli=[500 1800 1000 1000 1000]';
+% tripatchmod(Mesh,veli(Mesh.cellattr));
+map=(1:length(veli))';map(:,2)=veli;
+save sinkhole1.map map -ascii
+system('ttmod -vEV -t1e-3 -a sinkhole1.map -p sinkhole1.bms -o sinkhole1.sgt synth1.sgt');
+Res=readunishot('sinkhole1.sgt');
+plotshot(Res);minmax(getva(Res))
+if 1, return; end
+%%
+% [W,t]=waymatrix(Mesh,Res,veli(Mesh.cellattr));
+% plotshot(Res,t);
+%% inversion
+MM=Mesh;MM.cellattr(:)=2;savemesh(MM,'para1.bms');
+system('ttinv -vG -z0.3 -p para1.bms -t 1e-3 sinkhole1.sgt');
+vel=load('velocity.vec');tripatchmod(MM,vel)