UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Some new approaches to gravity modelling Green, William Robert 1973

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1973_A6_7 G74_3.pdf [ 6.36MB ]
Metadata
JSON: 831-1.0052851.json
JSON-LD: 831-1.0052851-ld.json
RDF/XML (Pretty): 831-1.0052851-rdf.xml
RDF/JSON: 831-1.0052851-rdf.json
Turtle: 831-1.0052851-turtle.txt
N-Triples: 831-1.0052851-rdf-ntriples.txt
Original Record: 831-1.0052851-source.json
Full Text
831-1.0052851-fulltext.txt
Citation
831-1.0052851.ris

Full Text

C ' SOME NEW APPROACHES TO GRAVITY MODELLING by WILLIAM ROBERT GREEN B.A.Sc, Un i v e r s i t y of B r i t i s h Columbia, 1970 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE' i n the Department of GEOPHYSICS We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA May, 1973 In presenting t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r reference and study. I f u r t h e r agree that permission f o r extensive copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s r e p r e s e n t a t i v e s . I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my w r i t t e n permission. Department of Geophysics The U n i v e r s i t y of B r i t i s h Columbia Vancouver 8, Canada Date Mav 28. 1973 i ABSTRACT Geophysical i n v e r s i o n methods are most e f f e c t i v e when applied to l i n e a r f u n c t i o n a l s ; i t i s therefore advantageous to employ l i n e a r models fo r geophysical data. A two-dimensional l i n e a r model co n s i s t i n g of many h o r i z o n t a l prisms has been developed f o r i n t e r p r e t a t i o n of g r a v i t y p r o f i l e s . A Backus-Gilbert i n v e r s i o n which finds the acceptable model "nearest" to an i n i t i a l estimate can be r a p i d l y computed; i t e r a t i v e a p p l i c a t i o n of the technique allows a s i n g l e - d e n s i t y model to be developed at modest expense of computer time. Gravity data from the Guichon Creek batholi'th" were invert e d as a t e s t of the method, with r e s u l t s comparable to a standard polygon model. The entropy of these l i n e a r models i s a useful property which can be minimized to f i n d an optimum "structured" or "compact" model. Since a numerical optimization i s used, computations become p r o h i b i t i v e l y long f o r any large number of parameters. Several simple models have been found by minimizing the entropy of an i n t i a l model under the constraints imposed by a known gravity p r o f i l e . S i m i l a r r e s u l t s can be obtained by using a simpler obje c t i v e function. I t i s also p o s s i b l e to maximize model entropy; t h i s procedure tends to evenly d i s t r i b u t e the anomalous mass beneath the p r o f i l e , while minimum entropy t r i e s to concentrate mass i n as few prisms of the model as p o s s i b l e . i i TABLE OF CONTENTS Page ABSTRACT i LIST OF TABLES i v LIST OF FIGURES v ACKNOWLEDGEMENTS v i i CHAPTER 1 INTRODUCTION TO GRAVITY MODELLING - 1 CHAPTER 2 GRAVITY DATA AND INVERSION TECHNIQUES .• - " • 6 2.1 Common Methods of Gravity Modelling 6 2.2 Non-Uniqueness i n Gravity Inversion 8 2.3 New Approaches to Geophysical Inversion Problems 11 a. The Backus-Gilbert Method 11 b. The Generalized Inverse Approach 16 c. Monte Carlo Modelling 18 CHAPTER 3 LINEAR MODELS FOR GRAVITY DATA 21 CHAPTER 4 BACKUS-GILBERT APPROACHES TO GRAVITY MODELLING 35 4.1 Models "Cl o s e s t " to an I n i t i a l Estimate 36 4.2' A Method of "Weighted-Distance" Mini mi zatLon 50 CHAPTER 5 NUMERICAL OPTIMIZATION OF GRAVITY MODELS USING ENTROPY 70 5.1 Entropy of Gr a v i t y Models 70 5.2 Numerical Optimization of Model Entropy 75 CHAPTER 6 CONCLUSION 92 i i i REFERENCES 95 APPENDIX A MATHEMATICAL CONCEPTS 100 APPENDIX B UNITS IN GRAVITY MODELLING 102 APPENDIX C DERIVATIONS FOR BACKUS-GILBERT INVERSIONS 104 APPENDIX D NUMERICAL OPTIMIZATION TECHNIQUES 109 APPENDIX E THE COMPLEX METHOD 114 APPENDIX F COMPUTER PROGRAM 118 i v LIST OF TABLES TABLE I. Gravity e f f e c t s of prisms ( i n the standard model) and 30 equivalent c y l i n d e r s . TABLE I I . Gravity e f f e c t s of prisms ( i n the large-block models) 31 and equivalent c y l i n d e r s . TABLE I I I . Correction factors used i n computing Frechet k e r n e l s . 32 TABLE IV. Computation time f o r numerical optimizations. 90 V LIST OF FIGURES Figure 1. Linear models f o r gravity anomalies. 25 Figure 2. Geometric q u a n t i t i e s required f o r the gravity e f f e c t 27 of a prism ( a f t e r Parasnis, 1962, p. 43). Figure 3. An a r t i f i c i a l g r a v i t y p r o f i l e . 38 Figure 4. An i n v e r s i o n of the a r t i f i c i a l p r o f i l e using a zero 41 i n i t i a l model. Figure 5. A more r e s t r i c t e d model f o r Figure 3. 42 Figure 6. Approximate model derived from Figure 5. 43 Figure 7. An unsuccessful i n v e r s i o n using a large-block, model. 45 Figure 8. Simple bodies composed of l a r g e prisms. 46 Figure 9. Inversion of data from Figure 8Ca). 48 Figure 10. Inversion of data from Figure 8(b). 49 Figure 11. A r t i f i c i a l g r a v i t y p r o f i l e s used with the Weight 54 program. Figure 12. Inversion of the p r o f i l e i n Figure 11(a). 56 Figure 13. A f i n a l model f o r the p r o f i l e of Figure 11(a). 58 Figure 14. Inversion of the p r o f i l e i n Figure 11(b). 59 Figure 15. F i n a l model f o r the data of Figure 11(b). 62 Figure 16. Residual anomaly over the Guichon Creek b a t h o l i t h . 64 Figure 17. The f i r s t i n v e r s i o n of the b a t h o l i t h anomaly. 66 Figure 18. A model f o r the Guichon Creek b a t h o l i t h . 67 Figure 19. A b a t h o l i t h cross-section obtained by standard 68 modelling methods (from Ager et a l , 1972). Figure 20. Entropy of h o r i z o n t a l c y l i n d e r s of equal mass. 74 Figure 21. A simple entropy minimization. 80 Figure 22. Minimum entropy from a d i f f u s e s t a r t i n g model. 81 v i Figure 23. An. a r t i f i c i a l g r a v i t y p r o f i l e f o r entropy minimization. 82 Figure 24. A minimum entropy model from the a r t i f i c i a l p r o f i l e . 83 Figure 25. The functions x ln(x) and x(x-A). 85 Figure 26. An optimum model f o r Figure 23, using x(x-A). 86 Figure 27. A maximum entropy model corresponding to Figure 21. 88 Figure 28. A comparison of maximum and minimum entropy models fo r Figure 23. 89 v i i ACKNOWLEDGEMENTS I would l i k e to thank my supervisor, Dr. G.K.C. Clarke, f o r d i r e c t i n g my a t t e n t i o n to the problems discussed i n the t h e s i s ; f o r many h e l p f u l comments during the progress of the research; and f o r c r i t i c a l evaluation of the manuscript. I am g r a t e f u l to Chuck Ager f o r providing some r e a l g r a v i t y data, and f o r some i n s i g h t i n t o the p r a c t i c a l aspects of g r a v i t y work. A p l o t t i n g program provided by Gary J a r v i s was i n s t r u -mental i n preparing many of the f i g u r e s . I am indebted to my wife f o r general encouragement and f o r an e x c e l l e n t job of typing. F i n a n c i a l support during this work was provided by a N a t i o n a l Research Council Postgraduate Scholarship, a Chevron Standard Scholarship, and from Dr. Clarke's research funds. 1 CHAPTER 1 Introduction To Gravity Modelling There are three major steps necessary to obtain u s e f u l information from geophysical explorations; data a c q u i s i t i o n , a n a l y s i s , and i n t e r p r e t a t i o n . Data must f i r s t be recorded i n the f i e l d with s u f f i c i e n t accuracy to allow drawing g e o l o g i c a l conclusions. Secondly, the data must be analyzed and put i n t o a form which enable them to be r e l a t e d to properties of the explo-r a t i o n o b j e c t i v e . F i n a l l y , by applying t h e o r e t i c a l knowledge of the r e l a -tionships between p h y s i c a l properties of the earth and the type of data under consideration, one can estimate those p h y s i c a l properties i n the explored region. This thesis w i l l consider some new variants of t h i s t h i r d stage i n the treatment of g r a v i t y data. Before any geophysical data can be used to estimate the p h y s i c a l parameters of the earth, i t i s necessary to solve the " d i r e c t problem" f o r the p a r t i c u l a r type of data. The d i r e c t Cor forward) problem consists of p r e d i c t i n g values of a geophysical function from known parameters of the earth. For example, i f free o s c i l l a t i o n data are to y i e l d information about the earth's density, the r e l a t i o n s h i p between the periods of the o s c i l -l a t i o n s and a known density function must be known; i . e . one must be able to accurately p r e d i c t the o s c i l l a t i o n periods f o r a hypo t h e t i c a l earth of known density. Fortunately the d i r e c t problem of f i n d i n g the g r a v i t a t i o n a l a t t r a c t i o n of a given density d i s t r i b u t i o n (e.g. a body of known shape, l o c a t i o n , and density) can be solved using p o t e n t i a l theory (MacMillan 1930). I f the forward problem can be solved, i t may be pos s i b l e to f i n d 2 solutions to the "inverse problem", which attempts to find physical para-meters of the earth from observed values of certain geophysical functions. Geophysical inversion frequently is a problem in modelling: one seeks parameters for a model whose properties correspond to observations of the real earth, by using the solution to the direct problem for the model para-meters and type of data available. The parameters may be determined in two ways: by an analytic method based on the functional relationship of the parameters to the data; or by a trial-and-error process in which an inter-preter adjusts model parameters in an attempt to improve the f i t to the data (in this case, he is continually solving the direct problem, rather than finding a formal solution to the inverse problem). A simple example of inversion is refraction seismology, where the objectives are the velocities and thicknesses of crustal layers which produce the observed travel times at the surface. In many cases, inversion is a subjective process without exact solution. To begin with, the infinite detail of the real earth cannot be uniquely determined from the finite number of observations available. It may also be true that several models can satisfy the same set of data, i f those data are not in themselves sufficient to determine a unique model. This is true of gravity data, as i t is well known that many quite different density distributions may produce identical gravity anomalies (Skeels, 1947). The non-uniqueness of inverse problems has been treated in different ways. Many techniques impose enough constraints on the model being sought to ensure that i t is unique. A similar, but conceptually simpler, viewpoint is to seek particular members of the set of models which satisfy the data, usually 3 by optimizing some property of the model, f o r example, f i n d i n g the "nearest" acceptable model to a s p e c i f i e d i n i t i a l model which does not s a t i s f y the data. Several d i f f e r e n t types of models may be used with the same data, and the i n v e r s i o n procedure w i l l n a t u r a l l y depend on the type of parameters being used. For example, the parameters of a gravity model might be the external shape and constant density of an anomalous body, or perhaps the density at various points i n a given subsurface region. In some problems, the speed and e f f i c i e n c y of the i n v e r s i o n may be improved by d e f i n i n g a new type of model whose parameters are more simply r e l a t e d to the observed data. Linear models are often preferable, since they obviate the need f o r i t e r a t i v e procedures required i n f i t t i n g non-linear models to the data. By a l i n e a r model, we mean a model whose parameters are l i n e a r l y r e l a t e d to the observed geophysical data. The inverse problem i n g r a v i t y e x p l o r a t i o n has usually been ap-proached by considering r e s t r i c t e d classes of models with only a few para-meters, and adjusting an i n i t i a l model to f i n d a " b e s t - f i t " s o l u t i o n ; t h i s has generally meant non-linear models and hence i t e r a t i v e methods. With the viewpoint of f i n d i n g p a r t i c u l a r models i n a large s e t , the l i n e a r i t y of the density-gravity r e l a t i o n s h i p can be e x p l o i t e d to develop new modelling techniques. Before considering methods, a b r i e f d e s c r i p t i o n of the objec-t i v e s and requirements of the i n v e r s i o n process i s i n order. Several types of anomalies may be found i n the earth's gravity 4 f i e l d , by applying d i f f e r e n t corrections to the t o t a l f i e l d (e.g. f r e e - a i r anomaly, Bouguer anomaly, i s o s t a t i c anomaly: see Garland, 1965, Chapter 4); d i f f e r e n t i n t e r p r e t a t i o n methods can be developed for d i f f e r e n t types of anomalies. Modelling techniques f o r exploration data usually attempt to define a s i n g l e subsurface formation to account f o r an i s o l a t e d anomaly (Grant and West, 1965, p. 268, term t h i s " q u a n t i t a t i v e i n t e r p r e t a t i o n " ) . The anomaly i s defined by applying a l l necessary corrections to the raw data (e.g. l a t i t u d e , f r e e - a i r , Bouguer, ter r a i n ) to f i n d the Bouguer anomaly (Grant and West, 1965, pp. 236-243), and then removing regional v a r i a t i o n s from the Bouguer gr a v i t y to locate i s o l a t e d , l o c a l features. The anomaly should then be due s o l e l y to some i s o l a t e d body whose density i s d i f f e r e n t from the c r u s t a l average. The separation of l o c a l anomalies from the r e g i o n a l f i e l d requires great care, since the shape of r e s i d u a l anomalies can be a l t e r e d by the separation process (Skeels, 1967; Ulrych, 1968). Any a n a l y t i c i n v e r s i o n procedure can only be s u c c e s s f u l i f the anomalous g r a v i t y corresponds to a r e a l feature of the earth, and i s not even p a r t l y a r e s u l t of a f i l t e r i n g process. I f the regional separation i s s u c c e s s f u l , a model can be formu-l a t e d s t r i c t l y i n terms of the density contrast between an anomalous feature and the average c r u s t a l rocks. In t h i s case, "density" can be taken to imply "density contrast." The corrections applied i n f i n d i n g r e s i d u a l anomalies contain many un c e r t a i n t i e s , and thus g r a v i t y models need not produce an exact f i t to the observations (and usually cannot do so). For the purpose of studying 5 modelling techniques, we w i l l assume that the anomalous g r a v i t y i s exactly known, but frequently w i l l accept a close f i t rather than attempting to produce p e r f e c t agreement with the data. The remainder of the thesis i s concerned with f i n d i n g simple models from r e s i d u a l anomalies. 6 CHAPTER 2 Gravity Data And Inversion Techniques 2.1 Common Methods of Gravity Modelling There are sev e r a l ways to b u i l d models f o r the c r u s t a l structures which produce g r a v i t y anomalies. A model could be formed by s p e c i f y i n g : the shape and p o s i t i o n of a body of constant density; the l o c a t i o n of a number of simple bodies; the density (contrast) at various points i n the subsurface; or other parameters. The basis of the i n v e r s i o n process i s to e s t a b l i s h the r e l a t i o n s h i p between the model parameters and the surface g r a v i t y ( i . e . to solve the forward problem); and then apply the r e l a t i o n s h i p to compute a set of model parameters corresponding to the observed g r a v i t y . The computational method w i l l depend on the nature of the parameters, and i n some cases i n merely t r i a l and e r r o r adjustment to improve an i n i t i a l model. Perhaps the most common methods of determining simple models from g r a v i t y anomalies stem from the polygon methods of Taiwan! et a l (1959) and Talwani and Ewing (1960), developed to solve the d i r e c t problem f o r two and three-dimensional features, r e s p e c t i v e l y . The two-dimensional method assumes the body to be an N-sided polygon of known, constant density; and uses the p o s i t i o n s of the v e r t i c e s to compute the g r a v i t y e f f e c t of the body. By s p e c i f y i n g the l o c a t i o n s of each vertex i n an i n i t i a l polygon, and comparing the computed g r a v i t y to an observed surface g r a v i t y p r o f i l e , one can make i t e r a t i v e changes to the coordinates of each vertex u n t i l there i s accept-able agreement with the r e a l data. A three-dimensional body can be modelled 7 with a s e r i e s of h o r i z o n t a l lamina, each of polygonal o u t l i n e . Again the gr a v i t y e f f e c t can be computed from the coordinates of each vertex and the given density, and i t e r a t i v e changes can be made to an i n i t i a l model to produce good agreement with the measured gr a v i t y on the two-dimensional surface. Rather than using a t r i a l and e r r o r basis f o r successive changes to the model, Corbato (1965) suggested a least-square e r r o r approach to improve a two-dimensional polygon, and thus accelerated convergence to the observed anomaly. Another frequently used approach i s to consider models consis-t i n g of many bodies of simple shape, f o r example, spheres or rectangular prisms. Tanner (.1967) developed an i t e r a t i v e procedure to develop two-dimensional models c o n s i s t i n g of constant density rectangular blocks. He assumed that the density, width, and depth to top of each block are known, so that the unknown parameters to be obtained are the depths to the bottom of each block. Dyrelius and Vogel (19 72) also used t h i s approach. Negi and Garde (1969) and Agarwal (1971) developed s i m i l a r models, but allowed each v e r t i c a l column to be subdivided i n t o units of d i f f e r e n t density. T h e i r methods were aimed at f i n d i n g g r a v i t y models of sedimentary basins, where several layers of d i f f e r e n t density might be expected. Nagy (1966) used models co n s i s t i n g of three-dimensional prisms of various s i z e s and d e n s i t i e s , and then used t r i a l and e r r o r adjustment to improve the f i t of an i n i t i a l model to a g r a v i t y map. C o r d e l l and Henderson (1968) found three-dimensional models made up of v e r t i c a l prisms, using an automatic i t e r a t i v e method to solve f o r the depths of the prisms beneath a gridded two-dimensional surface. 8 Modelling i s by no means the only way to u s e f u l l y i n t e r p r e t gra-v i t y data. Many other approaches have been used, some of which are des-cribed by Garland (1965) and Grant and West (1965). One example i s to examine the c h a r a c t e r i s t i c s of p o t e n t i a l f i e l d s i n the s p a t i a l frequency domain (e.g. Odegard and Berg, 1965; Berezhnaya and Telepin, 1968). Another widely used method i s downward continuation, which attempts to f i n d the "topography" of a density i n t e r f a c e (Grant and West, 1965). However, as th i s thesis i s intended to give i n s i g h t i n t o modelling techniques, such other methods w i l l not be discussed f u r t h e r . 2.2 Non-Uniqueness i n Gravity Inversion The greatest problem to be overcome i n f i n d i n g the source of a gr a v i t y anomaly i s the non-uniqueness associated with p o t e n t i a l f i e l d i n v e r -s i o n s ; i . e . d i f f e r e n t bodies or c r u s t a l structures can produce the same gr a v i t y e f f e c t at the surface, and there i s no way to d i s t i n g u i s h which of the acceptable models corresponds to the r e a l earth, from the gra v i t y data alone. Skeels (1947) presented s e v e r a l examples of quite d i f f e r e n t forma-tions which produce the same surface g r a v i t y . He also developed an a n a l y t i c proof (using Green's theorem to integrate the v e r t i c a l g r a v i t y e f f e c t around the body) to show that a la y e r of v a r i a b l e density can have the same gra v i t y p r o f i l e as a point mass at greater depth. Skeels noted that the ambiguity i n i n t e r p r e t a t i o n can be considerably reduced i f other g e o l o g i c a l or geo-p h y s i c a l data are a v a i l a b l e . In a fu r t h e r analysis of non-uniqueness, Roy (1962) studied s i t u a t i o n s i n which a unique model could be developed; f o r example, i f the 9 anomalous density i s confined to a plane at known depth, there i s only one density function which can produce a given g r a v i t y anomaly. He noted that this assumption i s i m p l i c i t i n s t r u c t u r a l determination by downward c o n t i -nuation. Another unique class of models are those of constant density and known external shape; models of this type usually t r e a t s i z e , p o s i t i o n , and o r i e n t a t i o n of the anomalous body as unknown parameters. Al-Chalabi (1971) examined methods f o r f i n d i n g models which pro-vide unique s o l u t i o n s to the inverse problem. He discussed other ( i . e . non-potential) sources of non-uniqueness, i n c l u d i n g incomplete knowledge of the anomaly; observational e r r o r s ; and the use of simple models to describe the complicated r e a l earth. He studied the r e s u l t s of modelling with polygons, using a r t i f i c i a l g r a v i t y p r o f i l e s from polygonal bodies. Unique solutions were p o s s i b l e i f the model polygon had as many sides as the ' r e a l ' body; however, the uniqueness could be destroyed by inadequate p r o f i l e length or other f a c t o r s . I f the model polygon had fewer sides than the r e a l body, any s o l u t i o n would not be unique. Al-Chalabi's con-clusions were: a s a t i s f a c t o r y s o l u t i o n can be obtained by s p e c i f y i n g only one model parameter, provided none of the factors contributing to ambiguity are too great; and, i n cases where there are strong sources of ambiguity, i t i s desirable to produce a number of s o l u t i o n s by examining i n t e r v a l s i n the parameter 'hyperspace' corresponding to acceptable agreement with the observed anomaly. Most techniques f o r i n v e r t i n g g r a v i t y data avoid the problem of non-uniqueness by using r e s t r i c t e d classes of models; s o l v i n g f o r a very 10 l i m i t e d number of unknown parameters may ensure that there i s only one model of that type which acceptably f i t s the data. The usual method i s to a l t e r chosen parameters of an i n i t i a l model to obtain agreement with the observed g r a v i t y , and thus the f i n a l model may l a r g e l y depend on the nature of the i n i t i a l model. In such methods, the "uniqueness" of the s o l u t i o n may l i e i n being the acceptable model which most c l o s e l y resembles the i n i t i a l estimate. I t i s p o s s i b l e that an inverse method thought to be unique i s not; f o r example, Parker (1973) found that the standard models f o r magnetization of the oceanic crust (which are very s i m i l a r to g r a v i t y modelling using a constant thickness l a y e r ) , do not provide a unique inverse f o r the magnetic anomalies observed at the ocean surface, although t h i s i s commonly b e l i e v e d to be the case. One d i f f i c u l t y with many methods i s that the r e s t r i c t i o n s on the model are incorporated i n t o the numerical techniques used, and the i n t e r -p reter may not be f u l l y aware of them. For example, the polygon methods produce a model of a given number of s i d e s , but do not guarantee that a d i f f e r e n t polygon w i l l not also be acceptable. The Cordell-Henderson method requires a reference plane to mark the top, midpoint, or bottom of the prisms which comprise the model. Such r e s t r i c t i o n s often necessitate an i n i t i a l model which f i t s the data reasonably w e l l . For example, C o r d e l l and Henderson found t h e i r r e s u l t s to be p h y s i c a l l y p l a u s i b l e only f o r res-t r i c t e d choices of reference plane; to obtain a reasonable s p h e r i c a l model, the reference plane had to be set through the center of the sphere. The non-uniqueness of g r a v i t y inversions i s then not an i n s u r -11 mountable obstacle. I t i s usually p o s s i b l e to produce a simple model (with a l i m i t e d number of parameters) to account f o r any g r a v i t y anomaly; however i t w i l l not n e c e s s a r i l y be the only acceptable simple model. Before making a g e o l o g i c a l i n t e r p r e t a t i o n , i t i s e s s e n t i a l to examine the i n v e r s i o n method to understand the l i m i t a t i o n s of the s o l u t i o n . The a v a i l a b i l i t y of other geophysical or g e o l o g i c a l data can be of great help i n e s t a b l i s h i n g an i n i t i a l model which allows a reasonable s o l u t i o n to the inverse problem. In s i t u a t i o n s where there i s l i t t l e p r i o r information, the choice of i n i t i a l model may become d i f f i c u l t ; i n t h i s case, a method r e q u i r i n g few assumptions about the form of the anomalous body would be advantageous. Some new techniques f o r g r a v i t y data w i l l be deve-loped l a t e r i n hopes of f i n d i n g inversions without r e q u i r i n g a d e t a i l e d i n i t i a l model. 2.3 New Approaches to Geophysical Inverse Problems a. The Backus-Gilbert Method In recent years, a strong t h e o r e t i c a l basis f o r geophysical inverse problems has been e s t a b l i s h e d , p a r t i c u l a r l y i n the work of Backus and G i l b e r t (1967, 1968, 1970). A major innovation i s that the non-uniqueness of inversions i s e x p l o i t e d by seeking p a r t i c u l a r models which * s a t i s f y the data from the H i l b e r t space of a l l p o s s i b l e models. This can usually be achieved by optimizing some property of the model subject to * the constraints imposed by the observations; f o r example, the "distance" of the model from an i n i t i a l model ( i n a parameter space) might be minimized, * See Appendix A. 12 or the average properties of adjacent models might be obtained by finding the "smoothest" model. Some of the older inversion methods may produce similar results, but the Backus-Gilbert approach emphasizes the true nature of the solution; in addition, their formalism is much more flexible in that similar algorithms can be employed to produce models of different properties. Parker (1970) has given a good summary of their techniques. A central concept in the Backus-Gilbert method is that models which satisfy the data are points in a Hilbert space which includes a l l possible models. In one sense, such models can not be unique for any geo-physical data, since only a f i n i t e number of data are available, but the earth's properties can be i n f i n i t e l y variable i f viewed on a sufficiently fine scale. Particular data, for example potential f i e l d observations, may also contain other sources of non-uniqueness. A model with N parameters Cor one parameter evaluated at N spatial positions) may be considered to l i e in a N-dimensional subspace. If the geophysical functions under consideration can be expressed * as inner products defined on the Hilbert space, the Backus-Gilbert formalism can exploit the properties of the inner product to develop a model from the observations. A linear functional can always be written as an inner product CHoffman and Kunze, 1961, p. 235); non-linearity can also be handled, so long as the functional is Frechet-differentiable and can be linearized in the region near an i n i t i a l model. In geophysical inversions, this r e s t r i c -tion i s usually no problem, since most earth data have been shown to be * See Appendix A. 13 F r e c h e t - d i f f e r e n t i a b l e ; Backus and G i l b e r t (1967) c i t e s e v e r a l examples, i n c l u d i n g mass, moment of i n e r t i a , t r a v e l times of seismic waves, and normal mode o s c i l l a t i o n frequencies. The b a s i c notation of the Backus-Gilbert method i s as follows: a) E." o r y, are the observed data j J b) E.(M) i s the data f u n c t i o n a l of the model M J Frechet d i f f e r e n t i a b i l i t y implies c) Ej (M + AM) = Ej (M) + (F , AM) + £ (AM). d) Fj i s the Frechet k e r n e l , and (F.., AM) i s an inner product e) A = E a i ^ i i s 3 1 1 averaging k e r n e l f o r the model M The i n v e r s i o n s t a r t s with computation of the Frechet kernels f o r the p a r t i c u l a r data functionals and model parameters. D i f f e r e n t approaches are then a v a i l a b l e , depending on the nature of the desired s o l u t i o n . For example, the "distance" of an acceptable model from an i n i t i a l guess can be minimized, subject to the constraints that the data functionals have c e r t a i n known values. This reduces to a c l a s s i c a l calculus of v a r i a t i o n s problem, e a s i l y solved v i a the Lagrangian m u l t i p l i e r technique ( d e t a i l s are given i n Appendix C). I f the functionals are non-linear, the s o l u t i o n must be i t e r a t e d i n small l i n e a r steps from the previous model. . Another technique i s to f i n d average properties of acceptable models, using an averaging k e r n e l which i s a l i n e a r combination of the Frechet kernels. The averaging k e r n e l i s given by 14 A = Z a.F. i 1 1 The 1968 paper demonstrates that the average model <m> can then be expressed as a l i n e a r combination of the observations, i . e . so long as the functi o n a l s are l i n e a r , or can be l i n e a r i z e d i n the neigh-bourhood of the average model. E s t a b l i s h i n g a c r i t e r i o n f o r the q u a l i t y of the l o c a l averages leads to computational methods f o r determining the c o e f f i c i e n t s a^ of the averaging kernel. Backus and G i l b e r t (1970) minimized the "spread" of A, given by to determine averaging kernels f o r whole earth models. The b a s i c aim of the process i s to f i n d an average with the shortest r e s o l u t i o n length pos-s i b l e . The form of the averaging kernels gives an estimate of the r e s o l v i n g power of the i n v e r s i o n , i d e a l l y , A(r) should be a d e l t a function (hence the term "<5-ness c r i t e r i a " i n the 1968 paper) . i f the model i s l i n e a r ; i n add i t i o n , a l l models which are "near" to each other ( i n the sense that (m) - Ej(m*) i s l i n e a r i n (m-m'), where m, m' denote d i f f e r e n t models) have the same average p r o p e r t i e s . The averaging process thus obtains unique solutions from a f i n i t e s e r i e s of observations; <m> = Am dV = E a y One advantage of th i s approach i s that the averages are unique 15 however i t does not remove other sources of non-uniqueness. When the data contain errors, the model averages can no longer be precisely determined. The problem that then arises i s that the standard error in <m> increases as the resolving length decreases, i.e. there is a tradeoff between the resolution of the averaging kernel and the standard error in the resulting averages. In their 1970 paper, Backus and Gilbert were concerned with computing tradeoff curves for noisy data, and with strategies which lead to a reasonable compromise. In some cases, the resolving power of the data may be so poor that different models cannot effectively be distinguished. This i s certainly true for gravity data. Since the gravity-density relationship i s linear (.see Chapter 3) , a l l possible density models are "near" to one another i n the Backus-Gilbert sense. It has been noted earlier (Section 2) that quite different models can produce the same gravity effect; the only properties common to a l l models are the total (anomalous) mass, and the surface position of the center of mass (which can easily be found from the data alone, e.g. Grant and West, 1965,'pp. 227-230). Backus and Gilbert applied their methods to the problem of deter-mining the density structure of the earth, and have been able to estimate the resolution limits imposed by the f i n i t e data set available. Parker (19 70, 1972) successfully adopted their technique to model the conductivity structure of the mantle, and to make gross estimates of the core densities of the outer planets. Der and Landisman (1972) used the same basic approach to produce crustal models from surface wave observations, and examined the a b i l i t y of the data to resolve density and shear velocity. 16 b. The Generalized Inverse Approach Other authors i n recent years have attempted to solve inverse problems using the generalized inverse of a matrix (Penrose, 1955). I f the a v a i l a b l e data can be r e l a t e d to model parameters by a matrix equation, the generalized inverse w i l l give a p a r t i c u l a r model which s a t i s f i e s the data, the s o l u t i o n being the one of l e a s t Euclidean "length" (Smith and F r a n k l i n , 1969) . Here "length" implies distance from the o r i g i n i n parameter space (see Appendix A). The usual a p p l i c a t i o n i s to solve f o r corrections to i n i t i a l model parameters, which i s equivalent to the Backus-Gilbert method of minimizing "distance" between models. The b a s i c formulation of a generalized inverse problem i s as follows (Wiggins (1972) and Jackson (1972) give more complete d e t a i l s ) . One seeks N model parameters P. knowing M observations 0. of the r e a l earth. Changes i n model parameters are r e l a t e d to changes i n the data f u n c t i o n a l by a matrix A' of " v a r i a t i o n a l parameters". I A ' J [AP'J * {AC'J (1) where l A ^ ' J = I3C i/9P j] C^ = F^(Pj) i s the l i n e a r i z e d f u n c t i o n a l corresponding to observation i A C . ' = 0 i - C i Solving t h i s system using the generalized inverse of A y i e l d s a set of 2 2 parameter corrections A P ^ , such that both | A P| and e = |AAP - Ac| are minimized (Wiggins, 1972). I f model parameters have d i f f e r e n t dimensions, 17 scaling i s necessary to make the minimization reasonable. Error s t a t i s t i c s may also be incorporated in weighting. One then solves the modified system A AP = AC (2) , . _-l/2 t Tl/2 where A = S A' W AP = W"1/2 AP' AC = S - 1/ 2 AC' S = covariance matrix of the observations W = covariance matrix of the parameters The inversion i s performed by finding the eigenvalues and eigenvectors of A. Three further benefits of this analysis are: near-zero eigenvalues can be rejected to remove potential i n s t a b i l i t y of the solution caused by noisy data; the eigenvectors corresponding to columns of A are a measure of the resolution of the parameters; and the row eigenvectors can indicate the information distribution among the observations. The generalized inverse method i s becoming popular, since i t can easily be implemented using standard matrix techniques. Jordan and Franklin (19 71) used a variation of the basic method, and found earth density models by considering them to be outputs of a linear f i l t e r ; this formulation allows rejection of models which are not "smooth". Jackson (19 72) examined the theoretical performance of a generalized inverse for underdetermined and overconstrained systems. Wiggins (1972) also studied resolution in deducing density models from surface wave and free oscillation data. Like the Backus-18 G i l b e r t method, the generalized inverse requires l i n e a r or l i n e a r i z e d f u n c t i o n a l s ; non-linear models can be found by i t e r a t i v e s o l u t i o n of a l i n e a r i z e d system. B r a i l e et a l (1973) used the method to solve f o r the d e n s i t i e s of multi-prism gravity models. The number, s i z e , and p o s i t i o n of prisms was l e f t to the d i s c r e t i o n of the i n t e r p r e t e r , allowing considerable use of other data. The generalized inverse found corrections to d e n s i t i e s of an i n i t i a l model, with reasonable success i n such a p p l i c a t i o n as a c r u s t a l model f o r a p r o f i l e across Texas. However, t h e i r success depended on being able to construct a good i n i t i a l model comprised of prisms of various s i z e s , and t h e i r method might not be p r a c t i c a l i n areas where l i t t l e other information i s a v a i l a b l e . c. Monte Carlo Modelling Monte Carlo techniques f o r geophysical modelling have become p r a c t i c a l with the advent of large, f a s t d i g i t a l computers. In essence, these methods construct models whose parameters are randomly d i s t r i b u t e d over s p e c i f i e d i n t e r v a l s , and then t e s t these random models to f i n d those whose properties agree with observations of the r e a l earth. A l l the accep-table models are examined to estimate the average values of model parameters, t h e i r standard deviations, and the r e s o l v i n g power of the data (Wiggins, 1972). Monte Carlo methods i n geophysics are usually attempts to randomly sample the space of a l l models which s a t i s f y the data, i n hopes of f i n d i n g the bounds of a c c e p t a b i l i t y . The main d i f f i c u l t y i s that almost a l l of the 19 random models are rejected; hence the model b u i l d i n g process must be res-t r i c t e d to ensure f i n d i n g some acceptable models, and the t e s t i n g s i m p l i f i e d to reduce t o t a l computing time. The former i s accomplished by generating parameters randomly wi t h i n the r e s t r i c t i o n s of present knowledge; f o r example density values are picked i n a small i n t e r v a l about the accepted earth models, and usually made to increase monotonically below the upper mantle. The t e s t i n g of models against the observations usually employs l i n e a r i z e d v a r i a -t i o n a l parameters rather than exact computations of non-linear func t i o n a l s (the v a r i a t i o n a l parameters of Wiggins (1968) are frequently used to compute free o s c i l l a t i o n periods of earth models). Monte Carlo methods are often used i n problems where many d i f f e r e n t data are to be inverted. Press (1968, 1970) considered models whose para-meters were density, shear v e l o c i t y , and compressional v e l o c i t y at 88 r a d i i w i t h i n the earth; the values were generated randomly at 23 points, and i n t e r p o l a t e d elsewhere. Models were tested against 9 7 eigenperiods, various t r a v e l time data, and the earth's mass and moment of i n e r t i a . Press found 11 acceptable models from a t o t a l of 5 m i l l i o n , r e q u i r i n g 20 hours of com-puter time (1968). Refinements to the procedure enabled him to f i n d 11 successes i n one hour of computation (1970). Other applications have been discussed by Ke i l i s - B o r o k and Yanovskaya (1967) and Wiggins (1969) i n inversions of body-wave data, and Anderssen (1970), who sought bounds on the conductivity of the lower mantle. A f i n a l point of discussion i s the question of whether the method provides adequate sampling of the model space. Backus and G i l b e r t (1970» 20 p. 126) suggest that i t cannot; however Press (1970) demonstrated that a small number of random models (25) can e f f e c t i v e l y span the parameter space i n which they were constructed. Anderssen et a l (19 72) discussed s e v e r a l points of contention, and concluded that the various methods of i n v e r s i o n should be able to provide equivalent information about the earth. 21 CHAPTER 3 Linear Models For Gravity Data As noted earlier (Section 2.1), the commonly-used methods for modelling gravity data involve non-linear functionals; i.e. the parameters of the model are not linearly related to the surface gravity. For example, the polygon methods use the subsurface location of vertices as parameters; and the unknowns in Tanner's models are depths to bottom of blocks. Linear density models can be constructed, since surface gravity is a linear function of density; the inversion techniques to be developed w i l l apply to linear models (densities at different locations are the only parameters). A linear approach was inspired by the simplified methods developed for linear (or linearized) functionals, which were discussed in Chapter 2. Since gravity data may be measured along a line or over a surface grid, models of the subsurface density distribution may be two-dimensional or three-dimensional, respectively. In either case, analytic relationships can be established to compute the gravity effect of the model; however the two-dimensional case requires the assumption that the model has i n f i n i t e extent in the third spatial coordinate. The two-dimensional equations are thus always an approximation, but are generally acceptable when the length of an anomalous body is greater than about five times i t s width (Grant and West, 1965). Only two-dimensional models w i l l be considered here, for they involve fewer parameters and are therefore more practical for testing new techniques. The approach would be essentially the same for three-dimensional models, but computational tests would be much more expensive. The primary 22 purpose i s to t e s t new modelling methods, so i t w i l l be assumed that the g r a v i t y data are exactly known; i n any case, non-uniqueness precludes d e t a i l e d consideration of data e r r o r s . The g r a v i t a t i o n a l p o t e n t i a l of any mass d i s t r i b u t i o n i n a Cartesian system i s * 00 * 00 » 00 U(x,y,z) - T P ( x 0 » y o » Z o > dx dy dz m J —co J —oo J _oo -r 2 2 2 2 where r = [(x-x ) + (y-y ) + (z-z ) J i \ Q ' o o and y = the g r a v i t a t i o n a l constant I f a body i s e s s e n t i a l l y two-dimensional, one assumes that p does not vary i n the y Q d i r e c t i o n , and performs one i n t e g r a t i o n to obtain . OD . OO U(x,z) = 2y P ( X 0 ' Z O ) l n ( R ) d x o d z o ( 2 ) J —00 i —o 2 2 2 where R = ( X - X Q ) + ( Z - Z Q ) (Grant and West, 1965,p.230) Gravity surveys measure the v e r t i c a l component of g r a v i t a t i o n a l a t t r a c t i o n at the earth's surface (z = 0); the v e r t i c a l g r a v i t y i s r e l a t e d to the p o t e n t i a l by g (x,0) = 2 -3U(x,z) I 3z 'z = 0 23 hence _oo J -oo (x-X ) + Z O O As noted earlier, i s a linear function of p. Linearity requires F(ax + y) = aF(x) + F(y) (Hoffman and Kunze, 1961, p. 91) From (3) OO » CO LI Z o ( a p l + P 2 } ^ o ^ o g z(ap 1 + p 2) = 2 Y , , 2 (x-x ) + z o o or 0 0 /• 0 0 z p, dx dz r °° r 0 0 z p„dx dz " 1 ' ' o 2 o o j.oo ^ 0° z n ax az r°°r[ | 0 1 2 ° 2+2W f * —OO * —CO f V — " V ^ 4" »7 ^ —00 * -8 z C a p l + P 2 } = 2 Y a J J , ,ZM 2 • - J j , , 2 ^ 2 —00 —00 (x-x ) + z — 0 0 —°° (x-x ) + z o o o o = ag (p ) + g (P 9) Q.E.D. z 1 Z Z Knowing that (3) i s linear in p, the Frechet kernel for g is immediately z seen to be 2yz G(x,x , z J = o' o' , . 2 , 2 (x-x ) + z o' o and the ver t i c a l gravity may be written as an inner product 0 CO # 00 g z(x) = (G(x,x o,z o), P(x Q,z o)) = J j G(x,x z ) p(x o,z )dx d z o (5) 24 G(x,x ,z ) may also be considered to be the Green's function which generates o o the v e r t i c a l g r a v i t y from the source term p ( x Q , z o ) . Considering g as a d e r i v a t i v e of the p o t e n t i a l , the t o t a l ano-z malous mass which causes g can be determined by an a p p l i c a t i o n of Gauss' z theorem; f o r a two-dimensional body the r e s u l t i s f M = 7—- g (x)dx (Grant and West, 1965,p.273) (6) 4Y J_oo Z A l l models which f i t a given g r a v i t y p r o f i l e must have the same excess mass, since the i n t e g r a l i s not taken over the source region, but i s r e l a t e d s o l e l y to the data. Equation (3) cannot be applied d i r e c t l y , since an exact i n t e g r a t i o n would require a knowledge of p at a l l values ( X Q , Z O ) ; the equations (3), (4), and (5) must be adapted to consider models which spec i f y p at only s e l e c t e d values of x ,z . The models used her e a f t e r w i l l have the form shown i n o" o Figure 1. The subsurface region underlying the g r a v i t y p r o f i l e i s divided i n t o rectangular c e l l s , centered at ( X O , Z Q ) . The model parameters are con-stant d e n s i t i e s p ( x Q , z o ) assigned to each c e l l ; the surface gravity i s then the sum of the g r a v i t y e f f e c t s of each c e l l (again because of l i n e a r i t y ) . The i n t e g r a l (4) i s now w r i t t e n as a summation ; ( x ± ) = Z Z G(x i,x^,z k) p(x^.,xfc) j k (7) provided that the Frechet kernel G ( X , X Q , Z O ) i s modified from the point mass Gravity stations: = i*dx dx ' , l , L Prism centers: x, = i*dx j z, = k*dz k Ca) Standard Model Gravity stations: xx^ = i*dx dx d Prism centers: x^ = C2j - 1/2)*dx = 2(k - 1/2)*dz (b) Large-block Model Fig. 1. Linear models for gravity anomalies. 26 expression (4) to a form corresponding to a c e l l of dimensions dx, dz. M o t t l and Mottlova (1972) used such a regular g r i d i n t h e i r simple two-dimensional models (using 20 to 30 c e l l s to model a 5 point g r a v i t y p r o f i l e ) , but r e s t r i c t e d density to values of 0 or 1 only, and thus required an i t e r a t i v e method of i n v e r s i o n (numerical optimization of a "shape pre-ference" f u n c t i o n ) . The models of B r a i l e et a l (1973) also used prism den-s i t i e s as parameters to e x p l o i t l i n e a r i t y , but reduced the t o t a l number of parameters by using blocks of d i f f e r e n t s i z e s , u s u a l l y much greater than the surface s t a t i o n spacing. The models suggested here d i f f e r i n that a l l blocks w i l l be of the same dimensions, generally equal to the s t a t i o n spacing; and the subsurface geometry w i l l not be a l t e r e d f o r each set of data to be inverted. Modelling the earth with these d i s c r e t e elements i s a sampling process, and thus the s p a t i a l frequencies which can be represented have an upper l i m i t . However, since the g r a v i t y data are also sampled, high frequen-cies of the r e a l earth density d i s t r i b u t i o n w i l l be a l i a s e d , and the model should be able to represent a l l r e a l frequencies present i n the data. I f prisms of dimension twice the surface spacing are used, the model cannot represent a l l p o s s i b l e frequencies i n the data; such models may be adequate i f large features are expected (which i s usually true i n g r a v i t y e x p l o r a t i o n ) . In e i t h e r case, one must hope that the s t a t i o n spacing used i n measuring the data i s small enough to prevent a serious a l i a s i n g problem. There are two ways to consider the nature of these models. F i r s t , they can be viewed as models of a s i n g l e parameter (density), which i s 27 Gravity s t a t i o n x • l d 1 3 D 2 4 I b <J>1> $2> Q^y fa a r e angles to respective corners r ^ , r ^ , r ^ , r ^ are distances to respective corners F i g . 2. Geometric q u a n t i t i e s required f o r the g r a v i t y e f f e c t of a prism ( a f t e r Parasnis, 1962, p.43). s p e c i f i e d at regular i n t e r v a l s i n the subsurface plane. A l t e r n a t e l y the density of each c e l l might be considered as an independent parameter, and i n t h i s case we consider properties of the N-dimensional parameter space, where N i s the number of prisms i n the model. gular prism must be calculated. They are obtained from (3) by i n t e g r a t i n g only over the rectangular prism, where the density i s assumed constant. Parasnis (1962, p. 44) gives the following r e l a t i o n s h i p f o r the v e r t i c a l g r a v i t y Before the models can be used, the Frechet kernels f o r a rectan-g (x) = 2yp[xln + b i n The various geometric q u a n t i t i e s are shown i n Figure 2. The Frechet k e r n e l 28 follows immediately from (8) as g.Cx) GCx-Vzo> = "V" ( 9 ) where the geometric quantities in (8) would correspond to a prism centered at ( X Q , Z o ) . Equations (8) and (9) give an exact expression for the Frechet kernel, but computations can be made much simpler by approximating each c e l l as a point mass located at i t s center (in the two-dimensional case, a "point" mass is of course a line mass). From (3) the gravity effect of this source is 2yz (mass) = f \ l + 2 <10> (x-x ) + z o o where the mass i s equivalent to that of the prism, i.e. (mass) = p(dx)(dz) (11) which i s expressed i n gm/length for the two-dimensional model. To ensure that the approximation i s valid, the results of Equation (8) and (10) must be compared for different cells in the model. Since we w i l l always employ the same geometry of c e l l location and have dx = dz, a simple correction can be made to (10) to give the same results as (8), i f the agreement i s not acceptable. We consider distances measured 3 in kilometers, densities in gm/cm , and gravity in mill i g a l s ; hence the numerical value of y must be 6.67 (see Appendix B). The comparison was made for prisms at various depths, and for prism 29 dimension 1 km or 2 km. The gravity was computed at 1 km surface inter-vals, with x=0 being the ground position directly above the center of the c e l l . The station spacing (and related prism dimension) is unimportant, as we require only the ratio of the gravity effects. The results are displayed in Table I and II. The approximation is almost always accurate within 0.5%, however for gravity stations near the surface c e l l s , the d i f f e -rence i s great enough to warrant a correction; as should be expected from the inverse-square nature of gravity. The Frechet kernels for the linear model are now obtained from (10), applying the correction factor i f the approximation i s off by more than 0.5%. The kernel is written G(i,j,k) - 2T <* d* f (12) ( x i - X j > + «k and the gravity effect at surface position x^ i s g (i) = Z Z G(i,j,k) p(j,k) (13) j k where p(j,k) is the density of the prism centered at (x^.,^). The correction factor " f " i s usually 1, the other values that were used are shown in Table III. Using this formulation, the gravity effect of different models i s obtained simply by evaluating the inner product (13) with different p(j,k); the Frechet kernels need be computed only once and then retained for later use. In this procedure, the time saved by computing the approximate kernel IN THE FOLLOWING TABLE "X" IS THE DISTANCE FROM THE GROUND POSITION DIRECTLY ABOVE THE CENTER OF PLATE OR CYLINDER THE DENSITY OF THE PLATE IS 1.00 THE WIDTH OF THE PLATE IS 1.00 THE HEIGHT OF THE PLATE IS 1.00 THE DEPTH TO THE CENTER IS 0 .5 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.0 23.1051 26.6800 -3.5749 0.86601 1.0 5. 2370 5.3360 -0.0990 0.98145 2.0 1.5638 1.5694 -0.0056 0.99643 3.0 0.7204 0.7211 -0.0006 0.99912 4.0 0.4104 0.4105 -0.0001 0.99979 14.0 0.0340 0.0340 -0.0000 0.99916 15.0 0.0296 0.0296 -0.0000 0.99936 THE DEPTH TO THE CENTER IS 1 .5 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.0 8. 8645 8.8933 -0.0288 0.99676 1.0 6.1684 6.1569 0.0115 1.00187 2.0 3.2018 3.2016 0.0002 1.00005 3.0 1.7783 1.7787 -0.0004 0.99978 THE DEPTH TO THE CENTER IS 2 .5 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.0 5.3337 5.3360 -0.0023 0.99957 1.0 4.6005 4.6000 0.0005 1.00011 2.0 3.2543 3.2537 0.0006 1.00020 THE DEPTH TO THE CENTER IS 9 .5 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.0 1.4043 1.4042 0.000 1 1.00005 1.0 1.3888 1.3888 0.0000 1.00001 2.0 1.3445 1.3446 -0.0001 0.99992 11.0 0.4428 0.4427 0.0000 1.00007 15.0 0.4018 0.4020 -0.0002 0.99959 TABLE I. Gravity- effects of prisms (in the standard model) and equivalent cylinders. 31 THE DENSITY OF THE PLATE IS 1.00 THE WIDTH OF THE PLATE IS 2.00 THE HEIGHT OF THE PLATE IS 2.00 THE DEPTH TO THE CENTER IS 1.0 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.5 43.3750 42.6880 0.6871 1.01609 1.5 16.1710 16.4184 -0.2475 0.98493 2.5 7.2643 7.3600 -0.0957 0.98700 3.5 4.0053 4.0272 -0.0218 0.99458 4.5 2.5050 2.5111 -0.0061 0.99757 5.5 1.7055 1.7075 -0.0020 0.99880 6.5 1.2329 1.2338 -0.0008 0.99934 7.5 0.9317 0.9321 -0.0003 0.99964 8.5 0.7283 0.7285 -0.0002 0.99973 9.5 0.5847 0.5848 -0.0001 0.99982 10.5 0.4796 0.4796 -0.000 1 0.99986 THE DEPTH TO THE CENTER IS 3.0 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.5 17.2690 17.3059 -0.0369 0.99787 1.5 14.2519 14.2293 0.0226 1.00159 2.5 10.5119 10.4970 0.0149 1.00142 3.5 7.5358 7.5332 0.0026 1 .00035 4.5* 5.4722 5.4728 -0.0006 0.99988 5.5 4.0776 4.0785 -0.0009 0.99978 6.5 3.1229 3.1235 -0.0007 0.99979 7.5 2.4529 c 2.4533 -0.0004 0.99983 THE DEPTH TO THE CENTER IS 5.0 KM GRAVITY EFFECT (MGAL) X (KM) PLATE CYLINDER DIFFERENCE CORRECTION FACTOR 0.5 10.5624 10.5663 -0.0039 0.99963 1.5 9.7904 9.7908 -0.0004 0.99996 2.5 8. 5394 8.5376 0.0018 1.00021 3.5, 7.1640 7. 1624 0.0016 1.00023 4.5 5.8970 5.8961 0.0009 1.00014 5.5 4.8292 4.8290 0.0003 1 .00005 6.5 3.9673 3.9673 0.0000 1.00000 7.5 3.2836 3.2837 -0.0001 0.99997 TABLE i i . Gravity effects of prisms (in the large-block models) and equivalent cylinders. 32 Standard model: z ,x = subsurface p o s i t i o n n m x. = surface p o s i t i o n 3 n 1 I 1 0 0.86601 1 1 " 0.98145 2 0 0.99676 Large-block model: xx^ = surface p o s i t i o n define XA = |x -xx.|/dx ' m 2 n XA f 1 0.5 1.01609 1 1.5 0.98493 1 2.5 0.98700 1 3.5 0.99458 2 0.5 0.99787 Note: These values are v a l i d only f o r dx = dz i n Figure 1. TABLE I I I . Correction factors used i n computing Frechet kernels. may not be important, p a r t i c u l a r l y i f a large computer i s a v a i l a b l e . On the IBM 360/67, about 5 seconds would be used i n computing the exact kernels f o r a 300 c e l l model and a 30 s t a t i o n p r o f i l e ; the approximate kernels are obtained i n about 1/25 of th i s time (0.2 s e c ) . The c a l c u l a t i o n s i n Tables I and II also demonstrate that i s o l a t e d anomalies need not be modelled to a great depth ( r e l a t i v e to the p r o f i l e 33 length). I f a r e s i d u a l anomaly i s to be nearly zero at the ends of the p r o f i l e , prisms at a depth greater than about one-third the p r o f i l e length cannot make a s i g n i f i c a n t contribution to the gr a v i t y . For example, a un i t c e l l at depth 9.5 km. has a gr a v i t y e f f e c t of 1.4 mgal at x = 0 km., and 0.4 mgal at x = 15 km., a r a t i o of only 3.5; we conclude that a c e l l at that depth cannot make a large c o n t r i b u t i o n at the center of the p r o f i l e , i f the t o t a l g r a v i t y at the end of the p r o f i l e i s small. For this reason, the models studied l a t e r w i l l have a maximum c e l l depth of z = 9.5 (to center) f o r a 30-unit p r o f i l e length. From s i m i l a r considerations, i n many cases c e l l s near the ends of the p r o f i l e need not be considered. We can now summarize the advantages of the suggested l i n e a r density models. L i n e a r i t y i s the main b e n e f i t ; i t allows easy computation of the g r a v i t y of any model from the Frechet kernels. The point mass ( a c t u a l l y l i n e mass) approximation s i m p l i f i e s the computation of the kernels. If an i n i t i a l model i s used, a model which f i t s the data can be obtained i n one step. Unfortunately, there are two major disadvantages. Since the model consists of blocks, only rough approximations to complex shapes are po s s i b l e ; however most g r a v i t y methods share t h i s d i f f i c u l t y , since the non-uniqueness problem prevents precise g r a v i t y i n t e r p r e t a t i o n . The v a r i a b l e density between blocks of the model may obscure g e o l o g i c a l i n t e r p r e t a t i o n , as the in v e r s i o n w i l l not n e c e s s a r i l y i n d i c a t e an anomalous str u c t u r e of only one or two d e n s i t i e s . In most uses of gr a v i t y exploration, a si n g l e - d e n s i t y model i s desired, as the exploration o b j e c t i v e i s usually a body of a s p e c i -f i c mineral ore, whose density should be e s s e n t i a l l y constant. The in v e r s i o n procedures tested here w i l l therefore attempt to f i n d models i n which the 34 density contrast of i n d i v i d u a l c e l l s i s e i t h e r zero or a constant value, i n hopes of def i n i n g a simple anomalous body. 35 CHAPTER 4 Backus-Gilbert Approaches To Gravity Modelling The techniques developed by Backus and G i l b e r t have mainly been used to construct models of the whole earth, but there i s no inherent l i m i t a t i o n against using them f o r more r e s t r i c t e d examples. In p a r t i c u l a r , l i n e a r density models of the type proposed i n Chapter 3 are e a s i l y obtained from Backus-Gilbert i n v e r s i o n s . I t was previously observed (Section 2.3.a) that the averaging procedures used by Backus and G i l b e r t are not of much b e n e f i t to the suggested g r a v i t y models, since a l l l i n e a r models have the same average p r o p e r t i e s , and the only unique properties obtainable from gravity data are the t o t a l mass and center of mass coordinates. For th i s reason, attention w i l l be confined to the technique of f i n d i n g p a r t i c u l a r models by optimizing some property of the model under the constraints imposed by the observations ( i n the present case, values of the v e r t i c a l g r a v i t y component at s p e c i f i c points along a surface p r o f i l e ) . The inherent lack of r e s o l u t i o n of gr a v i t y data also prevents any extensive analysis of the e f f e c t s of data errors on the models obtained. The approach used here i s merely to seek approximate density models which s a t i s f y a g r a v i t y anomaly wit h i n a s p e c i f i e d e r r o r . The generalized inverse gives r e s u l t s equivalent to the Backus-G i l b e r t "distance" minimization. For the proposed gravity models, the Backus-Gilbert formulation i s more compact, since taking inner products r e -duces the dimension of matrices which need to be inverted. The generalized 36 inverse uses a matrix of dimension N, where N i s the number of parameters; while the inner product matrix of a Backus-Gilbert inverse has dimension K, where K i s the number of data. In addi t i o n , the e r r o r and r e s o l u t i o n analyses of the generalized inverse cannot be used to f u l l advantage. The generalized inverse was thus not pursued f u r t h e r . 4.1 Models "Clo s e s t " to an I n i t i a l Estimate The f i r s t method (following Backus and G i l b e r t , 1967) i s to f i n d the model c l o s e s t to an i n i t i a l guess. The s o l u t i o n i s obtained by mini-mizing the distance of an acceptable model from a given i n t i a l model; the constraints imposed by the observed g r a v i t y are incorporated v i a Lagrange m u l t i p l i e r s . We then seek a stat i o n a r y point of Z .IgjCMg) - + (G^M-Mg)] where M, M denote the f i n a l and i n i t i a l models G J i s the Frechet k e r n e l for g r a v i t y at s t a t i o n j Aj = gra v i t y observed at s t a t i o n j g.(Mp) = gr a v i t y of i n i t i a l model 3 (G jM-M-) denotes an inner product J « J = number of gr a v i t y data used a are Lagrange m u l t i p l i e r s and the f a c t o r — w i l l cancel out l a t e r There are two b a s i c steps i n the i n v e r s i o n ; f i r s t , solve f o r the Lagrange m u l t i p l i e r s a. i n 37 For t h e . l i n e a r models of Figure l a , the inner product (G.,G ) i s 3 k (G.,G.) = n G, G. j k inm knm J n m J where {n,mj defines the subsurface p o s i t i o n of each prism. The model M which s a t i s f i e s the data i s now given by M = M G + ^ a k G k C 2 ) To f i n d the d e n s i t i e s f o r each block i n the model, (2) i s computed f o r each subsurface p o s i t i o n x m» Z q. A complete d e r i v a t i o n of these equations i s given i n Appendix C. The system of equations (1) and (2) i s e a s i l y programmed f o r d i g i -t a l computer a p p l i c a t i o n . The spacing of gr a v i t y s t a t i o n s i s s p e c i f i e d to e s t a b l i s h the s i z e and p o s i t i o n of c e l l s i n the model, and to compute the Frechet kernels. Given regularly-spaced g r a v i t y observations, and a set of i n i t i a l d e n s i t i e s f o r each prism, an exact inverse i s obtained immediately. Since the models are l i n e a r , any i n i t i a l model w i l l s t i l l y i e l d a s o l u t i o n which f i t s the data. To t e s t the program, a g r a v i t y p r o f i l e was computed f o r the a r t i -f i c i a l body shown i n Figure 3. I t i s a f a i r l y complex object, and should provide a reasonable example f o r inversions which cannot i n d i c a t e the precise 38 o I o I • I o VO # # «-) U> S3 O o as • H O X *— 6-" H > £K O O O t a © W 00 o « w to 0 O o • 1 o I I • * o o • in o o o o in o o o CM o o • in (N DISTANCE ALONG PROFILE IN KB. o o o p = 1.0 I 1 1 z = 6.0 km (a) Shape of the anomalous body. F i g . 3. An a r t i f i c i a l g r a v i t y p r o f i l e . I N I T I A L MODEL (KM) / Z= 0.50 1.50 2.50 3.50 4.50 5.50 6.50 8.00 : : 0.0 0.0 0.0 0.0 0.0 0.0 0.0 9.00 : 0.0 0.0 0.0 0.0 0.0 0. 0 0.0 10.00 • 0.0 0.0 0.0 0.0 0.0 0.0 0.0 11.00 : 0.0 1.000 1.000 1.000 0.0 0.0 0.0 12.00 . : 0.0 1.000 1.000 1.000 0.0 0.0 0.0 13.00 j : 0.0 1.000 1.000 1.000 0.0 0.0 0.0 14.00 . : 0.0 1.000 1.000 1.000 2.000 2.000 2.000 15.00 : 0.0 1.000 1.000 1.000 2.000 2. 000 2. 000 16.00 j : 1.000 1.000 1.000 1.000 2.000 2.000 2.000 17.00 • : 1.000 1.000 1.000 1.000 0.0 0.0 0.0 18.00 : 1.000 1.000 1.000 1.000 0.0 0.0 0.0 19.00 : 1.000 1.000 1.000 1.000 0.0 0.0 0.0 20.00 « 0.0 0.0 0.0 1.000 0.0 0.0 0.0 21.00 : : 0.0 0.0 0.0 1.000 1.000 1.000 1.000 22.00 : : 0.0 0.0 0.0 1.000 1.000 1.000 1.000 23.00 ' : 0.0 0.0 0.0 0.0 0.0 0.0 0.0 (b) Computer representation of the body. 40 d e t a i l of the earth. Figure 4 shows the model obtained by using a "zero" i n i t i a l model. The inverse s o l u t i o n cannot be represented as an anomalous body, since the computed d e n s i t i e s are d i f f e r e n t f o r each subsurface prism. A r e l a t e d problem i s that no r e s t r i c t i o n s on density are b u i l t i n t o the program, and both negative and p o s i t i v e d e n s i t i e s are obtained, p a r t i c u l a r l y on the edges of the model region. In most exploration s i t u a t i o n s , i t would be reasonable to assume that a body of constant density contrast i s causing the g r a v i t y anomaly. To reduce the ambiguity of the s o l u t i o n s , d i f f e r e n t s t r a t e g i e s were employed to seek, constant density models. A surface g r a v i t y p r o f i l e should extend beyond the l i m i t s of an anomalous body, i f the g r a v i t y anomaly i s to be completely defined. I t i s therefore reasonable to consider models which span a region only under the center of the p r o f i l e . This i s simply achieved by s p e c i f y i n g l i m i t s on x^ and z^, which reduces the number of parameters and thus should help remove negative d e n s i t i e s ( i n the case of p o s i t i v e anomalies), since the constraint of known mass must also be s a t i s f i e d . Figure 5 shows the i n v e r s i o n obtained by using a model of 105 prisms, rather than 300 as i n Figure 4. Variable density i s s t i l l present, but sharp changes i n density are evident, p a r t i -c u l a r l y i n the near-surface region of the model. To transform these models i n t o a rough p i c t u r e of a s i n g l e - d e n s i t y body, a program was w r i t t e n to replace the exact d e n s i t i e s by d i s c r e t e values at a s p e c i f i e d increment, simultaneously r e j e c t i n g any negative values. The r e s u l t i n g approximate models do not s a t i s f y the observed gra-v i t y values p r e c i s e l y ; however they can i n d i c a t e the extent and general shape COMPUTED DEN S I T I E S X (KM) / Z= 0. 50 1.50 1.00 : -0.127 -0.063 2.00 . : -0.096 -0.072 3.00 : -0.104 -0.076 ' 4.00 . -0.118 -0.081 5.00 : -0.135 -0.086 • 6.00 : -0.153 -0.086 7.00 : -0.167 -0.075 8.00 « -0.164 -0.040 9.00 . : -0.111 0.040 10.00 : 0.089 0. 187 11.00 : 0.455 0. 378 12.00 j 0.684 0.542 13.00 : 0.776 0.651 14.00 : 0.791 0.727 15.00 : 0.720 0.846 16.00 : 1 .582 1.050 17.00 : 1.452 1.099 18.00 j 1.366 1.013 19.00 : 1.255 0.784 20.00 j 0.061 0.405 21.00 : -0.010 0. 181 22.00 : -0.022 0.085 23.00 : -0.048 0.029 24.00 : -0.078 -0.010 25.00 : -0.099 -0.037 26.00 : -0.106 -0.053 27.00 : -0.105 -0.061 28.00 : -0.100 -0.064 29.00 : -0.097 -0.065 30.00 : -0.128 -0.058 2 .50 3.50 0.036 - 0 . 016 0. 043 - 0 . 019 0. 046 - 0 . 018 0. 045 - 0 . 0 13 0. 041 - 0 . 003 0. 030 0. 0 15 0.005 0. 044 0. 0 39 0. 089 0. 114 0. 152 0. 221 0. 233 0. 346 0. 324 0. 466 0. 414 0. 565 0. 495 0. 650 0. 565 0. 739 0. 626 0. 820 0.666 0. 839 0. 667 0. 772 0. 6 18 0. 623 0. 523 0. 430 0. 405 0. 268 0. 292 0. 162 0. 20 3 0. 093 0. 136 0. 046 0. 088 0. 013 0. 053 0. 008 0. 028 0. 022 0. 011 0. 029 0. 000 0. 032 - 0 . 005 0. 028 - 0 . 005 4 .50 5 . 50 0. 001 0. 016 0. 002 0. 019 0. 006 0. 025 0. 0 14 0. 036 0. 027 0. 051 0. 048 0. 072 0. 078 0. 100 0. 118 0. 136 0. 171 0. 180 0. 234 0. 230 0.303 0. 283 0. 372 0. 336 0. 435 0. 384 0. 489 0. 425 0. 531 0. 4 54 0. 553 0. 467 0. 548 0. 462 0. 512 0. 436 0. 449 0. 3 92 0. 370 0. 337 0. 290 0. 278 0. 220 0. 223 0. 162 0. 175 0. 117 0. 135 0. 082 0. 103 0. 056 0. 078 0. 03 8 0. 059 0.025 0. 045 0. 017 0. 035 0. 013 0. 029 6 .50 7.50 0. 028 0 .039 0. 03 3 0 .045 0. 04 1 0 .054 0. 05 3 0 .066 0. 068 0 .00 1 0. 088 0 . 100 0. 114 0 . 123 0. 145 0 . 149 0. 182 0 . 180 0. 222 0 .212 0. 264 0 . 245 0. 305 0 .277 0. 342 0 . 306 0. 372 0 . 329 0. 393 0 . 344 0. 402 0 . 350 0. 396 0 . 346 0. 377 0 .33 1 0. 346 0 .308 0. 306 0 .278 0. 262 0 .245 0. 219 0 .212 0. 180 0 . 180 0. 146 0 . 150 0. 117 0 . 125 0. 093 . 0 . 104 0. 075 0 .086 0. 060 0 .072 0. 049 0 .06 1 0. 041 0 .052 8.50 9.50 0.047 0.054 0.055 0.062 0.064 0.071 0.075 0.082 0.090 0.096 0.107 0.111 0. 127 0. 129 0.150 0.148 0. 175 0. 169 0.202 0.191 0.228 0.212 0.253 0.232 0.275 0.250 0.293 0.263 0.304 0.272 0.309 0.275 0.305 0.272 0.294 0.264 0.277 0.251 0.254 0.233 0.229 0.213 0.202 0.192 0. 176 0. 171 0.152 0.150 0. 1 30 0. 1 31 0.110 0.114 0.094 0.099 0.081 0.087 0.069 0.076 0.060 0.067 F i g . 4. An i n v e r s i o n of the a r t i f i c i a l p r o f i l e u s i n g a zero i n i t i a l model. -OHPUTED DENSITIES (KH) / Z= 0.50 1.50 2.50 3.50 4 .50 5 .50 6.50 8.00 : 0.011 -0.043 0.240 -0.105 -0. 319 -0. 210 0. 105 9.00 : 0.044 0.010 0.032 -0.014 -0. 029 0. 083 0.311 10. 00 : 0.131 -0.078 -0.033 0.029 0. 115 0. 255 0.456 11.00 : 0.257 0.723 0.910 1.053 0. 196 0. 362 0.559 12.00 : 0.369 0.654 0.885 1.074 0. 252 0. 437 0. 636 13.00 : 0.435 0.570 0.899 1. 105 0. 299 0. 495 0.698 14.00 « : 0.441 0.706 0.941 1.146 0. 344 0. 544 0.751 15.00 : 0.364 0.807 1.010 1.194 0. 385 0. 586 0. 796' 16.00 : : 1.217 1.003 1.084 1.233 0. 415 0. 618 0. 835 17.00 : 1.082 1.049 1. 104 1.24 1 0. 426 0. 637 0. 864 18.00 : : 0.998 0.967 1.042 1.204 0. 412 0. 641 0. 883 19.00 : 0.910 0.7'46 0.895 1. 125 0. 380 0. 634 0. 891 20.00 : -0.212 0.366 0.674 1.025 0. 354 0. 625 0. 884 21.00 : -0.155 0.073 0.350 0.970 0. 398 0. 626 0. 846 22.00 : : -0.034 -0. 110 -0.511 1.275 0. 654 0. 610 0.717 I n i t i a l model: density 1.0 i n 48 prisms ( 11 « x « 22, 0.5 < z < 3.5 ) Fi g . 5. A more r e s t r i c t e d model f o r Figure 3. APPROXIMATE DENSITIES X / Z= 0.50 1.50 2 . 50 3. 50 4. 50 5. 50 6 .50 8.00 : : 0.0 0.0 3. D 0. 0 0.0 0. 0 0 .0 9.00 : 0.0 0.0 0 .0 0. 0 0.0 0. 0 0 .0 10.00 : : 0.0 0.0 0 .0 0. 0 0.0 0. 0 0 .0 11.00 : 0.0 1.00 1 .00 1. 00 0.0 0. 0 1 .00 12.00 , : 0.0 1.00 1 .00 1. 00 0.0 0. 0 1 .00 13.00 : 0.0 1.00 1 .00 1. 00 0.0 0. 0 1 .00 14.00 : . 0.0 1.00 1 .00 1. 00 0.0 1. 00 1 .00 15.00 j 0.0 1.00 1 .00 1. 00 0.0 1. 00 1 .00 16.00 : 1.00 1.00 1 .00 1. 00 0.0 1. 00 1 .00 17.00 : 1.00 1.00 1 .00 1. 00 0.0 1. 00 1 .00 18.00 : 1.00 1.00 1 .00 1. 00 0.0 1. 00 1 .00 19.00 : 1.00 1.00 1 .00 1. 00 0.0 1. 00 1 .00 20.00 « : 0.0 0.0 1 .00 1. 00 0.0 1. 00 1 .00 21.00 : 0.0 0.0 0 .0 1. 00 0.0 1. 00 1 .00 22.00 i 0.0 0.0 0 .0 1. 00 1.00 1. 00 1 .00 EXACT GRAVITY - APPROX GRAVITY - ERROR 1 12.1100 11.7917 0.3183 2 13.8700 13.4755 0.3945 3 16.0400 15.5472 0.4928 4 18.7600 18. 1360 0.6239 5 22.2500 21.4309 0.8191 6 26.8200 25.7176 1. 1024 7 32.9700 31.4474 1 .5226 8 41.5300 39.3640 2.1660 9 53.8500 50.7185 3.1315 10 71.8000 67.2621 4.5379 11 94. 1600 87.7264 6.4336 12 113.2600 104.4841 8.7759 13 127.2900 116.2606 11.0294 14 137.4700 125.0552 12.4148 15 146.6200 134.5364 12.0836 16 '168.8500 159. 1159 9.7341 17 168.5000 162.6880 5.8120 18 159.3200 158.0511 1. 2688 19 140.4500 143.3264 -2.8764 20 99.3600 104.7375 -5.3775 21 78.6500 84. 1919 -5.5419 22 64.7900 69.3524 -4.5624 23 53.4800 57.0884 -3.6084 24 43.9100 46.8110 -2.9010 25 36.0500 38.4215 -2.3715 26 29.7800 31.7336 -1.9536 27 24.8300 26.4540 -1.6240 28 20.9300 22.2833 -1.3533 29 17.8300 18.9673 -1.1373 30 15.3400 16.3065 -0.9665 F i g . 6. Approximate model derived from Figure 5. 4 4 of an anomalous body. In some cases, the approximate gravity profile f i t s the real values within 10%, a degree of accuracy which may be satisfactory in modelling exploration data (C. Ager, personal communication, 19 73). Approximate models (for the inversion shown in Figure 5) and their gravity effects are given by Figure 6. These approximations are not always geolo-gically r e a l i s t i c ; in the example shown, the deeper section of the model i s not connected to the near-surface features. The distance-minimization algorithm of Equations (1) and (2) can also be used for the simpler models consisting of larger prisms (dimension two station spacings), i l l u s t r a t e d i n Figure 1. These models are of course a poorer means of representing a complex geological formation, but could be useful in estimating gross features. The computer program is very similar, with only the adaptations necessary for different geometry. Equations (1) and (2) have the identical form here; the differences in size and position of blocks are incorporated in the Frechet kernels. As before, approximate solutions are constructed as ah aid i n interpretation. Attempts to invert the previously-used a r t i f i c i a l data with the simpler models were not successful. The usual results were very unrealistic models, with large positive and negative densities. One such example i s shown in Figure 7. To test the vali d i t y of the method, a r t i f i c i a l data were generated from simple bodies which can be exactly described by a number of these large prisms; two examples.are shown in Figure 8. Inversions of these data were 45 COMPUTED DENSITIES X (KM) / Z = 1.00 3.00 5.00 7.00 1.50 : : - 0 . 238 4. 101 -2 . 729 -2.922 3.50 : - 0 . 024 -0 . 942 -0 . 884 -1.495 5.50 j - 0 . 226 3. 184 0. 120 -0.647 7.50 : - 0 . 026 -1 . 962 - 0 . 148 -0.227 9.50 : : - 0 . 228 1. 947 0. 402 0. 135 11.50 : 0. 028 -0. 767 0. 511 0.450 13.50 : : - 0 . 168 3. 642 1. 049 0.772 15.50 • : - 0 . 267 -4. 920 0. 973 1.110 17.50 : : - 1 . 562 14. 342 3. 280 1.438 19.50 : - 0 . 279 -4. 824 0. 602 0.985 21.50 j : - 0 . 010 -0 . 092 0. 100 0.608 23.50 : - 0 . 026 -0. 123 0. 080 0.482 25.50 : : - 0 . 017 - 0 . 074 - 0 . 250 0.400 27.50 : - 0 . 019 0. 432 - 1 . 832 0.510 BVED GRAVITY - MODEL GRAVITY - ERROR 1 12.1100 12.1081 - 0 . 0019 2 13.8700 13.8557 - 0 . 0143 3 16.0400 16.0697 0. 0297 4 18.7600 18.7783 0. 0183 5 22.2500 22.1984 - 0 . 0516 6 26.8200 26.8416 0. 0216 7 32.9700 33.0330 0. 0630 8 41.5300 41.5146 - 0 . 0154 9 53.8500 53.8385 - 0 . 0115 10 71.8000 71.8266 0.0266 11 94.1600 94.1672 0. 0072 12 113.2600 113.2385 - 0 . 0215 13 127.2900 127.2491 - 0 . 0409 14 137.4700 137.4915 0. 0215 15 146.6200 146.6283 0. 0083 16 168.8500 168.8568 0. 0068 17 168.5000 168.4854 - 0 . 0146 18 159.3200 159.3201 0. 0001 19 140.4500 140.5591 0. 1091 20 99.3600 99.4003 0. 0403 21 78.6500 78.6886 0. 0386 22 64.7900 64.8147 0. 0247 23 53.4800 53.5144 0. 0344 24 43.9100 43.9765 0. 0665 25 36.0500 36.1273 0. 0773 26 29.7800 29.8009 0. 0209 27 24.8300 24.9238 0. 0938 28 20.9300 20.9478 0. 0178 29 17.8300 17.9088 0. 0788 30 15.3400 15.4146 0. 0746 Fig. 7. An unsuccessful inversion using a large-block model. o o t o o I CM I ^ o -J O <: t o o E \£> X E-< H > O •oi O OS • U O (N a *-> « w VJ o ca o O • © 00 I I I o o t o o o * * o o « in o o o o t -in o o t o CN O o t in CM DISTANCE ALONG PROFILE IN KH. P = 1.0 z = 5.0 (a) Fig. 8. Simple bodies composed of large prisms. 47 — o <-A © •4) . o o E O H H > O «< o o i • o m r -a w > » W t/1 O CO o O . o m o o t m * * o o o o o o O o o o o o o • • ' t • • 1 in o in o in o *— *— CN CN cn DISTANCE ALONG PROFILE IN KH« 1 L , L _l _ z = 5.Q p = 1.0 (b) COMPUTED DENSITIES (KM) / Z= 1.00 3.00 5.00 7.0 3.50 : - 0 . 002 0. 049 -0.222 -0. 239 5.50 : 0. 000 - 0 . 107 0.057 0. 006 7.50 : : 0. 999 0. 851 0.410 0. 223 9.50 : 0. 999 0. 921 0.603 0. 375 11.50 : 0. 998 0. 992 0.679 0. 450 13.50 : 1. 001 0. 937 0.676 0. 457 15.50 : 0. 998 0. 937 0.609 0. 401 17.50 : 1. 002 0. 757 0.421 0. 280 19.50 : - 0 . 001 - 0 . 138 0. 116 0. 120 21.50 . : 0. 000 - 0 . 061 -0.041 - 0 . 025 23.50 : - 0 . 002 0. 058 -0.203 -0. 163 I n i t i a l model: a l l zero. APPROXIMATE DENSITIES X (KM) / Z= 1.00 3.00 5.00 7.00 3.50 : 0.0 0.0 0.0 0.0 5.50 : : 0.0 0.0 0.0 0.0 7.50 : 1.00 1.00 0.0 0.0 9.50 ; 1.00 1.00 1.00 0.0 11.50 : 1.00 1.00 1.00 0.0 13.50 : 1.00 1.00 1.00 0.0 15.50 : 1.00 1.00 1.00 0.0 17.50 : 1.00 1.00 0.0 0.0 19.50 : 0.0 0.0 0.0 0.0 21.50 : 0.0 0.0 0.0 0.0 23.50 : 0.0 0.0 o.o 0.0 9. I n v e r s i o n of data from F i g u r e 8Ca). I N I T I A L MODEL X (KM) / Z= 1.00 3.00 5.00 7.00 COHPUTED DENSITIES X (KM) / Z= 1.50 3.50 5.50 7. 50 9.50 11.50 13.50 15.50 17.50 19.50 21.50 23.50 25.50 27.50 1.00 3.00 5.00 7.0 0. 005 0.036 -0.026 0.019 0. 002 -0.051 0.038 0.089 0. 004 -0.005 0. 127 0. 183 0. 001 0.011 0.324 0.319 0. 004 1.311 0.700 0.471 0.001 1.407 0.796 0.533 0. 003 0.565 0.610 0.489 0. 002 0.562 0.486 0.412 0. 003 0.489 0.373 0. 326 0. 001 0. 137 0.221 0.235 0. 001 0.028 0.116 0. 160 0. 000 -0.036 0.058 0. 111 0. 000 -0.020 0.037 0.086 0. 001 -0.037 0.032 0.082 F i r s t r e s u l t from zero i n i t i a l model. 3.50 : 0.0 0.0 0.0 0.0 5.50 : : 0.0 • 0.0 0.0 0.0 7.50 : 0.0 0.0 0.0 0.0 9.50 : 0.0 1.000 1. 000 1.000 11.50 -: 0.0 1.000 1.000 1.000 13.50 : 0.0 0.0 0.0 1.000 15.50 : 0.0 0.0 0.0 1.000 17.50 : 0.0 0.0 0.0 1.000 19.50 : 0.0 0.0 0.0 0.0 21.50 : 0.0 0.0 0.0 0.0 23.50 : 0.0 0.0 0.0 0.0 APPROXIMATE DENSITIES X (KM) / Z= 1.00 3.00 5.00 7.00 3.50 : 0.0 0.0 0.0 0.0 5.50 ; 0.0 0.0 0.0 0.0 7.50 : 0.0 0.0 0.0 0.10 9. 50 : 0.0 1.00 1. 10 1. 10 11.50 -: 0.0 1.00 1.10 1.10 13.50 : 0.0 0. 30 0. 20 1.20 15.50 : : 0.0 0.40 0.20 1.20 17.50 : 0.0 0. 30 0. 20 1.20 19.50 : 0.0 0. 10 0.10 0.20 21.50 . : 0.0 0.0 0.0 0. 10 23.50 : 0.0 0.0 0.0 0.0 F i n a l r e s u l t . VO F i g . 10. Inversion of data from Figure 8(b). 50 very s u c c e s s f u l . Figure 9 shows the i n v e r s i o n of the data from Figure 8a; the approximate s o l u t i o n , s t a r t i n g from a zero i n i t i a l model, corresponds exactly to the body which produced the data. Figure 10 i l l u s t r a t e s the i n -version of data from Figure 8b. S t a r t i n g from a zero i n i t i a l model, the process was repeated three times, using p r i o r r e s u l t s to improve the i n i t i a l model i n each case. A shape corresponding to Figure 8b can be i n f e r r e d from the f i n a l r e s u l t . The more immediate success of the f i r s t example i s pro-bably a r e s u l t of density being concentrated near the surface. The surface blocks n a t u r a l l y have the greatest u n i t contribution to the g r a v i t y ; and thus t h e i r d e n s i t i e s are most c r i t i c a l i n f i t t i n g the data. In Figures 9 and 10, the computed den s i t i e s of surface prisms are usually very close to the d e n s i t i e s of the " r e a l " body, regardless of the i n i t i a l model. Two b a s i c conclusions follow from these r e s u l t s . F i r s t l y , the use of large-block models y i e l d s a more nearly unique inverse, since approximate models which p r e c i s e l y f i t the data can sometimes be obtained i n one step. However, i n i n v e r t i n g data from more complex s t r u c t u r e s , the simple model inverse i s often u n r e a l i s t i c ; t h i s i s l i k e l y a manifestation of an a l i a s i n g problem. I t appears that when the data contain frequencies which cannot be represented by the model, a p h y s i c a l l y u n r e a l i s t i c s o l u t i o n w i l l r e s u l t . The problem a r i s e s from the lower cutoff wavenumber of the large-block models; i . e . the model cannot sample the subsurface density as f i n e l y as the obser-vations sample the g r a v i t y f i e l d . 4.2 A Method of "Weighted-Distance" Minimization In seeking a model to represent a constant density body, one hopes 51 to f i n d compact configurations of the subsurface density i n the l i n e a r model. Using models with r e s t r i c t e d s p a t i a l range was one attempt i n t h i s d i r e c t i o n ; another i s to e s t a b l i s h some c r i t e r i o n f o r compactness, and then devise a scheme to s e l e c t a p a r t i c u l a r model s a t i s f y i n g that c r i t e r i o n . The Backus-G i l b e r t method lends i t s e l f to t h i s approach, i f a simple property of the model can be optimized. A simple view of compactness i s that anomalous density should be confined as much as p o s s i b l e to a r e s t r i c t e d subsurface region. A simple method i s then to minimize = f p 2 R 2 dV (3) J v where R i s a v a r i a b l e weighting f a c t o r , which i s smallest i n the most favoured subsurface regions. Following an a n l y t i c procedure s i m i l a r to that of Section 4.1, the following system i s obtained. The Lagrange m u l t i p l i e r s are the s o l u t i o n of J G1 °k and the acceptable model i s given by J G, P = ^ \ T (5) 1 * RZ Equation (5) i s evaluated at each x^, z^ to give a l l the prism d e n s i t i e s . A more complete d e r i v a t i o n i s shown i n Appendix C. 52 In the form (4) and (5), i t i s d i f f i c u l t to use the system repe-t i v e l y to improve the f i r s t s o l u t i o n ; so i t i s modified to allow use of i n i t i a l models as w e l l . Any R(x , z ) can be used i n (4) and (5), so long m n as i t i s not a function of p(x ,z ). To allow i n t e r a c t i v e use of the method, m n we consider R as an a r b i t r a r y weighting f a c t o r applied to each c e l l i n the model, and minimize the weighted "distance" of an acceptable model from an i n i t i a l guess. This use of R i s somewhat akin to the penalty functions of some optimization methods (see Appendix D); one poss i b l e a p p l i c a t i o n i s to maintain c e r t a i n d e n s i t i e s at t h e i r i n i t i a l value. We now obtain the Lagrange m u l t i p l i e r s by s o l v i n g J G, G, k With the resultant model being J V P = P G + I \ ^ 2 ( 7 ) The "weighted-distance" method i s more e f f e c t i v e than the previous examples, since one can improve an i n t i a l model while d i s c r i m i n a t i n g against those prisms which do not appear to contribute strongly to the anomaly. I t i s also a convenient way of keeping c e r t a i n blocks at a known density. The computer program, l i s t e d i n Appendix F as WEIGHT and SMOOTH, i s of course s i m i l a r to the previous ones, but re q u i r i n g a d d i t i o n a l input to e s t a b l i s h the weighting f a c t o r s . Rather than s p e c i f y i n g the numerical weight f o r each c e l l , a choice of four factors i s allowed. A weight of 1.0, WT, WT2, or WT3, i s assigned, depending on whether the i n d i c a t o r REG(N,M) i s 0.0, 1.0, 2.0, 53 or 3.0 r e s p e c t i v e l y . This method allows a simple display of v a r i o u s l y weighted regions, and f a c i l i t a t e s the use of d i f f e r e n t numerical factors i n the same region of the model. The program was used to i n v e r t a r t i f i c i a l g r avity p r o f i l e s , two of which were generated by the simple bodies i n Figure 11. Four or f i v e r e p e t i t i o n s of the in v e r s i o n were usually s u f f i c i e n t to construct a reaso-nable si n g l e - d e n s i t y model. Improvement i n the f i n a l stages was often a t r i a l - a n d - e r r o r i n t e r p r e t i v e process; however the i n i t i a l i n v e r s i o n from very simple models was the key to f i n d i n g an acceptable s o l u t i o n . In a l l cases, a p p l i c a t i o n was di r e c t e d towards f i n d i n g a body of a s i n g l e known density; i n any gra v i t y exploration, some density estimate should be av a i -l a b l e , even i f i t i s j u s t a guess to be tested. Inversions of the a r t i f i c i a l data are i l l u s t r a t e d i n Figures 12, 13, 14, and 15. The s t a r t i n g models were very simple; zero density f o r the data of Figure 11a; a large block of density 0.5 gm/cc f o r l i b . The f i r s t inversions were improved, p a r t l y by i n t e r p r e t i v e judgement, and p a r t l y through subsequent use of the program. In l a t e r stages, weighting factors are applied to keep c e r t a i n prisms at zero density, f o r example, those on the boundaries of the subsurface region. Figure 12 shows the i n i t i a l i n v e r -sion and an intermediate stage f o r the data of Figure 11a, the f i n a l model a f t e r f i v e r e p e t i t i o n s of the program i s displayed i n Figure 13. The model's gr a v i t y and the " r e a l " data usually agree w i t h i n 1.0 mgal; the worst e r r o r i s l e s s than 3.0 mgal. Figure 14 and 15 show the corresponding steps i n i n v e r t i n g the 54 1 o 1 o 1 • o vO >l a o ca >-) o «* o M « U o H =r C>4 w E-t OS <: o CM o O • o CN H > CC o O 1 O 1 1 * # • * * * * o o • in o o I o o o in o o o o o • in CN O o * o DISTANCE ALONG PROFILE IN KM. z = 5.0 km (a) Fig. 11. A r t i f i c i a l gravity profiles used with the WEIGHT program. GRAVITY OF A R T I F I C I A L BODY 0 . 0 30.00 60 .00 90.00 120.00 150.00 5.00 10.00 15.00 20.00 25.00 30.00 * * * 56 COMPUTED DENSITIES (KM) / Z= 0.50 1 .50 2.50 3.50 4 .50 5.50 7.00 : -0.004 - 0 . 005 -0.054 -0.008 -0.046 -0. 1 18 8.00 • -0.016 -0. 005 -0.002 0.004 -0.010 -0.039 9.00 : -0.029 0. 013 0.036 0.045 0.043 0.034 10.00 -0.022 0. 052 0.086 0.099 0. 103 0. 105 11.00 : 0.048 0. 129 0. 154 0.162 0.167 0. 173 12.00 : 0.242 0. 250 0. 236 0. 228 0. 228 0. 236 13.00 : 0.543 0. 376 C.310 0.285 0.281 0.292 14.00 : 0.613 0. 424 0. 348 0.321 0. 3.20 0. 337 15.00 : 0.438 0. 382 0.343 0.333 0.343 0.370 16.00 : 0.304 0. 308 0. 310 0.323 0. 350 0. 391 17.00 : 0.194 0. 232 0.264 0.299 0.344 0.400 18.00 : 0.098 0. 162 0. 214 0. 266 0. 325 0.397 19.00 : 0.026 0. 102 0. 163 0.004 0.006 0.008 20.00 : -0.015 0. 056 0. 113 0. 172 0. 244 0. 343 21.00 : -0.025 0. 022 0.061 0.106 0.168 0.272 22.00 : -0.015 -0. 008 -0.001 0.031 0.045 0. 141 23.00 -: -0.003 - 0 . 003 -0.086 0.010 -0. 182 -0.117 Ca) F i r s t result from a zero i n i t i a l model. I N I T I A L MODEL (KM) / Z= 0.50 1.50 2. 50 3. 50 4.50 5.50 10.00 : 0.0 0.0 0.0 0.0 0.0 0.0 11 .00 : 0.0 0.0 0.0 0.0 0.0 0.0 12.00 : : 0.0 0.0 0.0 0.0 0.0 0.0 13.00 : 0.0 0.800 0.800 0. 800 0.0 0.0 14.00 j 0.0 0. 800 0.800 0.800 0.0 0.0 15.00 : : 0.0 0.0 0. 800 0. 800 0. 800 0.0 16.00 : 0.0 0.0 0.0 0.800 0.800 0.0 17.00 : 0.0 0.0 0.0 0. 800 0.800 0. 800 18.00 : 0.0 0.0 0.0 0.0 0.800 0.800 19.00 • : 0.0 0.0 0.0 0.0 0. 800 0. 0 20.00 : 0.0 0.0 0.0 0.0 0.0 0.0 21.00 • : 0.0 0.0 0.0 0.0 0.0 0.0 (b) A subsequent inversion's i n i t i a l model. Fig. 1 2 . Inversion of the profile in Figure 1 1(a). APPROXIMATE DENSITIES X (KM) / Z= 0.50 1.50 2.50 3.50 4.50 5.50 10.00 : 0.0 0.0 0.0 0.0 0.0 0.0 11.00 : 0.0 0.0 0.0 0.0 0.0 0.0 12.00 : 0.0 0.0 0.0 0.0 0.0 0.0 13.00 : 0.0 1.00 1.00 1.00 0.0 0.0 14.00 : 0.0 1.00 1.00 1 .00 0.0 0.0 15.00 : 0.0 0.0 1.00 1.00 1.00 0.0 16.00 : 0.0 0.0 0.0 1 .00 1.00 0.0 17.00 : 0.0 0.0 0.0 1.00 1.00 1.00 18.00 ; 0.0 0.0 0.0 0.0 1.00 1.00 19.00 : 0.0 0.0 0.0 0.0 1.00 0.0 20.00 : 0.0 0.0 0.0 0.0 0.0 0.0 21.00 : 0.0 0.0 0.0 1.00 1.00 0.0 EXACT GRAVITY - APPROX GRAVITY - ERROR 1 3.1800 3. 8327 -0. 6527 2 3.6600 4. 3786 - 0 . 7186 3 4.2600 5. 0502 -0. 7902 4 5.0200 5. 8894 -o. 8694 5 6.0100 6. 9565 -o. 9465 6 7.3200 8. 3422 - 1 . 0222 7 9.1000 10. 1855 - 1 . 0855 8 11.6200 12. 7081 - 1 . 0881 9 15.3100 16. 2755 -0. 9655 10 20.8900 21. 5070 -o. 6171 11 29.3400 29. 4415 -0. 1015 12 41.2800 41. 3973 -o. 1173 13 54.0400 55. 1 149 -1 . 0749 14 59.1600 60. 9357 - 1 . 7757 15 55.5800 57. 1730 -1. 5930 16 49.5500 51. 1622 - 1 . 6123 17 42.5000 45. 5201 -3. 0 201 18 35.0000 40. 2489 - 5 . 2489 19 27.9100 35. 3689 -7. 4589 20 21.8000 30. 9025 - 9 . 1025 21 16.9500 26. 6177 -9. 6677 22 1 3. 2700 22. 3314 - 9 . 0614 23 10.5400 18. 2792 -7. 7392 24 8.5100 14. 8068 - 6 . 2968 25 6.9800 12. 0270 -5. 0470 26 5. 8100 9. 8659 - 4 . 0559 27 4.9100 8. 1947 -3. 2847 28 4. 1900 6. 89 3 5 - 2 . 7035 29 3.6200 5. 8688 -2. 2488 30 3. 1500 5. 051 1 - 1 . 9011 Cc) The approximate model obtained from Ch) • 58 o o o M o w a o & o o • o c o o CN • # # * * •M- * O O in o o o o • m o o o CN O O in CN O o o ro DISTANCE ALONG PROFILE IN KM. Ftg. 13. A f i n a l model for the profile of Figure 11(a). COMPUTED DENSITIES 59 (KM) / Z= 0 .50 1 .50 2. 50 3. 50 4.50 5.50 6.00 : -0. 005 -0. 159 0.047 -0.305 -0.254 -0.019 7.00 : - 0 . 073 0. 021 0.058 0.000 0.017 0. 125 8.00 : -0. 089 0. 221 0. 216 0. 184 0. 188 0.240 9.00 : 1. 072 0. 557 0.390 0.321 0.305 0. 327 10.00 : 1. 235 0. 673 0. 486 0. 409 0. 384 0. 389 1 1.00 : 0. 382 0. 537 0.490 0.451 0.433 0.433 12.00 : 0. 422 0. 455 0. 473 0. 470 0. 464 0.464 13.00 : 0. 312 0. 433 0.476 0.488 0.489 0.488 14.00 : 0. 371 0. 473 0. 508 0. 517 0. 514 0.509 15.00 : 0. 591 0. 562 0.563 0.556 0.542 0. 529 16.00 : 0. 662 0. 647 0.631 0. 600 0.570 0.545 17.00 : 0. 592 0. 771 0.711 0.643 0.591 0.554 18.00 . : 1. 478 0. 982 0.785 0.668 0. 596 0. 552 19.00 : 1. 413 1. 034 0.788 0.652 0.575 0.532 20.00 . : 1. 4 18 0. 887 0. 688 0. 580 0. 521 0.494 21.00 ' : 0. 388 0. 555 0.513 0.463 0.437 0.438 22.00 : 0. 186 0. 304 0. 333 0. 316 0. 323 0.368 23.00 : 0. 049 0. 159 0.210 0. 129 0. 163 0.292 24.00 : 0. 010 -0. 005 0. 294 -0.251 -0.096 0.237 APPROXIMATE DENSITIES X (KM) / Z= 0.50 1.50 2.50 3.50 4.50 5.50 6.00 : 0.0 0.0 0.0 0.0 0.0 0.0 7.00 : 0.0 0.0 0.0 0.0 0.0 0.0 8.00 : 0.0 0.0 0.0 0.0 0.0 0.0 9.00 : 1.00 1.00 0.0 0.0 0.0 0.0 10.00 • : 1.00 1.00 0.0 0.0 0.0 0.0 11 .00 : 0.0 1.00 0.0 0.0 0.0 0.0 12.00 : : 0.0 0.0 0.0 0.0 0.0 0.0 13.00 : 0.0 0.0 0.0 0.0 0.0 0.0 14.00 : 0.0 0.0 1.00 1 .00 1.00 1.00 15.00 . : 1.00 1.00 1.00 1.00 1.00 1.00 16.00 : 1.00 1.00 1.00 1 .00 1.00 1.00 17.00 : : 1.00 1.00 1.00 1.00 1.00 1.00 18.00 : 1.00 1.00 1.00 1 .00 1.00 1.00 19.00 j 1.00 1.00 1.00 1.00 1.00 1.00 20.00 : 1.00 1.00 1.00 1 .00 1.00 0.0 21.00 : 0.0 1.00 1.00 0.0 0.0 0.0 22.00 : 0.0 0.0 0.0 0.0 0.0 0.0 23.00 : 0.0 0.0 0.0 0.0 0.0 0.0 24.00 : 0.0 0.0 0.0 0.0 0.0 0.0 (a) F i r s t t r i a l : i n i t i a l model: p = 0.5 for a l l cells, uniform weighting. Fig. 14. Inversion of the profile in Figure 11(b). 60. I N I T I A L MODEL X (KM) / Z = 0.50 1.50 9.00 : 1.000 0.500 10.00 : 1.000 1.000 11.00 -: 0.0 0.500 12.00 : 0.0 0.500 13.00 : 0.0 0.0 14.00 : 0.0 0. 500 15.00 : 0.0 0.500 16.00 : 0.0 0. 500 17.00 : 0.500 0.500 18.00 : 1.000 1.000 19.00 : 1.000 1.000 20.00 : 1.000 1.000 21.00 : 0.0 0.500 22.00 : : 0.0 0. 0 2.50 0.500 0. 500 0.500 0. 500 0.500 0.500 0.500 0.500 0.500 1.000 1 .000 0. 500 0.500 0.0 3.50 0.0 0.0 0.500 0. 500 0.500 0.500 0.500 0. 500 0.500 0. 500 0.500 0. 500 0.500 0.0 4.50 0.0 0.0 0.0 0.0 0.500 0. 500 0.500 0. 500 0.500 0.500 0.500 0. 500 0.0 0.0 5.50 0.0 0.0 0.0 0.0 0.0 0.0 0.500 0.500 0.500 0.0 0.0 0.0 0.0 0.0 APPROXIMATE DENSITIES X (KM) / Z= 0.50 1.50 2.50 3.50 4.50 5.50 9.00 . : 1.00 1.00 0.0 0.0 0.0 0.0 10.00 : 1.00 1.00 1.00 0.0 0.0 0.0 11.00 : 0.0 1.00 1.00 1.00 0.0 0.0 12.00 • : 0.0 1.00 1.00 1.00 0.0 0.0 13.00 : : 0.0 0.0 1.00 1 .00 1.00 0.0 14.00 • : 0.0 1.00 1.00 1.00 1.00 0.0 15.00 : 0.0 1.00 1.00 1 .00 1.00 1.00 16.00 : : 0.0 1.00 1.00 1.00 1.00 1.00 17.00 : 0.0 1.00 1 .00 1 .00 1.00 1.00 18.00 : : 1.00 1.00 1.00 1.00 1.00 0.0 19.00 : 1.00 1.00 1.00 1 .00 1.00 0.0 20.00 : : 1.00 1.00 1.00 1.00 1.00 0.0 21.00 : 0.0 1.00 1 .00 1.00 0.0 0.0 22.00 : : 0.0 0.0 1.00 1.00 0.0 0.0 (b) A subsequent inversion with an improved i n i t i a l model. WEIGHTED REGIONS X (KM) / Z= 0.50 1.50 2.50 3. 50 4.50 5.50 9.00 : 0.0 0.0 0.0 1.000 1.000 1.000 10.00 0.0 0.0 0.0 0.0 1.000 • 1.000 11.00 1.000 0.0 0.0 0.0 0.0 1.000 12.00 : 1.000 0.0 0.0 0.0 0.0 0.0 13.00 : 1.000 0.0 0.0 0.0 0.0 0.0 14.00 : • 1.000 0.0 0.0 0.0 0.0 0.0 15.00 0.0 0.0 0.0 0.0 0.0 0. 0 16.00 : 0.0 0.0 0.0 0.0 0.0 0.0 17.00 0.0 0.0 0.0 0.0 0. 0 0.0 18.00 : 0.0 0.0 0.0 0.0 0.0 0.0 19.00 0.0 0. 0 0.0 0.0 0.0 0. 0 20.00 • 0.0 0.0 0.0 0.0 0.0 0.0 21.00 1.000 0.0 0.0 0.0 0.0 0.0 22.00 : 1.000 1.000 0.0 0.0 0.0 0.0 THE WEIGHTING FACTOR IN CELLS MARKED 1. 0 IS 50.0 r HE WEIGHTING FACTOR IN CEL L S MARKED 2. 0 I S 1.0 THE WEIGHTING FACTOR IN CELLS MARKED 3. 0 IS 1.0 Or HER REGIONS HAVE UNIT WEIGHT I N THE MINIMIZATION (c) The weighting f a c t o r s used i n (b ) . 62 o o • o in l l l © I o t o CN «-«JJ O s © » o H • H O t> cn ** as O w o a © o • « o o o o m o o o o t in o © t o o o m o o t o CN O O • in CN O O © DISTANCE ALONG PROFILE IN KM. Lz = 5.0 F i g . 15. F i n a l model f o r the data of F i g u r e 11(b). 63 data of Figure l i b . In t h i s case, the f i n a l model was obtained i n four i t e r a t i o n s , and has a maximum e r r o r i n gr a v i t y of about 6%. These examples show that the weighted-distance method can quite r e a d i l y produce acceptable models from a r t i f i c i a l g r a v i t y data. The l i n e a r . formulation allows several inversions to be made without r e q u i r i n g excessive computation time; the e n t i r e sequence of operations to produce the model i n Figure 15 consumed about 40 seconds of CPU time on the IBM 360/67. The ultimate aim of any i n v e r s i o n method i s to produce models of the r e a l earth, and thus the f i n a l t e s t of the WEIGHT program was an i n v e r -sion of gr a v i t y data taken over the Guichon Creek b a t h o l i t h i n south c e n t r a l B r i t i s h Columbia. A d e t a i l e d g r a v i t y survey was taken i n 19 71; other data were used as an a i d i n three-dimensional modelling (using a polygon method); and cross-sections were constructed from the complete model. The d e t a i l s of t h i s work have been discussed by Ager (1972) and Ager e t a l (1972). The b a t h o l i t h i s a long, narrow e l l i p t i c a l body, and thus should be amenable to two-dimensional modelling. The data treated by the new method were obtained from the Bouguer anomalies on one p r o f i l e across the center of the b a t h o l i t h ; a constant value could be substracted to s a t i s f a c t o r i l y remove the regional f i e l d (C. Ager, personal communication, 1973). The r e s u l t i n g g r a v i t y p r o f i l e i s shown i n Figure 16. Surface g e o l o g i c a l i n v e s t i g a t i o n i n d i c a t e s that the density con-3 t r a s t of the b a t h o l i t h i s -0.15 gm/cm , and i t s surface outcrop extends f o r approximately 13 miles along the p r o f i l e . From t h i s information, the i n i t i a l 64 o o o «* I o s 35 H O O >< • o «*: CN = i o ss a H (/> w 03 O o • o ro I o o o I o o CO o o VO o o CN O o CN m DISTANCE ALONG PROFILE IN KM. X GZO 1. 6 - 2 . 1000 3 . 2 - 4 . 1000 a . 8 - 6 . 9000 6. 4 - 9 . 8000 8 . 0 - 1 2 . 7000 9 . 6 - 1 5 . 5000 1 1 . 2 - 1 8 . 7000 12.8 - 2 2 . 9000 14. 4 - 2 7 . 5000 16. 0 - 3 0 . 5000 17. 6 - 3 1 . 6000 19. 2 - 3 2 . 2000 2 0 . 8 - 3 2 . 7000 22. 4 - 3 2 . 0000 24. 0 - 2 9 . 4000 25. 6 - 2 5 . 3000 27 . 2 - 2 1 . 2000 28. 8 - 1 7 . 6000 30. 4 - 1 3 . 7000 32. 0 - 9 . 2000 3 3 . 6 - 5 . 0000 3 5 . 2 - 1 . 9000 These values were obtained by subtracting -116 mgal from the Bouguer anomalies supplied by C. Ager. (Note: h i s r e s i d u a l (Figure 19) was derived v i a a wavelength f i l t e r i n g method.) T i g . 16. Residual anomaly over the Guichon Creek b a t h o l i t h . 65 model was chosen as a surface layer of density 0.15 i n 13 c e l l s (since the s t a t i o n spacing was one mi l e ) . The dens i t i e s and the gra v i t y data were con-sidered to be p o s i t i v e , since the SMOOTH routine s e l e c t s only p o s i t i v e den-s i t i e s i n making approximate solutions ( i t could e a s i l y be adapted to treat e i t h e r p o s i t i v e or negative density; the in v e r s i o n i t s e l f does not depend on s i g n ) . To insure that the s o l u t i o n also had a surface density of 0.15, a weighting f a c t o r of 100.0 was applied to a l l the surface c e l l s . The approximate model obtained from the i n v e r s i o n i s shown i n Figure 17. There i s a reasonable f i t to the data; however the model i s composed of three unconnected blocks, but a s i n g l e anomalous body i s expected. Three furt h e r a p p l i c a t i o n s of the program developed the single-body model shown i n Figure 18. The shape i n d i c a t e d by the WEIGHT program inversions compares favourably with the cross-section obtained from a standard model-l i n g method (Figure 19). I t i s evident that the new method can be success-f u l l y applied to r e a l data. One must be somewhat c a r e f u l i n applying the method to r e a l data, f o r those data w i l l not be p e r f e c t l y accurate. In working with the batho-l i t h data, u n r e a l i s t i c inverses with large p o s i t i v e and negative d e n s i t i e s often r e s u l t e d from using models which di d not span the e n t i r e region beneath the p r o f i l e . The problem can be overcome by using a f u l l - r a n g e model ( i . e . the model's subsurface width equals the length of the p r o f i l e ) , since the program can then assign small d e n s i t i e s near any st a t i o n s with erroneous data; rather than f i t t i n g the data by u n r e a l i s t i c density configu-rations beneath the center of the p r o f i l e . APPROXIMATE DENSITIES 66 X (KM) / Z = 0. 80 2.40 1. 60 : 0.0 0.0 3.20 : 0.0 0.0 4. 80 : 0.0 0.0 6.40 : 0.0 0. 15 8.00 : : 0.0 0.30 9.60 : 0.15 0.0 11.20 . : 0. 15 0.0 12.80 : 0.15 0.0 14.40 , : 0. 15 0.15 16.00 : 0. 15 0. 15 17. 60 : 0. 15 0.15 19.20 : 0.15 0. 15 20.80 : 0. 15 0.15 22.40 : : 0.15 0. 15 24.00 : : 0. 15 0.15 25.60 : 0.15 0. 15 27. 20 : 0. 15 0.0 28.80 « : 0. 15 0.0 30.40 ; 0.0 0.30 32.00 : 0.0 0. 15 33.60 : 0. 0 0.0 35.20 : 0.0 0.0 4.00 0.0 0. 0 0.0 0. 0 0.0 0. 0 0.0 0. 0 0. 15 0. 15 0.15 0. 15 0. 15 0. 15 0.15 0. 0 0.0 0.0 0.15 0.0 0.0 0.0 5.60 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.15 0. 15 0.15 0. 15 0.15 0. 15 0.15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.20 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0. 15 0.15 0. 15 0.15 0. 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8.80 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 EXACT GRAVITY - APPROX GRAVITY - ERROR 1 2.1000 3. 2526 -1. 1526 2 4. 1000 4. 4243 - 0 . 3243 3 6.9000 6. 5287 0. 3713 4 9.8000 9. 7191 0. 0809 5 12.7000 12. 6765 0. 0235 6 15. 5000 17. 7263 - 2 . 2263 7 18.7000 19. 5345 -0. 8345 8 22.9000 22. 5727 0. 3273 9 27.5000 26. 7013 0. 7987 10 30.5000 29. 9065 0. 5935 11 31.6000 31. 7112 -0. 1112 12 32.2000 32. 3741 - 0 . 1741 13 32.7000 32. 0374 0. 6626 14 32.0000 30. 7138 1. 2862 15 29.4000 28. 3698 1. 0302 16 25.3000 25. 0967 0. 2033 17 21.2000 21. 6822 -0. 4822 18 17.6000 19. 5061 - 1 . 9061 19 13.7000 14. 3055 -0. 6055 20 9.2000 11. 0188 - 1 . 8188 21 5.0000 7. 4272 -2. 4272 22 1.9000 5. 0234 - 3 . 1234 Fig. 17. The f i r s t inversion of the batholith anomaly. \ \ o o I 1-3 < O U O sc • CN X I M > 03 O a o • O I o s / /* o o o i o o CO o o VO o o cr CN O O CN CO DISTANCE ALONG PROFILE IN KH. The dashed l i n e i n d i c a t e s the r e s i d u a l anomaly of Figure Fig. 18. A model for the Guichon Creek batholith. On o cn E > o C P -icH -20- A g ( regional value derived by filtering the complete Bouguer anomaly ) — F i g . 19. A b a t h o l i t h cross-section obtained by standard modelling methods (from Ager et a l , 19 72). 69 In conclusion, i t i s suggested that the models produced by a "weighted-distance" minimization are a reasonable means of i n t e r p r e t i n g gravity data. A r t i f i c i a l data f o r two simple bodies, as w e l l as a r e a l g r a v i t y p r o f i l e , have been inverted by successive applications of the WEIGHT program; i n each case a simple model of one density was obtained. The process i s not completely o b j e c t i v e , since improvement of i n i t i a l models i s at the d i s c r e t i o n of the i n t e r p r e t e r . However, the l i n e a r functionals allow rapid computation of each inverse, so the method can be applied s e v e r a l times at a modest cost i n computer time. 70 CHAPTER 5 Numerical Optimization Of Gravity Models Using Entropy 5.1 Entropy of Gravity Models It was noted earlier that gravity interpretation often has the objective of defining the shape of a body of anomalous density. If the i n -version technique uses a combination of simple bodies, a method to find a compact grouping of those units i s desirable. To produce such single-density models from Backus-Gilbert inversions, i t i s usually necessary to apply considerable interpretive judgement. In hopes of developing a more objective method for finding compact models, we can investigate c r i t e r i a for selecting an acceptable model of the desired form. In particular, the con-cept of an entropy assigned to a density model w i l l be applied as a measure of "order" or "structure" of the model, following a suggestion of G.K.C. Clarke (personal communication, 1971). Entropy concepts are used i n many areas of physics and mathematics, and the success of entropy methods in other fields led to the present invest-igation. In s t a t i s t i c a l mechanics, the entropy of a closed system i s H = - E p s log p s where p g is the probability that the system is in state s. In general, maxi-mum entropy corresponds to the equilibrium state ( K i t t e l , 1958). Information theory uses entropy as a measure of average information (or average uncer-tainty) of a message source (Lathi, 1968). In spectral analysis, the maxi-mum entropy method developed by Burg (1967) has been quite successful, par-ticularly in resolving frequency components from short records (Lacoss, 1971; Ulrych, 1972) 71 Clarke's suggestion f o r using entropy of a density model d i f f e r s from the usual a p p l i c a t i o n s , since minimum entropy w i l l be sought as a c r i t e r i o n of maximum order of the model. In th i s sense, the idea i s another method f o r s e l e c t i n g a p a r t i c u l a r model from a l l those which s a t i s f y the gr a v i t y data. Minimum entropy i s an apparent c o n f l i c t with the p r i n c i p l e s of s t a t i s t i c a l mechanics and thermodynamics, since the e q u i l i b r i u m state of any system should have maximum entropy. Nevertheless, minimum entropy can be a us e f u l property, i f i t can produce the most structured model compatible with the data, since structures or concentrated bodies are the usual target of geophysical exploration. In ad d i t i o n , the maximum entropy p r i n c i p l e i n p h y s i c a l systems i s not i n v i o l a b l e . Fast (1968) observed that i n some s u b s t i t u t i o n a l a l l o y s , the tendencies toward entropy maximization and energy minimization c o n f l i c t . Prigogine et a l (1972) examined development of b i o l o g i c a l ordered systems, and concluded that f l u c t u a t i o n s i n i r r e v e r s i b l e systems could produce low entropy configurations, which need not be the lowest p r o b a b i l i t y state of the system. In l i g h t of these examples, i t i s reasonable to look f o r minimum entropy models, without being concerned that they may be p h y s i c a l l y improbable i f the c r u s t a l density d i s t r i b u t i o n i s a random process. The knowledge that compact ore bodies do e x i s t i n the earth's crust i n i t s e l f supports the argument. By analogy with other a p p l i c a t i o n s , the entropy of a density model i s defined by 72 H - - |v£mgdv ( i ) The normalized v a r i a b l e £ i s used to give the appearance of a p r o b a b i l i t y M density function, which requires M dV = 1 0 <s< 1 M dV i s monotonically in c r e a s i n g For the l i n e a r density models, (1) becomes P P u v v urn - rnm H = - I E —— In n m M M (2) where M = Z Z p ( c e l l dimensions cancel i n (2)) nm n m To f i n d the desired compact or structured model, (2) must be minimized under the constraints imposed by the known gravity anomaly. In the present formulation,.these are Z Z G. p = g. m X (3) inm nm j j n m J In addition, bounds can be imposed on density to produce g e o l o g i c a l l y reason-73 able models. We also require p to have the same sign everywhere, i f ^  i s to have the properties of probability. The additional constraints are now expressed as 0 p m $ A (4) nm (or opposite signs for a negative anomaly). To find a method to solve this problem, we f i r s t examine the properties of entropy in density models. A decrease in entropy as models become more compact i s confirmed by considering horizontal cylinders of different density but equal mass. Such cylinders w i l l produce the same gravity pro f i l e i f centered at the same depth. To compute the entropy of these simple models, (1) is simplified by considering density to be constant inside the cylinder. Since the cylinder mass M i s constant ^ = 77 i s M V constant for a given cylinder (V i s the volume). Then E - - j (y) In c|) dV = In V = - In £ (5) Since 11m x ln(X) = 0, x-*0 integration i s done only over the cylinder. For the two-dimensional models, V i s actually the cross-sectional area. The variation of entropy with the density and radius of the cylinder i s shown in Figure 20; we consider c y l i n -ders of different density as acceptable models for the same gravity profile. The entropy minimum w i l l correspond to a point mass (radius -*• 0) , and the maximum to a cylinder of i n f i n i t e radius (p -> 0) . 74 — MODEL ENTROPY — 15.00 17.50 20.00 22.50 25.00 27.50 * D I * E * N . * S 1.00. * I . * T * X - * * 0 2.00. * F . * * C * Y * L 3.00. * 1 * N • * D * E . * R 4.00. * — MODEL ENTROPY — 15.00 17.50 20.00 22.50 25.00 27.50 • * a . * k * D 1.00. * I . * a * s * 0 2.00. * F * m * C . * 1 • . * L 3.00. * I . * M • * D * E • * B 4.00. * Fig. 20. Entropy of horizontal cylinders of equal mass. 75 Figure 20 demonstrates that entropy has an absolute minimum only i n the l i m i t i n g sense of the model becoming a point mass. We require bounds on density (as i n equation ( 4 ) ) , and thus entropy minimization must be done i n a bounded density i n t e r v a l . A major problem i s that the minimum i s not 8H a stationary point ( i . e . we do not have -r— = 0) and thus the c l a s s i c a l mini-dp mization procedures cannot be used. The upper bound on density, A, should be chosen as the desired density of the f i n a l model, since (5) implies that the minimum entropy model w i l l have the l e a s t p o s s i b l e volume. Minimizing entropy of the l i n e a r models should then assign a density of zero or A to each prism; hence a constant density model would be produced. Another consequence of equation (5) i s that model entropy i s a u s e f u l property only i f density i s a v a r i a b l e parameter. Any model which assumes constant density also has constant entropy, since the mass i s a known constant; i . e . by s p e c i f y i n g the density of a model, one i s also s p e c i f y i n g i t s entropy. 5.2 Numerical Optimizations of Model Entropy Minimizing entropy under the constraints (3) and (4) i s a d i f -f i c u l t non-linear optimization problem, which cannot be solved by the tech-niques of v a r i a t i o n a l c alculus, since the minimum i s not a s t a t i o n a r y point. As a r e s u l t , a somewhat d i f f e r e n t viewpoint of the model i s developed; we consider the d e n s i t i e s of N prisms as N independent parameters, rather than viewing density as a s i n g l e parameter which varies with the s p a t i a l coor-dinates. The problem i s then an optimization i n an N-dimensional space. There are many methods f o r optimization i n general parameter spaces, but 76 most are designed f o r s p e c i a l i z e d problems, so that a p a r t i c u l a r example i s often t r a c t a b l e with only one or two procedures. Beveridge and Schechter (1970) have analyzed many numerical methods and discussed t h e i r f e a s i b i l i t y f o r various problems. On t h e i r recommendation, and that of the UBC Computing Center (C. B i r d , personal communication, 1973), the COMPLEX search method of Box (1965; Box et a l 1969) was adapted to minimize entropy. The major advantage of the COMPLEX method f o r the entropy problem i s that i t makes no use of de r i v a t i v e s or gradients. Derivatives of model entropy cannot be used, since they become i n f i n i t e f o r those c e l l s whose density approaches zero, and the desired compact models w i l l have prisms of zero density. E = " Z M 1 1 1 M = " 1 P' 1 0 P' Jjjrr- CP' In p') - k p ' + l (6) and l i m In p' = - °° p '-*<>; This d i f f i c u l t y l e d to the r e j e c t i o n of most popular optimization methods. A b r i e f d e s c r i p t i o n of the general ideas behind these techniques i s given i n Appendix D. The COMPLEX method follows an i n i t i a l strategy somewhat s i m i l a r to Monte Carlo modelling. S t a r t i n g with one acceptable model (considered to be a point i n the N-dimensional parameter space), a set of "po i n t s " with 77 random coordinates are generated. The random points can be made to define new acceptable models (pr " f e a s i b l e points") by moving them towards the centroid of f e a s i b l e points u n t i l a l l constraints are s a t i s f i e d . The f e a s i b l e points define the v e r t i c e s of a fi g u r e termed a "complex". Searching f o r the optimum of an objective function consists of evaluating the function at each vertex, r e j e c t i n g the vertex of "worst" value, and re p l a c i n g i t by r e f l e c t i o n through the centroid of remaining points. At each step, tests are required to make sure each new point i s f e a s i b l e ; i n ad d i t i o n , a d d i t i o n a l s t r a t e g i e s can be employed to keep the search moving towards the optimum point. The d e t a i l s of the method, along with modifications developed f o r the present use, are described i n Appendix E. i s quite straightforward. The coordinates of the parameter space are the den s i t i e s of the N prisms which comprise the model. A l l d e n s i t i e s are assumed to l i e i n the i n t e r v a l [0,AJ, where A i s the expected constant density of an anomalous body. The i m p l i c i t constraints are a known gr a v i t y anomaly, Xy which must be s a t i s f i e d w i t h i n a s p e c i f i e d l i m i t <5; thus the constraints are w r i t t e n The a p p l i c a t i o n of a COMPLEX search to the minimum entropy problem I A, - E E G. p 1^6 1 j inm nm 1 J n m J or -<5 < ( A . - E E G . jnm nm ) $ 6 (7) along with the bounds 0 •$ P nm S A (8) 78 Equations (7) and (8) are linear in Pnm> and thus the feasible region is convex, since i t can be shown that any region whose boundaries are defined by linear constraints is convex (Beveridge and Schechter, 1970, p. 113; Cooper and Steinberg, 1970, p. 68). It i s convenient to write the necessary computer programs in several parts; the programs used here were as follows. The START program computes a l l necessary arrays (such as the Frechet kernels) and passes them, along with an acceptable model and the gravity data, to the searching routines. POINTS generates the i n i t i a l complex from the starting model. The TEST subroutine ensures that each point satisfies the implicit constraints (7). SEARCH does the rest of the job, computing entropies, finding "worst" points, and generating new ones. The optional routine RESTAR (very similar to POINTS) w i l l create a new complex from the current best point i f : (a) the centroid appears to be unfeasible; (b) the search is stalled at one vertex; or (c) a specified number of iterations have been completed. Several subtle d i f f i c u l t i e s can arise i n obtaining reasonable results from these programs. Since the search region is convex, the centroid of feasible points must always l i e within i t ; however the f i n i t e precision of the computer may cause an apparently unfeasible centroid when the search i s near a boundary. This problem can usually be overcome by generating a new complex from the current best point. To maintain a reasonable rate of progress, a f a i r l y generous error should be permitted in applying the con-straints; a general rule of thumb i s to set <S at 5% to 10% of the maximum gravity. The limits of model gravity allow acceptable models to have a range 79 of masses; frquently the search w i l l i n i t i a l l y only reduce the mass before concentrating i t i n as few cells as possible. In these cases, the f i n a l model tends to have a smaller mass than is indicated by the data; the den-si t i e s could be multiplied by an appropriate scale factor to restore the true mass. The f i r s t test of the method was a simple example to show that the minimum entropy model w i l l i n fact be the most compact. A simple 9-prism model was used, with gravity constraints generated by a density of 0.9 in the center c e l l . The search started from a model with mass 0.5 in the central block, and Q.4 distributed around i t , as shown i n Figure 21. As expected, entropy minimization concentrated a l l the mass i n the center block. A more d i f f i c u l t problem i s to start from a very diffuse i n i t i a l model; in the previous example, the center c e l l was clearly dominant in the i n i t i a l model. This time a 12-parameter model was used, with gravity data generated by a unit-density c e l l i n the third block of the middle column. The i n i t i a l model had density 0.11 in the nine cells grouped around the expected source of the anomaly. This search did not successfully converge, as the cells whose mass i n i t i a l l y increased were not in the central location. The model obtained i n 600 iterations i s shown in Figure 22. It seems li k e l y that a very diffuse i n i t i a l model may be brought to an intermediate stage of compactness which cannot be further improved. To test the method on a possible exploration objective, the gravity p r o f i l e of the body shown in Figure 23 was used. A starting point was ob-80 INITIAL MODEL Z (KM)/ X= 14.00 15.00 16.00 1.50 2.50 3.50 0.050 0.050 0. 050 0.050 0.500 0.050 0.050 0.050 0. 050 TOTAL MASS IS MODEL - ENTROPY 0.9000 1 1.61116 15 FEASIBLE POINTS WILL BE PICKED 400 IS THE MAXIMUM NUMBER OF ITERATIONS ALLOWED A NEW COMPLEX WILL BE GENERATED EVERY 125 ITERATIONS NO RESTARTING OF THE COMPLEX AFTER 5 ITERATIONS THE GRAVITY WILL AGREE WITHIN 0.400 MGAL UPPER DENSITY LIMIT IS 1.00 GM/CC CONVERGENCE IF CENTROID ENTROPY AGREES WITHIN 0.10E-THE RANDOM NUMBER GENERATOR IS INITIALIZED AT 0.0 THE EXPANSION PARAMETER OF THE COMPLEX IS 1.40 CONVERGENCE WITHIN SPECIFIED LIMIT; FEASIBLE POINT NO. 9 Z (KM) / X= 14.00 15.00 1.50 : 2.50 : 3.50 : 0.002 0.000 0.001 0.001 0.951 0.001 16.00 0.002 0.001 0. 000 CENTROID ENTR0PY= 0.0001397 CENTROID MASS= 0.9590 Fig. 21. A simple entropy minimization. 0.0 I T E B A T I 0 N C 0 0 N T 81.00. 181.00. 281.00. 381.00. 481.00. 581.00. ENTROPY OF MODEL 0.50 1.00 1.50 2.00 2.50 * * * * INITIAL MODEL (KM)/ X= 14.00 15.00 16.00 0.50 : 0.0 0.0 0.0 1.50 : 0.110 0. 110 0.110 2.50 : 0.110 0. 110 0. 110 3.50 : 0.110 0. 1 10 0. 110 THE BEST POINT IS SHOWN BELOW FEASIBLE POINT NO. 2 (KM) / X= 14.00 15.00 16.00 0.50 • « 0.006 0.000 0.000 1.50 • 0.000 0.406 0.012 2.50 • * 0.1 64 0.006 0.000 3.50 • 0.000 0.350 0.001 Fig. 22. Minimum entropy from a diffuse starting model. co 82 o 1 X o 1 a . • o o 1 a vo 1 1 1 «c 1 H 1 u \ H O 1 CM O 1 H . • H O 1 1 1 1 1 O 1 X o 1 EH O 1 M • 1 > O 1 <C CN 1 05 1 o 1 o o t * # * * * * * * * * * * * * o o • vn o o o o o o o CN O O in CN * # * * O O o ro DISTANCE ALONG PROFILE IN KM. I z = 5.0 Computer Display Z (KM)/ X= 13.00 14.03 15.00 16.00 0. 50 1. 50 2.50 0.0 0.0 -0. 0 1.000 1.000 1.000 1.000 1.000 0.0 0.0 0.0 1.000 Fig. 23. An a r t i f i c i a l gravity prof i l e for entropy minimization. INITIAL MODEL Z (KM) / X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 0.50 • 0.0 0. 100 0.200 0.200 1.000 0. 100 0.0 1.50 • • 0.100 0.400 0.400 0.500 0.800 0. 100 0.200 2.50 * • 0.0 0.400 0.700 0.700 0.300 0.0 0.0 TOTAL MASS IS 6.2000 MODEL - ENTROPY 1 2. 52968 THE BEST POINT IS SHOWN BELOW • FEASIBLE POINT NO. 18 Z (KM) / X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 0.50 • 0.000 0.002 0.000 0.000 0.994 0. 002 0.000 1.50 • 0.000 0.000 0.675 0.766 0.988 0.000 0.001 2.50 • 0.002 0.586 0.998 0.942 0.002 0.000 0.017 CENTROID ENTROPY 1. 9187727 CENTROID HASS= 5. 9904 Model gravity accurate within 3.0 mgal. Fig. 24. A minimum entropy model from the ar t i f i c i a l profile. 84 tained by using a "zero" i n i t i a l model in the WEIGHT program of Chapter 4. The results of this example are shown in Figure 24; the allowed error has permitted a solution different from the "real" body. The f i n a l solution also seems to reflect the starting model, since the largest densities of that model were a l l increased. Nonetheless, the method i s pa r t i a l l y . f u l f i l l i n g our basic purpose; to provide an objective means for producing a compact model from the diffuse inversions produced by the methods of Chapter 4. One advantage of the COMPLEX method i s that i t can easily opti-mize different objective functions. A simple function which may also produce a f i n a l model of density A is F(p) - - E E G> - A) (9) nm nm n m Each term in the summation is zero at p=0 and p=A, hence minimizing F(p) should set each p^ to a value of 0 or A. The similarity to entropy is ill u s t r a t e d i n Figure 25. This is to be expected, since the f i r s t term in a Taylor series for ln(x) is (x-1) (Abramowitz and Stegun, 1965, p. 68), and thus the f i r s t term in (x) ln(x) i s x(x-l). The main difference i s that FCp) is not normalized, which can be an advantage i n the minimization, particularly i f the mass of the model is greater'than about 3.0 (here mass refers to the sum of a l l prism densities, as in Equation (2)). In such cases, the contribution to entropy from each c e l l i s restricted to the interval x<0.3 in Figure 25; and an increase in density for any c e l l w i l l not improve the objective function. In minimizing the unnormalized function (9), hringing an individual density towards A w i l l s t i l l improve the value p R I S M D E N S I T X 0.20 0.40 0.60 0.80 0. 10 ENTROPY = X*LN (X) 0.20 0.30 ,0.40 F=X (X-A) -- OBJECTIVE FUNCTION 0.0 . 0.10 0.20 0.30 1.00* F i g . 25. The f u n c t i o n s x l n ( x ) and x(x-A). oo INITIAL MODEL 2 (KM)/ X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 ................. ... ... ...... ^  . 0.50 : 0.0 0. 100 0.200 0.200 1.000 0. 100 0.0 1.50 : 0.100 0.400 0.400 0.500 0.800 0.100 0.200 2.50 : 0.0 ,0.400 0.700 0.700 0.300 0.0 0.0 TOTAL MASS IS 6.2000 OBJECTIVE FUNCTION 2.59999 THE BEST POINT IS SHOWN BELOW Z (KH) / X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 0.50 : 0.0 0.0 0.001 0.0 0.989 0.0 0.0 1.50 : , 0.004 0.997 0.0 0.996 0.996 0.0 0.0 2.50 : 0.0 0.845 0.998 0.138 0.001 0.029 0.001 MASS IS 5.9950 CENTROID FUNCTION= 0.3207084 Model gravity accurate within 3.0 mgal Fig. 26. An optimum model for Figure 23, using x(x-A) . 87 of the objective function. The new function can be optimized with the same basic programs, except that the SEARCH routine no longer computes entropies. The example of Figure 24 was repeated, with results shown in Figure 26. In the same number of iterations, F(p) has been more successful than entropy i n putting a l l the mass in cells of density A, which is probably a result of entropy being normalized by a mass M=6.0. Once again the error limit has prevented exact duplication of the body in Figure 23. This example suggests that F(p) is a reasonable alternative to entropy, particularly for multi-block models; however the nature of the search w i l l l i k e l y prohibit any dramatic saving of computer time. The optimization programs can also be easily adapted to find maximum entropy models, which may be of interest i f maximum entropy does indicate the most probable state (although i t may be suggested that the most probable density configurations preclude the existence of gravity anomalies). By analogy to the arguments suggesting minimum entropy as a measure of compactness, maximum entropy should produce the most diffuse model compatible with the observations. Maximum entropy models were computed for two of the earl i e r examples for direct comparison to minimum entropy. Figure 27 is the result of 400 iterations of the simple 9-parameter model; the maximization i s clearly trying to evenly distribute the mass. Figure 28 is the comparison for the gravity profile of Figure 23. Maximizing entropy has not made any 88 I N I T I A L MODEL Z (KM)/ X= 14.00 15.00 16.00 1. 50 2.50 3.50 0.050 0.050 0.050 0.050 0.500 0.050 0.050 0.050 0.050 TOTAL MASS IS 0.9000 MODEL - ENTROPY 1 1.61116 THE BEST POINT IS SHOWN BELOW FEAS I B L E POINT NO. 13 Z (KM)/ X= 14.00 15.00 1.50 : 2.50 : 3.50 : 16.00 0.100 0.040 0.077 0.029 0.228 0.091 0.143 0.232 0.103 CENTROID ENTRCPY= 2.1702032 CENTROID MASS= 1.0422 Model gravity accurate within 0.4 mgal. Fig. 27. A maximum entropy model corresponding to Figure 21. Z (KM)/ X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 0.50 : 0.000 1.50 : 0.000 2.50 : 0.002 (a) Minimum entropy. 0.002 0.000 0.000 ' 0.675 0.586 0.998 0.000 0.994 0.766 0.988 0.942 0.002 0.002 0.000 0.000 0.001 0.000 0.017 Z (KM)/ X= 12.00 13.00 14.00 15.00 16.00 17.00 18.00 0.50 : 0.001 0.128 0.217 0.261 0.946 0.135 0.017 1.50 : 0.187 0.473 0.401 0.591 0.561 0.149 0.239 2.50 : 0.281 0.493 0.674 0.532 0.258 0.125 0.109 (b) Maximum entropy. F i g . 28. A comparison of maximum and minimum entropy models f o r F i g u r e 23. co vo 90 N K ITMAX CPU TIME COMMENTS _ _ (sec.)  9 15 400 34 converged 12 15 600 58 partly converged 18 21 400 61 partly converged 21 25 500 82 partly.converged 60 62 1000 501 no progress N = Number of parameters K = Number of points in the complex ITMAX = Number of iterations TABLE IV. Computation time for numerical optimizations. dramatic change in the model, while minimum entropy has had some success i n concentrating the mass. Our expectations regarding maximum entropy models seem to be upheld. The numerical optimizations were not successful in developing any really complex models, largely because they require very extensive computa-tions. Table IV shows the CPU time consumed for various t r i a l s of the mini-mum entropy program. The 60-parameter example was an attempt to improve a simple model for the a r t i f i c i a l p rofile of Figure 11a, Chapter 4 (the i n i t i a l model was obtained using the WEIGHT program); there was no discernible improvement towards compactness in 1000 iterations. We conclude that the numerical method i s not practical for any excessive number of parameters (i.e. prisms In the model). 91 The r e s u l t s of numerical optimizations support the suggestion that minimizing entropy w i l l produce a compact model, i . e . a model of maxi-mum allowable density. Simpler objective functions can also be adopted towards the same goal. Unfortunately, the straightforward a p p l i c a t i o n of these c r i t e r i a i s very time-consuming. I t appears that these optimum models, as presently derived, are not p r a c t i c a l f o r any complex system r e q u i r i n g many parameters. However, t h i s does not n e c e s s a r i l y mean that a p r a c t i c a l minimum entropy method cannot be developed; since there may be many other approaches to the problem, f o r example transformation i n t o a space where entropy i s a more simply-behaved function. 92 CHAPTER 6 Conclusion Geophysicists are continually seeking new ways to obtain accurate data and extract information from them. In many applications, modelling plays an integral role i n interpretation. Gravity exploration data lend themselves to modelling techniques, and i t i s hoped that some of the new developments for gravity problems presented here can be useful in exploration situations. The linear models proposed in Chapter 3 have proven quite useful for inverting gravity profiles via a Backus-Gilbert approach (developed in Chapter 4). I n i t i a l experience with the "weighted-distance" method indicates that i t i s flexible, particularly i n iterative use, when certain densities from the previous solution can be held constant. Approximate models (composed of prisms of a single density) can frequently give an adequate f i t to a given profile after 4 or 5 repetitions, at quite modest expense in terms of computer time. Reasonable care must be taken with real data to ensure an acceptable solution. It is perhaps best to use models spanning the subsurface region beneath the pr o f i l e , so that the inversion can more easily account for noise in the data. If a model has too few prisms to adequately represent the data, an oscillatory solution (i.e. large positive and negative densities) may result. There are many po s s i b i l i t i e s for further improvement of the methods of Chapter 4. The present work certainly j u s t i f i e s more application to real 93 data. In some situations, i t appears that the requirement of exact agreement with the data i s detrimental; the inversion produces an exceptionally close f i t by using unrealistic densities. In these examples, i t may be helpful to relax the constraints, by demanding agreement only within a specified error. Gilbert (19 71) and Wiggins (1972) have discussed ways in which this might be achieved. The large-block models may also prove useful for real data, i f the problem of bad solutions for high-frequency data can be overcome. It may be possible to find acceptable inverses in the spatial frequency domain, and f i t the model only to the lower frequency components of the data. The large blocks could then be subdivided for further improvement using the complete data set. The investigation of entropy in Chapter 5 demonstrates that optimi-. zing a compactness property of a density model i s a feasible idea, although not necessarily a practical one. Other approaches to the entropy problem may well be more ef f i c i e n t . Numerical optimization was of some help i n solving a d i f f i c u l t problem, and might conceivably be useful with many other types of models. The COMPLEX method described i n Chapter 5 might be easily adapted to other problems i f three simple changes are made: (1) allow a variable number of constraints (the present program requires 30); (2) incorporate objective function evaluations as a separate routine; and (3) add an i n d i -cator to choose maximization or minimization. It appears that the COMPLEX method w i l l work best with a small number of parameters, so the large-block models might be of some benefit here. 94 A l l the i n v e s t i g a t i o n s herein were pursued i n hopes of gaining some b a s i c i n s i g h t i n t o the nature of the modelling process. Perhaps the greatest b e n e f i t from these studies i s the knowledge that new methods from other areas of geophysics and applied mathematics can e a s i l y be adapted to the problems of exploration geophysics. 95 REFERENCES Abramowitz, M., and I. Stegun, 1965, Handbook of Mathematical Functions, Dover, New York. Agarwal, B.N.P., 1971, Direct gravity interpretation of sedimentary basin using d i g i t a l computer, Part I and II, Pure Appl. Geophys. v. 86, pp. 5-17. Ager, C.A., 1972, A gravity model for the Guichon Creek batholith, unpub-lished M.Sc. thesis, University of B.C. Ager, C.A., W.J. MacMillan, and T.J. Ulrych, 1972, Gravity, magnetics and geology of the Guichon Creek batholith, Bulletin 62, B.C. Dep't of Mines and Petroleum Resources. Al-Chalabi, M., 1971, Some studies relating to nonuniqueness i n gravity and magnetic inverse problems, Geophysics, v. 36, pp. 835-855. , 19 72, Interpretation of gravity anomalies by non-linear optimi-sation, Geophys. Prosp., v. 20, pp. 1-16. Anderssen, R.S., 1970, The character of non-uniqueness in the conductivity modelling problem for the earth, Pure Appl. Geophys., v. 80, pp. 238-259. Anderssen, R.S., M.H. Worthington, and J.R. Cleary, 1972, Density modelling by Monte Carlo inversion, I - Methodology, II - Comparison of recent earth models, Geophys. J., v. 29, pp. 433-444; 445-457. Backus, G.E., and J.F. Gilbert, 1967, Numerical application of a formalism for geophysical inverse problems, Geophys. J., v. 13, pp. 247-276. , 1968, The resolving power of gross earth data, Geophys. J., v. 16, pp. 169-205. , 1970, Uniqueness in the inversion of inaccurate gross earth data, Ph i l . Trans. Roy. Soc. London, A, v. 266, pp. 123-192. Berezhnaya, L.T., and M.A. Telepin, 1968, Determination of density from gravimetric data, Exploration Geophys., v. 47, pp. 107-113. Beveridge, G.S.G., and R.S. Schechter, 1970, Optimization: Theory and Practice, 773 pp., McGraw-Hill, New York. Box, M.J., 1965, A new method of constrained optimization and a comparison with other methods, Computer J., v. 8, pp. 42-52. , 1966, A comparison of several current optimization methods, and the use of trans-formations in constrained problems, Computer J., v. 9 pp. 67-77. 96 Braile, L.W., G.R. Keller, and W.J. Peebles, 1973, Inversion of gravity profiles for two-dimensional density distributions, paper presented at Utah Inversion Symposium, Salt Lake City, March 13, 1973. Burg, J.P., 1967, Maximum entropy spectral analysis, presented at 37th Annual Meeting, SEG, Oklahoma City, October 31, 1967. Cooley, W.W. and P.R. Lohnes, 1971, Multivariate Data Analysis, 364 pp., Wiley, New York. Cooper, L., and D. Steinberg, 1970, Introduction to Methods of Optimization 381 pp., W.B. Saunders, Philadelphia. Corbato, C.E., 1965, A least-squares procedure for gravity interpretation, Geophysics, v. 30, pp. 228-233. Cordell, L., and R.G. Henderson, 1968, Iterative three-dimensional solution of gravity anomaly data using a d i g i t a l computer, Geophysics, v. 33, pp. 596-601. Der, Z.A., and M.Landisraan, 1972, Theory for errors, resolution, and sepa-ration of unknown variables i n inverse problems, with application to the mantle and crust i n Southern Africa and Scandinavia, Geophys. J., v. 27, pp. 137-178. Dunford, N., and J.T. Schwartz, 1958, Linear Operators Part I, Pure and Applied Math., v. 7, 858 pp., Interscience, New York. Dyrelius, D., and A. Vogel, 1972, Improvement of convergence in iterative gravity interpretation, Geophys. J., v. 27, pp. 195-205. Fast, J.D., 1968, Entropy, 332 pp., Philips Technical Library, Eindhoven, Holland. Garland, G.D., 1965, The Earth's Shape and Gravity, 183 pp., Pergamon, London. Gilbert, F., 1971, Ranking and winnowing gross earth data for inversion and resolution, Geophys. J., v. 23, pp. 125-128. Grant, F.S., and G.F. West, 1965, Interpretation Theory in Applied Geophysi 584 pp., McGraw-Hill, New York. Guin, J.A., 1968, Modification of the Complex method of constrained optimi-zation, Computer J., v. 10, pp. 416-417. Hoffman, K., and R. Kunze, 1961, Linear Algebra, 332 pp., Prentice-Hall, Englewood C l i f f s , N.J.., Jackson, D.D., 1972, Interpretation of inaccurate, insufficient, and incon-sistent data, Geophys. J., v. 28, pp. 97-109. 97 Jackson, D.D., 1973, The EDGEHOG procedure, paper presented at Utah Inversion Symposium, Salt Lake City, March 12, 19 73. Jordan, T.H., and J.N. Franklin, 1971, Optimal solutions to a linear inverse problem in geophysics, Proc. Natl. Acad. Sci. U.S., v. 68, pp. 291-293. Keilis-Borok, V.I., and T.B. Yanovskaya, 1967, Inverse problems of seismology (structural review), Geophys. J., v. 13, pp. 223-234. K i t t e l , C , 1958, Elementary S t a t i s t i c a l Physics, 228 pp., Wiley, New York. Lacoss, R.T., 1971, Data adaptive spectral analysis methods, Geophysics, v. 36, pp. 661-675. Lathi, B.P., 1968, An Introduction to Random Signals and Communication Theory, 488 pp., International Textbook Co., Scranton, Pa. MacMillan, W.D., 1930, Theoretical Mechanics: The Theory of the Potential, 469 pp., McGraw-Hill, New York, (reprinted 1958, Dover, New York). Mottl, J., and L. Mottlova, 1972, Solution of the inverse gravimetric problem with the aid of integer linear programming, Geoexploration, v. 10, pp. 53-62. Nagy, D., 1966, The gravitational attraction of a right rectangular prism, Geophysics, v. 31, pp. 362-371. Negi, J.G., and S.C. Garde, 1969, Symmetric matrix method for rapid gravity interpretation, J. Geophys. Res., v. 74, pp. 3804-3807. Odegard, M.E., and J.W. Berg, 1965, Gravity interpretation using the Fourier integral, Geophysics, v. 30, pp. 424-438. O'Reilly, D., 1971, Introduction to non-linear function optimization, Infor-mal Documentation, Computing Centre, University of B.C. Parasnis, D.S., 1962, Principles of Applied Geophysics, 176 pp., Meuthen, London. Parker, R.L., 1970, The inverse problem of e l e c t r i c a l conductivity i n the mantle, Geophys. J., v. 22, pp. 121-138. , 1972, Inverse theory with grossly inadequate data, Geophys. J., v. 29, pp. 123-138. , 1973, A new method of modelling gravity and magnetic anomalies;, presented at Utah Inversion Symposium, Salt Lake City, March 13, 1973. Penrose, R., 1955, A generalized inverse for matrices, Proc. Cambridge Phil. Soc, v. 51, pp. 406-413. Pierre, D..A., 1969, Optimization Theory with. Applications, 612 pp., Wiley, New York. 98 Press, F., 1968, Earth models obtained by Monte Carlo inversion, J. Geophys. Res., v. 73, pp. 5223-5234. , 1970, Earth models consistent with geophysical data, Phys. Earth Planet. Interiors, v. 3, pp. 3-22. Prigogine, I., G. Nicolis, and A. Babloyantz, 1972, Thermodynamics of evo-lution, Parts I and II, Physics Today, v. 25, no. 11, pp. 23-28; v. 25, no. 12, pp. 38-44. Rastrigin, L.A., 1965, Solution of inverse problems by s t a t i s t i c a l optimi-zation methods, Rev. Geophys., v. 3, pp. 111-114. Roy, A., 1962, Ambiguity in geophysical interpretation, Geophysics, v. 27, pp. 90-99. Skeels, D.C, 1947, Ambiguity in gravity interpretation, Geophysics, v. 12, pp. 43-56. , 1967, What is residual gravity?, Geophysics, v. 32, pp. 872-876. Smith, M.L., and J.N. Franklin, 1969, Geophysical applications of generalized inverse theory, J. Geophys. Res., v. 74, pp. 2783-2785. Sokolnikoff, I.S., and R.M. Redheffer, 1966, Mathematics of Physics and Modern Engineering, 752 pp., McGraw-Hill, New York. Talwani, M., J.L. Worzel, and M. Landisman, 1959, Rapid gravity computations for two-dimensional bodies with application to the Mendocino submarine fracture zone, J. Geophys. Res., v. 64, pp. 49-59. Talwani, M., and M. Ewing, 1960, Rapid computation of gravitational attrac-tion of three-dimensional bodies of arbitrary shape, Geophysics, v. 25, pp. 203-225. Tanner, J.G., 1967, An automated method of gravity interpretation, Geophys. J., v. 13, pp. 339-347. Ulrych, T.J., 1968, Effect of wavelength f i l t e r i n g on the shape of the r e s i -dual anomaly, Geophysics, v. 33, pp. 1015-1018. • , 19 72, Maximum entropy power spectrum of truncated sinusoids, J. Geophys. Res., v. 77, pp. 1396-1400. Wiggins, R.A., 1968, Terrestrial variational tables for the periods and attenuation of the free oscillations, Phys. Earth Planet. Interiors, v. 1, pp. 201-266. , 1969, Monte Carlo inversions of body-wave observations, J. Geophys. Res., v. 74, pp. 3171-3181. 99, Wiggins, R.A., 1972, The general l i n e a r inverse problem: i m p l i c a t i o n of surface waves and free o s c i l l a t i o n s f o r earth s t r u c t u r e , Rev. Geophys. and Space Phys., v. 10, pp. 251-285. 100 APPENDIX A Mathematical Concepts Centroid: The coordinates of the centroid of N points in M-dimensional space are: 1 N m. = — E x^. j = l,m (Cooley and Lohnes, 1971, p.30) ^ i = i x3 Conjugate directions: Any two directions (i.e. unit vectors) r^, are conjugate for a given matrix A i f r^ A r 2 = 0 (Pierre, 1969, p.314) Convex region: A region i s convex i f the line segment between any two points of the region is contained entirely within the region. (Beveridge and Schechter, 1970, p.113) Distance in an inner product space is given by the norm of X^ - X 2» i.e. D = || X x-X 2 || = ( X r X 2 , X x - X 2 ) 1 / 2 where X^, X 2 denote two points i n the space. Depending on the type of inner product defined, this might be D 2 = | (x x -J v 2 2 N 2 x,) dV or D = Z ^ ( x u - x 2 1 ) Z Frechet d i f f e r e n t i a l : A function is Frechet differentiable i f 101 E(m + dm) = E(m) + (F,dm) + e(dm) where (F,dm) i s an inner product and e(dm) approaches zero faster than ||dm|| F is called the Frechet derivative or Frechet kernel of E. (Backus and Gilbert, 1967, p. 249) Hilbert space: a infinite-dimensional linear vector space, over the f i e l d of complex numbers, with an associated inner product whose properties are: a) (x,x) = 0 i f and only i f x = 0 b) (x,x) > 0 c) (x+y,z) = (x,z) + (y,z) d) (ax,y) = a(x,y) it e) (x,y) = (y,x) * = complex conjugate (Dunford and Schwartz, 1958, p.242) Inner product: A scalar valued function of a pair of vectors (x,y) of the properties shown above. Examples are: (A,B) = £ E A...B * (f,g) = f(t) g*(t) dt i j J 0 (Hoffman and Kunze, 1961, pp. 221-222) Norm: The norm on an inner product space i s ||x|| = (x.x) 1^ 2 (Dunford and Schwartz, 1958) 102 APPENDIX B Units In Gravity Modelling Before applying the equations for gravitational attraction of the cells comprising the linear model, they are adapted to allow convenient units for the density and distance parameters. The units to be used are: 3 2 [g ] = m i l l i g a l =10 x cm/sec z [p] = gram/cubic centimeter [z], Ix] = kilometer The ver t i c a l gravity due to a horizontal cylinder, which we take to appro-ximate the gravity of a prism, i s = (2.) CY) ( Z q ) (mass) ( 1 )  Z 2 (x-x r . 2 o + z o —8 2 2 The gravitational constant i s y = 6.67 x 10 dyne-cm /gm = 6.67 x 10 ^  cm~Vsec2-gm (mass) here i s mass per unit length, equivalent to a horizontal prism dx by dz, of density p. i.e. (mass) = p(dx)(dz) The appropriate units are 2 3 (mass) = gm-km /cm Thus g z in (1) has units z I g J = IYJ (mass) I ° z J L' J v ' L . ,2 . 2 (x-x ) + z o o 3 2 _ cm x gm-km x 1_ 2 3 km sec -gm cm = km/ sec 2 103 To transform to m i l l i g a l , we multiply by 10 , since 1 km = 10 cm and -3 2 1 mgal = 10 cm/sec . We are thus using a gravitational constant Y = 6.67 in calculating the Frechet kernels. 104 APPENDIX C Derivations For Backus-Gilbert Inversions I. We wish to find a model M which satisfies the observed gravity data, i.e. g (M) =. X j = l,k (1) and i s the "closest" such model to an i n i t i a l estimate M . The model para-meters are the densities of various rectangular prisms; linearity of g with density enables us to write gjCM) = CG^M) (2) th where G^  i s the Frechet kernel for gravity at the j station and CGj,M) denotes an inner product In parameter space, the distance of M from M,., is given by Or [|M-MG||2 = CM-MG, M-M6) = | (P-P G) 2 dV (3) where p and p are the densities of the models, and V denotes the subsurface region containing the models. The gravity effect of the i n i t i a l model i s g jCM G).= CG.,MG) which i s combined with CI) and C2) to give Xj - g jCM G) = CGy M-MG) = | GjCp-Pg) dV j = 1,K (4) The problem i s then to minimize C3), subject to the constraints (4); which can be solved using Lagrangian multipliers CSokolnikoff and Redheffer, 1966, p.346). The minimum is located by finding a stationary point of 105 F(P) = \ ( V ( P _ P G ) 2 d V " I " j ^ j ^ + J v G j ( P " P G } d V " X j ] ( 5 ) 3F Taking -r— = 0, and observing that dp i 1 av , . [ v ( P-p G) dV ap j„ GJ ( p- pG> d v - ( G. J v J dV We have dV = 0 For the solution to be valid at a l l points of V, we require K P - P G + E 0 j G j for a l l points of V, which i s written as K M = M- + E a. G,, G 1 2 j (6) To evaluate the , substitute (6) into (4) \ A - g 4(M 0) = G4(p-p„) dV JV 1 G, G, dV or 106 K (7) in the inner product notation. Equation (7) is a linear system which can be solved for ct^  by standard matrix methods, since the inner products (G,,G.) form a square K J matrix of dimension K. The steps in the solution are: a) Compute the Frechet kernels and their inner products b) Compute gravity of i n i t i a l model c) Solve (7) for the Lagrangian multipliers ct^  d) Use the ct^ in (6) to find the f i n a l model M. II. Similar systems can be devised to optimize other properties of the model. One possibility i s to minimize ||MR|| , where R i s a spatial weight-ing factor, which allows us to discriminate against subregions of V i n finding a model. In this case the constraints analogous to (4) are simply Aj - 8 j 0 0 and (5) becomes F(P) = 7 f P 2 R 2 dV - Z <x.[ f G. p dV - X . ] •J V 1 J J V J - J 8F Taking = 0 as before, the solution i s K CL G P = I (8) 1 R Again, we find the ct^ by applying the constraints to get 107 K X. = | 3 1 K K dV or • R 2 K G G. in the inner product notation. This is also a linear system, which we solve for ctj^  to use in (8) . III. The most flexible approach i s to combine the methods of I and II, allowing use of an i n i t i a l model as well as spatial weighting factors. We now want to minimize | | (M-M )R| I 2 = f (p- 2 2 P G r R dV Equation (9) becomes K G G and the f i n a l model is given by K o G p = p + Z - L i (11) 8 1 R Z Note. A l l these derivations have written inner products as integrals, e.g. •I, g j - C P dV For the multi-prism models, the inner products are really summations of the 108 form g = E E G. p j inm Knm J n m J where j denotes surface position and [n,m] indicates subsurface location of the prism Thus in applying these results, the integrals are replaced by the appro-priate summations. The derivations could be repeated using summations, but this would not change the results. Notation. M denotes a complete model; i.e. a set of densities. p denotes density as a spatial variable, or density of individual prisms. ! 109 APPENDIX D Numerical Optimization Techniques Finding a model to f i t any geophysical data very often requires optimization processes. When the data functionals are non-linear, one fre-quently seeks unique model parameters by minimizing errors of the model predictions compared to the observed data. If a non-unique class of models is being considered, a particular model can be selected by optimizing some property of the model with the constraints that the predictions of the model agree with the observations (within some l i m i t ) . The nature of the model, the number of parameters, and the type of "objective function" to be opti-mized often make the analytical methods of variational calculus unsuitable. In these cases, the basic approach to the problem i s to seek an optimal point in parameter space; i.e. the N model parameters are each considered as coordinates of an N-dimensional space, and the optimum model i s then a point i n that space whose location i s determined by the optimum value of the objective function. Parameter optimization methods have seen some geophysical applica-tions i n recent years. Rastrigin (1965) developed a theoretical procedure to minimize error between model predictions and observed data, and dis-cussed some simple methods of "searching" for the optimum point. Al-Chalabi (1971, 19 72) used searches i n parameter space to find polygonal models (following the Talwani procedure) which gave optimum f i t s to gravity profiles. The HEDGEHOG and EDGEHOG methods (Jackson, 1973) employ searches to produce density and velocity models of the upper mantle from seismic data (EDGE-110 HOG is a refinement aimed at defining the envelope of a l l models which f i t the data). A considerable literature has now developed on the subject of numerical optimization. Cooper and Steinberg (1970) gave a good review of basic mathematical concepts, along with descriptions of several commonly used procedures. Box et al (1969) b r i e f l y reviewed several methods with a view towards choosing methods for a particular problem. Pierre (1969) described the commonly used techniques, with some discussion of computer application. An excellent book by Beveridge and Schechter (1970) examined the subject in detail, recommending methods proven in practice, and discussing the procedures best suited to particular problems. In general, numerical optimization methods attempt to locate an optimum in-a series of moves starting from a given point. The basis of any method i s the means of selecting the direction and length of each move. Methods can be roughly classified using two basic c r i t e r i a : (1) whether or not constraints on the coordinates (parameters) can be incorporated; and (2) whether or not gradients are used. The usual terminology is that a "direct" search does not employ gradients. The efficiency and speed of search differs between these groups; however the type of problem to be sol-ved may dictate which class of procedure to use. As a general rule, un-constrained gradient methods are fastest. An evaluation of methods f a l l i n g into each class can be found in the following sections of Beveridge and Schechter: 1) Unconstrained direct search pp. 363-406 I l l 2) Unconstrained gradient search pp. 407-433 3) Constrained direct search pp. 448-456 4) . Constrained gradient search pp. 435-448; 457-483 The following brief review follows their descriptions. Unconstrained optimization i s naturally much faster than a con-strained problem, since motion in parameter space need not be confined to particular regions. Direct search methods may either move in prescribed directions, or use values of the objective function to determine the search direction. The basic systematic search method is univariant search, i n which an optimum is found along each coordinate axis from the i n i t i a l point; unfortunately this i s usually inef f i c i e n t . Rosenbrock's method (Beveridge and Schechter, 1970, p.404) is an effective and reliable adaptation of the basic idea, allowing acceleration in distance and change of direction. The Sequential Simplex method (p.372) determines moves from objective function values; basically one evaluates the objective function at the vertices of an n-sided "simplex" in parameter space, and replaces the vertex of worst function value by reflecting i t through the centroid of the figure. Faster convergence i s often possible i f gradients of the objec-tive function can be used. The simplest application i s to move a given distance i n the direction of steepest ascent or descent (i.e. the gradient direction) in hopes of eventually finding a stationary value of the objective function. The method i s not usually recommended, since the gradient i s often only a local property which does not "point" to the optimum; in addi-tion, oscillations about the optimum can result from the fixed distance of 112 each step. An efficient method which w i l l converge i s searching along con-jugate directions (See Appendix A for definition) un t i l the gradient direction becomes.perpendicular, whereupon new conjugate directions are selected. Several techniques have been suggested for selecting conjugate directions; the frequently recommended Davidon-Fletcher-Powell method (Box et a l , p.30; Pierre, p.320; Beveridge and Schechter, p.426) employs the gradient of the objective function at the i n i t i a l point. Constrained optimization i s a more d i f f i c u l t task, and a common approach i s to avoid constraints as much as possible. In some cases, this can be achieved by transformation of variables (Box, 1966). Another method is to create a new objective function by adding the constraints (multiplied by a suitably large constant) and then use an unconstrained method (Beveridge and Schechter, pp.443-448; 477-482). Incorporating constraints as "penalty functions" i n this way keeps the search out of unfeasible regions (where the constraints are not satisfied and the modified objective function has a value far from the optimum). The penalty function methods are often im-practical, since the objective function is distorted near boundaries, and since the choice of weighting factors to apply to violated constraints may be c r i t i c a l to the success of the optimization. Many techniques have been developed to incorporate constraints directly. In general, a l l constraints may be considered to be inequalities, since an equality can be expressed as two inequalities, and since equalities involving real data usually have some allowable error. The most widely used methods are probably those of linear programming (Beveridge and Schechter, 113 pp. 287-322), which deal with problems where a l l constraints and the object-ive function are linear functions of the parameters; the feasible region then has linear boundaries. Since the objective function i s also linear, the optimum must l i e on a corner of the feasible region, and ef f i c i e n t matrix methods are available to find this optimum corner. More general problems may be solved with variations of the unconstrained methods. The usually recommended direct method i s the Complex method (Box, 1965), which stems from the Sequential Simplex method. Details of the method are given in Appendix E. Gradient methods usually require linearization of constraints to be most effective, and usually incorporate constraints into the choice of search direction to avoid leaving the feasible region. The hemstitching technique (Beveridge and Schechter, p.456) allows any method of choosing the search direction, but follows the gradient of a constraint to move back into the feasible region whenever any constraint i s violated; unfortunately this often means that the search i s not moving towards the optimum point. Rosen's gradient projection (Beveridge and Schechter, p.469), or the method of riding constraints (Box et a l , p.42) always use the constraints in choosing the search direction to ensure that the search stays in the feasible region. Linear constraints can be added to the Davidon-Fletcher-Powell method of choosing conjugate directions (Box et a l , p.47), and thus this e f f i c i e n t method can be applied to some constrained problems. Beveridge and Schechter suggest that most of these methods have the basic problems of steepest ascent searches, and recommend the Complex method for most non-linear constrained problems, particularly i f the optimum is expected near a boundary of the region. 114 APPENDIX E The Complex Method The Complex method of numerical optimization (Box, 1965) can be applied to many different problems, since i t does not require any special type of objective function. The search procedure is similar to the Sequen-t i a l Simplex method (see Appendix D), but incorporates constraints on the parameters of the objective function. The only problem in implementation is that an i n i t i a l feasible point (i.e. one which satisfies a l l the con-straints) is required before the search begins. The allowable constraints are bounds on a l l the model parameters, A t < x ± < B ± i = 1,N (1) plus implicit constraints of the form 7i « j = 1,M (2) where y^ are functions of the x^. With a known i n i t i a l point, the other vertices of the "complex" are found as follows: the coordinates x. of each x point are generated randomly in the intervals (A^,B^); the quantities y^ are calculated for the point; i f any of the implicit constraints are not sa t i s -fied, the point i s moved halfway towards the centroid of feasible points as many times as are necessary to make the t r i a l point satisfy a l l the con-straints. This method ensures that other feasible points can be found, so long as the region defined by the constraints is convex (see Appendix A). The process i s repeated u n t i l a complex of k vertices (k>N+l) i s established (each vertex is a point in the "feasible" region of the parameter space, i.e. the region where a l l the constraints are satisfied). 115 Once the complex has been set, the objective function i s evaluated at each vertex, and the vertex of worst value (smallest f o r maximization; l a r g e s t . f o r minimization) i s replaced by o v e r - r e f l e c t i o n through the centroid of remaining points. The coordinates of the new vertex are X'. _ = a(X. _ - C.) + C. (3) i,R i,R 1 1 where X. _ = coordinates of worst vertex i,R X 1. = coordinates of replacement i,K C. = coordinates of the centroid of a l l other v e r t i c e s , l given by C i = K=T ^ = 1 X i , j - Xi,R> W and a i s an expansion parameter (Box suggests a=1.3) Of course, the replacement point must s a t i s f y a l l the constraints (equations (1) and ( 2 ) ) . I f any coordinate X 1 does not l i e i n the required i n t e r v a l (A^,B^), i t i s given a value j u s t i n s i d e the appropriate bound. I f any im-p l i c i t c o n s t r a i n t (equation 2) i s v i o l a t e d , the t r i a l point i s moved halfway towards the centroid of other points and tested again; the new point i s then X " i , R = 0 ' 5 < X , i , R + V <5> Once an acceptable new point i s found, the r e j e c t i o n process i s repeated with the current worst po i n t . The search could s t a l l i f the replacement point s t i l l has the worst function value. To avoid t h i s problem, Box 116 suggested moving such a point halfway toward the centroid, rather than reflecting i t back again. Guin (1968) suggested some modifications to the method, which allow i t to be used in non-convex regions. In a non-convex region, the cen-troid may eventually f a l l outside the region, making i t impossible to find new feasible points by moving towards the centroid. In this eventuality, Guin suggested creating a new complex using the "best" point as the i n i t i a l point, and creating other points by random generation in the interval bet-ween the best point and the old centroid. The coordinates of a new point are then X. . - X. + R (C. - X. ) (6) i , j i,o l i,o' where X. are the coordinates of the best point C^ are coordinates of the old centroid R is a random number (0<R<1) Each point i s tested for f e a s i b i l i t y and moved towards the current centroid i f necessary. Guin's other ideas were to reject the second worst point i f one point is continually the worst, and to move a t r i a l point halfway towards the centroid i f a bound on a variable (i.e. coordinate) i s exceeded. An additional feature has been implemented by the University of B.C. Computing Centre. They suggest periodically restarting the complex (following Box's method) to accelerate searching when far from the optimum 117 (O'Reilly, 1971); a new complex is created from the best point after a spe-c i f i e d number of iterations, but not after a maximum number has been passed. The actual search procedure used for the minimum entropy problem follows Box's basic method, incorporating some of the other suggestions and some new variations. Restarting is optional after a specified number of iterations; the new complex may be produced by random number generation i n the complete interval (A^,B^), or in the region near the best point and the old centroid. In the latter case, Guin's method has been slightly altered by using a random number (R in equation (6)) evenly distributed between -1 and +1, thus the coordinates of a new point can f a l l on both sides of the best point relative to the old centroid. Guin's suggestion for generating a new complex i f the centroid becomes an unfeasible point was adopted. When one point i s consistently the worst, a new complex is generated from the current best best point, since this problem was found to be an indication that the search was not progressing very well. This strategy proved more successful than Box's suggestion, since repeated moves toward the centroid can lead to apparent convergence of the search, i f the centroid i t s e l f has the worst function value. When the optimum has been reached, a l l the vertices of the complex should approach the centroid; hence the convergence test i s to evaluate the objective function at the centroid, and stop the search when several suc-cessive evaluations agree within a specified limit. In addition, a specified maximum number of iterations w i l l stop searching, to avoid using excessive computer time when progress is slow. 118 APPENDIX F Computer Program The required inputs are as follows. Ml, M2, NI, N2 define the locations of blocks to be used in the model. JA = 1 defines the f i r s t gravity station CM i s measured from JA). JB = number of data i n the profile. DX = DZ = station spacing = c e l l dimension. WT, WT2, WT3 are weighting factors to apply to each c e l l i n the model. REGCN,M) an indicator C0.0, 1.0, 2.0, 3.0) of which weighting factor to apply Cl-0, WT, WT2, WT3). RH0CN,M) = c e l l densities for the i n i t i a l model. GZOCJ) = the gravity pr o f i l e . Other routines mentioned. FSLE: a UBC Computing Centre routine which solves IPR0DCJ,K)J IYCJ)J = JBCK)]. HIST: an optional line-printer plot of the data. In SMOOTH, C i s the desired increment of density to be displayed. C "WEIGHT" C ** C A PROGRAM TO COMPOTE MODELS FROM WEIGHTED REGIONS OF THE C SUBSURFACE GRID REAL Z (10) ,G (30,10,30) ,PROD (30,30) ,GZO (30) , E (30) , DENS (10, 30) , 1X (30) ,GZM (30) ,T (30, 30) ,RR (10,30) ,Y (30) ,MASS,REG (10,30) , 2RHO (10,30) ,B (30) ,GZ (30) ,YNAM (100) ,XNAM (100) INTEGER IPERM(100) C EACH CELL HAS AN ASSIGNED WEIGHT OF "WT" OR 1.0,DEPENDING ON C THE SPEC I F I E D VALUE OF REG(N,M) : THE MODEL IS FOUND BY MINIMIZING C THE SUM OF (|DENS-RHO|*RR) SQUARED READ(5,3) M1,M2,N1,N2,JA,JB,DX,DZ,WT,WT2,WT3 DO 5 1=1,30 Z (I)= (I-.5) *DZ 5 X(I)=I*DX 3 FORMAT (615, 5F10. 3) C THE I N I T I A L MODEL AND REGIONS TO BE WEIGHTED ARE GIVEN READ (5,4) ( (REG (N,M) ,N=1, 10) ,M=M1,M2) 4 FORMAT (10F5.2) READ (5,4) ( (RHO (N,M) ,N=1, 10) ,M=M1,M2) C DISPLAY THE I N I T I A L MODEL WRITE(6,16) (Z (N) , N=N 1 , N 2) 16 FORMAT (//IX,'INITIAL MODEL',/,1X,• X (KM) / Z= «,10F8.2) WRITE(6,31) DO 17 M=M1,M2 17 WRITE(6,40) X (M) , (RHO (N , M) , N=N 1, N2) C ** C THE FRECHET KERNELS FOR GZ ARE COMPUTED; "C" CORRECTS THE CYLINDER C EXPRESSION TO THAT OF A RECTANGULAR PRISM OF EQUAL MASS DO 10 J=JA,JB DO 10 M=M1,M2 DO 10 N=N1,N2 IF (J.EQ.M.AND.N.EQ.1)C=.86601 I F (IABS (J-M) . EQ. 1. AND.N.EQ. 1) C=.98145 IF (J.EQ.M.AND.N.EQ.2) C=0.99676 I F (IABS (J-M) . GT. 1.0R.N.GT. 1) C=1.0 RR (N,M) =1 .0 IF (REG(N,M) .EQ. 1.0) RR(N,M)=WT IF (REG (N,M).EQ.2.0) RR(N,M)=WT2 I F (REG (N,M) . EQ. 3.0) RR(N,M)=WT3 10 G(J,N,M) =2.0*C*Z (N) *DX*DZ*6.67/ ( (X (J ) - X (M)) **2+Z (N) *Z (N) ) C THE WEIGHTED REGIONS ARE DISPLAYED WRITE (6,27) (Z (N) ,N=N1 ,N2) 27 FORMAT(//1X,'WEIGHTED REGIONS',/,1X, 1' X (KM) / Z= «,10F8.2) WRITE (6, 31) DO 28 M=M1,M2 28 WRITE(6,40) X (M) , (REG(N,M) ,N=N1,N2) WRITE(6,29) WT,WT2,WT3 29 FORMAT(/1X,'THE WEIGHTING FACTOR IN CELLS MARKED 1.0 IS',F9.1 2/,IX,•THE WEIGHTING FACTOR IN CELLS MARKED 2.0 I S ' , F 9 . 1 , 3/,1X,'THE WEIGHTING FACTOR IN CELLS MARKED 3.0 I S ' , F 9 . 1 , 1/,IX,'OTHER REGIONS HAVE UNIT WEIGHT IN THE MINIMIZATION *) C *** C THE INNER PRODUCTS OF THE FRECHET KERNELS ARE GENERATED DO 15 J=JA,JB DO 15 K=JA,JB PROD (J , K) =0.0 DO 15 M=M1,M2 DO 15 N-N1,N2 15 PROD(J,K)=PROD(J,K) +G(J,N,M) *G(K,N,K)/RR (N,M) C * * * C THE LAGRANGE MULTIPLIERS ARE FOUND BY SOLVING THE MATRIX PROD; C " F S L E " IS A COMPUTING CENTRE ROUTINE FOR SIMULTANEOUS EQUATIONS READ (5, 20) (GZO (J) , J=JA, JB) 20 FORMAT (6F12.3) WRITE(6,18) 18 FORMAT(//1X,« X GZ (MODEL) GZ (REAL) DIFFERENCE'/) DO 22 J=JA,JB GZ (J)=0.0 DO 21 M-=M1,M2 DO 21 N=N1,N2 21 GZ (J) = GZ (J) + G(J,N, M) *R HO (N , M) B (J)=GZO (J ) - G Z (J) 22 WRITE(6,23) J , GZ (J) , GZO (J) , B (J) 23 FORMAT (1X,14,3F10.4) JLN=JB-JA*1 LENA = 3 0 LENBX=30 LENT=30 NS0L=1 CALL FSLE(JLN,LENA,PROD,NSOL,LENBX,B,Y,IPERM,LENT,T,DET,J EXP) C • * * C THE FINAL MODEL IS COMPUTED USING THE LAGRANGE MULTIPLIERS Y C * * * DO 26 M=M1,M2 DO 26 N=N1,N2 DENS(N rM)=0.0 DO 25 J=JA,JB 25 DENS (N,M)=DENS (N,M) + Y (J ) * G (J,N,M) 26 DENS (N, M) = DENS (N, M) /RR (N , M) +RHO (N,M) WRITE (6,30) (Z (N),N=N1,N2) 30 FORMAT (1X,//,1X,'COMPUTED DENSITIES',//,1X,• X (KM) / Z= •, 1 10F8.2) WRITE (6,31) 31 FORMAT (1X, * •) DO 35 M=M1,M2 35 WRITE (6,40) X (M) , (DENS (N,M) ,N=N1,N2) 40 FORMAT (1X,F8.2, • : «,10F8.3) C * * * C THE GRAVITY EFFECT AND MASS OF THE FINAL MODEL ARE COMPUTED C AND DISPLAYED WRITE (6,42) 42 FORMAT (1X,//,1X,'OBSERVED GRAVITY - MODEL GRAVITY - ERROR'//) DO 50 J=JA,JB GZM(J)=0.0 MASS=0.0 DO 45 M=M1,M2 DO 45 N=N1,N2 MASS=MASS+DENS (N,M) 45 GZM (J) =GZM (J) +DENS (N , M) *G(J,N,M) E (J) =GZO ( J )-GZM (J) 50 WRITE(6,55) J , GZO (J) , GZM (J) , E (J) 55 FORMAT (1X,15,2F11.4 ,F 10.4) WRITE (6,60) MASS 60 FORMAT ( 1 X , / / / , 1 X T O T A L MASS IS ',F10.4//) C PLOT THE GRAVITY PROFILE READ (6,65) XMN,XMX 65 FORMAT (2F10.4) READ (5,70) (YNAM (I) ,1=1 ,80) READ (5,70) (XNAM (I) , 1=1, 80) 70 FORMAT (80A1) CALL HIST(GZO,JLN,DX,YNAM,XNAM,XMN,XMX,1) C COMPOTE AN APPROXIMATE MODEL CALL SMOOTH(DENS,G,M1,N1,M2,N2,X,Z,MASS,JA,JB,GZO) STOP END SUBROUTINE SMOOTH(DENS,G,M1,N1,M2,N2,X,Z,MASS,JA,JB,GZO) REAL ACENS (10,30) ,GZO (30) ,GA (30) ,MASS,AMASS, ER (30) ,C,X(30) 1Z (30) , DENS (10, 30) ,G (30, 10,30) 9 READ (5,10) C 10 FORMAT (F10.4) IF (C.EQ.0.0) GO TO 80 12 DO 13 M=M1,M2 DO 13 N=N1,N2 ADENS (N, M) =C*FLOAT (I FIX ( (DENS (N,M) +0.5*C) /C) ) IF (ADENS (N,M) .LE.0.0) ADENS (N,M)=0.0 13 CONTINUE C NOW COMPUTE AND COMPARE GRAVITIES DO 15 J=JA,JB GA(J)=0.0 AMASS=0.0 DO 14 M=M1,M2 DO 14 N=N1,N2 GA (J)=GA (J)+ADENS (N,M)*G (J,N,M) 14 AMASS=AMASS + ADENS(N,M) 15 ER (J)=GZO (J)-GA (J) WRITE (6,30) (Z (N) ,N=N 1,N2) 30 FORMAT(//1X,'APPROXIMATE DENSITIES ',//,1X, 1'X (KM) / Z = « , 10F8.2) WRITE (6,33) 33 FORMAT (1X, • ') DO 35 M=M1,M2 35 WRITE (6,37) X (M) , (ADENS (N , K) , N=N 1 , N2) 37 FORMAT(1X,F8.2,' : »,10F8.2) WRITE (6,42) 42 FORMAT (//1X,'EXACT GRAVITY - APPKOX GRAVITY - ERROR'/) DO 50 J=JA,JB 50 WRITE (6,55) J,GZO (J) ,GA (J) , ER (J) 55 FORMAT (1X, 15, 2F11. 4, F10. 4) WRITE (6,60) MASS, AMASS 60 FORMAT (1X,///, 'TRUE MASS IS ',F10.4,/,1X, 1'APPROX. MASS IS ',F10.4) GO TO 9 80 RETURN END Co 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052851/manifest

Comment

Related Items