GRAV3D manual: 
Introduction

The Program and its purpose

The program library GRAV3D is a suit of algorithms for inverting gravity data gathered over a three dimensional earth. Version 1.0 of this facility is described in detail in Li and Oldenburg (1998), and subsequent versions added greater flexibility and efficiency to the basic forward modelling and inversion algorithms. In particular it allowed larger problems to be solved through the use of wavelet transforms, and it allowed geophysical constraints, in the form of upper and lower bounds on the density of each cell, to be included. Notes about Version 3.0 are below.

The problem addressed by GRAV3D involves gravimetric data gathered anywhere at or above the surface of the Earth. These data are the vertical component of the gravity field caused by a three dimensional distribution of density contrast within the volume of ground directly beneath the survey area. This subsurface volume (with optional topography) is modelled as a set of rectangular cells each with constant density contrast. For forward calculations the anomalous density in each cell is known and the data that be measured over this known Earth model are calculated. The inverse problem involves estimating the density contrasts of all the cells based upon measurements gathered during a field survey. In the remainder of this introduction we point out some highlights, and we emphasize the importance of understanding what kind of models the program can recover, and how the program operates. Changes implemented for GRAV3D Version 3.0 are outlined at the end of this chapter.

Summary of sections on theoretical background

The introduction provides a basic understanding of how gravity data relate to the Earth, and the goal of inverting such data.

The solution to the forward problem is described, and an example of a forward calculation is provided.

The method of solving the gravity inverse problem is outlined. In summary, the inversion is solved as an optimization problem with the simultaneous goals of i) minimizing an objective function on the model and ii) generating synthetic data that match observations to within a degree of misfit consistent with the statistics of those data.

• By minimizing the model objective function, distributions of subsurface density contrast are found that are both close to a reference and smooth in three dimensions. The degree to which either of these two goals dominates is controlled by defining length scales for smoothness. This is a crucial step and allows the user to incorporate a priori geophysical or geological information into the inversion. Explicit prior information may also take the form of upper and lower bounds on the density contrast in any cell.
• The relative importance of the objective function the misfit term is controlled by a regularization parameter. This parameter is determined in one of three ways, and depends upon how much is known about errors in the measured data. Section 1.6 discusses this problem, and it is important to understand how the choice of options affects the outcome from the inversion program.

Potential fields data have no inherent information about the distance between source and measurement, therefore the incorporation of a depth or distance weighting term in the formulation is critical. Section 1.4 describes this issue, and explains the options available for controlling how cells in the model enter into the solution regardless of their depth or distance from the measurements.

The large size of useful 3D inversion problems is mitigated by the use of wavelet compression. Parameters controlling the implementation of this compression are available for advanced users, and section 1.5 provides some details on how wavelet compression is applied.

Summary of program elements, instructions for use, and examples

Files and file formats are similar to those of other UBC-GIF codes, although there are some aspects specific to GRAV3D V3.0. Highlights are noted below, but you must read Chapter 2 of the manual for details.

The mesh file defines how the 3D subsurface volume of ground is discretized. Specify cell sizes that will image your target with adequate detail without resulting in a too many model cells. Then add padding cells outside the region of investigation. Recall that the program must invert a sensitivity matrix that has a size proportional to N x M, where N is the number of data points and M is the total number of cells in the model. This sensitivity matrix should reside within the computer’s memory for efficient execution. Problems with more than 10,000 to 20,000 model cells, and/or more than a few thousand data points would
be considered large, and may be expected to require a considerable amount of computing time.

Also, it is considered prudent to upward continue data that are gathered very close to the ground so that measurements appear as if collected at an elevation roughly equal to half a cell width. This is especially true in the presence of severe topography.

Files that define topography, observation locations and observations (measured data) are self explanatory. The manual includes details on how to ensure consistency regarding elevations and topography.

There are separate model files defining density contrasts for the initial model, reference model, upper and lower bounds on density contrast, and final output models, all with similar structures. These models are defined via the mesh file, and each cell has a constant density contrast. Also, you should be clear on how cells above topography are managed.

Forward modelling, sensitivity calculation and inversion programs are run from the command line with parameters specified in separate files. In the descriptions of these parameter files look for details regarding depth or distance weighting, wavelet compression, and how to specify model objective function parameters. There are also important notes on inversion "mode". Your choice of mode determines what method the program uses to find the regularization parameter. This choice is based upon what type of information you have on data errors, and error statistics. Output files are also specified, and you should take particular note of the information provided in the log file. This file summarizes the progress of inversion, and careful examination of its information is a critical aspect of quality control. The log file contents depend upon which mode was use for inversion; see the final section of chapter 3 for details.

The manual concludes with two synthetic examples. The first illustrates the newly added facility for specifying upper and lower bounds on recovered cell density contrasts. The second example illustrates inversion of a larger problem.

Conclusions

Successful application of inversion results to geological problems demands understanding of the inversion process. While models will be as smooth and as close to a reference model as data misfit will allow, the actual results are controlled by careful selection of length scales, weighting functions, and constraining parameters. The size of the problem is reduced using wavelet compression and the management of misfit is controlled by the choice of mode and associated parameters. This degree of flexibility makes it imperative that the user has a good understanding of general inversion theory and the specifics of its implementation for GRAV3D. Finally, please be sure to read the notes regarding changes implemented for Version 3.0.

NOTES: GRAV3D Ver 3.0 (June 2005) - changes to the code and the manual

As might be expected, more recent UBC-GIF codes have features that have been found important for solving practical problems but these features were not included in earlier program libraries. The upgrades described below address this issue. The revised codes are more uniform in capabilities and are more computational efficient.

Improvements since version 2.0

• A new preconditioner for solving the Gauss-Newton system of equations results in significantly improved performance.
• All values, except for the stored sensitivity matrix, are now in double precision. This results in more accurate calculations.
• When entering UTM coordinates, which have a lot of significant digits, the user no longer has to subtract a constant.
• A file sensitivity.txt is output after running gzsen3d.exe. It contains the average sensitivity for each cell. This file can be used for depth of investigation analysis or for use in designing special model objective function weighting.
• A file gzinv3d_nopos.den is output during the first part of the inversion. It contains the densities without the bounds constraints imposed.
• A file gzinv3d_XX.den is output after each beta iteration.
• In gzinv3d.exe the user can enter either alpha values or length scales.
• The reference model is now included in the calculation of the model norm.
• The reference model is now scaled by the depth weighting before starting the no-positivity iterations.
• Sensitivity calculations carried out by gzsen3d.exe are now more efficient.

Changes to run-time files for GRAV3D Version 3.0

None.

Changes to this manual for GRAV3D Version 3.0

This introductory note, and a few additions to "execution".

Notes on computation speed

Run times for GRAV3D Ver2.0 and Ver3.0 are compared in the figures above for a moderately complicated problem. This was a real example involving 1839 data points and 920856 cells in the model, and moderate topography.

• The Version 3.0 code is significantly faster in all cases.
• Speed also depends strongly on the computer. Times when using two computers are shown below. Both were Pentium IV processors with 1 Gbyte of RAM. The significant difference however concerns the size of the CPU’s cache memory. Most good computers sold since roughly 2005 have 1.0Mb of cache memory and this results in a significant improvement in time to complete computation-intensive jobs.
• The complexity of the problem also affects computation times. In the table, compare times when bounds are included to times when no bounds are included. All other aspects of the inversions are identical.
• These results are presented for illustration only. The time to compute any given problem is strongly dependent upon the number of data points, the size of the mesh, and how how you set up all parameters for the inversion, including data, constraints, regularization, compression, etc.