Straight Ray Tomography Example

OBJECTIVE: Solve Traveltime Tomography Ls=t by s=[L^TL]^-1L^Tt.

BACKGROUND: Straight ray tomography does not take into account ray bending but can provide a good quick first velocity model. Code below is quite simple and does a transmission tomography experiment where sources are along top boundary and receivers along bottom boundary.

PROCEDURE:

1. Download all matlab codes straight_ray.zip . It includes "timetomo.m", "coord.m", "raymat.m", "merge.m", and "raypath.m" .

2. Traveltime Matrix SLL: Drag and drop the following code that builds up the raypath matrix SLL for a crosswell transmission traveltime geomtery.

   clear all;hold off;
   nx=20;nz=nx;dx=10;dxx=dx;firstsx=0;ds=10;ns=20;firstgx=0;dg=10;ng=20;
   [SLL,LL]=raypath(nx,nz,dx,firstsx,ds,ns,firstgx,dg,ng);alpha=.01;
   %
   % Create a Vertical Layer Slowness Model
   %
   
   slow=ones(nx,nz)/2000;
   slow(round(nx/2):nx,:)=slow(round(nx/2):nx,:)*.5;
   nxm=nx-1;slow1(1:nxm,1:nxm)=slow(1:nxm,1:nxm);
   %slow1=slow1';%                                  Model becomes 2 layer horizontal model
   %nst=round(nx/2);nen=round(nst*1.3);%          Add a slowness anomaly
   %for i=nst:nen;slow1(i,nst:nen)=1/2000;end%     Add a slowness anomaly
   subplot(211);imagesc([1:nx-1]*dxx,[1:nx-1]*dxx,1./slow1');
   colorbar;ylabel('Z (m)');title('Velocity Model');
   text(dxx*nx,-1.5*dxx,'Velocity (km/s)')
   

   Which index is associated with raypath number and which index is associated with slowness inknown? What is dimension of SLL (Use size(SLL)). Type SLL(3,:) to see the raypath segment lengths of 3rd ray. Is this an overdetermined system of equations? Is the associated normal equation matrix well conditioned? What is the conditon number of the nromal equations?

3. Regularization: Now solve the system of equations by dragging and dropping the following code:

   
   Solve for least square solution of Ls=t;
   %
   TIME=SLL*slow1(:);LTL=SLL'*SLL;LTt=SLL'*TIME;
   d=diag(LTL);D=diag(1./d,0);ID=eye(length(d))*alpha;   % alpha is damping parameter
   appslow=inv(D*LTL+ID)*D*LTt;
   

   The resulting tomogram can be plotted by typing

   %
   % Plot Tomogram
   %
   appslow=reshape(appslow,nx-1,nx-1);subplot(212);
   imagesc([1:nx-1]*dxx,[1:nx-1]*dxx,1./appslow');
   colorbar;xlabel('X (m)');ylabel('Z (m)');title('Slowness Tomogram');
   text(dxx*nx,-dxx,'Velocity (m/s)')
   
   

   To understand the importance of damping repeat above except set the damping paramater alpha =0. Explin this new result using the word condition number.

4. Model Covariance Matrix: Type timetomo, and it will first give you a movie of rays, where shot is along top surface and receivers is along bottom boundary. The model is a vertical layer model. Extract the diagonal components of the [L^TL]^-1 matrix by typing "di=diag(inv(SLL'*SLL));". Plot them up by typing imagesc(reshape(di,nx-1,nx-1)');colorbar. Make sure that you add a diagnoal matrix with damping parameter of .1 to LTL matrix before you invert. This will lessen singularity of matrix so covariance values are better displayed. This is the plot of the variance terms of the slowness model. Do the high variance values correspond to poorly resolved parts of model? Why? Which cells have rays with widest variability of ray angles? Are they best resolved according to the covariance matrix?

5. Models in Null Space: Adjust calculations by changing model to a horizontal layer model. Model can be changed to a 2-layer model by uncommenting the 2-layer model line in timetomo.m.

   Does the new covariance matrix vlaues predict areas of poor resolution? Make sure that you add a diagnoal matrix with damping parameter of .1 to LTL matrix before you invert. This will lessen singularity of matrix so covariance values are better displayed.

6. Adjust calculations for an anomaly model. An anomaly model can be included by uncommenting the 2 anomaly lines in timetomo.m.

7. Explain why vertical layer model is well resolved but horizontal layer model is not. Comment about resolvability of anomaly model.

   (Extra Credit) Try different damping parameters and constraints to see if the result improves.

8. (Extra Credit) Determine optimal Marquardt damping parameter by creating a plot similar to that in Figure 2.4 in here.