I assume here that you have advanced in your GEOP 205 to  where you are familiar with some of the basic solutions of LIP (Linear Inverse Problem). In this page I present some simple but important examples of LIP that illustrate some of the important features of the different inverse problems encountered in geophysical interpretations work.

1. The Dampened Least Squares Problem:

The objective of this example is to illustrate the computational mechanics of this important type pf solution to the inverse problem. Recall that the damped least squares solution is derived under a combination of criteria:

(1) Minimum prediction error

(2) Minimum solution length:

The objective function to minimize, therefore is:

Which yields the solution:

We generate a set of synthetic data from a quadratic model:  y = a + b x + c x²  .

Here is the MATLAB^® script to do so:

Figure (1) to the right shows our model together with the noisy observations.  Our objective is to attempt to recover estimates of the model parameters m = [10,  3.6,  -2.89] from the noisy observational data d^obs .

Here is how we may form the design matrix G:

G = [X.^{^}0    X     X.^{^}2];

The first 4-rows of G:

We want to solve the problem for several values of the dampening parameter m and save the results of each solution for examination and comparison. We do so via a for-loop as follows:

The results are shown in Figure (2),  which illustrates the trade-off between minimum prediction error E and minimum length solution L as the damping factor m increase towards unity.

2. Delay Time in Seismic Refraction Surveys   

The concept of delay time is important in the interpretation of refraction seismic data, especially in the  Hagedoorn plus—minus method of interpretation. In this method, it is assumed that the layers present are homogenous, that there is a large velocity contrast between the layers, and that the angle of dip of the refractor is less than 10^o. The method uses intercept times and delay times in the calculation of the depth to a refractor below any geophone location. Referring to Figure (1), the delay time (Dt) is the difference in time between:

1. the time taken for a ray to travel along a critically refracted path from the source via the sub-surface media and back to the surface (T_ABCD along the path ABCD);

2. the time taken for a ray to travel the equivalent distance of the source-geophone offset (x) along the refractor surface (i.e. T_PQ along the projection PQ of the refracted ray path onto the refractor surface).

The total delay time Dt is effectively the sum of the ‘shot-point delay time’ (Dt_s) and the ‘geophone delay time Dt_g. In cases where the assumptions on which the method is based are realistic (especially the refractor dip is not too large), it is sufficiently accurate to consider the distance PQ to be approximately the same as the source-geophone offset X. The mathematical equations relating the various delay times and velocities are then as follows:

The total delay time is:

Consider now the refraction delay time experiment shown in the figure to the right, where two shots are received at three stations each. Assuming a two-layer Earth model with velocities V1 > V2 for the top and bottom layers respectively, then the travel-time between shot (i) and receiver (j) is:

Where S2 is the slowness of the refractor (= 1/V2) and   are the delay times associated with the shot points and receivers respectively.  The matrix formulation of this problem is of the form:

d = G m

Note the special structure of the G-matrix of this problem. Many of its elements are zero and ones. This is called a sparce matrix and is of special importance in geophysical computations.

3. Travel Time Seismic Tomography   
Tomography is the process of obtaining a 2 and 3D images of a slice through the interior of a 3-dimensional object by passing rays of some sort through the object and observing how the rays are affected (attenuated or delayed in time) by passage through the material.

There are applications of tomographic imaging in many different fields. In medicine, methods include X-ray tomography, ultrasonic tomography, magnetic resonance imaging, and positron emission tomography (PET). In the earth sciences, tomography is widely used  In seismic and EM imaging. In cross-well seismic tomography, for instance, acoustic waves are generated within one well, and received by sensors in a second well. In cross well radar tomography, ground penetrating radar waves (at frequencies of hundreds of MHZ) are used. In seismic tomography of the deep earth, seismic waves from naturally occurring earth quakes are sensed at stations on the surface of the earth. Tomography techniques are also used in the analysis of controlled source seismic data generated by reflection of seismic waves from geological structures.

