INTRODUCTION

In the last class we talked about several ways to interpret resistivity data. These included empirical correlation, master curves, trial-and-error computer modeling and least-squares parameter inversion. These methods all have strengths and weaknesses. Today we’ll look at some more complex methods that avoid some of the problems of the simple methods.

SMOOTHEST MODELS

The basic concept of this method of inversion comes from William of Occam (1300-1349) who said "it is vain to do with more what can be done with fewer". This principle is called "Occam’s Razor", the "principle of parsimony" or the "principle of simplicity".

Constable et al. (1987) interpret this to mean that the simplest possible model for a sounding is the one that is as close to uniform resistivity with depth as possible. By this reasoning an infinite half-space is the simplest model in mathematical terms because it can be described by just one parameter, the resistivity. Thus the goal is to find a model that fits the field data within experimental error and that has the minimum possible variation of resistivity with depth. In other words the aim is to find the "smoothest" model. Recall that the least-squares method finds the "roughest" or "least smooth" model.

The great advantage of smoothest models is that they test whether the parameter values for some layer are necessary or not. They may demonstrate that the extreme layer-to-layer jumps in some models are simply not necessary. Thus they are very useful in giving the simplest alternative to other models we may have.

From a philosophical point of view I think we can argue that geological simplicity trumps mathematical simplicity when dealing with real data. For example, here in New Hampshire we expect to find big resistivity jumps between marine clay, glacial till and bedrock. For sediments and sedimentary strata elsewhere we don’t see thick layers of uniform resistivity. To the contrary, well logs show large resistivity variations on all scales. Uniform thick layers also conflict with ubiquitous stratification brought about by climate change, tectonism, water depth and so on.

Nowadays it is common to constrain least-squares inversions with a smoothness constraint. These programs combine the simplicity of Occam’s models with the precise data matching of least-squares models.

MONTO CARLO METHOD

With Monte Carlo methods our aim is to find the smallest and largest resistivity values possible at each depth (Sternberg, 1979). In other words we want to find "extremal bounds" of the distribution of resistivity versus depth. We don’t try to find a single "best" model but rather we seek features essential to all models that match the field data acceptably.

In practice we simply let the computer try models at random. To avoid excessively simple models we might select only models with 100 layers, all equal in thickness on a log scale. We might also limit resistivities to those typical of earth materials (i. e. no pure metal layers).

Within these limits the computer simply calculates models with random parameters. The models either fit the data and are saved or don’t fit and are rejected. After we obtain thousands of acceptable models we plot resistivity versus depth. This plot clearly shows the extremal bounds at each depth. Thus (approximately) all possible models must fit within the bounds. Any model with resistivity outside the bounds would not fit the field data. The program can be written to pick out the least-squares and smoothest models as a by-product.

The great advantages of the Monte Carlo method is that the results are independent of the starting model (there is none!) and that common features of all good models are explicitly displayed.

The main disadvantage is that huge amounts of computer time are needed for models with more that a few layers or data points. If we use too few layers the bounds will be unrealistically narrow. If we use too many layers the models will be unnecessarily rough. In addition it may be difficult to include geological information.

ZOHDY METHOD

This method is fundamentally a least-squares method that has been widely used by the United States Geological Survey (Zohdy, 1989). The key features are that the number of layers is equal to the number of data points and that the layer thicknesses are equal on a log scale.

Given these constraints the algorithm readily finds the least-squares solution. In the original program it is difficult to include geologic constraints. The models may change radically if a data point is added or subtracted.

GENETIC ALGORITHMS

Rather than seeking a "best" model of some class, genetic algorithms try to explore the range of parameters for all "good" models; that is, ones that fit the field data within experimental error. For example we might want to get histograms of layer resistivities and layer thicknesses for all good four-layer models. To the extent that we are looking for attributes of all good models, this method resembles the Monte Carlo method. It seeks, however, to do so with far less computation.

We start with an initial "population" (set of models) of, say, 100 three-layer models. In the simplest case all five parameters (three resistivities plus two thicknesses) are selected at random by the computer. We may, of course, impose some geological reality to constrain the range of random numbers allowed.

Then the computer tests the models against the data using some measure of misfit such as least-squares. The worst 50 models are discarded whereas the best 50 pass on to the next "generation" of models. To bring the second generation population back up to 100, we add 50 "mutants". The mutants have parameters picked at random from all the models (good and bad) of the preceding generation. The parameters of the mutants can then be randomly varied within, say, plus or minus ten percent of the values of the parents.

Now the second generation is tested against the field data. Again the best and worst 50 per cent are found and modified as before.

The process continues for a set number of generations (100 in my program). Generation by generation the population becomes more fit. In other words the total or average misfit of the whole population decreases.

At this point almost all the models are acceptable and we can output graphs of resistivity versus depth (approximate extremal bounds) and histograms for all the parameters. In addition we can search the final population for approximations to the least-squares model and the smoothest model.

I have used essentially the same program to model resistivity soundings, gravity and magnetic profiles and seismic refraction profiles.

SIMULATED ANNEALING

This method aims to combine the pure, unbiased randomness of the Monte Carlo method with the directedness of the least-squares method of steepest descent. We want to speed up the Monte Carlo method and avoid the problems of local minima and initial model of the least-squares method.

We start with a reasonable initial model. This model may not fit very well but hopefully is in the right ball-park. Now the computer perturbs the model parameters at random and recalculates the misfit. For a while we repeat this procedure allowing large random parameter changes and allowing the misfit to decrease or increase for iteration to iteration. The large parameter jumps allow the model to avoid getting stuck in local misfit minima. Eventually the size of the allowable random perturbations is reduced and only better fitting models are accepted. In this manner the program homes in on the least-squares model.

Given that simulated annealing models are essentially least-squares models, they share the same pros and cons. The algorithm may be more efficient and less prone to difficulties with local minima and initial models but they are still roughest models.

SUMMARY

We have seen that soundings can be interpreted using empirical correlations, modeled by trial-and-error and inverted by several different schemes. The same basic methods are used throughout geophysics and, indeed, whenever we have to deal with data and quantitative relationships.

Despite the sophistication of these computer methods we should not lose sight of the limitations of the models and data. The world has three dimensions, not just one. Uniform layers are probably the exception, not the rule. Our sounding data are always incomplete and subject to experimental error. Don’t let slick computer packages dazzle you.

REFERENCES

Constable et al., 1978, Geophysics, 52, 289-

Sternberg, 1979, Journal of Geophysical Research,84, 212-

Zohdy, 1989, Geophysics, 54, 245-