Cross-well travel time tomography

From the article
``Geophysical Inverse Problems''

in
Mathematica in Education
Vol.3 No.2
Spring 1994
(C) TELOS/Springer-Verlag

by
John Scales
Department of Geophysics
Colorado School of Mines
Golden, Colorado 80401-1887
jscales@dix.mines.colorado.edu

Bug fix: 2.0 -> 2.2

If you are running Version 2.2 of Mathematica, you will have to evaluate the following code in order for the Transpose function to work properly with the multi-dimensional lists in this notebook.

Unprotect[Transpose];
Transpose[v_?VectorQ] := v
Protect[Transpose]

Set up basic arrays for sources/receivers

nsources = 25;
nreceivers = 25;
nx = 20;
nz = 20;
xmax = 10.;
zmax = 10.;
ncells = nx nz;
nrays = nsources nreceivers;

dr = zmax/(nreceivers - 1);
ds = zmax/(nsources - 1);
zsrc = Table[(i-1) ds,{i,nsources}];
zrec = Table[(i-1) dr,{i,nreceivers}];
xsrc = Table[0, {i,nsources}];
xrec = Table[xmax, {i,nreceivers}];

nrays = nsources nreceivers;
ncells = nx nz;
dx = xmax/nx;
dz = zmax/nz;
gridx = Table[(i-1) dx,{i,nx+1}];
gridz = Table[(i-1) dz,{i,nz+1}];

izsrc = Table[Min[Floor[zsrc[[i]]/dz]+1,nz],{i,nsources}];
ixsrc = Table[Min[Floor[xsrc[[i]]/dx]+1,nx],{i,nsources}];
izrec = Table[Min[Floor[zrec[[i]]/dz]+1,nz],{i,nreceivers}];
ixrec = Table[Min[Floor[xrec[[i]]/dx]+1,nx],{i,nreceivers}];

Compute the Jacobian matrix.

This form of the Jacobian calculation is due to Lydia Deng

raynum = 1;
Jac = Table[0.,{i,nrays},{j,ncells}];

Do[ (* loop over sources *)
 ixs = ixsrc[[i]]; izs = izsrc[[i]];
 Do[ (* loop over receivers *)
         delz = zrec[[j]]-zsrc[[i]];
     delx = xrec[[j]]-xsrc[[i]];
     tang=delz/delx;
     ix = ixs; iz = izs;
     xbeg = xsrc[[i]]; zbeg = zsrc[[i]];
                While[ix<=nx && iz<=nz && iz>0,
                cellnum = ix + (iz -1) nx;
                        (*Print["ix=",ix,"  iz=",iz," ncell=",cellnum];
                        Print["nray=",raynum," ncell= ",cellnum];*)
                        ix=ix+1; xend=gridx[[ix]];
                        zend = zbeg + tang (xend-xbeg);
                        Which [delz>0,
                                Which[
                                 zend>gridz[[iz+1]],
                                 zend=gridz[[iz+1]];
                                 xend = xbeg + (zend-zbeg)/tang;
                                 ix=ix-1; iz=iz+1,
                                 zend == gridz[[iz+1]],
                                 iz=iz+1;
                                ],
                                delz<0,
                                Which[ zend<gridz[[iz]],
                                 zend=gridz[[iz]];
                                 xend = xbeg + (zend-zbeg)/tang;
                                 ix=ix-1; iz=iz-1,
                                 zend == gridz[[iz]],
                                 iz=iz-1;
                                ];
                        ];
                        unit = Sqrt[(xend-xbeg)^2+(zbeg-zend)^2];
                        xbeg = xend; zbeg = zend;
                        Jac[[raynum,cellnum]]=Jac[[raynum,cellnum]]+unit;
                ];
        raynum = raynum + 1;
  ,{j,nreceivers}];
,{i,nsources}];

This is the Jacobian matrix (Frechet derivative).  It mainly consists of zeros.