Tomography depends on the assumption that energy from the transmitter follows ray paths through the 3-dimensional body to receivers. Under this assumption we can evaluate path integrals along ray paths to determine travel times and attenuations. Of course, computing the travel time or attenuation from the known physical parameters of the 3-dimensional problem is a forward problem. What we want is to determine the physical parameters which result in the observed travel times or attenuations.

In the simplest case, ray paths follow straight lines with essentially no refraction (bending rays) or diffraction (scattering from small features). The problem becomes nonlinear and much more difficult as we consider refraction effects which result in curved ray paths, and especially when we consider diffraction effects. The straight ray tomography problem, however, can be formulated as a linear inverse problem that can be addressed using the techniques that we introduce in GEOP 205.

which leads to a set of linear equations of the form G m = d.

Now consider a simple tomographic experiment in which we discretize the object into a 4 x 4 blocks. If we use 4 E-W and 4 N-S rays,  we get a rough image of the target; more details may be obtained if additional diagonal scans are performed. Here is how the equations for the first row and first column scans look like:

For the diagonal scans, the equations are not as easy to code; and things get even messier when considering radial scans. Here is the equation for the main diagonal scan:

I am sure that with this equation and the figure above, you can write down the equations for all diagonal scans of the image.

Now that we know the basics of tomography, let us setup a simple tomography problem and try to work out its solution. The objective of this exercise is to illustrate again that important Geophysical Inverse Problems are not easily solvable. They have innate characteristics of being under-determined,  ill-posed and unstable. Hence, they require special algorithms to handle them properly.

We select for our problem the target image of the checkerboard shown in the figure to the right.  This is the matrix of slow nesses of the various blocks that make up the object. We have intentionally selected a situation of high velocity contrast (~ 100 m/sec) between adjacent blocks just to make the problem somewhat simpler to solve.

Here are the steps to follow to generate synthetic data and set up the problem:

1. Construct the design matrix G;

2. Using a set of 16 velocity values, convert to slowness s_t (call this the true slowness vector);

3. Compute synthetic arrival-times data d_t (true data) from:      d_t = G s_t

Here is how your [22 x 16] G-matrix  should look like assuming you have followed the scan sequence:

Row => Column => NE-SW => NW-SE. This sequence determines the arrangements of your equations and consequently the structure of the G-matrix.

Colored matrix elements (pixels) contain values of 1 or √(2); all others (white background) contain zeros. This indeed is a special type of matrix that we may expect to give us some problems when attempting to invert it. How would we know the quality of G^-1 ? Find the condition number, or rank or the SVD of G. The appropriate MATLAB commands return:

cond(G) = 9.7892e+015                    rank(G) = 15  .. .. ooops! i.e. rank deficient (= singular).

Let us look at the SVD of G; the command  [U S V] = svd(G) returns:

diag(S) = [...., 244.95,  244.95,  218.84,  200.00,  200.00,  200.00, 141.42,  141.42,  116.33,  000.00 ]

So, the 16^th singular value is ZERO. This says that we can not expected to recover the image from direct inversion or simple Least Square solution. We have to try other techniques that will allow inversion to proceed without destabilizing the computation to a disastrous end! There several such techniques developed specially for tomographic imaging of which are the so called ART, KAC and SIRT methods. These inversion strategies were developed specially for seismic tomography and subsurface imaging. In the figure below are the solutions obtained by direct (un-regularized)  and by Damped Least Squares. You can see for yourself what works best!

Top of Page   

If you find this problem interesting and/or challenging, here is a collection of MATLAB  functions TOMO that will generate the G-matrix for the 4 x 4 image and the horizontal, vertical and diagonal scans.  It also includes solver functions of all the three methods mentioned above (ART, KAC and SIRT) and an additional method (not discussed here) called TOTAL Least Squares. Experiment with them under different conditions of random errors, scan sequence, etc.

Have FUN