sergey@sep.stanford.edu, jon@sep.stanford.edu

As the former time made light the land of Zebulun and the land of Naphtali, so the latter hath honoured the way of the sea, beyond the Jordan, Galilee of the nations.Isaiah 9:1

## ABSTRACTWe applied the inverse linear interpolation method to process a bottom sounding survey data set from the Sea of Galilee in Israel. Non-Gaussian behavior of the noise led us to employ a version of the iteratively reweighted least squares (IRLS) technique. The IRLS enhancement of the method was able to remove the image artifacts caused by the noise at the cost of a loss in the image resolution. Untested alternatives leave room for further research. |

Figure 1

The Sea of Galilee, also known in Israel as Lake Kinneret (Figure galgol), is a unique natural object because of both its ancient history and its geography. It is a rare case of a freshwater lake below sea level. From a geological point of view, the location of the lake is connected to a great rift crossing east Africa. Our project addressed processing of the bottom sounding survey data collected by Zvi Ben Abraham at Tel-Aviv University .

The original data include over 100,000 triples (*x*_{i},*y*_{i},*z*_{i}) . The range of
*x*'s is about 12 km, and the range of *y*'s is about 20 km.
The depth coordinates
range from 211 m to 257 m below sea level.
Consecutive measurements are
close to each other and cover the area of the lake with a fairly even grid.
However, a number of empty bins remains even after sparse binning (Figure
galsea). Using primitive preprocessing editing, we got rid of several
measurements containing zeroes or evident mispositioning errors.

Figure 2

The major goal of this project was to interpolate the data to a regular grid. Properly interpolated data can enable searching for the local features of interest, including geological structures as well as sunken ships and other exotic archeological objects. Our attention was focused on developing a data interpolation technique that might be useful in a wide variety of geophysical applications.

The general method we chose for the problem is linear inverse interpolation with a known filter, as described by . This method did produce the desired result, but a non-Gaussian noise distribution caused serious problems in its implementation. To cope with these problems, a version of the iteratively reweighted least squares (IRLS) technique (, ) was applied.

METHOD: INVERSE LINEAR INTERPOLATION WITH KNOWN FILTER

The theory of inverse linear interpolation with a known filter is
given in *Applications of Three-Dimensional Filtering*
(), section 2.6, and can be easily extended to
the 2-D case. The algorithm is based on the following procedure. To invert
the data vector given on an irregular grid for a regularly sampled
model ,
we run the conjugate-gradient solver on the system of equations

(1) |

(2) |

(3) |

FIRST RESULTS

To test the inverse linear interpolation method, we cut a 4 by 4 km patch from the initial data plane and posed the problem of interpolating the data to a regular mesh. Figure galpat shows the patch after sparse, medium, and dense simple binning.

Figure 3

The result of the inverse linear interpolation on the dense grid after 200 iterations is shown in Figure gallap. For a better display, we convolve the interpolated data with a set of first-order derivatives, taken in different directions, such as:

.

Figure 4

The result looks encouraging, since the interpolation produced a more detailed image than simple binning. However, this image is not satisfactory. Studying the sharp dipping discontinuity in the upper-right quarter of the patch revealed that it was produced not by a real geological structure but by an inconsistent track present in the data. The results of applying the interpolation technique to the whole lake bottom (the left side of Figure galsen) show more examples of that kind of noise and demonstrate that the problem of data inconsistency among different tracks is common for most parts of the lake.

Such an inconsistency is easily explained by taking into account the fact that the measurements were made on different days with different weather and human conditions. To allow for these circumstances, we had to reformulate the initial least squares optimization approach applied to inversion.

Figure 5

ENHANCEMENTS OF THE METHOD: IRLS

A minor improvement in the method was the change from the Laplacian filter laplacian in equation regression to the gradient filter consisting of the two orthogonal first-order derivatives

(4) |

The major enhancement is the introduction of a weighting operator into equation linear. Thus system linear - regression takes the form

(5) |

(6) |

(7) |

- 1.
- Bad (inconsistent) data points are indicated by large
values of the residual
*r*=*d*-*Lm*left after conventional least squares inversion. - 2.
- Abnormally large residuals attract most of the conjugate gradient solver's efforts, directing it the wrong way. The residual should be whitened to distribute the solver's attention equally among all the data points and emphasize the role of the ``consistent majority.''

Figure galdat demonstrates changes in the residual space caused by the weighting. Each plot contains one long vector cut into several traces for better viewing. The top left plot is the original data from the patch used in the first experiments. The top right plot is the residual after 200 conventional conjugate gradient iterations. Note two major features of the residual distribution: statically shifted tracks (the zero-frequency component) and local bursts of energy (mostly in the upper part of the plot). To take away the zero-frequency component, we convolved the residual with the first-derivative operator. The result of differentiation appears on the bottom left plot. Finally, to take into account both large spikes created by differentiation and local bursts of energy, we constructed the following weighting function:

(8) |

Figure 6

Thus the weighting operator in equation irls includes two
components: the first-derivative filter taken in the data space
along the recording tracks and the weighting function *w*. Since *w*
depends upon the residual, the inversion problem becomes nonlinear,
and the algorithm approaches the IRLS class ().

The result of the IRLS inversion is shown in Figure galrep. While the noisy portions of the model disappeared as we had expected, a large-scale structure (valley) crossing the lake from north to south remained on the image. Comparison of Figures galavl and galavn demonstrates that the price for that improvement is a certain loss in the image resolution. This type of loss occurs because we rely on the model smoothing to remove the noise influence. The role of the smoothing increases when we apply IRLS in the manner of piece-wise linearization. The first step of the piece-wise linear algorithm is the conventional least squares minimization. The next step consists of reweighted least squares iterations made in several cycles with reweighting applied only at the beginning of each cycle. Thus the problem is linearized within a cycle, and the conjugate gradient solver opens each cycle with the steepest descent iteration. What happens when the nonlinear reweighting is applied on each iteration is shown in Figure galres. Round ``villages'' on the bottom of the lake were identified as the consequence of occasional erroneous data points.

Figure 7

Figure 8

DISCUSSION

The method of inverse linear interpolation with IRLS enhancement was able to produce a fairly clear image of the Sea of Galilee bottom, showing large-scale, presumably geological structures. The lessons learned from the Sea of Galilee project include the following:

- If a model (interpolation result) contains sharp edges, it can be preferable to use a gradient filter instead of a Laplacian one for smoothing purposes.
- The iteratively reweighted least-squares approach is capable of removing non-Gaussian noise influence on the result of inversion.
- The conjugate-gradient solver can behave differently in nonlinear and piece-wise linear implementations of nonlinear problems.

The main cause of the regular noise in the Galilee data set is the inconsistency among repeated measurements. An alternative approach could be applied to the problem of the data track inconsistencies. A well-known example of such an approach is the seismic mis-tie resolution technique (). Combining a technique of this kind with the inverse linear interpolation is a possible direction for future research.

Another interesting untested opportunity is preconditioning. As Bill Harlan has pointed out, the two equations in the system linear - regression contradict each other, since the first one aims to build details in the model, while the second tries to smooth them out. Applying the model preconditioning method could change the contradictory nature of the algorithm and speed up its convergence. The inversion problem takes the form

(9) |

(10) |

Figure 9

Figure 10

[MISC,SEP,SEGCON,GEOPHYSICS,paper]

5/9/2001