Algorithms for Highly Nonlinear Inverse Problems

Principal Investigator

Malcolm Sambridge

Research School of Earth Sciences

An inverse problem arises whenever data provide only indirect constraints on some (Earth) property or parameter of interest. They occur in many fields of the physical sciences. An example in seismology arises when one attempts to determine information on the seismic structure of the Earth at depth from seismograms (following an earthquake) recorded at the earth's surface. Inversion is the process by which indirect information is retrieved from data. Inversion techniques have been applied to a wide range of problems in the Earth sciences. The most difficult type of inverse problem is one where the relationship between the observations and unknowns is highly complex, or non-linear. Examples include inversion of body wave or high frequency seismic waveforms, used in studies of crustal structure and by the exploration industry. The past two decades have seen a rapid growth in both the development and application of non-linear inversion techniques for this type of problem. At present the most popular approaches are based on stochastic or randomized (Monte Carlo) searching of a parameter space of unknowns, where the dimension of the space is equal to the number of unknowns in the problem.

The objective of the present work is to design and implement a new class of non-linear inversion algorithm based on stochastic sampling of a high dimensional parameter space. One which will improve upon previous approaches in being more efficient at solving the inherent global optimization problem and also to extract robust 'error' information on the sought after unknowns.

Projects

r58 - PC


What are the results to date and the future of the work?

The current phase of the project is to investigate the use of the concept of natural neighbours (taken from the field of computational geometry) for general inverse problems involving 10's to 100's of parameters. This means we have to work in parameter spaces of 10-100 or more dimensions and the computational demands are considerable. The primary aim is to apply the concepts of irregular meshing to non-linear inverse problems, in particular the non-linear optimization and numerical integration in high dimensional spaces. Much of the development



work is currently being carried out on Sun workstations. I plan to fully port the code to the SGI and then the VPP for testing in the near future. At present a new stochastic non-linear search algorithm has been developed for the global optimization problem. This method finds 'many models of good data fit' rather than a single model. Preliminary results seem to indicate that its robustness and stability compares favourably with a genetic algorithm method that is currently being used for the same problem in the seismology group of RSES.

What computational techniques are used?

Ultimately the aim is to solve the global optimization problem (i.e. minimizing the misfit between observables and predictions from a set of model parameters) efficiently while generating a large collection of models (set of unknowns) in the region of the global and local minima. These models will then provide the basis of the 'error analysis' stage which will require multi-dimensional integration techniques.

The main approaches used are iterative re-sampling techniques in a high dimensional parameter space. This method was designed specifically for the project in hand and is based on performing random walks according to high dimensional probability density functions. To achieve this I've designed an algorithm that requires the solution of a 'nearest neighbour' problem many times over (~107 or more). A reasonably efficient approach has been developed for this problem, although this is the subject of further study. The next stage of the project will require multi-dimensional integration in high dimensional spaces. I hope to be able to use or adapt 'canned' routines for this task, although it is not yet clear whether this will be possible.

- Appendix A