ListDensityPlot[Jac, 
	Mesh -> False,
	ColorFunction :> (Hue[# + .67]&)]

Show[%, ColorFunction :> (RGBColor[0,0,#]&),
	AspectRatio -> 532/416,
	Axes -> False,
	Frame -> False];

Sum the jacobian over rays for each cell.  This tells us how many rays pierce each cell.  This is called the "hit count" and is a crude measure of resolution.  Clearly cells which have only a few rays in there are unlikely to be very well resolved.  ListPlot3D is a nice option.  But it's much slower if you have lots of pixels.

hitcount = Table[Sum[Jac[[i,j]],{i,nrays}],{j,ncells}];
ListDensityPlot[Partition[hitcount,nx],PlotRange->All]

A simple example of travel time tomography

Now we'll generate a veocity model whose features we'll recover from travel times.  The background model we'll choose is identical but for the cross shaped perturbation in the center.

back = Table[1.,{i,nx},{j,nz}];
tmp = back;
Do[
	Do[
		If[((i-10)^2 + (j-10)^2) < 5 , tmp[[i,j]] = 2];
			,{i,nx}];
				,{j,nz}];
perturbed = Flatten[tmp];
backgr = Flatten[back];

ListPlot3D[tmp]

Compute the travel times through the true and background model.  The background model is supposed to be our first guess in an iterative reconstruction.  Since this problem is, by construction, completely linear, we only do one iteration.

ttback = Table[ Sum[Jac[[i,k]]/backgr[[k]],{k,ncells} ],{i,nrays}];
tttrue = Table[ Sum[Jac[[i,k]]/perturbed[[k]],{k,ncells} ],{i,nrays}];
ttdelta = tttrue - ttback;

ListPlot[tttrue,PlotJoined->True]

Look in the singular value reconstruction of the solution in detail

First, do the SVD, then reconstruct the solution piece by piece.

{v,mm,u} = SingularValues[Jac];

svalues = 10;
modelpert =  Sum[u[[i]]( (v[[i]].ttdelta) /mm[[i]]),
		{i,svalues}];
updatedmodel = 1./((1./backgr) + modelpert);
ListPlot3D[Partition[updatedmodel,nx],
	PlotRange->All]

svalues = 50;
modelpert =  Sum[u[[i]]( (v[[i]].ttdelta) /mm[[i]]),
		{i,svalues}];
updatedmodel = 1./((1./backgr) + modelpert);
ListPlot3D[Partition[updatedmodel,nx],
		PlotRange->All]

svalues = Length[mm];
modelpert =  Sum[u[[i]]( (v[[i]].ttdelta) /mm[[i]]),
		{i,svalues}];
updatedmodelsvd = 1./((1./backgr) + modelpert);
ListPlot3D[Partition[updatedmodelsvd,nx],
		PlotRange->All]

We could plot each model singular vector a list

ListPlot[u[[1]],PlotJoined->True]

But each of these singular vectors is a legitimate model,
sort of an elemental model out of which all others are built,
like the sines and cosines of a Fourier series, or the normal 
modes of a mechanical system.

Plot the two most well resolved and poorly resolved eigenvectors

nd = Length[mm];
svalplot = ListPlot[mm,PlotJoined->True]
dataplot = ListPlot[datacomp,PlotJoined->True,PlotRange->All]
ListPlot3D[Partition[u[[1]],nx],PlotRange->All]
ListPlot3D[Partition[u[[2]],nx],PlotRange->All]
ListPlot3D[Partition[u[[nd - 1]],nx],PlotRange->All]
ListPlot3D[Partition[u[[nd]],nx],PlotRange->All]

Next we will compare the results of the direct SVD
calculation with an iterative method based on conjugate
gradients.

(* This one took 4748.28 seconds on a NeXT turbo with 48 MB ram. *)

nits = 150;
x = Table[0,{i,ncells}];
s = ttdelta;
p = Transpose[Jac].s;
rold = p;
q = Jac.p;
Do[
	alpha = (rold.rold)/(q.q);
	x = x + alpha p;
	s = s - alpha q;
	r = Transpose[Jac].s;
	beta = (r.r)/(rold.rold);
	rold = r;
	p = r + beta p;
	q = Jac.p;
,{i,nits}];

updatedmodelcg = 1./((1./backgr) + x);
ListPlot3D[Partition[updatedmodelcg,nx],PlotRange->All]

Now look at the difference between the SVD and the CG solutions. 
We could make this smaller by doing more CG iterations, but note
that this residual looks a lot like the singular vector associated
with the small singular value.

ListPlot3D[Partition[ (updatedmodelsvd - updatedmodelcg),nx]]

Finally, we can look at the scaling factors that weight each model
eigenvector in the SVD reconstruction.  Each factor is proportional
to the component of the data vector in each data eigenvector, and
inversely proportional to the associated singular value.

nd = Length[mm];
datacomp = Table[v[[i]].ttdelta, {i,nd}];
svalplot = ListPlot[mm,PlotJoined->True]
dataplot = ListPlot[datacomp,PlotJoined->True,
	PlotRange->All]
ListPlot[datacomp/mm,PlotJoined->True,PlotRange->All]

Scales Figure 13

Here is a simple example with a 2 by 2 diagonal matrix.  The condition number is 10.  First, steepest descent.  This implementation is due to Herman Jaramillo.  It's a little more complicated than  required by the basic iteration, but that's so we can keep track of the path of the solution for the pictures.

nits = 19;
R = Table[{0,0},{i,nits+1}];
X = Table[{0,0},{i,nits+1}];
lam1 = 10.;
lam2 = 1.;
A = {{lam1,0},{0,lam2}};
x = {0.,0.};
h = {1.,-1.};

quad[x1_,x2_] := {x1,x2} . (A . {x1,x2})/2 - h . {x1,x2};
cont = {quad[x[[1]],x[[2]]]};
r = h - A . x;
R[[1]] = r;
X[[1]] = x;
err = Sqrt[r.r];
its = 0;

While[its < nits,
	{
		a = (r . r)/(r . (A . r));
		x = x + a r;
		cont = Join[cont,{quad[x[[1]],x[[2]]]}];
		X[[its + 2]] = x;
		r = h - A . x;
		R[[its + 2]] = r;
		err = Sqrt[r.r];
		its = its + 1;}
	];

plot1 = 
ListPlot[X, PlotJoined -> True,
	AspectRatio -> Automatic,
	Axes -> False,
	DisplayFunction -> Identity]

plot2 = 
ContourPlot[quad[x1,x2], 
	{x1,-.25,.45}, {x2,-2.2,0.2},
	ContourShading -> False,
	AspectRatio -> Automatic, Contours -> cont,
	PlotPoints -> 50, FrameLabel -> None,
	DisplayFunction -> Identity]

Show[plot1, plot2,
	DisplayFunction -> $DisplayFunction]

You can readily convince yourself that the only way steepest  descent can terminate in a finite number of steps is if the starting point lies along one of the principle axes of the quadratic form.

Scales Figure 15

Now let's see what conjugate gradient does.

nits = 3;
R = Table[{0,0},{i,nits+1}];
X = Table[{0,0},{i,nits+1}];
lam1 = 10.;
lam2 = 1.;
A = {{lam1,0},{0,lam2}};
x = {0.,0.};
h = {1.,-1.};

quad[x1_,x2_] := {x1,x2} . (A . {x1,x2})/2 - h . {x1,x2};
cont = {quad[x[[1]],x[[2]]]};
r = h - A . x;
R[[1]] = r;
X[[1]] = x;
err = Sqrt[r.r];
its = 0;

p = r;
rold = r;
While[its < nits,
        {
                a = (rold . rold)/(p . (A . p));
                x = x + a p;
                cont = Join[cont,{quad[x[[1]],x[[2]]]}];
                X[[its + 2]] = x;
                r = rold - a (A . p);
                b = (r . r)/(rold . rold);
                R[[its + 2]] = r;
                rold = r;
                p = r + b p;
                err = Sqrt[r.r];
                its = its + 1;}
        ];

plot3 =
ListPlot[X, PlotJoined -> True,
	AspectRatio -> Automatic,
	Axes -> False,
	DisplayFunction -> Identity]

Show[plot3, plot2,
	DisplayFunction -> $DisplayFunction]

Finally, ask yourself:  what is implied by finite termination?  In other words, under what circumstances is the residual at the i-th step equal to zero for CG?