Home

Complete - Earth, Ocean and Atmospheric Sciences Dept. UBC

Contents

1. Background theory Elements of the library Executing the programs Examples Graphical user interface UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 39 00 MAG3D manual Background theory This page Introduction Forward modelling borehole Inversion Depth weighting Wavelet compression Regularization Example 1 1 Introduction This manual presents theoretical background numerical examples and explanation for implementing the program library MAG3D This suite of algorithms developed at the UBC Geophysical Inversion Facility is needed to invert magnetic responses over a three dimensional susceptibility distribution The manual is designed so that a geophysicist who is familiar with the magnetic experiment but who is not necessarily versed in the details of inverse theory can use the codes and invert his or her data A magnetic experiment involves measuring the anomalous magnetic field produced by magnetically susceptible materials beneath the surface which have been magnetized by the earth s main magnetic field The material with susceptibility rfr y zis magnetized when the earth s main field with flux intensity Bo impinges upon the subsurface formation The magnetized material gives rise to a magnetic field which is superimposed on the inducing field to produce a total or resultant field By measuring the resultant field and removing the inducing field from the
2. comment npt E1 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 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 to npt The topographic data need not be supplied on a regular grid MAG3D 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 0 0 50 0 0 0 1000 0 50 0 1000 0 0 0 50 0 1000 0 1000 0 50 0 NOTE 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 MAGI NV3D also includes the cells that are excluded from the model but these cells will have a value of 1 0 as identifier FI LE obs loc
3. Inversion mode MAG3D Version 3 x can determine the tradeoff parameter in either of three ways 1 Select chifact and enter a value for chifact The program chooses the tradeoff parameter by carrying out a line search so that a target value of data misfit is achieved This target value is di chifact X N where N is the number of data points This is a common option if errors on data are assumed to be Gaussian and un correlated 2 Select user and enter a fixed value for the tradeoff parameter This option is rarely used in practice 3 Select GCV to ask the program to calculate a tradeoff parameter by applying the GCV analysis to the inversion without positivity This option is computing intensive See the user manual section 1 6 for details regarding these three methods of choosing regularization Setting bounding constraints This section allows for specification of upper and or lower bounding values for susceptibilities of each cell There are three ways of specifying bounds 1 Default If the user chooses to not define upper and lower bounds the program employs default bounds for susceptibilities in every cell of 1 0 and 0 0 respectively S I units While it is true that some rocks have susceptibility greater than 1 0 MAG3D Version 4 0 still assumes small susceptibilities as this code has done since the original version When there are very high susceptibilities the relation between incident and induced magnetization is no
4. This file is used to specify the inducing fieldparameters anomaly type and observation locations The following is the file structure of obs loc comments Incl Decl geomag Aincl Adecl idir ndat El N1 Elev1 aincll adecll E2 N2 Elev2 aincl2 adecl2 Endat Nndat Elevndat ainclndat adeclndat Parameters are top lines beginning with are comments comments Incl jnclination and declination of the inducing magnetic field The declination is specified east with Decl respect to the northing used in the mesh obs loc and obs mag files geomag strength of the inducing fieldin nT Aincl jnclination and declination of the anomaly projection Adecl idir Q multi component data set Observations have different inclinations and declinations aincln and adecin specifying the projection direction of the anomaly In this case Aincl and Adecl should be equal to Incl and Decl respectively 1 single component data set All observations have the same inclination and declination of the anomaly projection Aincl Adecl If idir is missing it is assumed to be equal to 1 ndat number of observations When single component data are specified the number of observations is equal to the number of data locations When multi component data are specified the number of observations will exceed the number of data locations For example if three component data are specifie at N locations the number of observations is 3N En Nn eas
5. Z ww m 1000 The depth labelled in each section is 500 500 referenced to the surface elevation of 125 m E 0 D sog 1100 fl 2000 2500 343000 o 500 1100 lid ECS The large block at depth and the Peete Fastinn im Small block near the surface have a susceptibility of 0 08 SI unit and the other blocks have a susceptibility of 0 05 Click this button to see the true synthetic model UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 39 00 MAGS3D Ver 4 0 Help Graphical User Interface One page version UBC GIF i This page contains brief instructions on E Untitled MAG3D User Interface aE using the MAG3D graphical user File view inversion Heip interface shown to the right It is i a N a provided as a printable version of the g HTML point and click version Mag observations He available by clicking here Browse ViewDala ypc Panes ape Depth weight Hesh Menus E eae ee F dph E detak Browse The four menus provide essentially the C tance Com P 20 Create mesh same functions as the tool bar buttons Wavelet compression METEEN i pa ie C described below except for k SEO Fas FA Flt File i oor f relative reconstruction enor o e The File menu Open option na C pE permits opening of a parameter file i if l d ist Mode Bours inp if one already exists E chiad m E Delmi I bwar ef e The File menu New option C constant wadeni a GO basically
6. PDF format has complete details on every parameter and all file Structures
7. R g 92 93 i lt 6 91 92 93 R h lo Is The rotation matrix R therefore allows measured components in local coordinates to be rotated into global coordinate or the components of the regional field to be rotated into local coordinates for use in regional removal Numerical Implementation We divide the region of interest into a set of 3D cuboidal cells by using a 3D orthogonal mesh and assume a constant susceptibility within each cell By eq 1 we have an uniform magnetization within each cell and its field anomaly can be calculated using eqs 3 and 6 The actual anomaly that would be measured at an observation point is the sum of field produced by all cells having a non zero susceptibility value The calculation involves the evaluation of eq 3 in a 3D rectangular domain define by each cell The program that performs this calculation is MAGFORS3D As input parameters the coordinates of the observation points and the inclination and declination of the anomaly direction must be specifie for each datum For generality each component in a multi component data set is specified as a separate datum with its own location and direction of projection As an illustration of the forward modelling program we calculate the total field anomaly above the surface and three component anomaly in boreholes produced by a synthetic model The model consists of two vertical prisms buried at different depths Figure 1 displays one cross section and
8. clears the interface to l k eni iraia model Reference model eS eee G Delak C Value Si 0 001 Delak C Vake Sl 0 yas rel 8 8 8 Save button be ie l led Alpha Length arabes The parameter set in the GUI window C Ae Ae AN AR mn fF 7 must be saved before an inversion can be SAMAR y Le Ln Lz run The results of inversion will be UBC Geophysical Inversion Facility placed in the directory where this input file is saved Important note about directories Save each inversion run into a new directory because all output files get saved using default names Help button This basically gives the program name and version number and some notes about it s development Please be aware of the copyright notice Run and kill buttons The inversion can be run only after the parameters have been saved using the save button Use the kill run button to stop an inversion before it ends by itself Do not end an inversion prematurly by closing the DOS window This button is inactive unless an inversion is underway View resulting log data and model files Four buttons here provide access to inversion results and to log files summarizing the inversion program s performance e Clicking logs pops up the sensitivity log file into a text editor window This is only useful for debugging e Clicking logi pops up the inversion log file into a text editor window This contains details about each step and iteration of the inversion
9. for one logarithmic solution as described 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 MAGPRE3D This utility multiplies a model by the sensitivity matrix in maginv3d 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 maginv3d mtx can utilize the available sensitivity matrix to carry out model studies Command line usage magpre3d maginv3d mtx obs loc model sus Input files maginv3d mtx the sensitivi
10. 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 and 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 specifie 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 MAG3D 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 O 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
11. maginv3d sus Recovered susceptibility model maginv3d pre The predicted data maginv3d 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 maginv3d kap Temporary file containing the susceptibility model produced at different stages of the inversion This is used only for the purpose of restarting the inversion maginv3d_nopos sus This file is output during the first part of the inversion It contains the susceptibilities without the bounds constraints imposed Added for mag3d version 4 0 maginv3d_XX sus These files are output after each beta iteration Added for mag3d version 4 0 Example of maginv3d 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 1rest 1 mode 1 0 O par tolc obs nois obsf a E Ps maginv3d mtx mtx file 1 bal i i L Md pd 0 001 initial model 0 0 reference model 0 5 0 001 Lower and uppper bounds SI null lower upper bounds 0 0001 1 1 1 alphaS alphaE alphaN alphaz null 3D weighting 0 idisk Log file explanation This explanation comes from the GRAV3D manual the MAG3D log file is basically identical The log file maginv3d log contains more detailed information about the convergen
12. one plan section of the model The inducing field has inclination 65 and declination D 25 The surface anomaly is calculated at an interval of 25m over seven east west lines spaced 100 m apart The borehole data are calculated in three vertical holes at an interval of 25m Figure 2 shows the contour map of the surface data and the depth profiles of the borehole data B 188 mm E epp sae Aap i 380 200 100 S 188 2800 328 ai 488 0 085 0 080 e pE mg T 108 i a 0 015 t5 g o10 0 05 S 10B Figure 1 _3BB Two sections of a susceptibility model consisting of two prisms buried in a nonsusceptible background 388 The collar positions of three vertical boreholes are 300 2 0 a8 180 288 2 indicated on the plan section EASTING m Jaa Hole A Hole B o meg JAA aan 20 100 E TE Ad 7 160 150 F T H os 128 E 200 E4 LU ZETT PA D 250 40 a00 20a u 5 4i 35 300 400 306 2BB if g 16a 288 fe 200 200 W0 0 200 200 0 200 DATA nT DATA nT DATA nT Figure 2 The panel on the left is the total field anomaly on the surface produced by the model in Figure 1 under an inducing field with declination and inclination noted above The panels on the right are the three component anomalies in three vertical boreholes Different components are identifie by the labels beside the curves EASTING m 1 3 Inversion Methodology for Ver 3 1 Notes regarding modifications for Ver 4 0 are
13. provided below Let the set of extracted anomaly data be d lay dy tand the susceptibility of cells in the model be ESE aes Sway The two are related by the sensitivity matrix d Gk 7 The matrix has elements gij which quantify the contribution to the tn datum due to a unit susceptibility in the ie cell see Section 1 2 The program MAGSENSD performs the calculation of the sensitivity matrix which is to be used by the Subsequent inversion The sensitivity matrix gives the forward mapping from a model to data during the entire inverse process We will discuss its efficient representation via the wavelet transform in a separate section The inverse problem is formulated as an optimization problem where an objective function of the model is minimized Subject to the constraints in eq 7 For magnetic inversion the first question that arises concerns definition of the model We choose amp as the model since the anomalous field is directly proportional to the susceptibility This is the choice for the inversion program MAGINV3D For generality we introduce the generic symbol m for the model element Having define a model we next construct an objective function which when minimized produces a model that is geophysically interpretable The details of the objective function are problem dependent but generally we need the flexibility to be close to a reference model mo and also require that the model be relatively smooth in three spatial dir
14. the file structure of model sus Susi 1 1 SUS1 1 2 SUS1 1 NV SUS1 2 1 SUSi 4 k SUSNN NE NV Each susj j k is the susceptibility 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 1 0 to avoid confusion with the other model elements 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 421 1 1 W ZNN 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 sus however they can be
15. the interpreter has high confidenc in the reference model at a particular region he can specify ws 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 To perform a numerical solution we discretize the model objective function in eq 8 using a finite difference approximation on the mesh defining the susceptibility model This yields m m m mo wiw WIW Wy W WIW m mo m ma WEW m m mo 9 W m mp where m and mo are M length vectors The individual matrices Ws Wx Wy and Wz are straight forwardly calculated once the model mesh and the weighting functions w r and ws Wx Wy Wz are defined The cumulative matrix WoW is then formed The next step in setting up the inversion is to define a misfit measure Here we use the 2 norm measure oby g r W Ga d 10 For the work here we assume that the contaminating noise on the data is independent and Gaussian with zero mean Specifying Wg to be a diagonal matrix whose i element is 1 7 where f is the standard deviation of the t datum makes gq a chi squar
16. wavelet compression Parameters controlling the implementation of this compression are available for advanced users The research 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 Services Kennecott Exploration Company Newmont Gold Company Noranda Exploration Placer Dome and WMC Since then improvements have been implemented as time and resources permit especially in the context of the ICIS consortium project 2005 2006 which supported development of version 4 0 The theoretical framework for MAG3D is provided in the following papers see the UBC GIF website publications page for details Li Y and Oldenburg D W 1996 3 D inversion of magnetic data Geophysics 61 no 02 394 408 Li Yaoguo and Oldenburg Douglas W 1998 Separation of regional and residual magnetic field data Geophysics 63 no 02 431 439 Li Y and Oldenburg D W 2000 Joint inversion of surface and three component borehole magnetic data Geophysics Geophysics 65 2 pp540 552 Two short papers including examples of applying MAG3D in mineral exploration contexts are Cost effectiveness of g
17. will be set to 1 automatically 1 sensitivity matrix will be accessed from disk when needed UPDATES FOR Ver4 0 IN THE APPENDIX Input files Input files can be any file name If there are spaces in the path or file name you MUST use quotes around the entire path filename obs mag input data file The file must specify the standard deviations of the error By definition these are greater than zero maginv3d mtx sensitivity matrix and depth weighting functionfile calculated by MAGSEN3D ini sus jnitial model stored in the same way as model sus If null is entered the default value of 0 001 is used For a constant initial model enter a value ref sus 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 sus Version 4 0 only Three options are possible here either a file name two constants or NULL bounds sus file name for a two column file with bounds for each cell 2 numbers two constants representing lower and upper bounds respectively NULL default susceptibility bounds of 0 0 and 1 0 will be used 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 maginv3d log The log file containing more detailed information for each iteration and summary of the inversion
18. 108 282 380 388 228 Lh bc NORTHING m E i i i amp a Le La 308 J00 200 10E 8 ide 2a ode EASTING m O36 ae 028 Mat 020 1g Diz O04 total field surface data and three component borehole data The distance weighting function is used in this inversion This model clearly define both prisms and whose true positions are indicated by the solid whire lines 1 8 References Bhattacharyya B K 1964 Magnetic anomalies due to prism shaped bodies with arbitrary magnetization Geophysics 29 517 531 Li Y and Oldenburg D W 1996 3 D inversion of magnetic data Geophysics 61 394 408 Li Y and Oldenburg D W 2000 Joint inversion of surface and three component borehole magnetic data Geophysics 65 pp540 552 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 490 493 Sharma P V 1966 Rapid computation of magnetic anomalies and demagnetization effects caused by bodies of arbitrary shape Pure Appl Geophys 64 89 109 UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 39 00 MAG3D manual Elements of the program 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 ob
19. 188 288 DEPTH mi JAE age aigi 308 200 i00 A i188 28H 388 oe 388 Ov tae o CEL 288 alee a iag p 010 Ra D Dig E g D one S 18E Hi 208 30E 3BF 208 1BB B 188 200 300 EASTING im three component data shown above This simulates a more realistic situation since in many practical applications only the total field anomaly can be extracted accurately from borehole measurements The recovered model images both prisms The success of this inversion demonstrates that single component borehole data can provide the information that is complementary to surface data when multicomponent data are not Example 2 The true model is illustrated to the right The figure showsthe relative sizes and positions of seven susceptible bodies buried in a non susceptible background The six smaller bodies are placed at shallow S depths to simulate small scale anomalies and the large block is placed at a greater depth to generate a broad anomaly over which the small anomalies are superimposed mi g 1000 1500 The surface topography above this model is shown in the next figure The elevation of the surface ranges mostly between 0 and 125 m with a few points reaching 150 m available 150 100 i i TORA F PREP HER Simulated aeromagnetic data calculated at a constant terrain clearance of 75 m are shown next The data are located on a grid with a 50 m spacing in both directions This
20. Furthermore the inversion program assumes that the anomalies are produced by a positive susceptibility contrast distribution It is crucial that the data be prepared as such Examples of obs mag file The following two example data files The firs example file specifie a set of total fieldanomaly data and the second example file provides a set of multi component borehole data File 1 single component data surface data 45 0 45 0 50000 0 incl decl geomag 45630 45 0 aincl adecl direction of anomaly idir 441 of data 0 00 0 00 1 0 0 174732E 02 0 598678E 01 0 00 50 00 1 0 0 265550E 02 0 613890E 01 0 00 100 00 1 0 0 311366E 02 0 629117E 01 1000 00 900 00 1 0 0 109595E 01 0 530682E 01 1000 00 950 00 1 0 0 902209E 01 0 523738E 01 1000 00 1000 00 1 0 0 397501E 01 0 518496E 01 File 2 multi component data borehole data l 65 00 25 00 50000 00 incl decl geomag 65 00 25 00 0 aincl adecl anomaly idir 144 of data 12 50 137 50 12 50 0 0 0 0 0 134759E 03 0 200000E 01 12 50 137 50 37 50 0 0 0 0 0 162606E 03 0 200000E 01 12 50 137 50 62 50 0 0 0 0 0 165957E 03 0 200000E 01 237 50 12 50 337 50 90 0 0 0 0 662445E 02 0 200000E 01 237 50 12 50 362 50 90 0 0 0 0 693134E 02 0 200000E 01 237 50 12 50 387 50 90 0 0 0 0 608605E 02 0 200000E 01 FILE model sus This file contains the cell values of the susceptibility model The susceptibility must have values in SI units The following is
21. It s format depends upon which mode was run chifact constant tradeoff or GCV and information is more useful to those who understand in detail the theory of the processing algorithms e Clicking the model button opens the inversion model in a new MeshTools3D window MeshTools3D has it s own user manual e Clicking the predicted data button opens a new data viewing window with the observed data at the top and the predicted data on the bottom Predicted data were calculated by the inversion program on the final model You can choose to see errors or misfit instead of predicted data See the gm data viewer exe help notes for more details Data file Specify the input data file here See the MAG3D user manual for data file format You can enter the file name including full path or drag and drop from Windows Explorer The View Data button opens the data set in a new window by running the gm data viewer exe program Mesh The mesh which defines the model is specified here See the user manual section 2 2 FILE mesh for details If a mesh has not been generated a default mesh can be created using the Create a mesh x Create mesh button A data file must already be specified This dialogue appears peun Save Cancel Specify the size of square cells in the core region under the data set by Eltvs oriet the top of the mesh entering their width in units of metres Cell depths will be half this width Specify th
22. MAG3D A Program Library for Forward Modelling and Inversion of Magnetic Data over 3D Structures VERSION 4 0 Developed under the consortium research project JOINT COOPERATIVE INVERSION OF GEOPHYSICAL AND GEOLOGICAL DATA UBC Geophysical Inversion Facility Department of Geophysics and Astronomy University of British Columbia Vancouver British Columbia May 2005 UBC Geophysical Inversion Facility 1998 2005 MAG3D manual Program library On this page Description Package contents Licensing Installation Description MAGS3D is a program library version 4 0 as of August 2005 for carrying out forward modelling and inversion of surface airborne and or borehole magnetic data in the presence of a three dimensional Earth The program library carries out the following functions Forward modelling of the magnetic field anomaly response to a 3D volume of susceptibility contrast Data are assumed to be the anomalous magnetic response to buried susceptible material not including Earth s ambient field The model is specified using a mesh of rectangular cells each with a constant value of susceptibility and topography is included The magnetic 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 Assumptions O This code assumes susceptibilities are small This means results will b
23. 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 39 00 MAG3D manual Introduction MAGS3D is a program library for forward modelling and inversion of magnetic data over 3d structures It was developed under the consortium research project OINT COOPERATIVE INVERSION OF GEOPHYSICAL AND GEOLOGICAL DATA by the UBC Geophysical Inversion Facility MAGS3D Ver 4 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 3 1 The user can input lower and upper bounds This feature can be important in helping control the range of acceptable models 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 A file sensitivity txt is output after running magsen3d exe It contains the average sensitivity for each cell This file can be used for depth of investigation analysi
24. ained next Depth weighting for borehole data For data sets that contain borehole measurements the sensitivities do not have a predominant decay direction therefore a weighting function that varies in three dimensions is needed We generalize the depth weighting used in surface data inversion to form such a 3D depth weighting function called distance weighting ay l4 N 1 l cay wF 5 aa l j Hle M 14 WAVE ayy iiy Ro where 3 3 0 and Vj is the volume of jt cell Ri is the distance between a point in and the t observation and Ro is a Small constant used to ensure that the integral is well define 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 borehole data it is necessary to use this more general weighting Note This weighting function is also advantageous if surface data with highly variable observation heights are inverted The weighting function is directly incorporated in the sensitivity file generated by program MAGSEN3D This program allows user to specify whether to use the depth weighting or the distance weighting for surface data When borehole data are present however distance weighting must be used 1 5 Wavelet Compression of Sensitivity Matrix The two major obstacles to the solution of large scale magnetic inversion problem are the large amount of memory required for storing th
25. all on one line or broken up over several lines Since the weights for a derivative term are applied 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 FILE bounds sus This file contains the cell values of the lower and upper bounds on the sought model It is only required optionally by maginv3d The bounds have the same dimension as the susceptibility contrast The following is the file structure of bounds sus lk 1 1 uki 1 1 Let 2 URL 4 2 1k1 1 NvV uk1 1 NV lk 2 1 uki 2 1 1ki 4 k uki 4 k LKNN NE NV UXKNN NE NV Parameters are Iki j k is the lower bound on cell ij k Uki j k is the upper bound on cell i j k 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 or
26. and nested inside each barrier iteration is a 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
27. arrier method the bounds constraints are implemented as a logarithmic barrier term The new objective function is given by A g gah 2 AY In x a t In c 12b gel where 4 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 A and as approaches zero the sequence of solutions approach the solution of eq 11b 1 4 Depth Weighting Depth weighting for surface or airborne data It is a well known fact that static magnetic data have no inherent depth resolution A numerical consequence of this is that when an inversion is performedwhich minimizes mF dv subject to fitting the data the constructed susceptibility is concentrated close to the observation locations This is a direct manifestation of the kernel s decay with the distance between the cell and observation locations Because of the rapidly diminishing amplitude the kernels of magnetic data are not sufficient to generate a function that possess significant structure at locations that are far away from observations In order to overcome this the inversion needs to introduce a weighting to counteract this natural decay Intuitively such a weighting will approximately cancel
28. at surface topography optional 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 and or borehole beta znot parameters defining the depth weighting function When iwt 1 beta and znot are used as 7 and Zo to define the depth weighting according to eq 13 background When iwt 2 beta and znot are used as 7 and Zo to define the distanc e weighting according to eq 14 background If null is entered on this line line 5 then the program sets beta 3 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 3 0 Note this is different than the value used by GRAV3D 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 i
29. ce of the inversion Depending on how the inversion 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
30. d may be deleted once the work is completed 2 Afile sensitivity txt is output after running magsen3d 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 Added for mag3d version 4 0 3 4 MAGINV3D This program performs 3D magnetic inversion Command line usage is maginv3d maginv3d inp For a sample input file type maginv3d inp Format of the control file maginv3d inp is as follows irest mode par tolc obs mag maginv3d mtx ini sus ref sus bounds sus Ys a CF oF w dat idisk Control parameters irest 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 mode 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 par tolc 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 p
31. dered 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 negative value e g 1 0 to avoid confusion UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 39 00 MAG3D manual Running the programs Contents of this page 4 general files MAGFORS3D performs forward modelling MAGSENSD calculates sensitivity and the depth weighting function MAGI NV3D performs 3D magnetic inversion MAGPRE3D 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 In a MS WindowsXX operating environment codes are best run using the Graphical User Interface GUI See the separate GUI instructions 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 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
32. e minimization from crossing over to the infeasible region The method solves a sequence of nonlinear minimizations with decreasing A and as approaches zero the sequence of solutions approach the solution of eq 11 This methodology provides a basic framework for solving 3D magnetic inversion 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 tradeoff parameter that ultimately determines how well the data will be fit and the logarithmic barrier method to obtain the solution with positivity 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 tradeoff 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 Notes regarding MAG3D Ver 4 0 In MAG3D versions 3 1 and earlier the condition of positivity in the model was imposed by implementing a modification to the model objective function In Version 4 0 a more general condition can be imposed such that lower and upper bounds on values of susceptibility in the model can be defined This is useful because often there are well defined bounds on the susceptibility contrast based on direct sampling or other geological information The procedure for i
33. e different inclinations and declinations aincln and adecin specifying the projection direction of the anomaly In this case Aincl and Adecl should be equal to Incl and Decl respectively 1 single component data set All observations have the same inclination and declination of the anomaly projection Aincl Adecl If idir is missing it is assumed to be equal to 1 ndat number of observations When single component data are specified the number of observations is equal to the number of data locations When multi component data are specified the number of observations will exceed the number of data locations For example if three component data are specifie at N locations the number of observations is 3N En Nn easting northing and elevation of the observation measured in meters Elevation should be above Elevn the topography for surface data and below the topography for borehole data The observation locations can be listed in any order aincin adecin inclination and declination of the anomaly projection for observation n This is used only when idir 0 The brackets indicate that these two field are optional depending on the value of idir Magn magnetic anomaly data measured in nT Errn standard deviation of Magn This represents the absolute error It CANNOT be zero or negative NOTE It should be noted that the data Magn are extracted anomalies which are derived by removing the regional from the fieldmeasurements
34. e elevation in metres for the top of the mesh This must be above the topograpy There are remarks on topograpy and meshes in several places in the manual including section 2 2 FILE mesh and section 2 2 FILE topo dat The default mesh created by this dialogue will include 3 padding cells around the core region These cells increase in size at greater distances from the core Wavelet compression parameters The compression used to speed up dense matrix multiplication is specified here e Use default to use parameters that have worked well in nearly all problems This is the recommended mode e Specify user to set your own compression parameters Details are in the user manual section 3 3 The value of daub specifies the type of wavelet as per section 3 3 Selecting relative reconstruction error is equivalent to setting itol 1 as per section 3 3 Selecting relative threshold is equivalent to setting itol 2 as per section 3 3 The value is set based upon which of these two choices is made e Selecting none will result in no compression and probably rather long computation times Topography file No topography file is necessary if the Earth s surface under the data set is assumed flat If topography is significant then it should be specified as per the user manual section 2 2 FILE topo dat There are remarks on topograpy and meshes in several places in the manual including section 2 2 FILE mesh and section 2 2 FILE topo dat
35. e sensitivity matrix and the CPU time required for the application of the sensitivity matrix to model vectors The MAG3D program library overcomes these difficultie 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 MAG3D library Each row of the sensitivity matrix in a 3D magnetic 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 coefficient are nearly or identically zero When the coefficient with small magnitude are discarded the process of thresholding the remaining coefficient 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 so
36. e wrong when susceptibilties are high enough to cause self demagnetization O There is no method for incorporating remanent magnetization in this code Inversion of surface airborne and or borehole magnetic data to generate 3D models of susceptibility contrast O 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 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 susceptibility 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 susceptibility contrast in any cell as of version 4 0 O 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
37. ections Here we adopt a right handed Cartesian coordinate system with positive north and positive down Let the model objective function be Guw Fi m r mo fa GT Omi cs f wawi Tiir moly du F ay wad a foe E fw Do lV 7 r i oy 1 Oz Y where the functions ws Wx Wy and wz are spatially dependent while ts x My and tz are coefficient which affect the relative importance of different components in the objective function Here the function wr is a generalized depth weighting function The purpose of this function is to counteract the geometrical decay of the sensitivity with the distance from the observation location so that the recovered susceptibility is not concentrated near the observation locations The details of the depth weighting function will be discussed in the next section The objective function in eq 8 has the flexibility to construct many different models The reference modelmo may be a general background model that is estimated from previous investigations or it could be the zero model The reference model would generally be included in the first component of the objective function but it can be removed if desired from the remaining terms often we are more confident in specifying the value of the model at a particular point than in supplying an estimate of the gradient The relative closeness of the final model to the reference model at any location is controlled by the function ws For example if
38. ed variable distributed with N degrees of freedom Accordingly Ely provides a target misfit for the inversion The inverse problem is solved by finding a model m which minimizes m and misfits the data by a pre determined amount Since the susceptibility is positive by definition we also need to impose the constraint that all model elements be positive Thus the solution is obtained by the following minimization problem minimize a py 11 subject to m gt 0 where is a tradeoff parameter that controls the relative importance of the model norm and data misfit When the standard deviations of data errors are known the acceptable misfit is given by the expected value 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 tradeoff 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 positivity constraint is implemented as a logarithmic barrier term The new objective function is given by Ad eh o d lem 245 Infamy 12 j l where 4 is the barrier parameter and the tradeoff parameter is fixed during the minimization As the name suggests the logarithmic barrier term forms a barrier along the boundary of the domain of feasible solutions and prevents th
39. eophysical 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 magnetic surveys The MAG3D library WindowsXxX or Linux platforms consists of three major programs and one utility O MAGFOR3D performs forward modelling MAGSENS3D calculates sensitivity and the depth weighting function O MAGINV3D performs 3D magnetic inversion MAGPRE3D multiplies the sensitivityfile by the model to get the predicted data eA oe user interface is supplied for the WindowsXxX platforms only Facilities include 2 gt MAG3D GUI EXE a primary interface for setting up 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 O MESHTOOLS3d a utility for displaying resulting 3D models as volume renderings Susceptibility volumes can be sliced in any direct
40. fy the file name if a susceptibility model defined on the mesh in use for inversion 1s available for use in the model objective function as a reference Model objective function coefficients These are the parameters that control how much emphasis is placed on each part of the model objective function e Select default to specify default alpha parameters of alphas 0 0001 alphax alphay alphaz 1 0 e Alternatively specify your own values using either alpha parameters or length scales and following these guidelines Influence of Alphas on results x Ey 2 i Consider the ratios and a Larger ratios result in smoother models As a rule of thumb s tt E keep this table in mind O Q Fe l Structure is penalized Constant reference models has little effect a a a O Q Na I gt Smallest term dominates so models are rough a a a To estimate values for the s for a specific case start by defining two length scale terms as follows L J L and L a G Ya Then consider these rules of thumb to help choose the s In general keep Ly and y approximately the same as the shortest array separation and keep Ly and y approximately equal to Lz Also the following relations are useful LY or y gt mesh cell width Lz gt mesh cell thickness LY or y lt total width of the mesh Lz lt total depth of the mesh Top of MAG3D GUI help The user manual for MAG3D
41. hows the presence of the shallow prism at the correct location but it does not give a clear indication of a separate deeper prism There is only a broad zone of low susceptibility extending from the high susceptibility zone The vertical extent of the anomaly is not well defined Click the button to see the susceptibility model recovered by inverting the three component borehole data The model shows two regions of high susceptibility at locations corresponding to the true prisms and the recovered depths agree well with the true depths However the amplitude of the shallow prism is small and there is no clear separation from the deep prism Despite this the model provides a good result considering that there are only three widely separated boreholes and that the inversion has no explicit information regarding where to place the magnetic material Click to see the susceptibility model obtained from the joint inversion of the surface and three component borehole data This model combines the merits of the models from individual inversions Both target prisms are well define in horizontal and vertical locations and their amplitudes are comparable to those of the true model This model provides the best representation of the true model D This is the total field Anomaly Inversion It shows the model recovered by inverting the surface and borehole total fieldanomaly The total field anomaly in the boreholes is first computed from the
42. ion 65 and D 25 is assumed The E 20a total field anomaly was calculated on the surface at an interval of 25 pi m along east west lines spaced 100 m apart resulting in a total of 30e 175 data Three component data in easting northing and vertical jae Daath directions are calculated in the boreholes There are 16 observation 200 200 174 A 12A 24g 3pg Mao locations spaced 25 m vertically in each hole and the total number of 200 ee data in three holes is 144 DEM 208 ae Gaussian noise having a standard deviation of 2 nT has been added to m D fe all data and the resulting simulated observations are shown in the g 128 l k O 010 next figure where the surface data are shown as gray scale contours l a p 010 and the three component borehole data are plotted as functions of E k depth For inversion we use a m region and divide it into 24x24x16 B_iagg 2 00a cubic cells of 25 m on a side This yields a total of 9216 cells For all sI the inversions presented below the distance weighting is used 20A 30a 380 200 188 A 1B 200 3p8 EASTING mi 200 aB ed ia T aa oni x 160 pi i a 120 5 Bo 108 40 288 aie 4 308 2AB 2p 106 F 12A Z200 3g 200 O 200 200 O 200 200 O 200 EASTING im DATA nT DATA nT DATA nT This figure shows the susceptibility model recovered from inverting the surface data alone The true positions of the prisms are indicated by the white lines This model clearly s
43. ion or isosurface renderings can be generated Documentation is elsewhere via the menu to the left Example data sets and excercises are provided on the IAG 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 magnetic 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 versions 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 MAG3D version 4 0 educational version For users with a 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 mag3d 2 No further installation is necessary 3 Follow instructions in one of this CD ROM s exercises or refer to the program documentation 4
44. jective function is specifie by choosing gt 0 0001 t y y 7 1 anda zero reference model PE The 3D weighting functions are all set to unity The inversion is performed by setting the target misfit to the expected value of 3 600 and executing MAG3D with mode 1 The inversion uses 60 Mb of memory and lasts 146 minutes on the SUN Py f TT Morn Sparc20 workstation or 110 minutes 3000 on a 233 MHz MMX Pentium PC having 64 Mb of memory 1995 2500 Thus the entire procedure from the calculation and compression of the sensitivity to the inversion requires about 392 minutes on the SUN workstation or 285 minute on the PC 2000 Northing im The last figure displays the susceptibility obtained from the 2000 200 1000 000 00 s00 inversion The model is shown in six 3000 plan sections The different bodies in 2500 the true model are well imaged In particular the large block at depth is _ pg cleary visible The blank area in the section at Z 87 5 indicates the borthing m a0 00 4001 2500 2000 1500 1000 500 1006 1500 2000 2500 3000 000 2000 2000 1500 1000 200 io 1500 Z000 2500 3000 2000 2300 2000 1500 a aC H 1000 1000 1500 1300 20170 25010 3000 2000 2300 3000 0 068 0 067 0 054 0 047 0 040 0 034 0 027 0 020 region above the topographic surface
45. ld Hy J KHy 1 where Ho Bo to and 0 is the free space magnetic permeability This essentially ignores the self demagnetization effect by which the secondary field reduces the total inducing field within the susceptible region and results in a weaker magnetization than that given by eq 1 The anomalous field produced by the distribution of magnetization 7 is given by the following integral equation with a dyadic Green s function F 2t Bal F i YFF F dJdv 2 where is the position of the observation point V represents the volume of magnetization The above equation is valid for observation locations above the earth s surface It is also valid in the boreholes provided we assume that the magnetic permeability is H When the susceptibility is constant within a volume of source region the above equation can be written in matrix form as Tu Tiz Tia Ji Gq 19 1322 199 ERG Iy Ty Faq j 3 nou T Ho i where Ti is given by dr f Oa Aaj r r i where x X2 and x3 represent x y and z respectively The expressions of Tij for a cuboidal source volume can be found in Bhattacharyya 1964 and Sharma 1966 Here we assume that the effect of borehole cavity can be neglected Since T is symmetric and its trace is equal to 1 when the observation is inside the cell and is 0 when the observation is outside the cell only five independent elements need to be calculated Once T is formed the mag
46. longer linear and the problem becomes more complicated Inverting data in the presence of very high susceptibilities is still a topic of research and MAG3D Version 4 0 does not allow for high susceptibilities in the solution 2 Upper and lower values Use this option if you want to set a single upper bound value and a single lower bound for all cells 3 Using a file Each cells can have individual upper and lower bounds if a bounds sus file is specified See the MAG3D manual Chapter 2 Initial model This is the first susceptibility model the program works with The outcome of the inversion should not depend much on this model but if the initial model is close to the final one then convergence can be expected much more quickly e Selecting default causes the program to use a default initial model of zero susceptibility e Select Value and enter a susceptibility if the host rocks are not non magnetic e Select File and specify the file name if a susceptibility model defined on the mesh in use for inversion is available Reference model This is the reference susceptibility model that is part of the model objective function The inversion will try to minimize the difference between the final model and this reference e Selecting default causes the program to use a default reference model of zero susceptibility e Select Value and enter a susceptibility if you know the background susceptibility of host rocks e Select File and speci
47. lvable 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 G G 15 where G is called the transformed matrix The thresholding is applied to individual rows of by the following rule to form the sparse representation gije gig Oy i o 1e N 16 0 gig b where fiis the threshold level and 4 and j are the elements of G and respectively The threshold level i are determined according to the allowable error of the reconstructed sensitivity which is measured by the ratio of norm of the error in each row to the norm of that row ri i It can be evaluated directly in the wavelet domain by the following expression 43 gt gj lala yg J 1 N 17 Here the numerator is the norm of the discarded coefficients For each row we choose jj such that ri i r where r is the prescribed reconstruction accuracy However this is a costly process Instead we choose a representative row lo and calculate the threshold level fio This threshold is then used to define a relative threshold Big ac Gig i The absolute threshold level for each row is obtained by 6 e max g E51 N 18 J The program that implements this com
48. measurements through numerical processing one obtains the distribution of the anomalous field due to the susceptible material Very often the susceptible materials underground possess a certain amount of natural remanent magnetization In this program library however we make the assumption that no remanent magnetization is present and restrict our attention to induced magnetization The data from a typical magnetic survey is a set of magnetic field measurements acquired over a 2D grid above the surface or along a number of boreholes within the volume of interest These data are first processed to yield an estimate of the anomalous field due to the susceptible material in the area The goal of the magnetic inversion is to obtain from the extracted anomaly data quantitative information about the distribution of the magnetic susceptibility in the ground Thus it is assumed that the input data to the inversion program is the extracted residual anomaly and the programs in the library are developed accordingly 1 2 Forward Modelling General Formulation For a given inducing field By the magnetization T depends upon the susceptibility through a differential equation However to the first order approximation when the actual susceptibility is very small as is most often the case with material encountered in mineral explorations the magnetization 7 IS proportional to the susceptibility and is given by the product of susceptibility with inducing magnetic fie
49. mes Also the filename extensions are not important Many prefer to use the txt filename convention so that files are more easily read and edited in the Windows environment The files are Pe a Nou FILE mesh 3D mesh defining the discretization of the 3D model region topo dat specifies the surface topography obs loc specifies the inducing field parameters anomaly type and observation locations obs mag specifies the inducing field parameters anomaly type observation locations and the observed magnetic anomalies with estimated standard deviation model sus susceptibility model file w dat contains the 3D weighting functions bounds sus optional file containing values for upper and lower bounds new for version4 0O 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 AENE Ani Ano nnn Ayi Ave Ayvnv 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 mag and in topo dat see the relevantfile for description Aen Cell widths in the easting direction from W to E nn Cell widths i
50. mplementing upper and lower bounds on model values is the same as that used in GRAV3D Version 2 0 and later In MAG3D the positivity constraint is no longer necessary if upper and lower bounds are defined However if the user chooses to not define upper and lower bounds the program employs default bounds for susceptibilities in every cell of 1 0 and 0 0 respectively S I units While it is true that some rocks have susceptible greater than 1 0 MAG3D Version 4 0 still assumes small susceptibilities as this code has done since the original version When there are very high susceptibilities the relation between incident and induced magnetization is no longer linear and the problem becomes more complicated Inverting data in the presence of very high susceptibilities is still a topic of research and MAG3D Version 4 0 does not allow for high susceptibilities in the solution As a result of implementing upper and lower bounds equations 11 and 12 above are different for MAG3D Version 4 0 Instead of imposing positivity the solution is obtained by the following minimization problem which replaces equation 11 minimize a Hm a O 11b subject to Emn Ama where imin and max are vectors containing the lower and upper bounds on the model values and is the vector containing model values As before we use the primal logarithmic barrier method with the conjugate gradient technique as the central solver In the logarithmic b
51. mponent data borehole data 65 00 25 00 50000 00 incl decl geomag 65 00 25 00 0O It aincl adecl anomaly idir 144 of data 12 50 1I37 50 12 50 040 040 12 50 137 50 37 50 0 0 0 0 12 50 137 50 62 50 0 0 0 0 237 50 12 50 337 50 90 0 0 0 237 50 12 50 362 50 90 0 0 0 237 50 12 50 387 50 90 0 0 0 FILE obs mag This file is used to specify the inducing fieldparameters anomaly type observation locations and the observed magnetic anomalies with estimated standard deviation The values of parameters specifying the inducing field anomaly type and observation locations are identical to those in obs loc The output of the forward modelling program MAGFORSD has the Same structure except that the column of standard deviations for the error is omitted The following is the file structure of obs mag comments l Incl Decl geomag Aincl Adecl idir ndat El N1 Elev1 famine LL adecll Magl Erri E2 N2 Elev2 ainc1l12 adec12 Mag2 Err2 Endat Nndat Elevndat ainclndat adeclndat Magndat Errndat Parameters are top lines beginning with are comments comments Incl jnclination and declination of the inducing magnetic field The declination is specified east with Decl respect to the northing used in the mesh obs loc and obs mag files geomag strength of the inducing fieldin nT Aincl jnclination and declination of the anomaly projection Adecl idir Q multi component data set Observations hav
52. n 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 consists 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 no magnetic material 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
53. netic anomaly and its projection onto any direction of measurement is easily obtained by the inner product with the directional vector The projection of the field B onto different directions yields different anomalies commonly obtained in the magnetic survey For instance the vertical anomaly is simply a the vertical component of whereas the total field anomaly is to first order the projection of B onto the direction of the inducing fiels By Borehole Data In a borehole experiment the three components are measured in the directions of local coordinate axes 1 Y1 Z1 define according to the borehole orientation Assuming that the borehole dip fis measured downward from the horizontal surface and azimuth is measured eastward from the north a commonly used convention has thezj axis pointing downward along borehole x1 axis pointing perpendicular to the borehole in the direction of the azimuth Theyj axis completes the right handed coordinate system and is 90 clockwise from the azimuth and perpendicular to the borehole Based upon the above definition the rotation matrix that transforms three components of a vector in the global coordinate system to the components in the local coordinates is given by cosigsin sinwsn cos R amg cosh 0 5 cosacos sinweosf sinf If a vector is define in local coordinates as f fa In and in global coordinates as q1 92 gay then the following two relations hold 24 fg BY
54. ons using MAG3D Ver3 2 and Ver4 0 are compared in the figure below for a moderately complicated problem This was a real example involving 1564 data points and 417 582 cells in the model and moderate topography The Version 4 0 code is significantly faster in all cases Speed also depends strongly on the computer Times when using two computers are shown 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 More complex problems include topography complicated distributions of susceptibility weighting functions reference models bounds etc 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 you set up all parameters for the inversion including data constraints regularization compression etc MAG3D computation times on two computers for both versions 3 2 and 4 0 of data 1564 of cells 417 582 CPUT clock 3 0 Ghz cache 0 5Mb CPU clock 3 2 Ghz cache 1 0Mb Waid Vat Waid Vd 0 26 30 0 20 03 0 27 14 0 11 54 4 25 41 1 27 14 3 14 10 1 00 22 The remaining sections of the manual are on separate pages
55. pression procedure is MAGSEN3D 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 When both surface and borehole data are present two different relative threshold levels are calculated by choosing a representative row for surface data and another for borehole data For experienced users the program also allows the direct input of the relative threshold levele 1 6 Choice of Tradeoff Parameter The choice of the tradeoff parameter ultimately depends upon the magnitude of the error associated with the data The inversion of noisier data requires heavier regularization thus a greater value of ji is required In this section we discuss the various implementations for the choice of jt in the MAG3D library If the standard deviation associated with each datum is known then the data misfit define by eq 10 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 js Because of the positivity constraint our problem is nonlinear Thus for each jt a nonlinear solution using a logarithmic barrier method must be obtained This is computationally demanding and we therefore have developed the following strategy to red
56. ractical applications the estimate of data error is often begin not available Then the degree of regularization hence the value of j 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 involved and numerically PY i expensive to implement However applying GCV on the ines Mm P bq Be aa 3D magnetic inversion without positivity still produces a find po by line search ee iil reasonable estimate of the data error That error can serve as a Starting 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 jt obtained in this manner directly in the final inversion which has the positivity 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 MAGINV3D has implemented this approach as the third method for choosing the tradeoff parameter f dat hotm m by log barrier method Figure 3 illustrates the structure of the program MAGINV3D It has three options for determining the tradeoff parameters The controlling parameter is mode When mode 1 the line search based on kno
57. roduct 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 true 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 2 only par is used When mode 3 neither par nor tolc are used However the third line should always have two values s coefficient for the smallest model component Fe coefficient for the derivative in the easting direction fn coefficient for the derivative in the northing direction v coefficient for the derivative in the vertical direction If null is entered on the eighth line then the above four parameters take the following default values alphas 0 0001 alphae alphan alphav 1 0 None of the alpha s can be negative and they cannot be all equal to O at the same time idisk 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
58. s or for use in designing special model objective function weighting A file maginv3d_nopos sus is output during the first part of the inversion It contains the susceptibilities without the bounds constraints imposed A file maginv3d_ XX sus is output after each beta iteration In maginv3d exe the user can enter either alpha values or length scales The reference model is 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 magsen3d exe are now more efficient Changes to run time files for MAG3D Version 4 0 To implement the bound constraints the input file for maginv3d has an extra line where the user can choose one of three options a filename two column file with lower and upper bounds for each cell 2 numbers constant values for lower and upper bounds NULL default bounds of 0 and 1 There are several additional output files created see description above Changes to this manual for MAG3D Version 4 0 This introductory note Additions to section 1 4 Inversion Methodology are included in the Appendix Changes to the format for input file Maginv3d inp are described in the Appendix The optional bounds sus file that can be used to specify lower and upper bounds for susceptibility values in each model cell is described after page 22 Notes on computation speed Time to run computati
59. s the Haar wavelet and daub2 is the Daubechies 4 wavelet The Daubechies 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 Example of mesh obs nois null 2 null daub2 1 0 05 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 The detailed explanation of threshold level and reconstruction error can be found in the Background section 1 5 of this manual magsen3d inp control file mesh file data file topography iwt 1 depth 2 distance beta znot null wavelet type itol eps null Two output files 1 maginv3d 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 an
60. s3d loc contains the locations for the survey obs3d mag contains the observations model sus contains the cell values for the susceptibility mode w dat contains special weightings that alter the type of model produced in the inversions bounds sus specifies upper and lower bounds for susceptibility values 2 1 Introduction The MAG3D library consists of three major programs and one utility aa ale MAGFORS3D performs forward modelling MAGSENS3D calculates sensitivity and the depth weighting function MAGINV3D performs 3D magnetic inversion MAGPRE3D multiplies the sensitivityfile by the model to get the predicted data This rarely used utility multiplies a model by the sensitivity matrix in maginv3d 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 maginv3d 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 MAG3D Programs There are seven general files which are used in MAG3D version4 0 All are in ASCII text format I nput files can have any name you like Only program output files have restricted file na
61. simulates a normal data set that has been decimated it is NOT recommended that gridded data are used for inversion The inducing field is assumed to be in the direction 65 and D 25 and total field anomalies are calculated The data have been contaminated by independent Gaussian noise having a zero mean and a standard deviation of 5 nT The data show six anomalies due to the shallow blocks but they provide no indication about the presence of the deep block To invert these data a model region of 3 2 km by 3 2 km by 1 5 km is used The top of the mesh is placed at the elevation of 125 m The cell width is 50 m in both Northing mn 3000 200 2000 a al al 1500 1000 300 0 300 1000 1500 Easting m 2000 2000 3000 62 5 20 0 102 0 185 0 horizontal directions and the thickness varies from 25 m near the surface to 100 m at the bottom After the surface topography is discretized onto the mesh the resulting model contains a total of 110 000 cells The corresponding sensitivity matrix requires more than 1 5 Gb to store and that is beyond the memory limit of many workstations When compressed using the Daubechies 4 wavelet at a reconstruction accuracy of 5 a compression ratio of 76 is achieved The compressed sensitivity matrix requires on 43 5 Mb of storage The sensitivity calculation takes 245 minutes on a SUN Sparc20 workstation or 175 minutes on a 233 MHz MMX Pentium PC The model ob
62. the decay and give cells at different locations equal probability to enter into the solution with a non zero susceptibility For surface data the sensitivity decays predominantly in the depth direction Numerical experiments indicate that the function of the form z j zo 3 closely approximates the kernel s decay directly under the observation point provided that a reasonable value is chosen for Zo The value of 3 in the exponent is consistent with the fact that to first order a cubic cell acts like an dipole source whose magnetic field decays as inverse distance cubed The value ofzo can be obtained by matching the function 1 z zo 3 with the field produced at a observation point by a column of cells Thus we use a depth weighting function of the form 41 2 u 1 dz wry Gta l M_ 13 for the inversion of surface data where 7 3 0 and rj is used to identify the res cell and 4z 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 susceptibility model constructed by minimizing a model objective function in eq 8 subject to fitting the data places the recovered anomaly at approximately the correct depth Note that if the data set involves highly variable observation heights the normal depth weighting function might not be most suitable Distance weighting used for borehole data may be more appropriate as expl
63. 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 MAGFOR3D This program performs forward modelling Command line usage magfor3d mesh obs loc model sus 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 Details are in the Elements chapter mesh the 3D mesh obs loc the inducing field parameters anomaly type and observation locations model sus the susceptibility model topo dat surface topography optional If omitted the surface will be treated as being flat Output file This file is created by MAGFOR3D it s name can NOT be specified Details are in the Elements chapter magfor3d mag the computed magnetic anomaly data Since the data in this file are accurate the column of the Standard deviations for the error is not included 3 3 MAGSEN3D This program performs the sensitivity and depth weighting function calculations Command line usage magsen3d magsen3d inp For a sample input file type magsen3d inp Format of the control file magsen3d inp mesh obs mag topo dat iwt beta ZnO wvlet itol eps Input files and parameters mesh 3D mesh obs mag data file Contains the inducingfield parameters anomaly type observation locations and the observed magnetic anomalies with estimated standard deviation topo d
64. ting northing and elevation of the observation measured in meters Elevation should be above Elevn the topography for surface data and below the topography for borehole data The observation locations can be listed in any order aincln jnclination and declination of the anomaly projection for observation n This is used only when idir adecln Q The brackets indicate that these two field are optional depending on the value of idir The total fieldanomaly is calculated when Aincl equals Incl and Adecl equals Decl The vertical field anomaly is calculated when Aincl 90 and Adecl 0 The user can specify other Aincl Adecl pairs to calculate the anomaly component in those directions Easting northing and elevation information should be in the same coordinate system as define in the mesh Example of obs loc file We provide two example file below The first file is for calculating total field anomaly at 441 points The inducing field has an inclination of 45 and a declination of 45 The second file is for calculating multi component anomalies in boreholes and each datum is specified by its own inclination and declination of anomaly projection File 1 single component data surface data 45 0 45 0 50000 0 incl decl geomag 45 0 45 0 aincl adecl direction of anomaly idir 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 2 multi co
65. tyfile from MAGSEN3D obs loc the inducing field parameters anomaly type and observation locations model sus the susceptibility model Output file magpre3d mag predicted data UBC Earth and Ocean Sciences F Jones Tuesday September 05 2006 14 39 00 MAG3D manual Synthetic examples Ver 4 0 Introduction We present two synthetic examples to illustrate Version 3 0 of MAG3D Important functionalities of MAG3D are the ability to handle multicomponent borehole data and the wavelet compression for large scale problems The two synthetic examples are constructed to show these features The first model consists of two vertical prisms and both surface and borehole data are simulated We illustrate the inversion of individual data sets and joint inversion of surface and borehole data The second example is intended to be large and is composed of several prisms buried at different locations and depths beneath a topographic surface Aeromagnetic data are simulated for this example The size of this example is too large for the direct approach to handle We show that the wavelet compression allows the solution of such a large problem with little demand on computing resources Example 1 This figure displays the synthetic susceptibility model It consists of B two prisms buried in a non susceptible background The collar 1 positions of three vertical boreholes are indicated on the plan section H An inducing field in the direct
66. uce the cost It is observed that when plotted on a log log scale the misfit curves for 3D magnetic inversion with and without positivity 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 ajo that gives rise to 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 thatSo can be used to approximate the slope of the misfit curve when the positivity is imposed A rigorous line search incorporating positivity Starts with an initial guess of jt 0 5 o9 This usually yields a misfit that is very close to the target value If the misfit is not sufficient close to however a new guess for jt is obtained which makes use of the approximate slope The inversion with updated can be solved efficient if the logarithmic barrier algorithm is started from an initial model close to the fina 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 target 4 after testing two to four values of ji This strategy is implemented in MAGINV3D as the first method for choosing the tradeoff parameter In p
67. wn target value of data misfit is used Two stages as discussed above are used and several solutions for different values exit Figure 3 Flow chart illustrating the solution of the 3D l magnetic inversion by MAGI NV3D using different strategies of must be tested to obtain one that produces the target for choosing the tradeoff parameter misfit When mode 2 the user specifie 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 of I nversion We now invert the total field anomaly data in Figure 2 using program MAGI NV3D The model region is discretized into 24 by 24 by 16 cells The cells are all cubes of 25 m on a side A total of 344 data 175 surface data and 144 three component borehole data are inverted to recover the susceptibility in 9 216 cells For a model objective function we choose t 0 0001 amp x y z 1 0 and a distance weighting with Ro 6 2 m A zero reference model is used Figure 4 shows one cross section and one plan section of the recovered model The two target prisms are well define in both horizontal and vertical locations and their amplitudes are comparable to those of the true model Figure 4 The susceptibility model oe _ recovered from the joint inversion of E E DEPTH m hd ry 300 408 300 200 1090 6

Download Pdf Manuals

Related Search

Related Contents

CX5000 Y CX6000  Notebook User manual  Scarica - Farmalock  MC/R730 - 121ware.com  2012 Infiniti  事業契約書(案) (PDF形式, 638.63KB)  Kingston Technology ValueRAM Memory 4096MB 333MHz DDR ECC CL2.5 DIMM  BlackBerry• Wireless Handheld  Mackie MR3SWTK User's Manual  

Copyright © All rights reserved. Failed to retrieve file