Home

GRAV3D Manual - Earth, Ocean and Atmospheric Sciences Dept

1. JBA North 500 North 500 208 488 608 888 1008 208 488 608 S88 1808 aaa 288 408 gaa 1 aaa aun 0 m E ca E H 40 pag ch ona TOF Ana 0 800 Ana ana evant 0 640 agg Ogee 288 225 a a 8 100 288 408 an 1008 H 4 100 288 408 688 10208 H a ABA fa TEA 1AB a ABA 1AB 208 408 408 gaa 1008 3a Synthetic model 3b Recovered model using 3c Recovered model using constraining l constraining bounds of 0 4 bounds of 0 0 8 a a 188 188 2p 488 A apa ag B 208 OBE FEA IEEE 1 0 0 TR no A D TE 3 700 488 Ana o 800 M e 6 600 3 300 a D EE a 1 800 a 722 ADE ADE 1008 io _ a 2008 488 5008 gaa 1 0 100 a FASTING m 3d Variable bounds for result 3 Hee 1220 Example 2 1g B ABB ja E EASTING m 3e Recovered model using variable constraining bounds GBA jena 1209 1508 Easting Data with added noise Large Model Cross sections FARR 2500 a 258 Sg 758 A aan 1258 200 1880
2. The figure below left shows the gravity data simulated on the surface over 51 by 51 grid of 50 m grid spacing A total of 2601 data is generated Noisy observations were simulated by adding Gaussian random noise The data plot shows four of the five anomalies and the effect of the added noise 2588 4 4 06 a a 60 a 4 16 His g 2 70 E 1 1 80 1 1 36 500 a 2 90 a 0 46 a 530 1888 1588 2000 2580 SA 12009 1588 2000 2580 Easting Easting m Gravity data simulated over the large model Predicted data based on the recovered model To invert this data set we discretize the model region using 50 m cubes This results in 62500 cells and the corresponding sensitivity matrix requires over 600 Mb to store We use this example to illustrate the wavelet compression of the GRAV3D version 2 0 Using Daubechies 4 wavelet and a reconstruction accuracy of 5 a compression ratio of 30 was achieved The resulting matrix is stored in 42 Mb With the compressed sensitivity matrix the inversion was carried out readily without much demand on computer memory or CPU time The predicted data from the inversion are shown in the figure above right which is a smooth version of the observation The recovered density contrast model is shown via mouseover in the first Example 2 figure above The cutoff value is 0 17 g cm All five anomalous blocks are imaged The recovered model is shown in
3. ndat El N1 E2 N2 Elev1 Elev2 Endat Nndat Elevndat Parameters are En comments lines beginning with are comments Nn ndat number of observations Elevn easting northing and elevation of the observation measured in meters Elevation should be above the topography for surface data and below the topography for borehole data The observation locations can be listed in any order Easting northing and elevation information should be in the same coordinate system as define in the mesh Example of obs loc file We provide an example file below Example data file Test data 441 of data 0 00 0 00 1 0 0 00 50 00 1 0 0 00 100 00 1 0 1000 00 900 00 1 0 1000 00 950 00 1 0 1000 00 1000 00 1 0 FILE obs grv This file is used to specify the observation locations and the observed gravity anomalies with estimated standard deviation The values of parameters specifying the observation locations are identical to those in obs loc The output of the forward modelling program GZFORS3D has the same structure except that the column of standard deviations for the error is omitted The following is the file structure of obs grv comments ndat El E2 N1 N2 Elev1 Gravl Berl Elev2 Grav2 Err2 Endat Nndat Elevndat Gravndat Errndat Parameters are comments top lines beginning with are comments ndat number of observations En Nn Elevn easting northing and elevation of the
4. The ordering of the cells is the same as that for model cells i j k 1 1 1 is defined as the cell at the top south west corner of the model The total number of lines in this file should equal NN NE NV where NN is the number of cells in the North direction NE is the number of cells in the East direction and NV is the number of cells in the vertical direction The lines must be ordered so that k changes the quickest from 1 to NV followed by j from 1 to NE then followed by from 1 to NN If the surface topography topo dat file is supplied the bounds for cells above the surface will be ignored These values should be assigned a large negative value e g 100 0 to avoid confusion FILE w dat This file contains the values for a user supplied weighting function The following is the file structure for w dat W S1 1 1 W SNN NE NV W E1 1 1 W ENN NE 1 NV W N1 1 1 W NNN 1 NE NV W 21 1 1 W 4NN NE NV 1 Parameters are W Si j k cell weights for the smallest model W Ei j k cell weights for the interface perpendicular to the easting direction W Ni j k cell weights for the interface perpendicular to the northing direction W Zi j k cell weights for the interface perpendicular to the vertical direction Within each part the values are ordered in the same way as in model den however they can be all on one line or broken up over several lines Since the weights for a derivative term are applied
5. If null is entered the surface will be treated as being flat iwt an integer identifying the type of generalized depth weighting to use in the inversion 1 for depth weighting only for surface data 2 for distance weighting surface especially when topography is severe and or borehole beta znot parameters defining the depth weighting function When iwt 1 beta and znot are used as and Zo to define the depth weighting When iwt 2 beta and znot are used as 7 and Zo to define the distance weighting If null is entered on this line line 5 then the program sets beta 2 and calculates the value of Znot based upon the mesh and data location This is true for iwt 1 or 2 For most inversions however setting this input line to null is recommended The option for inputing beta and znot is provided for experienced users who would like to investigate the effect of the generalized depth weighting for special purposes The value of beta should normally be close to 2 0 Note the difference from the value used in MAG3D Smaller values of give rise to weaker weighting wvlet a five character string identifying the type of wavelet used to compress the sensitivity matrix The types of wavelets available are Daubechies wavelet with 1 to 6 vanishing moments daubl daub2 and so on and Symmlets with 4 to 6 vanishing moments symm4 symm5 symm6 Note that daubl is the Haar wavelet and daub2 is the Daubechies 4 wavelet The Daubechies
6. Z 75 i l N 12 0 0 where 4 is the threshold level and sj and g s the elements of and ts respectively The threshold level are determined according to the allowable error of the reconstructed sensitivity which is measured the ratio of norm of the error in each row to the norm of that row 0 It be evaluated directly in the wavelet domain by the following expression i 1 N 09 the numerator is the norm of the discarded coefficients For each row we choose such that ri 4 r where r is the prescribed reconstruction accuracy However this is a costly process Instead we choose a representative row and calculate the threshold level 4 0 This threshold is then used to define a relative threshold bin ar a igi The absolute threshold level for each row is obtained by e max gi PS ee N 14 1 The program that implements this compression procedure is GZSEN3D The user is asked to specify the relative errorr and the program will determine the relative threshold levelg Usually a value of a few percent is appropriate for r For experienced users the program also allows the direct input of the relative threshold level 1 6 Choice of Regularization Parameter The choice of the regularization parameter ultimately depends upon the magnitude of the error associated with the data The inversion of noisier data requires heavier regu
7. name of the executable program If the program is not in the current directory its path must be included also arg_n acommand line argument which is the name of a file It is usually one of those described in the preceding section or a control file containing input parameters 3 2 GZFOR3D This program performs forward modelling Command line usage gzfor3d mesh obs loc model den topo dat Input files All files are in ASCII text format they can be read with any text editor Input files can have any name the user specifies mesh the 3D mesh obs loc the observation locations model den the density model topo dat surface topography optional If omitted the surface will be treated as being flat Output file This file is created by GZFOR3D it s name can NOT be specified gzfor3d grv the computed gravity anomaly data Since the data in this file are accurate the column of the standard deviations for the error is not included 3 3 GZSEN3D This program performs the sensitivity and depth weighting function calculations Command line usage gzsen3d gzsen3d inp For a sample input file type gzsen3d inp Format of the control file gzsen3d inp mesh obs grv topo dat iwt beta wvlet itol eps I nput files and parameters mesh 3D mesh obs grv data file Contains the observation locations and the observed gravity anomalies with estimated standard deviation topo dat surface topography optional
8. 1588 SARE z300 258 7 3B 102p 1258 258 SpA 758 1 08 1258 fm DEPTH Rae 1 1500 2808 2520 EASTING rm 12809 1585 20800 z500 Easting m Predicted data Large Model Recovered Density Cross sections B Kag 1 1500 2008 2580 122a 1258 500 1800 1588 2800 25808 258 508 758 Ai agn 1250 1520 Z508 EABTING ra D Bk D ii 5 30 UBC Earth and Ocean Sciences Jones Tuesday September 05 2006 14 38 00 GRAV3D Graphical User I nterface The graphical user interface GUI for setting up and running 3D gravity inversions is identical to the GUI for MAG3D except for units and default values for depth weighting and reference initial or bounds models Therefore the instructions for using the MAG3D GUI should answer all questions about using the GRAV3D GUI Click here for the MAG3D GUI instructions UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual References Haaz B 1953 Relationship between the potential of the attraction of the mass contained in a finite rectangular prism and its first and second derivatives Geofizikai Kozlemenyek 11 No 7 Li Y and Oldenburg D W 1997 Fast inversion of large scale magnetic data using wavelets 67th Ann Internat Mtg Soc Expl Geophys Expanded Abstracts
9. 490 493 Li Y and Oldenburg D W 1998 3D inversion of gravity data Geophysics 63 109 119 Nagy D 1966 The gravitational attraction of a right rectangular prism Geophysics 31 361 371 UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00
10. GRAV3D manual Background theory This page Introduction Forward modelling Inversion Depth weighting Wavelet compression Regularization Example 1 1 Introduction This manual presents theoretical background numerical examples and explanation for implementing the program library GRAV3D This suite of algorithms developed at the UBCGeophysical Inversion Facility is needed to invert gravimetric responses over a dimensional distribution of density contrast or anomalous density henceforth referred to as density for simplicity The manual is designed so that geophysicists who are familiar with the gravity experiment but who are not necessarily versed in the details of inverse theory can use the codes and invert their data In the following we describe the basics of the algorithm but readers are referred to Li and Oldenburg 1997 1998 for an in depth discussion of various aspects of the algorithm Note that an understanding of these components is necessary for the user to have a global view of the algorithms and to use the program library correctly A gravity experiment involves measuring the vertical components of the gravity field produced by anomalous either excess or deficient mass beneath the surface A distribution of anomalous mass characterized by density x y 2 produces its own gravity field F s which is superimposed on the ambient gravity field By measuring the resultant field and removing the ambient f
11. contrast above in the upper region of the entire model to being very close to zero while allowing the density contrast to vary between 0 and be as high as 1 0 required by the data The recovered model is shown in Figure 3e The top surface of the dipping dyke is well imaged as expected because of the imposed bound constraints However constraints on the top surface of the dyke has greatly helped image the bottom Surface of the dyke Furthermore the lateral extent of the dyke is well imaged although we have not constrained it at all This demonstrates the advantages of imposing specific geologic information through variable bounds that help better define the targets It should be noted that it is possible to create nearly any model you like using bounding constraints They should not be employed unless geological information is very reliable 4 3 Example 2 The figure to the right displays a volume rendered image of the large test model Move your mouse over to see the recovered model It consists of five blocks of different density contrasts in a uniform background There is one large dipping dyke to the left that extends to a large depth Four smaller blocks of various shapes are located at Shallower depths to the right Three cross sections of this true model are shown via mouseover in last figure below This figure indicates the vertical lateral locations the five blocks The model occupies a volume of 2 5 km by 2 5 km by 1 25 km
12. rigorous line search incorporating positivity starts with an initial guess of jt O0 5 o9 This usually yields a misfit that is very close to the target value If the misfit is not sufficiently close to however new guess for is obtained which makes use of the approximate slope so The inversion with updated can be solved efficiently if the logarithmic barrier algorithm is started from an initial model close to the final solution That model is obtained by perturbing the solution corresponding to the previous away from the zero bound The line search using this strategy is often successful in reaching the targets after testing two to four values of jt This strategy is implemented in GZINV3D as the first method for choosing the tradeoff parameter In practical applications the estimate of data error is often begin not available Then the degree of regularization hence the value of jt needs to be determined based on other criteria A commonly used method in linear inverse problems is the generalized cross validation GCV technique The use of GCV in inverse problems with inequality constraints such as positivity is far more ro hot re involved and numerically expensive to implement 5 However applying GCV on the 3D gravity inversion without bounds still produces a reasonable estimate of the data error when the data are dominated by anomalies of the same sign That error can serve as a star
13. to the boundary between cells the weights have one fewer value in that direction For instance the weights for the derivative in easting direction has NE 1 NN NV values whereas the number of cells is NE NN NV If the surface topography topo dat file is supplied the cell weights above the surface will be ignored These weights Should be assigned a value of 1 0 to avoid confusion If null is entered instead of the file w dat then all of the cell weights will be set equal 1 0 UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual Running the programs Ver 3 0 Contents of this page 3 general files GZFOR3D performs forward modelling GZSENS3D calculates sensitivity and the depth weighting function GZINV3D performs 3D gravity inversion GZPRE3D multiplies the sensitivity file by the model to get the predicted data Log file Log file Explanation of the log file contents 3 1 Introduction a MS WindowsXX operating environment codes are best run using the Graphical User Interface GUI However all programs in the package also can be run by typing the program name followed by command line arguments With such a format they can be executed directly on the command line or in a shell script When a program is executed without any arguments it will print a simple message describing the usage The command format is described below Command format PROGRAM arg_l arg_2 arg_3 PROGRAM
14. when needed UPDATES FOR Ver4 0 IN THE APPENDIX Example of gzinv3d inp control file The inversion is started from scratch with a zero reference model The inversion will try to converge to a target misfit equal to the number of data The sensitivity matrix will be stored in memory 0 1 1 0 0 obs nois gzinv3d mtx null 0 0 0 0 1 0 100 100 10 null 0 Input files Input files can irest mode par tolc obsf mtx file initial model reference model lower and upper bounds alphaS alphaE alphaN alphadz 3D weighting idisk 0 be any file name If there are spaces in the path or file name you MUST use quotes around the entire path filename Details in the Elements chapter obs grv input data file The file must specify the standard deviations of the error By definition these are greater than zero gzinv3d mtx sensitivity matrix and depth weighting functionfile calculated by GZSEN3D ini den jnitial model stored in the same way as model den If null is entered the default value determined from the bounds is used For a constant initial model enter a value ref den reference model stored in the same way as model sus If null is entered the default value of 0 0 is used For a constant reference model enter a value bounds den lower and upper bounds on the recovered density value stored in the formats bounds den If null is entered on this line the default v
15. 1 av 4 j l Af 10 FAY d Rij Roy J i is the volume of jth cell Rij is the distance between a point AVj the ith observation and Ro is a small constant used to ensure that the integral is well defined chosen to be a quarter of the smallest cell dimension Similarly this weighting function is normalized to have a maximum value of unity For inversion of data sets acquired in areas with high topographic relief it is advantageous to use this more general form of weighting function The weighting function is directly incorporated in the sensitivity file generated by program GZSEN3D This program allows the user to specify whether to use the depth weighting or the distance weighting depending on the terrain of the observed data 1 5 Wavelet Compression of Sensitivity Matrix The two major obstacles to the solution of large scale gravity inversion problems are the large amount of memory required for storing the sensitivity matrix and the CPU time required for the application of the sensitivity matrix to model vectors The GRAV3D program library overcomes these difficulties by forming a sparse representation of the sensitivity matrix using a wavelet transform based on compactly supported orthonormal wavelets For more details the users are referred to Li and Oldenburg 1997 In the following we give a brief description of the method necessary for the use of the GRAV3D library Each row of the s
16. 4 wavelet should be used for most inversions while the others are provided for users experimentation If null is entered no compression is performed and the program generates a dense matrix in its original form itol eps an integer and a real number that specify how the wavelet threshold level is to be determined itol 1 program calculates the relative threshold and eps is the relative reconstruction error of the sensitivity A reconstruction error of 0 05 is usually adequate itol 2 the user define the threshold level and eps is the relative threshold to be used If null is entered on this line a default relative reconstruction error of 0 05 is used and the relative threshold level is calculated i e itol 1 eps 0 05 Example of gzsen3d inp control file mesh mesh file obs nois data file null topography 2 iwt 1 depth 2 distance null beta znot null daub2 wavelet type 1 0 05 itol eps null Two output files 1 gzinv3d mtx is the sensitivity matrix file to be used in the inversion This file contains the sensitivity matrix generalized depth weighting function mesh and discretized surface topography It is produced by the program and it s name is not adjustable It is very large and may be deleted once the work is completed 2 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 us
17. GRAV3D Version 3 0 A Program Library for Forward Modelling and Inversion of Gravity Data over 3D Structures UBC Geophysical Inversion Facility Department of Earth and Ocean Sciences University of British Columbia Vancouver British Columbia May 2005 UBC Geophysical Inversion Facility 2001 2005 GRAV3D manual Program library On this page Description Package contents Licensing Installation Description GRAV3D is a program library version 3 0 as of August 2005 for carrying out forward modelling and inversion of surface airborne and or borehole gravity data in three dimensions The program library carries out the following functions Forward modelling of the vertical component of the gravity response to a 3D volume of density contrast The model is specified using a mesh of rectangular cells each with a constant value of density contrast and topography is included The gravity response can be calculated anywhere within the model volume including above the topography simulating ground or airborne surveys and inside the ground simulating borehole surveys of surface airborne and or borehole gravity data to generate 3D models of density contrast The inversion is solved as an optimization problem with the simultaneous goals of i minimizing an objective function on the model and generating synthetic data that match observations to within a degree of misfit consistent with the statistics of those
18. a dipping structure The recovered density has a minimum of 0 as constrained by the lower bound of 0 The maximum value is slightly great than 1 0 g cm which is in fact very close to the true maximum density value in this case Overall we have a reasonably good inversion that delineates the essential structure of the true model Result 2 tighter upper bound constraint The second inversion result was obtained using a slightly tighter upper bound to 3 illustrate the use of simple upper bound This is useful when a reliable estimate North 500 of maximum density contrast is available Imposition of such a bound can often a 3200 408 50808 880 1008 improve the solution For this inversion we have set the lower and upper bounds 84 to be 0 0 and 0 8 The maximum value of density contrast in the true model is 1 0 The result is shown in Figure 3c This model is not very different from the 5 400 one recovered in the preceding inversion but the anomalous density appears to 9 800 be slightly wider This is to be expected since we now have smaller density contrast and the required anomalous body should have a large dimension to Da reproduce the same observed anomalies 0 400 0 000 Figure 3 True synthetic models and models recovered by various inversions i 228 Aa ba 1088 1 100 Synthetic model 3b Recovered model using constraining bounds of 0 4 _ Recov
19. alue of 2 0 and 2 0 are used for all cells For constant bounds enter two values for the lower and upper bounds in that order For example to construct a model that is bounded between 0 1 0 this line should be entered 0 1 0 w dat weighting function optional If null is entered the program assumes uniform weight of 1 0 Output files Output files are created by the programs They have fixed file names Detail file formats are in the Elements chapter gzinv3d log The log file containing more detailed information for each iteration and summary of the inversion gzinv3d den Recovered density contrast model gzinv3d pre The predicted data gzinv3d aux Auxiliary file listing the data misfit model norm and Lagrange multiplier at different stages of the inversion This is used only for the purpose of restarting the inversion gzinv3d rho Temporary file containing the density contrast model produced at different stages of the inversion This is used only for the purpose of restarting the inversion gzinv3d_nopos den This model file is output during the first part of the inversion It contains the densities without the bounds constraints imposed Added for grav3d version 3 0 gzinv3d_xx den These model files are output after each beta iteration Added for grav3d version 3 0 Log file explanation The log file gzinv3d log contains more detailed information about the convergence of the inversion Depending on how the i
20. ar region he can specifyw to have increased amplitude there compared to other regions of the model The weighting functions wx wy and wz can be designed to enhance or attenuate structures in various regions in the model domain If geology suggests a rapid transition zone in the model then a decreased weighting for flatness can be put there and the constructed model will exhibit higher gradients provided that this feature does not contradict the data The function w z in eq 4 is a depth weighting that depends upon the model discretization and observation location The weighting function is used to counteract the decay of the kernel functions Gij with depth Because of this decay an inversion which minimizes p oll subject to fitting the data will generate a density distribution that is concentrated near the surface We introduce a weighting of the form w z z zo 3 2 into the objective function to counteract this effect The values of and Zo are discussed in the next section They are important in defining the model objective function and thus the type of final model The next step in setting up the inversion is to define a measure of misfit Here we use the 2 norm measure 4 gw We assume that the noise contaminating the data is independent and Gaussian with zero mean SpecifyingWa to be a diagonal matrix whose i th element is 1 gi where jis the standard deviation of the i th datum makes a chi squared rand
21. 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 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 minimizing an objective fu
22. cations model den the density contrast model Output file gzpre3d mag predicted data UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual Synthetic examples 4 1 Introduction This page includes two synthetic examples to illustrate the GRAV3D Library The first model consists of a dipping slab of positive density contrast buried in a uniform background This is a relatively small model and we use it to illustrate both the basic operation of the library as well as the incorporation of variable bound constraints The second model consists of several blocks of different configurations buried in a uniform background This model has a relatively large number of observations and the number of cells in the inversion is large It is used to illustrate the utility of wavelet compression in Speeding up the inversion of large data sets For both examples we calculate the synthetic data from the true model and then add uncorrelated Gaussian noise to simulate noisy observations 4 2 Example 1 Figure 3a below shows one cross section and two plan sections of the true model This model consists of a dipping dyke in a uniform background The dyke has a westerly dip of 45 and width of 200 by 300 m in the easting and northing directions It extends from 50 m to 400 m in depth The density contrast is 1 0 The observations are simulated above the surface at an elevation of 1 m and over 21 by 21 gri
23. d can be expected to require a considerable amount of computing memory and time The vertical position of the mesh is specified in elevation This is to accommodate the inversion of a data set acquired over a topographic surface When there is strong topographic relief and one wishes to incorporate it into the inversion Special care should be taken to design the mesh A conceptually simple approach is first to design a rectangular mesh whose top specified by Vo is just below the highest elevation point and then to strip off cells that are above the topographic surface This is the approach taken in GRAV3D The number of cells to be stripped off in each column is determined by the user supplied topography file topo dat Only the remaining cells will be used in the forward modelling or included in the inversion as model parameters Example of the mesh file The following is a 10 x 10 x 5 mesh where each cell is 50m by 50m by 50m 10 10 5 0 0 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 50 0 FI LE topo dat This optional file is used to define the surface topography of the 3D model by the elevation at different locations topo dat has the following structure comment npt 1 N1 elevl E2 N2 elev2 Enpt Nnpt elevnpt Comments top lines beginning with are comments npt number of points specifying Ei Ni elevi Easting Northing and elevation of
24. d to within an error tolerance The details of the objective function are problem dependent and can vary according to the available a priori information but generally the objective function should have the flexibility of constructing a model that is close to a reference modelo and producing a model that is smooth in three spatial directions An objective function that accomplishes this is 2 po Dml O as ww 2 pg du ety 2920 du J T F j 4 og Gulez ag 2 J z where the functions Ws Wx Wy and wz are spatially dependent whereas 5 x are coefficients which affect the relative importance of the different components of the objective function The greater the ratio and so on the smoother the recovered model is in that axis direction A length scale of smoothness can be defined for each direction as B Of hg fy yf ir yf Oz fg and specifying greater values for the length scales will produce smoother models The objective function in eq 4 has the flexibility to allow many different models to be constructed The reference model fo may be a general background model that is estimated from previous investigations or it could be the zero model The relative closeness of the final model to the reference model at any location is controlled by the functionws For example if the interpreter has high confidence in the reference model at a particul
25. d with a grid interval of 50 m in both directions The standard deviation of noise added to the data is 2 plus 0 01 mGal The noisy data are displayed in Figure 1 It shows a peak that around 600 500 m and it decays steeply on the east side but slower on the west side due to the dipping structure of the dyke The small scale variations reflect the additive noise in the data 1800 1908 2 70 apa 2 40 apa 5 40 10 2 10 z epa 1 660 H agg 1 80 1 50 E 40a 1 20 498 1 20 0 90 Z on 0 60 apg 0 86 0 30 0 30 a d ee 288 468 188 Figure 1 Measured gravity anomaly mGal with Figure 2 from Resulk noise l Result 1 positive density contrasts Since the simulated data are entirely positive it is reasonable to invert for a model consisting of entirely positive density contrast We accomplish this by imposing a lower bound of 0 0 for positivity and a large upper bound of 4 0 The model objective function has the length scales set to 100 m in all three directions and the simple depth weighting with default parameters is used The predicted data from this inversion is shown in Figure 2 above it is a smooth representation of the noisy data The recovered model is shown in one cross section and two plan sections in Figure 3b The model is characterized by a broad density high at the location corresponding to the true dipping dyke and there is clear indication of
26. data To counteract the inherent lack of information about the distance between source and measurement the formulation incorporates a depth or distance weighting term By minimizing the model objective function distributions of subsurface density contrast are found that are both close to a reference model and smooth in three dimensions The degree to which either of these two goals dominates is controlled by the user by incorporating 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 regularization parameter controlling relative importance of objective function and misfit terms is determined in either of three ways depending upon how much is known about errors in the measured data 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 The research underlying this program library was funded principally by the mineral industry consortium Joint and Cooperative Inversion of Geophysical and Geological Data 1991 1997 which was sponsored by NSERC Canada s National Science and Engineering Research Council and the following 11 companies BHP Minerals CRA Exploration Cominco Exploration Falconbridge Hudson Bay Exploration and Development INCO Exploration amp Technical Serv
27. e txt filename convention so that files are more easily read and edited in the Windows environment The files are ee ee eS FILE mesh 3D mesh defining the discretization of the 3D model region topo dat specifies the surface topography obs loc specifies the observation locations obs grv specifies the observed gravity anomalies with estimated standard deviation model den density model file bounds den optional file containing values for upper and lower bounds w dat contains the 3D weighting functions mesh This file contains the 3D mesh which define the model region mesh has the following structure NE NN NV EO NO VO Aki Aro Ani Ano nnn Ave Parameter definitions NE Number of cells in the East direction NN Number of cells in the vertical direction NV Number of cells in the North direction E0 NO VO Coordinates in meters of the southwest top corner specifie in Easting Northing Elevation The elevation can be relative but it needs to be consistent with the elevation used to specify the observation location in obs loc or obs den and in topo dat see the relevant file for description Aen Cell widths in the easting direction from W to E nn Cell widths in the northing direction from S to N vn Cell depths top to bottom The mesh can be designed in accordance with the area of interest and the spacing of the data available in the area In general the mesh consis
28. e and target misifts is less than tolc Normally the value of par should be 1 0 if the correct standard deviation of error is assigned to each datum When 0 0 is entered for tolc the program assumes a default value of tolc 0 02 mode 2 par is the user input value of tradeoff parameter In this case tolc is not used by the program mode 3 none of the two input values are used by the program However this line of input still needs to be there NOTE When mode 1 both par and tolc are used When mode z2 only par is used When mode 3 neither par nor tolc are used However the third line should always have two values length scales in easting northing and depth directions respectively These parameters determine the weighting coefficients 5 fy ty fz in the model objective function The recommended value of the length scale is two to five cell widths in the corresponding direction If the input sets LE Ly Lz 0 0 the inversion will recover a smallest model If NULL or null is entered the length scale will be equal to two times the maximum cell width at the centre of the mesh For example if the cells are 50m x 50m x 25m at the centre of the mesh then the default values LE Ly Lz7 100m parameter which determines how the sensitivity matrix will be accessed 0 sensitivity matrix will be stored in memory If there is not enough memory idisk will be set to 1 automatically 1 sensitivity matrix will be accessed from disk
29. e in Section 1 2 above The model region is again discretized into 20 by 20 by 10 cells of 50 on side We invert the 441 noise contaminated data illustrated Fig 2 to recover the density contrasts within the 4000 cells For the model objective function we choose Lx Ly Lz 100 0 and unit 3D weighting functions for all components and the default parameters for w z A zero reference model is used We have also imposed a lower bound of 0 0 and an upper bound of 2 0 The final model is shown to the right in Fig 4 The top panel is a cross section through the centre of the model and the other two panels are horizontal slices at different depths The tabular shape of the anomaly and its dipping structure are very clear and the depth extent is reasonably recovered The dip angle inferred from the recovered model is close to the true value The amplitude of the recovered model is slightly less than the true value Over all however this recovered model compares favorably with the true model shown in Fig 1 Northing Northing im 125 UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual Elements of GRAV3D Ver 3 0 Contents of this page Introduction 7 general files mesh contains the finite difference mesh for the 3D modelling and inversions topo contains the topographic data of the earth s surface obs loc contains the locations for the survey obs
30. e in designing special model objective function weighting Added for grav3d version 3 0 3 4 GZI NV3D This program performs 3D gravity inversion Command line usage is gzinv3d gzinv3d inp For a sample input file type gzinv3d inp Format of the control file gzinv3d inp is as follows irest mode par tolc obs mag maginv3d mtx ini sus ref sus bounds sus Le Ly Lz w dat idisk Control parameters mode par tolc LE LN Lz idisk restarting flag 0 start inversion from scratch 1 restart inversion after it is interrupted Restart requires two file written out by MAGINV3D before the interruption maginv3d aux and maginv3d kap see below an integer specifying one of three choices for determining the tradeoff parameter see Figure 4 of background mode 1 the program chooses the tradeoff parameter by carrying out a line search so that the target value of data misfit is achieved mode 2 the user inputs the tradeoff parameter mode 3 the program calculates the tradeoff parameter by applying the GCV analysis to the inversion without positivity two real numbers that are used differently Their use depends upon the value of mode mode 1 the target misfit value is given by the product of par and the number of data N i e The second parameter tolc is the misfit tolerance The target misfit is considered to be achieved when the relative difference between the tru
31. ensitivity matrix in a 3D gravity inversion can be treated as a 3D image and a 3D wavelet transform can be applied to it By the properties of the wavelet transform most transform coefficients are nearly or identically zero When the coefficients with small magnitude are discarded the process of thresholding the remaining coefficients still contain much of the necessary information to reconstruct the sensitivity accurately These retained coefficients form a Sparse representation of the sensitivity in the wavelet domain The need to store only these large coefficients means that the memory requirement is reduced Further the multiplication of the sensitivity with a vector can be carried out by a Sparse multiplication in the wavelet domain This greatly reduces the CPU time Since the matrix vector multiplication constitutes the core computation of the inversion the CPU time for the inverse solution is reduced accordingly The use of this approach increases the size of solvable problems by nearly two orders of magnitude Let G be the sensitivity matrix and be the symbolic matrix representation of the 3D wavelet transform Then applying the transform to each row of G and forming a new matrix consisting of rows of transformed sensitivity is equivalent to the following operation 11 where cis called the transformed matrix The thresholding is applied to individual rows of by the following rule to form the sparse representation ts Gigs
32. er 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 J une 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
33. ered model using constraining bounds of 0 0 8 2 Variable bounds for result 3 2 Recovered model using variable constraining bounds NOTE All figures are shown on one page here suitable for printing Result 3 variable bounds 308 4ea 08 a a 1 ERASING im One of the new features in the GRAV3D Version 2 0 is the ability to impose variable bounds on the density contrast to be recovered in the inversion This provides users with one more tool for inputting geological information to improve the inversion For instance we may expect a lower density contrast in the region of an orebody than that in another region Similarly one region of the subsurface may have a negative contrast while the rest has a positive contrast In special cases imposing a lower and upper bounds that are very close to each other will effective fix the recovered density contrast to a known value during the inversion In this example we illustrate both the variable bounds as well as the use of tight bounds to fix the model values We invert the data introduced above by imposing a lower bound of zero constant throughout the model and a variable upper bounds shown in Figure 3d We assume that the top surface of the dipping dyke is known but we do not know the lateral extent of the dyke nor its depth extent Therefore we impose an upper bound of 0 01 above the dyke and an upper bound of 1 0 below that surface This effectively constrains the density
34. ering the solution with a non zero density For data acquired over a relatively flat surface the sensitivity decays predominantly in the depth direction Numerical experiments indicate that a function of the form 220 closely approximates the kernels decay directly under the observation point provided that a reasonable value is chosen for Zo Having the exponent to be 2 is consistent with the fact that the gravity field decays as inverse distance squared The value ofzo can be obtained by matching the function 2 20 with the field produced by a column of cells directly beneath the observation point In accordance with the discretization used in the inversion we use a depth eighting function of the form 1 2 1 az wF j a 1 9 Az z zp i I for the inversion and is used to identify the jth cell and 2 is its thickness This weighting function is first normalized so that the maximum value is unity Numerical tests indicate that when this weighting is used the density model constructed by minimizing a model objective function in eq 4 subject to fitting the data places the recovered anomaly at approximately the correct depth For data sets acquired over rugged terrain the simple depth decay does not describe the dominant variation in sensitivity Therefore a weighting function that varies in three dimensions is needed We generalize the depth weighting to form a distance weighting 1 4
35. gravity anomaly on the surface of the earth produced by the density model shown in Fig 1 The model consists of a dipping dyke with anomalous density of 1 0 g cm situated in a uniform background The model is represented by 4000 cells with 20 in each horizontal direction and 10 in the vertical direction The anomaly is calculated at a height of 0 5 m above the surface over a 21 by 21 grid with 50 m spacing in both directions A contour presentation of the data is shown in Fig 2 Figure 1 Synthetic model dipping dyke with anomalous density 1 g cm3 E mi 200 A cross section at north 525 Middle plan section at depth 125 a 400 ipog Bottom plan section at depth 225 1000 0 880 0776 Figure 2 The gravity anomaly in units of mGal 500 Uncorrelated Gaussian noise has been added to _ 0 607 the data The standard deviation of the additive noise 600 gs58 IS 0 01 mGal plus 5 of the magnitude of the datum E 400 0 444 1000 2 5000 lt 0 333 2 2800 200 35 Tala 0 222 800 2 0100 12 0 111 000 1 7300 0 000 00 1 4600 200 D B n 1 1800 600 400 0 9080 5 200 0 6320 200 0 3570 Depth 225 0 0 0817 g 0 200 400 a00 800 1000 200 400 600 800 1000 Easting m Easting im 1 3 Inversion Methodology The inverse problem is formulated as an optimization problem where an objective function of the density model is minimized subject to the constraints that the data be reproduce
36. grv contains the observations model den contains the cell values for the density mode bounds den specifies upper and lower bounds for density values w dat contains special weightings that alter the type of model produced in the inversions 2 1 Introduction The GRAV3D library consists of three major programs and one utility GZFOR3D performs forward modelling GZSEN3D calculates sensitivity and the depth weighting function GZINV3D performs 3D gravity inversion GZPRE3D multiplies the sensitivity file by the model to get the predicted data This rarely used utility multiplies a model by the sensitivity matrix in gZinv3d mtx to produce the predicted data This program is included so that users who are not familiar with the wavelet transform and the structure of gzinv3d mtx can utilize the available sensitivity matrix to carry out model studies Each of the above programs requires input files as well as the specification of parameters in order to run However some file are used by a number of programs Before detailing the procedures for running each of the above programs we first present information about these general files 2 2 General Files for GRAV3D Programs There are seven general files which are used in GRAV3D version 3 0 All in ASCII text format I nput files can have any name you like Only program output files have restricted file names Also the filename extensions are not important Many prefer to use th
37. ices Kennecott Exploration Company Newmont Gold Company Noranda Exploration Placer Dome WMC The theoretical framework for GRAV3D is provided in the following paper see the UBC GIF website publications page for details Li Y and Oldenburg D W 1998 3D inversion of gravity data Geophysics 63 No 1 109 119 Two short papers including examples of applying GRAV3D in mineral exploration contexts are Cost effectiveness of geophysical inversions in mineral exploration Applications at San Nicolas Nigel Phillips Doug Oldenburg and Jiuping Chen Yaoguo Li Partha Routh 2001 The Leading Edge Volume 20 Issue 12 p 1351 Applications of Geophysical Inversions in Mineral Exploration Problems Oldenburg D W Li Y Farquharson C G Kowalczyk P Aravanis T King A Zhang P and Watts A 1998 The Leading Edge 17 461 465 Software package contents The package that can be licensed includes the following components Executable programs for performing 3D forward modelling and inversion of gravity surveys The GRAVS3D library WindowsxXxX or Linux platforms consists of three major programs and one utility GZFOR3D performs forward modelling GZSEN3D calculates sensitivity and the depth weighting function GZINV3D performs 3D magnetic inversion GZPRE3D multiplies the sensitivityfile by the model to get the predicted data graphical user interface is supplied for the WindowsxXxX platforms o
38. ield from the measurements through numerical processing one obtains the field due to the anomalous mass The vertical component of the gravity field produced by the density x y z is given by 0 ol i da Se Fal where p o is the vector denoting the observation location is the source location V represents the volume of the anomalous mass is the gravitational constant Here we have adopted a Cartesian coordinate system having its origin on the earth s surface and the z axis pointing vertically downward The data from a typical gravity survey are a set of field measurements acquired on a 2D grid over the surface These data are first processed to yield an estimate of the anomalous field which is due to the excess or deficient mass below the data area The goal of the gravity inversion is to obtain from the extracted anomaly data quantitative information about the distribution of the anomalous density in the ground The GRAVS3D library consists of three programs 1 GZFOR3D calculates the surface gravity data produced by a given density model 2 GZSEN3D calculates the sensitivity matrix for use in the inversion 3 GZINV3D inverts the anomalous gravity field to generate a density model In the following we outline the basics of the forward and inverse procedures used by these programs 1 2 Forward Modelling Forward modelling of gravity data is a linear problem and can be ca
39. ile contains the cell values of the density contrast model model The density contrast must have values in The following is the file structure of model den deny 1 1 den1 1 2 den1 1 NV den1 2 1 deni dennn NE NV Each deni j k is the density contrast at location i j k ij k 1 1 1 is define as the cell at the top south west corner of the model The total number of lines in this file should equal NN x NE x NV where NN is the number of cells in the North direction NE is the number of cells in the East direction and NV is the number of cells in the vertical direction The lines must be ordered so that k changes the quickest from 1 to NV followed by j from 1 to NE then followed by from 1 to NN If the surface topography topo dat file is supplied the values above the surface will be ignored These values should be assigned a large negative value e g 100 0 to avoid confusion with the other model elements FILE bounds den This file contains the cell values of the lower and upper bounds on the sought density model It is only required optionally by gzinv3d The bounds have the same dimension as the density contrast The following is the file structure of bounds den ub1 1 1 lbi 1 ub1 1 2 1b1 1 NV ub1 1 NV 1 ub1 2 1 lbi 3 ubi 3 lbNN NE NV UDNN NE NV Parameters are Ibi j k is the lower bound on cell i j k ubj j k is the upper bound on cell ij k
40. in Mode 1 Mode 3 In this mode the program first performs a GCV estimate of data noise to obtain an approximate value of the regularization parameter and it then carries out a single logarithmic barrier minimization to obtain the inverse solution During the first stage a number of trial values of the regularization parameter is tested and two linear systems are solved for each value Thus the first part of the log file is organized according to the value of regularization parameter For each value the numbers of CG iterations data misfit model norm and GCV value are listed At the end of GCV search the log file lists the data misfit model norm GCV value according to the regularization parameter sorted in increasing order The regularization parameter corresponding to the lowest GCV value is used to obtain the final solution in the second stage The second part of the log file corresponds to a single logarithmic barrier minimization and it is identical to that when mode 2 is chosen 3 5 GZPRE3D This utility multiplies a model by the sensitivity matrix in gZinv3d mtx to produce the predicted data This program is included so that users who are not familiar with the wavelet transform and the structure of gzinv3d mtx can utilize the available sensitivity matrix to carry out model studies Command line usage gzpre3d gzinv3d mtx obs loc model den Input files gzinv3d mtx the sensitivityfile from GZSEN3D obs loc the observation lo
41. larization thus a greater value is required In this section we discuss the various implementations for the choice of jt in the GRAV3D library If the standard deviation associated with each datum is known then the data misfit defined by eq 5 has a known expected value which is equal to the number of data when the errors are assumed to be independent Gaussian noise with zero mean The value of should be such that the expected misfit is achieved This entails a line search based on the misfit curve as a function of jt Because of the positivity constraint our problem is nonlinear Thus for each nonlinear solution using a logarithmic barrier method must be obtained This is computationally demanding and we therefore have developed the following strategy to reduce the cost It is observed that when plotted log log scale the misfit curves for 3D inversion with and without bound constraints often parallel each other in the vicinity of the expected misfit The curve with positivity must lie above the curve without positivity Therefore we can first perform a line search without positivity to find that gives rise to pe This search also generates the slope So of misfit curve at jto This process is very efficient and the required CPU time is much smaller compared to the time required for the solution with positivity We next assume thatSg can be used to approximate the slope of the misfit curve when the positivity is imposed A
42. nction 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 inver
43. nd Thus the solution is obtained by the following constrained minimization problem 1 i nue 2 F rae 7 subject to Pumin lt lt Pimax where j min and max are vectors containing the lower and upper bounds the model values When the standard deviations of data errors are known the acceptable misfit is given by the expected values and we will search for the value of that produces the expected misfit Otherwise an estimated value of will be prescribed The details of various aspects of choosing a regularization parameter will be discussed in a following section We use a primal logarithmic barrier method with the conjugate gradient technique as the central solver In the logarithmic barrier method the bound constraints are implemented as a logarithmic barrier term The new objective function is given by Ad lA da 2 In 9 yn In pihl 8 j l where is the barrier parameter and the regularization parameter is fixed during the minimization As the name suggests the logarithmic barrier term forms a barrier along the boundary of the feasible domain and prevents the minimization from crossing over to the infeasible region The method solves a sequence of nonlinear minimizations with decreasing 4 and as approaches zero the sequence of solutions approaches the solution of eq 7 The above methodology provides a basic framework for solving 3D gravity in
44. nly Facilities include GRAV3D GUI EXE a primary interface for setting the inversion and monitoring the progress of calculations GM DATA VIEWER a utility for viewing raw surface or airborne data but not borehole data error distributions and for comparing observed to predicted data directly or as difference maps MESHTOOLS3d a utility for displaying resulting 3D models as volume renderings Susceptibility volumes can be sliced in any direction or isosurface renderings can be generated Documentation is elsewhere via the menu to the left Note there is currently no documentation for the GRAV3D GUI user interface However it is very similar to the MAG3D GUI and there is documentation for MAG3D GUI EXE Example data sets and excercises are provided on the AG CD ROM Licensing A constrained educational version of the program is avaliable with the IAG CD ROM The educational version is fully functional so that users can learn how to carry out effective and efficient 3Dinversions of gravity data However RESEARCH OR COMMERCIAL USE IS NOT POSSIBLE because the educational version will NOT work with more than 200 data points or 12 000 cells in the 3D mesh Licensing for an unconstrained academic version is available see the licensing policy document on the UBC GIF website NOTE all academic licenses will be time limited to one year You can re apply after that time This ensures that everyone is using the most recent ver
45. nversion is set up by the users the content of the log file is slightly different In general there are two stages In the first stage the program estimates an approximate regularization parameter In the second stage the program performs the inversion with bound constraints using a logarithmic barrier method of minimization which consists of outer iterations with the barrier parameter and inner iterations of conjugate gradient solution of an linear system within each barrier iteration The log file information is organized according to these two levels of iterations We can refer to the flow chart Figure 3 in the Background chapter to understand the output in the log file Below we briefly describe the content of the log file according to the parameter mode chosen for the inversion Mode 1 In this mode a user defined target misfit is supplied by specifying chifact The first stage of inversion performs a line search without the bound constraints to estimate an approximate value of the regularization parameter and the slope of the misfit curve The log file identifies this segment and lists summaries of each linear solution such as the number of conjugate gradient CG iterations data misfit and model norm The second stage of inversion carries out a number of logarithmic barrier solutions Each solution corresponding to a single regularization parameter is obtained by a sequence of barrier iterations and nested inside each barrier iteration is a
46. observation measured in meters Elevation should be above the topography for surface data and below the topography for borehole data The observation locations can be listed in any order Gravn gravity anomaly data measured in mGals Errn standard deviation of Gravn This represents the absolute error It CANNOT be zero or negative NOTE It should be noted that the data Gravn are extracted anomalies after standard data reduction and regional removal are applied A suitable data set is one in which there are just a few data values for each mesh cell at the Earth s surface Also it is prudent to upward continue data so they appear as if they were gathered at an elevation roughly equal to half a cell width above the Earth s surface NOTE Careful attention to the regional is particularly important for gravity data This is because any residual anomaly near mesh edges or which could be caused by features too large to be modelled with the given mesh must be accounted for in the inversion result even if the causative geologic materials are in reality outside the mesh The result will be erroneous density anomalies in the recovered model Example of obs grv file Examples file 441 of data 0 00 0 00 2 00 0 978189E 01 0 293238E 02 0 00 50 00 2 00 0 109359E 00 0 320770E 02 0 00 100 00 2 00 0 122059E 00 0 351328E 02 0 00 500 00 2 00 0 236872E 00 0 562092E 02 0 00 550 00 2 00 0 224701E 00 0 556866E 02 FI LE model den This f
47. om variable distributed with N degrees of freedom Accordingly E a N provides a target misfit for the inversion 9 5 The inverse problem is solved by finding density which minimizes misfits the data according to the noise level This is accomplished by minimizing 02 where ug 0 is a regularization parameter that controls the relative importance of data misfit and model objective function A value of jt is sought so that the data are neither fit too well nor too poorly To perform a numerical solution we first discretize the objective function in eq 4 according to the mesh that defines the density model This yields 0 Fo WIW WEW 2 7 Foy Win Wn 7 9 1102 where where and oare M length vectors The individual matrices Ws Wx Wy Wz are straightforwardly calculated once the model mesh the weighting functions 2 and Ws W x W w z defined The cumulative matrix W mWm is then formed The details of these quantities are specified by the user through the inputs to programs GZSEN3D and GZINV3D Since the density contrast is limited to a small range for any practical problems and often there are well defined bounds on the density contrast based on direct sampling or other geological information we also impose constraints to restrict the solution to lie between a lower and upper bou
48. rried out by performing the integration in eq 1 We divide the region of interest into a set of 3D prismatic cells by using a 3D orthogonal mesh and assume a constant density within each cell We discretize the density model in this manner since it is best suited for our inversion methodology Given such a discretization the gravity field at thei th location can be written as d F roi M 2 20 gt a J F Fl AV j l 2 id py Gij j l In eq 2 Pj Vj are the density and volume of the j th cell di is introduced as a generic symbol for the i th datum and Gij defined by the expression in brackets quantifies the contribution of thej th cell to the i th datum The solution for the integral in eq 2 can be found in Nagy 1966 and we have adopted the solution by Haaz 1953 here Expressed matrix notation the gravity data consisting of observations are given by G d 3 where 1 d1 dx is the data vector and 1 Pm is the vector containing the density values of the M cells The program GZFOR3D performs forward modelling calculation to generate the surface gravity data produced by a given density model Program GZSEN3D calculates the complete matrix G to be used in the inversion with optional wavelet compression applied To illustrate the forward modelling program and to introduce the example used in the inversion section of this manual we calculate the
49. set of CG iterations This portion of the log file begins with the value of regularization parameter initial values of data misfit and model norm It is then organized in segments corresponding to barrier iterations and lists the number of CG iterations data misfit model norm barrier parameter and the value of the barrier function at the conclusion of that iteration As barrier iterations progress the data misfit barrier parameter and barrier function should decrease monotonically The model objective function may increase or decrease depending on the nature of the initial model of the logarithmic barrier minimization However both data misfit and model objective function should plateau at the end One or more solutions may be obtained to complete the line search and to achieve the target misfit Each solution will have a distinct portion in the log file The trial values of the regularization are dynamically predicted during the inversion and they are not necessarily in any order Upon completion of line search and the inversion the log file will list a summary of the data misfit and model norm corresponding to different values of the regularization parameter sorted in increasing order Mode 2 In this mode the user specifies a value of regularization parameter and the program performs a single logarithmic barrier minimization to obtained the solution The log file mainly consists of the information for one logarithmic solution as described
50. sion 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
51. sions of codes Licensing for commercial use is managed by distributors not by the UBC GIF research group Details are in the licensing policy document For learning and documentation For links to documentation related utilities and examples see the menu to the left Installing GRAV3D version 3 0 educational version For users with copy of IAG only 1 Copy all files in this folder onto your computer Place them all together in a new folder such as c ubcgif grav3d 2 No further installation is necessary 3 Follow instructions in one of this CD ROM s exercises or refer to the program documentation 4 Recall that this is an educational version Codes will not work with more than 200 data points or 12 000 cells in the 3D mesh UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual Introduction The Program and its purpose The program library GRAVS3D 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
52. 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 GRAVSD computation times two computers for both versions 2 0 and 3 0 of data 1839 of cells 920 556 CPU clock 3 0 Ghz cache 0 5Mb CPU clock 3 2 Gh
53. the pn point The elevation in this file and Vo in the meshfile must be specified relative to a common reference The lines in this file can be in any order as long as the total number is equal tonpt The topographic data need not be supplied a regular grid GRAV3D assumes a set of scattered points for generality and uses triangulation based interpolation to determine the surface elevation above each column of cells To ensure the accurate discretization of the topography it is important that the topographic data be supplied over the entire area above the model and that the Supplied elevation data points are not too sparse If topo dat is not supplied the surface will be treated as being flat Example of topo dat file topographic data 4 0 0 1000 1000 NOTE O O O 1000 0 O 1000 O 50 0 O 50 0 O 50 0 50 0 Although the cells above the topographic surface are removed from the model they must still be included in the model file model sus as if they are a part of the model For input model files these cells can be assigned any value The recovered model produced by inversion program NV3D also includes the cells that are excluded from the model but these cells will have a value of 100 0 as identifier FI LE obs loc This file is used to specify the observation locations It is used by the forward modelling code only The following is the file structure of obs loc comments
54. three crosssections in the figure to the right Move your mouse over to see the true model The amplitudes of the recovered anomalous blocks are lower than the true value The depth extent of the large dipping block is also smaller This is expected since the area of the data is limited Over all this is a good representation of the true model and the inversion utilizing the wavelet compression has allowed the inversion to be carried out with very moderate demand on computational resources NOTE All figures are shown on one page here suitable for printing Large Model Recovered Density Cross sections 2 1588 2008 2520 1 20 1258 Neath 1200 a 1 1588 2588 2 750 A aan 1250 a 1880 7788 2800 2508 TASTING D kH 0 Bem g D E D Bek 2 4200 D UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 38 00 GRAV3D manual Examples figures only Example 1 1888 1888 atl apg 6 40 2 40 2 10 _ 2 10 sga 1 00 soa 1 80 1 50 1 60 1 20 48a 1 20 2 90 0 60 200 0 60 0 30 0 30 A A 22a 1806 g 22a 488 a apg 1800 Easting Easting m Figure 1 Measured gravity anomaly mGal with noise Figure 2 Predicted data from Result 1 Figure 3 models a a 1Ba 18 10 T Toa 288 200 aaa 5
55. ting point for further adjustment by the user based on his or her judgement Since no other information is assumed we have chosen to use the value of obtained in this manner directly in the final inversion which has bound constraints imposed In this case only one logarithmic barrier solution is needed Numerical tests have indicated that this Simplistic use of GCV is in fact surprisingly effective unless the data have a large negative bias or are distributed sparsely GZINV3D has implemented this approach as the third method for choosing the tradeoff parameter find pe by GCV b i m O update Ho by log barrier method The flowchart shown to the right illustrates the structure of the program GZINV3D It has three options for determining the tradeoff parameters The controlling exit parameter is mode When mode 1 the line search based on known target value of data misfit is used Two stages as discussed above are used and several solutions for different values of j must be tested to obtain one that produces the target misfit When mode 2 the user specifies a tradeoff parameter and a single solution is produced When mode 3 the program first performs GCV analysis on the inversion without positivity and then uses the resultant value of in the final inversion 1 7 Example Inversion We now illustrate results of using the program GZINV3D to invert the surface gravity data that was generated for the exampl
56. 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 chapt
57. ts of a core region which is directly beneath the area of available data and a padding zone Surrounding this core mesh Within the core mesh the size of the cells should be comparable with the spacing of the data There is no restriction on the relative position of data location and nodal points in horizontal direction The cell width in this area is usually uniform The maximum depth of the mesh used for inversion should be large enough so that the no density contrast below that depth would produce a noticeable anomaly with the length scale covered by the data area A rule of thumb is that the maximum depth should be at least half of the longest side of the data area Based upon the user s knowledge of the Survey area one may adjust the maximum depth as necessary The cell thickness in vertical direction usually increases Slightly with depth In the shallow region the ratio of thickness to width of about 0 5 is good especially when surface topography is present At depth a cell thickness close to the cell width is advisable Once this core mesh is designed it can be extended laterally by padding with a few cells possibly of variable width This padding is necessary when the extracted anomalies are close to the boundary of the core mesh or if there are influences from anomalies outside the area which cannot be easily removed Problems with more than 20 000 to 30 000 model cells and or more than a few thousand data points would be considered large an
58. version with arbitrary observation locations The basic components are the forward modelling a model objective function that incorporates a depth weighting a data misfit function a regularization parameter that ultimately determines how well the data will be fit and the logarithmic barrier method to obtain the solution with bound constraints Without getting into the algorithmic details we discuss three of these basic components in the next sections namely the depth weighting efficient forward mapping and the choice of the regularization parameter An understanding of these components is necessary for the user to have a global view of the algorithm and to use the program library correctly 1 4 Depth Weighting It is well known that gravity data have no inherent depth resolution As a result structures will concentrate near the surface when a simple smallest or flattest model is produced regardless of the true depth of the causative bodies In terms of model construction this is a direct manifestation of the kernels decay with depth Because of their rapidly diminishing amplitude at depth the kernels of surface data are not sufficient to generate a function that possess significant structure at depth In order to overcome this the inversion needs to introduce a depth weighting to counteract this natural decay Intuitively such a weighting will approximately cancel the natural decay and give cells at different depths equal probability of ent
59. z cache 1 0Mb 2 0 VSO 2 WoO sens 1 05 29 0 35 34 sens 0 57 23 0 16 48 mnj 3946 0 53 22 mnj 15607 0 38 35 bounded imwnj 6 59 12 1 30 39 bounded imn 2 20 40 1 05 20 invn hounded invn bounded invn inm 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 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 UBC Earth and Ocean Sciences F Jones Thursday October 12 2006 14 23 17

GRAV3D Manual - Earth, Ocean and Atmospheric Sciences Dept

Contents

Download Pdf Manuals

Related Search

Related Contents