INVERSION OFTHREE-DIMENSIONAL DIRECT CURRENT RESISTIVITY DATAbyYaoguo LiB.A.Sc. (Geophysics), Wuhan College of Geology, 1983A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDEPARTMENT OF GEOPHYSICS AND ASTRONOMYWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIAMarch 1992©YaoguoLi, 1992In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.Department of Geophysics & AstronomyThe University of British ColumbiaVancouver, CanadaDate April 22, 1992DE-6 (2/88)ABSTRACTA direct current (d.c.) resistivity experiment investigates subsurface geo-electrical structuresby measuring the electric field set up by introducing current into the earth. Information aboutgeo-electrical structures is extracted by inverting the observed data to generate an image of theconductivity or to construct a conductivity model. The goal of this thesis is to develop efficientinversion techniques for the interpretation of three-dimensional (3d) d.c. resistivity data. Thestudy assumes data consisting of pole-pole potentials measured over a regular grid on the surfacefor many current locations. The Born approximation is employed to linearize the inverse problem.The source of the electric field measured in the d.c. resistivity is the accumulated electriccharges. Different aspects of the charge accumulation are reviewed, enlarged with new insightsand presented in a unified notation. This provides the basis for understanding the fundamentalsof d.c. resistivity experiments. Two algorithms are developed to image simple 2d conductivities.The first constructs a structural image by combining the charge density images obtained byinverting multiple sets of common current potentials. The second constructs a conductivityimage directly. Processing and displaying the apparent conductivity, and constructing equivalentsources from secondary potentials are studied as the means of imaging. Assuming a multiplicativeperturbation to a uniform half-space, the potential anomaly of pole-pole arrays is expressed as adepth integral of the logarithmic perturbation convolved with a kernel function in the horizontaldirections. Applying the Fourier transform decomposes the data equation for a 3d problem into aset of id equations. A rapid approximate 3d inversion is developed based upon this decompositionby solving a sequence of id inversions in the wavenumber domain. The approximate 3d inversionis used to construct iterative inversion algorithms using the AIM (Approximate Inverse Mapping)formalism. The approximate inversion and an exact forward mapping are used to update themodel successively so that the final result reproduces the observed data. The AIM inversion isapplied to analyse a set of field data.11TABLE OF CONTENTSAbstract iiTable of Contents iiiList of Tables vList of Figures viAcknowledgment ixChapter 1 Introduction IChapter 2 Theoretical Basis of d.c. Resistivity Experiments 112.1 Basic Equations 132.2 Current Flow in d.c. Resistivity Problems 182.3 An Integral Equation for the Charge Density 202.4 Effect of Charge Density at the Earth’s Surface 232.5 An Analytic Example of Charge Accumulation 282.6 Image Method Explained by Charge Accumulation 302.7 Numerical Examples of Charge Distribution 352.8 Born Equation for d.c. Responses 402.9 Discussion 43Chapter 3 Imaging Subsurface Structures 463.1 Direct Imaging via Apparent Conductivity 473.2 Imaging via Equivalent Source from Secondary Potential 543.3 Charge Density Imaging of 2d Structures 583.4 Conductivity Imaging of 2d Structures 623.5 Discussion 66Chapter 4 approximate 3d Inversion 674.1 Data Equation 674.2 Data 721114.3 Kernel Functions .744.4 Inversion 804.4.1 Continuous Formulation 824.4.2 Discrete Formulation 874.5 Calculation of Background Conductivity 894.6 Synthetic Examples 924.7 Discussion 97Chapter 5 AIM Inversion 995.1 Review of the AIM Formalism 995.1.1 AIM-MS 1015.1.2AIM-DS 1035.2 AIM Inversion of 3d d.c. Resistivity Data 1055.3 Examples of the AIM Inversion 1095.3.1 The 5-prism model 1105.3.2 The geo-statistical model 1215.4 Discussion 132Chapter 6 Analysis of Field Data 1346.1 Geologic Objectives 1346.2 Data Acquisition 1356.3 Data Processing 1406.4 Inversion 1436.5 Discussion 153Chapter 7 Conclusions 155ivLIST OF TABLES4.1 Parameters of the 5-prism model.93VLIST OF FIGURES2.1 Basic geometry for current flow at a boundary 152.2 Change of the electric field by the accumulated charges 192.3 Geo-electrical model 222.4 A point source in a layer over a half-space 282.5 The charge density on a 2d topography 372.6 The charge density on a buried prism 392.7 The secondary potential above the buried prism 402.8 The comparison of the true and Born potentials 443.1 Varying depth mapping for a three-layered background 513.2 Comparison of the pole-pole and square array responses 533.3 Equivalent source imaging 573.4 Charge density imaging 613.5 Conductivity imaging 654.1 Geometry of the surface pole-pole arrays 694.2 Specification of the data map 734.3 Spatial domain kernel function 754.4 Amplitude change of the wavenumber domain kernel function 774.5 Wavenumber domain kernel function 794.6 A set of id kernels with different array separations 804.7 A set of id kernels at different wavenumbers 814.8 The perspective view of the 5-prism model 934.9 The comparison of the true and recovered 5-prism model at x = 450m 944.10 The comparison of the true and recovered 5-prism model at z = 30m 954.11 The comparison of the true and recovered 5-prism model at z = I 50m 96vi5.1 Comparison of the 5-prism model and the AIM-MS inversion result at x = 460 Iii5.2 Comparison of the 5-prism model and the AIM-MS inversion result at z =20m, I 50m .. 1125.3 Misfit reduction in the AIM-MS inversion for the 5-prism model 1135.4 AIM-DS inversion result for the 5-prism model 1155.5 Misfit reduction in AIM-DS inversion for the 5-prism model 1165.6 AIM-DS inversion result of noisy data from 5-prism model 1185.7 Misfit reduction in the AIM-DS inversion of noisy data from the 5-prism model 1205.8 Comparison of the geo-statisticai model and the AIM-MS inversion result in layer-i andlayer-3 1225.9 Comparison of the geo-statistical model and the AIM-MS inversion result in layer-6 andlayer-lO 1235.10 Comparison of the geo-statistical model and AIM-MS inversion result at x =420m .... 1255.11 Detailed comparison of the geo-statistical model and the AIM-MS inversion result at x =420mabovez=200m 1265.12 Misfit reduction in the AIM-MS inversion for the geo-statistical model 1275.13 AIM-DS inversion result for the geo-statistical model 1285.14 Misfit reduction in the AIM-DS inversion for the geo-statistical model 1295.15 The AIM-DS inversion result of noisy data from the geo-statistical model 1305.16 Misfit in the AIM-DS inversion of noisy data from the geo-statistical model 1316.1 Field acquisition procedure of the B-SCAN® experiment 1366.2 Grid layout of the field data set 1386.3 Four different patterns of data collection 1396.4 Composition of a common current gather 1426.5 Two apparent conductivity maps from the raw data 1446.6 Two processed apparent conductivity maps 1456.7 Two apparent conductivity pseudo-sections 146VII6.8 Misfit reduction in the AIM-DS inversion 1476.9 Misfit reduction in the AIM-DS inversion with a relaxation parameter 1486.10 Two predicted apparent conductivity maps 1506.11 Two predicted apparent conductivity pseudo-sections 1516.12 Two sections from the constructed conductivity model 152vii’ACKNOWLEDGMENTSTo my parents:—53fj!Jl.I express my deepest thanks to my supervisor, Dr. Doug Oldenburg, for his invaluable adviceand constant encouragement in every aspect of this work, and for his moral and financial supportthroughout my stay at UBC. He has made possible my attendance at several conferences. Hisgenerous help has made rough times so much easier for me to deal with.I thank Dr. Rob Ellis, Dr. Stan Dosso, David Lumley, Guy Cross and David Aldridge formany helpful and interesting discussions. I thank Dr. Rob Ellis, John Amor and Dr. GerryGrieve for their help and expertise with the computing facilities. I am grateful to Drs R. M.Ellis, R. J. Knight and R. D. Russell for acting as my committee members. I sincerely thankErik Blake for the help he has offered in my adjusting to life in Canada, and all my friends inthe Department of Geophysics and Astronomy who have made my studies here so enjoyable.The research was supported by a University of British Columbia Graduate Fellowship, NSERCgrant 5-84270, NSERC CRD grant 5-80141, and by E-SCAN Technologies Inc. and FMC GoldCompany. Finally, I sincerely thank my wife Yumei Hu, who has accompanied me across thePacific Ocean and supported me wholeheartedly in my pursuits.ixCHAPTER 1INTRODUCTIONThe direct current (d.c.) resistivity experiment is one of the oldest geophysical methodsdevised to investigate subsurface structures. An electric field in the earth is set up by puttingsteady current into the ground through a pair of electrodes. A separate pair of electrodes is usedto measure the field in the form of an electrical potential difference. When an electric currentfield is introduced into the earth, it sets up a distribution of accumulated electric charges bothon and beneath the earth’s surface. These charges exist in the region where there is a gradientof conductivity and a non-zero electric field parallel to it. It is these accumulated charges whichgive rise to the electric field whose potentials are being measured. The distribution of the chargeis determined by the conductivity structure. Thus different conductivity structures will result indifferent distribution of accumulated charges and, hence, different distributions of electric field.The measurement of the field, therefore, contains information about the conductivity.Many different configurations of electrode placement can be used. Electrodes can be placedon the ground or some can be in boreholes. This difference leads to the division between surfaceand borehole d.c. resistivity methods. The surface experiment has been the most commonly used,but recent years have seen a rapid development in borehole d.c. resistivity methods. Measurementsare generally one of the three types (Keller and Frischknecht 1966). The first is the potential inreference to a point at “infinity”. An example of this type is the data measured by a pole-polearray. The second is the potential gradient or the electric field intensity. The data from poledipole, Schiumberger, and Wenner arrays are of this type. The third is the curvature of the electricpotential generated by a point source. A dipole-dipole array collects such data. Among these,the potential data measured by a pole-pole array are the most basic since the other types can besynthesized from the pole-pole data by the principle of superposition.Due to the decay of the field away from the current source, the measurement will only besensitive to a limited region surrounding the active electrodes (those which are not placed atICHAPTER 1. INTRODUCTION 2“infinity”). Each measurement represents an integrated effect of the conductivity within thisregion. The extent of this effective region increases with both the separation of the currentelectrodes and the separation between source and potential electrode pairs and for this reason thed.c. resistivity experiment involves a series of spatially distributed measurements with differentarray separations. This series of measured data collectively carries the information about thesubsurface conductivity structure. Successful processing and inversion of the data yields suchinformation by generating an image of the conductivity structure or by constructing a conductivitymodel.The conductivity of the material in the earth’s crust depends upon a number of factors suchas the composition, structure and the physical processes involved. Within the regime of d.c.resistivity experiments, the conduction of electricity is mainly of two types, namely, electronicconduction and electrolytic conduction (Ward and Fraser 1971, Telford et at. 1976). The formeris by the free electrons in the material while the latter is by the ions in the electrolyte. Only asmall portion of the naturally occurring minerals is a good electrical conductor. Among theseare native metals, graphite, magnetite and most of the suiphides. However, the majority of therock building minerals are very poor conductors and the propagation of electricity in rocks ismainly electrolytic and depends upon the interstitial fluid. Thus the conductivity of the materialis influenced by its porosity, connectivity of the pores, fluid saturation and the concentration ofions. Since temperature changes the concentration and mobility of the ions in the fluids, theconductivity of the rock will increase with the temperature under most near-surface conditions.These dependencies make it possible for the d.c. resistivity experiment to have a wide range ofapplications.The first, and perhaps the most prominent area of application is in mineral exploration. Manymetallic mineral deposits such as massive suiphides have conductivities much higher than thesurrounding rock and form good targets to be located by d.c. resistivity experiments. More often,the local geological environment or the host of the ore deposit has anomalous conductivities.Identifying such geological settings serves to predict the localities of possible ore deposits. ForCHAPTER 1. INTRODUCTION 3instance, many hydrothermal deposits occur on the boundaries of intrusive structures, or theymay be associated with the vein formed in fault structures which are usually highly resistive.Another example is the mineralization resulting from alteration processes in the hydrothermaldeposits. The ascending hot fluid precipitates precious metals by altering the iron rich rock. Theprocess eventually forms a halo containing a significant amount of pyrite surrounding the oredeposits. The conductive halo thus forms a rather unique signature. Much of the literature ond.c. resistivity methods focus upon this application. Another application for the d.c. resistivitysurveys is in the exploration for geothermal resources which are structurally controlled and oftenare associated with conductive anomalies (e.g., Wright et al. 1985).The second area of application for d.c. surveys is in environmental, groundwater andgeotechnical problems. Recent years have witnessed the growing application of d.c. resistivitymethods in such problems. This is evidenced by extensive publications (e.g., Ward 1990). Mostapplications of this type are in the characterization and monitoring of waste disposal sites, theassessment of existing contamination sites, and the study of the transportation of groundwatercontaminants. Conductivity variations are diagnostic to answering specific questions in suchenvironments and d.c. resistivity experiments provide an indirect but effective means to detectthese changes. In contrast to mineral and geothermal exploration, these problems are usuallycarried out on a much smaller spatial scale and consequently the electrode separations are muchsmaller.The third area of application is in archeological problems and a considerable amount of workhas been published. An example is mapping man-made structures which often exhibit highlyresistive signatures in the more conductive sedimentary settings. Scollar et al. (1990) present arather good overview of the application of the d.c. resistivity method in archeology.The applications above are presented in the order of decreased spatial scale of the problem.However, there are also applications in the extreme ends of the scale. One experiment wasconducted to investigate the crustal structures in South Africa (Van Ziji and Joubert 1975),where the current electrode separation exceeded 400 km. On the small scale, techniques haveCHAPTER 1. INTRODUCTION 4been developed which use direct currents to perform a medical scan of the human body (e.g.,Brown 1986, Yorkey and Webster 1987). There are also techniques which use d.c. fieldsto reconstruct tomographically the conductivity of core-samples (Dines and Lytle 1981). Theelectrode separation in these applications can be as small as 1 cm.In exploration problems, the data of surface d.c. resistivity experiments are displayed in theform of apparent resistivity, which is the resistivity of a uniform half-space which produces thesame potential difference as that acquired by the field electrodes. In keeping with the recentlyadopted convention, I choose to work with the apparent conductivity throughout this thesis(conductivity is the reciprocal of resistivity). The display of data in such a form constitutesthe crudest way of extracting information about subsurface conductivity structures, since theapparent conductivity is only a convoluted representation of the true structure. Any quantitativeinformation can only be obtained by inverting the data to obtain conductivity models which canreproduce the observed data.The computational difficulty in inverting a data set increases with the number of the dimensions of the problem. The Id and 2d problems have been studied extensively. The methodsof automated Id inversion generally fall under the category of direct techniques or iterativeapproaches. Direct methods usually employ the properties of the Id field to construct directlya depth-conductivity profile from the given observations (e.g., Koefoed 1976, Szaraniec 1976,Szaraniec 1980, Coen and Yu 1981, Parker 1984). The most rigorous is the bilayer expansionmethod by Parker (1984). On the other hand, iterative approaches start from a initial model anditeratively update the model to produce a final solution based upon linearization of the mappingfrom conductivity to potential (e.g., Inman 1975, Jupp and Vozoff 1975, Oldenburg 1978) or uponsome ad hoc methods of matching the apparent resistivity curve (e.g., Zohdy 1974, Goldberg etal. 1982). The most rigorous solution is that of Oldenburg (1978), who derives an analytic formof Fréchet derivative and treats both model construction and appraisal problem using linearizedinverse theory of Backus and Gilbert. The 2d inversion is carried out mostly in parametric formand employs iterative techniques (e.g., Pelton et al. 1978, Smith and Vozoff 1984, Tripp al. 1984,CHAPTER 1. INTRODUCTION 5McGillivray and Oldenburg 1990). This usually involves, in each iteration, the generation ofthe sensitivity matrix of data with respect to the model parameters and the solution of the matrixequation by an optimization technique to obtain an updated model. Some work is also done toinvert for 2d interfaces (e.g., Lee 1972). This uses the asymptotic expansion of the Id kernelfunction in a 2d environment.In contrast, the work regarding the inversion of 3d d.c. resistivity data is rather limited.Alfano (1959) proposed an approach for a general 3d model based upon a rectangularly griddedmodel and the theory of charge accumulation. That work provides much insight into the natureof the 3d d.c. resistivity problems. Vozoff (1960) applied a similar formulation to simple 3dproblems and proposed an approximate algorithm to invert for a conductivity model consistingof a number of conductivity anomaly blocks in a uniform half-space. The geometry of the blockmust be specified in advance. Petrick et a!. (1988) devised a method which uses the concept ofalpha centres (Stefanescu and Stefanescu 1974) to invert for the central locations of the conductiveanomalies. However, this method is not capable of inverting for a general 3d conductivity model.Recently Ellis and Oldenburg (unpublished but presented at 10th workshop on EM induction,Ensenada 1990) formulated a linearized, iterative 3d algorithm. Park and Van (1991) presented asimilar iterative 3d inversion algorithm which handles a rather limited number of parameters.One of the reasons for the delayed development of 3d algorithms may be the lack of practicaldata for tackling 3d problems. In his study of the d.c. resistivity inverse problems via the integralequation approach, Stevenson (1934) concludes that a unique 3d conductivity solution exists intheory if the surface potential measurements form a complete 3d data set. One example of such adata set consists of the potential measured continuously over a 2d surface for all the current sourcescontinuously positioned along a curve across the surface. This conjecture is proven rigorouslyby Backus (see Vozoff 1960). Although this uniqueness theorem is not directly applicable topractical problems with incomplete and inaccurate data, it does however place a requirement onthe data. The number of dimensions of the data must be at least equivalent to, or no less than, thenumber of dimensions sought for the model. For 3d studies, such data have not been availableCHAPTER 1. INTRODUCTION 6until recent advancements in acquisition techniques. The E-SCAN® * data set seems to be thefirst of such data sets.E-SCAN® d.c. resistivity data are usually acquired on a regular grid over a target area. Thecurrent source is placed at each grid point in turn with the current sink at “infinity”. For eachsource location, the surface potential is measured at the remaining grid points in reference to asecond potential electrode placed at “infinity” (The details of the acquisition procedure will bediscussed in Chapter 6). Thus the data set consists of common source pole-pole potentials overthe grid for many different source locations. It provides the necessary data obtainable in practicefor attempting to construct a 3d conductivity model. It is the availability ofE-SCAN® data whichprompted the research work in this thesis.In principle, standard approaches used in 2d inversion can be applied to 3d problems. However, calculating and inverting a sensitivity matrix for a realistic 3d problem demands far morecomputing power than in the 2d cases. No matter what approach one takes to compute thesensitivity matrix, it always requires that the forward modelling be performed a number of times.For a 3d problem, this is very demanding computationally. The problem is further compoundedif the sensitivity matrix is actually inverted, since it is large and dense. Therefore, alternativemethods are needed to cope with the size of the problem in 3d cases.The final goal of the research in this thesis is to develop efficient techniques to facilitate theinterpretation of 3d d.c. resistivity data such as the B-SCAN® data. To put this in perspective,when the research was initiated in 1987, only 2d algorithms for inverting d.c. resistivity data(e.g., Pelton et al. 1978, Smith and Vozoff 1984) existed. Moreover these algorithms onlyinverted for a small number of parameters involving conductivities and dimensions of blocks.The inversions were formulated as an over-determined problem and generally suffered fromproblems of poor convergence. There was no algorithm for recovering a general 3d conductivitymodel characterized by many cells of unknown conductivity. My research, therefore, began withinterpreting 2d data and this led to the formulation of 2d imaging algorithms. I parameterized* E-SCAN® is the registered trade mark of Premier Geophysics Inc.CHAPTER 1. INTRODUCTION 7the problem with a large number of cells and obtained the solution by minimizing an objectivefunction comprised of the model norm and the data misfit. Thus the type of the model obtained wasdetermined by the norm of the model rather than a priori parameterization. As a natural extensionof the 2d imaging algorithms, I formulated a rapid approximate inversion algorithm for imagingcomplex 3d conductivity structures. It uses the data from fixed pole-pole arrays and solves the3d problem by decomposing it into a sequence of Id inversions in the wavenumber domain. It iscapable of inverting a large number of data to recover conductivities of a large number of cells.This algorithm is then used to construct iterative inversions that generate conductivity modelsreproducing the observed data.In Chapter 2, I present the theory of charge accumulation in the d.c. resistivity experiment.This is a different point of view from the more commonly used approach which focuses uponcurrent flow. In this chapter, different aspects of charge accumulation are reviewed, enlargedwith some new insight obtained through this research and presented in a unified notation. Thisprovides the basis for understanding the fundamentals of d.c. resistivity experiments. Thealgorithms developed later in the thesis all hinge upon the material in this chapter.The direct application of the theory developed in Chapter 2 results in the 2d imaging algorithmusing charge accumulation. This method generates a structural image of the subsurface 2dconductivity by constructing a sparse representation of the charge density accumulated on theboundaries of conductivity prisms. The final image is a composite of many charge density imagesobtained by inverting multiple sets of common source potential data. Thus the algorithm requirespotential data acquired along a traverse over the 2d structure. A variation of this approach is toconstruct a conductivity image directly, utilizing the relationship between the charge accumulationand the surface potential. This results in the second method for imaging 2d structures. Bothalgorithms depend upon an approximate (Born) equation for the surface potential. In the Bornapproximation an estimated background conductivity is assumed and the charge accumulation isgenerated only by the primary electric field and the internal interactions are neglected. TheseCHAPTER 1. INTRODUCTION 8two methods, together with brief discussions on data display and an equivalent source imagingmethod form Chapter 3.In Chapter 4, I present a rapid approximate inversion for the 3d problems. Using the Bornequation, I derive an approximate integral equation for the potential data of a surface pole-polearray. This equation states that the percentage potential anomaly of a surface pole-pole array,relative to a known uniform half-space, is equal to a depth integral of the logarithmic conductivityanomaly convolved with a kernel function in horizontal directions. Fourier transforming theequation, and applying the convolution theorem, decomposes it into a sequence of Id equationsin the wavenumber domain. By such a decomposition, each component of the surface potentialdata is given by a single depth integral of the corresponding component of conductivity anomalymultiplied by a kernel in the wavenumber domain. The number of discrete (complex) data foreach id equation is equal to the number of data maps for different pole-pole arrays. Solving theseId inversions yields the Fourier transform of the conductivity anomaly, and applying an inverseFourier transform produces the 3d conductivity model in the spatial domain. The solution of the3d inversion is, therefore, obtained by solving a sequence of id inversions in the wavenumberdomain. Various aspects of the algorithm are examined in this chapter. These include thepreparation of data, calculation of kernel functions, and the solution of each Id problem. The idinverse problem can be formulated either in continuous or in discrete form. Both approaches arepresented.The above solution for the conductivity model is only approximate and does not provide anacceptable model. Because of the way it is formulated, the algorithm itself cannot be furtherdeveloped to yield a full solution which reproduces the observed data within the statisticaltolerance of the error. However, it can be used as an approximate inverse mapping in AIM(Approximate Inverse Mapping) inversion (Oldenburg and Ellis 1991).The AIM inversion algorithm is a general formalism by which an approximate inversemapping operator can be used in conjunction with an exact forward mapping to successivelyimprove the model and yield a full solution reproducing the observed data. There are two waysCHAPTER 1. INTRODUCTION 9the AIM formalism can be applied. In the first, a perturbation in model space is sought so that theupdated model reproduces the observed data. This is the AIM-MS. In the second, a perturbation indata space is sought so that the updated data will, upon the application of the approximate inversemapping, produce an acceptable model. This is the AIM-DS. In Chapter 5, I construct an AIMinversion method for 3d d.c. resistivity problem using the approximate 3d inversion developedin the preceding chapter. I first present a brief review of the AIM formalism. I then outline theapplication to the d.c. resistivity problem of AIM formalism and the approximate 3d inversion.The AIM algorithm is demonstrated with accurate and noisy data sets from two synthetic models.The first model consists of 5 prismatic anomalies in a uniform half-space while the second oneis made of correlated random conductivity perturbations. The performance of the algorithm isexcellent in all cases.The combined application of the 3d approximate inversion and the AIM formalism results inthe iterative inversion scheme. It has several unique characteristics. First, the method is very fast.The algorithm constructs a 3d conductivity model by solving a sequence of Id inversion problems.Thus the computational effort is drastically reduced compared to solving the 3d problem directly.This approach of decomposing a large 3d problem into a sequence of smaller problems may proveto be advantageous in other geophysical problems. The method of decomposition, of course, willbe problem dependent. Utilizing the AIM formalism, the iterative model construction is veryefficient. At each iteration, only one application of the approximate inverse mapping and theforward mapping is used. The algorithm reduces the data misfit to an acceptable level withina very few iterations. This is in sharp contrast to linearized approach where many forwardmodellings have to be performed in each iteration, either in computing the sensitivity matrix orin search of a proper step-size for the model perturbation.Second, the method is capable of solving problems with a large number of unknown parameters. The storage requirement and the CPU time for the solution increase approximately linearlywith the number of parameters. This is important in practical applications, where even simplegeologies demand a great number of the model parameters to describe the possible variations soCHAPTER 1. INTRODUCTION 10that meaningful information can be drawn from the inversion. The algorithm is also capable ofinverting large data sets. This is closely related to the capability of dealing with large number ofparameters.Third, the method is stable and has less possibilities of generating spurious structures in themodel. By virtue of the decomposition, the solution is sought through a sequence of scales.Structures in the model on different scales horizontally are obtained by inverting the data component on the same scale. The variation is controlled by limiting the smallest scale componentof the data used in the inversion. At the same time, a minimum norm model is constructed in thevertical direction. Thus only the features required by the data appear in the model. Variationsbeyond the data resolution are excluded.In Chapter 6, a set of field data from an epithermal regime is analysed. An epithermal depositis typically formed in an area with well developed fracture and fault systems. It is usually closeto intrusive structures, which provide the hot fluid for the hydrothermal system. The faults andfractures form the pathway for the fluid, which alters the surrounding rocks and deposits theprecious metals. The deposits are often associated with silicified structures with deep “roots”,which usually exhibit strong resistive anomalies. The E-SCAN® experiment is employed to mapthe silicified zones and fault structures in an attempt to predict the localities of possible ore bodies.I present the detailed field procedure of the data acquisition and the format of the data set. I outlinethe different stages of processing applied to the data. First, the data set is corrected for the finitedistance of the “infinite” electrodes based upon an estimated regional conductivity value. Thepole-pole data maps are gathered from the corrected data set. Each map is then smoothed usinga generalized cross-validation (GCV) technique to suppress the outliers in the data. I then applythe 3d AIM-DS inversion to construct a conductivity model. The validity of the features in themodel is discussed, with the emphasis on the algorithmic performance.CHAPTER 2THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTSIn d.c. resistivity experiments, the information about the underground geo-electrical structureis contained in observed potential differences measured at the earth’s surface or at depth. Themeasured potential arises from two sources: (1) the potential due to a current source embeddedin a half-space of conductivity equal to that of the region immediately surrounding the currentelectrode; and (2) the potential that arises from the secondary charge distributions which are setup in the medium by the current input into the ground. The existence and distribution of thesecharges in the earth and on the earth’s surface are fundamental to d.c. potential signals and itfollows that a thorough understanding of charge accumulation should provide enhanced insightinto the forward modelling of d.c. potentials and in formulating new approaches to the inversionof d.c. data.Alpin (1947) was perhaps the first to point out that the source of potential in d.c. experiments is accumulated charge. That seminal paper presents the physical interpretation and basicmathematical description of the electric field in a conductive medium using charge accumulationand it also presents an integral equation for the surface charge density. Alfano (1959, 1960,1961) provides a comprehensive treatment of modelling and interpreting d.c. resistivity datausing charge accumulation. His presentation gives a basic understanding of the field behaviorthat arises in the presence of conductivity discontinuities in an orthogonally gridded 3d environment. Basic physical equations relating to the accumulation of electric charge in a conductivemedium and evaluation of electric potentials from boundary charges are lucidly developed inthe excellent work of Kaufman and Keller (1985, p 11-44). Kaufman (1985) presents a tutorialdissertation regarding the role of charge buildup in time varying problems; that work containsinsight and physical understanding about the d.c. limit which is of interest here. The use of chargedensity in correcting for the effects of topographic distortion is examined by Oppliger (1984).Jiracek (Lecture at 9th workshop on EM induction, Sochi, 1988) investigates the distortion of11CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 12electromagnetic measurements due to the charge accumulation resulting from topography andnear-surface conductivity variation. Finally, integral equations for the d.c. problem have beenderived and used in forward modelling by numerous authors (e.g., Dieter et al. 1969, Snyder1976, Okabe 1984). Integral equation forward modelling approaches have two steps. The firstcomputes the charge density on boundaries of blocks across which there is a conductivity contrast,and the second uses Coulomb’s law to compute the potential at the observer site that arises fromthis charge. This dissection tends to emphasize the importance and the physical role of chargeaccumulation.The goal of this chapter is to assimilate some of the diverse information on the subject intoa single document, enlarge it with new insight, and present results in a uniform notation. It ishoped that this will be of use to others who are involved with interpreting d.c. measurements.The chapter begins with the governing equations and boundary conditions for d.c. potentialsin a conducting medium; particular emphasis is paid to terms which involve charge density. Theapplication of an electric field to a polarizable medium produces a polarization charge. Althoughit is generally well known that this charge does not alter the observed potentials, the precise role ofthe polarization charge is sometimes not appreciated and hence is a source of confusion. I treat thisaspect explicitly. I next look at the relationship between refraction of currents at a boundary andthe change in the electric field caused by the surface charges. This emphasizes the fundamentalrole played by charge accumulation in current channelling problems. An integral equation forcharge density is developed and solved with the aid of the whole-space Green’s function. Thisillustrates how boundary element forward modelling results can be obtained and also suggestshow topographic problems can be handled. The effect of the charges accumulated on the earth’ssurface is studied with the aid of the integral equation and it is shown that formulating the problemin the half-space and in the whole-space are equivalent. Analytic expressions for charge densityresulting from a buried current source in a layer over a half-space is derived using the integralequations. Since the image method is often used to solve simple d.c. problems, I show howthis mathematical approach relates to the physical aspects of charge accumulation. NumericalCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 13examples are presented to quantify the charge density buildup on a simple 2d topography and ona 3d body buried in a half-space. Lastly, the Born equation of the d.c. potential is outlined; thisconstitutes the basis for imaging and inversion techniques in this thesis.2.1 Basic EquationsIn a d.c. resistivity experiment, current is put into the ground and potential differences aremeasured away from the source. The electric potential at any point in the medium is dependentupon the distribution of the conductivity within the earth and that potential may be evaluated byusing Maxwell’s equations, conservation laws, constitutive relations and the boundary conditions.In this section I outline the basic equations and boundary conditions required to understand andto solve the d.c. problem. Equations which involve charge density are emphasized. Additionalmaterial can be found in Stratton (1941), Grant and West (1965), Jeans (1966), Keller andFrischknecht (1966), Telford et al. (1976), Kaufman and Keller (1985, chapter 2), and Ward andHohmann (1988).For a steady state problem, only two of Maxwell’s equations are needed:VxE=Ô, (2.1)V.13=p, (2.2)where E is the electric field, ñ, the electric displacement and Pi’ the volumetric free chargedensity. V is the del operator and is defined in a cartesian coordinate system as V = L+j+k,where 1, j and are the unit vectors in the x-, y-, and z-directions respectively.The law of conservation of charge states thatV.y=—L, (2.3)where 5is the current density due to free charges. A positive current is defined to be related tothe flow of positive charge carriers. In steady state conditions V . = 0 is satisfied everywhereCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 14except at locations of sources or sinks of electric charges. If a current I is injected at a location-0“a,V .j= IS(r— p,). (2.4)In addition to field equations we also require constitutive relations. For the purpose of thisstudy, I assume that the medium is linear and isotropic. Thus13 =E, (2.5)(2.6)where o is the electrical conductivity and is the permittivity.Finally, at an interface separating media of different conductivities, the tangential componentof the electric field and the normal component of current density are continuous. ThusE1 =E2, (2.7)3m =32n, (2.8)where subscripts I and 2 refer to the respective media and subscripts t and n refer to tangentialand normal components. The unit normal vector i is chosen to point outward from medium-I atthe interface (see Fig. 2.1). The normal components of i3 and E satisfy— D1 Tf, (2.9)(2.10)COwhere T1 and Tt are respectively the surface densities of free and total charge. Thus the normalcomponents of 13 and . can be discontinuous if there is a surface charge distribution on theboundary. Using equation (2.1) and the vector identity V x V4 = 0, we can express the electric4CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS12 =Figure 2.1 Basic geometry for current flow at a boundary.field as the gradient of a scalar potential ,E=-Vci.15(2.11)The electric field is bounded away from sources. Consequently, the potential is continuous andwe have the boundary condition qfi = çb2.Any steady state conduction problem can be solved using the above equations and boundaryconditions. In the most general approach, substitution of (2.6) and (2.11) into (2.4) yieldsV (oVç5) = —IS(— it), (2.12)2 VT.Vqf I -Vq=— ——S(r—r5).0•(2.13)Equation (2.13) is recognized formally as Poisson’s equation. The two terms on the righthand side have units of p160 and therefore each can be thought of a charge density. For thepurpose of calculating the potential, the term involving the source current can be replaced by ann2= a1 E1which can be expanded to produceCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 16effective point charge of magnitude Q= 6oI/o, where o, is the conductivity at the location ofthe current source. The first term on the right hand side is of prime interest because it expressesthe charge buildup that exists as a result of changes in the electrical conductivity structure. It isnon-zero whenever there is a component of the electric field parallel to the conductivity gradient.Under these conditions there will be a physical buildup of electric charge with volumetric densityPt= VuVçf (2.14)In the limit that the conductivity gradient approaches infinity, i.e., when the medium suffers adiscontinuous change in conductivity, the volumetric charge density in (2.14) becomes a surfacecharge density confined to the boundary separating the two regions. It is this charge which createsa discontinuity in the normal components of D and E as in equations (2.9) and (2.10).By combining (2.7) and (2.10) and using Ohm’s law, the charge density can be written as= (1 — 2)E. (2.15)60Equations (2.14) and (2.15) show that the accumulated charge is negative when current flowsfrom a resistive into a conductive region. Conversely, positive charge accumulates when a currentflows from a conductive to a resistive region. This rule delineating the signs of the accumulatedcharge is very useful for predicting the general character of d.c. fields resulting from simplegeological structures.Equation (2.13) indicates that conductivity is the governing parameter for a stationary electricfield and that permittivity does not play a part. Yet, application of an electric field to a polarizablematerial creates polarization charges which produce an electric potential. This apparent contradiction is sometimes a source of confusion and I therefore address it here. It is best illustrated byexamining what is happening on the boundary between two media.An electric field applied to a polarizable medium generates an electric polarization or dipolemoment per unit volume. The polarization is proportional to the applied total electric field, andCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 17is given byP=XEoE, (2.16)where P is the polarization vector and x = (/o) — I is the electric susceptibility. At the boundaryof the medium there will be an polarization charge = P , where ii. is the normal vector. Thenet polarization charge density at a boundary separating media with permittivities e and 2 is= —(P2 — P1). (2.17)and the total charge is the sum of this polarization charge plus the free charge. Substituting (2.5)into (2.9), and (2.16) into (2.17) and summing yieldsTt =Tf + Th=€0(E2— E1).This result is precisely the same equation as (2.10). It is this total charge density on the boundarywhich is responsible for the continuity of the normal component of the current density and thischarge is not affected by variations in permittivity. The permittivity however, does determinehow much free charge has to be accumulated at the boundary so that boundary conditions aresatisfied.In general, when the physical properties are changing continuously, a volumetric chargedensity ofPt.. VcT.EEOis set up. Again, this charge density is determined only by the conductivity, but if the medium ispolarizable, then Pt can be explicitly formed as the sum of free and polarization charge densities.The expressions for these arepf=——Vo.E+E.V, (2.18)Vcr•E -p&=(e—eo) —E.V€. (2.19)CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 18The sum of P1 and Pb gives the total charge density which determines the final electric field.In summary, in a d.c. experiment, electric charges accumulate whenever there is a gradientof conductivity and a non-zero component of electric field parallel to it. The final electric fieldis produced by the primary source and by surface and volumetric charge distributions. Whenthe medium is polarizable, both free and polarization charges contribute to the total accumulatedcharge. However, the total accumulated charge is controlled only by electrical conductivity; thepolarization plays no part other than to supply a portion of the charges needed to satisfy theboundary conditions.2.2 Current Flow in d.c. Resistivity ProblemsIn this chapter I focus upon the accumulated electric charges in order to provide insightto the d.c. resistivity problem. A commonly used description of d.c. fields, however, isthrough the concept of current flow. The variations in conductivity structure alter the flow ofelectric charges and the final distribution of current is such that the energy loss due to ohmicdissipation is minimized. Physically, this results in current being channelled into the regionsof high conductivity and deflected away from resistive regions. But knowledge of the currentdistribution, or even its direction, does not imply immediate knowledge of the potentials. Inparticular, current direction at any point is determined by the gradient of the potential. Currentmagnitude also requires specification of an electrical conductivity.The relationship between current flow and charge accumulation is best illustrated by considering the refraction, or change in direction, of currents impinging upon a plane interface separatingtwo media. This refraction is a direct effect of surface charges which accumulate on the interface.This is quantified here but a similar discussion can be found in Kaufman and Keller (1985).Consider, as in Fig. 2.2(a), a point on the boundary and let Eb be the base electric field which isproduced by all sources away from this point. Continuity of the normal component of the currentdensity requires the existence of a charge density Tt (given by equation (2.10)). This chargeCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 19\ts22 02_______ _______CE) CE) CE) CE) CE) E)01_3J_,’\Is1 fl\0 \ 01<02(a) (b)Figure 2.2 (a) Electric charges accumulate on the boundary separating two conductivemedia and produce perturbation fields E,1 and E32 normal to the boundary. Thesefields are added to the base field Eb to produce a total field whose direction changes atthe boundary. The resultant current flow shown in (b) is equivalently refracted at theboundary.affects the normal components of the electric field so that-,E1 =Eb 71 — —,—(2.20)E2 =Eb fl + —.2e0An expression for the total accumulated charge, obtained by combining equation (2.20) and (2.8)using Ohm’s law, is=2°’ °2Eb•ñ. (2.21)O O+OThe normal components of the electric field in (2.20) are then2o2 -.= Eb 71,01+02 (2.22)2o —Eb•n.1 + 2The tangential component of the electric filed at the boundary is unaltered by the charge densityand is equal to Ebt. Using equation (2.22), we obtain‘i + 2 Ebttan81 =O2 ‘bnCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 20cr1 + 2tano2=—,2a Eand therefore_!_tan6i = !tan92 (2.23)This shows that a current line refracts as it crosses an interface; it bends towards the normalwhen entering a resistive medium and, conversely, it bends away from normal when entering aconductive medium. Equation (2.23) is identical to the usual refraction formula (e.g., Keller andFrischknecht 1966, Telford et at. 1976) which is derived directly from the continuity conditionsof the electric field and current density. The derivation here shows that the charge buildup onthe boundary causes a change in the normal component of the electric field so that the directionof current flow is altered as it passes into a medium with different conductivity. In general, theelectric field at any point within the medium is a vectorial sum of the primary field and the fieldproduced by the accumulated charge. It is this additional field which changes the direction ofcurrent flow in the medium and results in the channelling of the current into conductive regionsand the deflection away from resistive regions.2.3 An Integral Equation for the Charge DensityCharge accumulation plays a fundamental role in d.c. resistivity problems and we thereforeneed a technique to explicitly quantify that role and to illustrate how it can be applied in practice.An integral equation is a desirable tool because it enables us to examine quantitatively the chargethat is accumulated and it also provides an effective way to carry out the forward modelling of d.c.responses. This latter aspect has been studied by many authors (e.g., Dieter et a!. 1969, Snyder1976, Okabe 1984). In the following an integral equation is developed for charge densities whichclosely follows that presented by Snyder (1976).I begin with the initial differential equation (2.13) and Green’s second identity,Jj’J(V2G— GV2çf)dv =— G)ds, (2.24)CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 21where S is the surface bounding the volume V. Green’s identity is valid for any functions 4,and G that are continuous and have derivatives up to second order (Morse and Feshbach 1953).I choose 4, to be the potential function, satisfying (2.13), and G to be the whole-space Green’sfunctionG(,i’)= I’which satisfiesV2G(j = —4irS(— p’).Let V be the whole space. Then within V, both 4’ and C decrease as inverse distance from thesource and 94,/ôn and ÔG/ôn decrease at least as inverse distance squared. Therefore, if S istaken as the surface of a sphere with radius approaching infinity, the right hand side of (2.24)vanishes. Substituting (2.13) into (2.24) yields= G(,) + 1JJJVT(i’) .V4,(i’)G(?1)dV (2.25)4iro, 4irVwhere o is the conductivity at the current source location. The first term in (2.25) is recognizedas the potential due to the point source in a uniform background of c.onductivity o. The secondterm is the potential due to the accumulated charge distribution.Although equation (2.25) can be used directly, it is simpler to assume that the conductivitystructure is piecewise constant. Then Vo is zero everywhere except at boundaries betweenregions of different conductivities. The volume integral reduces to a set of surface integrals andthe integrands contain the surface charge density.Consider the geo-electrical model shown in Fig. 2.3. The background medium has a conductivity o. There are ii + I inhomogeneities embedded in this background; each has a constantconductivity Oj and boundary F1. The 0th inhomogeneity represents the volume of air above theearth and is bounded by the surface F0. All boundaries are assumed to be piecewise smooth andrnraEo)4Nasu r)•VG(,)+= ni0, — o eo 4ircr.,/7 r’L VG(’)ds,Coi=o j=O,...,n. (2.28)Here i E F, ñ = V operates on the field point and the integrals operate on the secondaryCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 22Figure 2.3 The geo-electrical model consisting of n + I bodies in a uniform background.n bodies are beneath topography and the air above constitutes the 0th body.ñ is the outward pointing normal vector. The point current source is located in the backgroundmedium.For this model, equation (2.25) becomesç(i) = G(r,r) + ff ri(r’)C(r.l)ds (2.26)4Tro 4ir._JJ Cwhere r2 is the charge density on the Lh boundary. Using equation (2.15) and writing the normalcomponent of the electric field in the background as E = —. V, I have= U9 — Vq!. (2.27)Co UJEquation (2.26) is valid everywhere, so I can substitute (2.26) into (2.27) to eliminate 4 andobtain the integral equation for ‘i-i:source points F’.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 23The integral over r’, is not straightforward to compute because the integrand in equation(2.28) is infinite when i’ —* . However, this singularity is removable. A circle of radius Scentred at divides F3 into two parts: let F denote the area inside the circle and let F be the areaoutside. This is always possible if S is sufficiently small since I assume that the boundaries areall piecewise smooth. As S —* 0, the charge density in T can be taken as constant and therefore,lim jJ ri(r’)ñ..VG(.r.I)ds = 6016Thus (2.28) becomes= IkVG(i)+CO 2iru3E2ir 60 (2.29)k3 ffr(’) ——n3.VG(r,r )ds,2rjj€owhere k3 = (o — o-)/(o, + oj) and is the j’ boundary with a small area around pointexcluded. All integrals in (2.29) are now proper.Equation (2.29) is a Fredholm equation of the second kind and is solved using standardtechniques to generate the charge density on all boundaries. The potential is then evaluated bysubstituting these charge densities into (2.26).2.4 Effect of Charge Density at the Earth’s SurfaceThe formulation thus far has been to solve for the potential in a whole-space. This generalityis advantageous because we ultimately wish to solve for the potential in the presence of surfacetopography. With the whole-space formulation the undulating surface may be considered to be theboundary of another body and the surface charge due to a point source can be evaluated directly.Most studies assume a half-space model with a flat upper surface. Under such an assumption,the charge density at the earth’s surface does not appear in the formulation and is equivalentlyCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 24taken into account by a zero flux boundary condition. Although an integral equation whichincorporates the zero vertical flux of current at the earth’s surface can be derived (Snyder 1976),further insight regarding the fundamental nature of the surface charges is obtained by usingequation (2.29) to compute the charge density on the surface of a flat earth and then to computeits secondary potential and its effect on the charges of buried bodies. I have carried out thisstudy and shown that the effect of the charge distribution at the earth’s surface can be completelyreproduced by a distribution of fictitious images.Adopt a right handed Cartesian coordinate system with z positive downward and z = 0denoting the earth’s surface. The coordinates of the current source are (x,, y, z,). From equation(2.29), the charge density on the earth’s surface isro(r) I+27O ((x— a,)2 + (y — y8)2 +‘E ffri(’) ds. (2.30)2W•JJ (( — + (y — yI)2 +The integral term over the earth’s surface (T0) does not appear in (2.30) because VG(i, r’) isperpendicular to the surface normal. The first term in (2.30) is the portion of the secondary chargeproduced by the primary field at the earth’s surface, which is equal to the surface charge densityin the absence of any underground inhomogeneities. The terms in the summation represent theportion due to the existence of the buried inhomogeneities as a result of interaction among thecharges on all interfaces.Denote the first part by i-i,rg(i) = I za1/2 (2.31)€27fl7.g ((x — x)2 + (y—y3)2 +We see that r(x, y) attains its maximum at a point directly over top of the current source andthat its amplitude falls off as.As the current source moves closer to the surface, theCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 25charge density becomes more concentrated. In the limit,1 zaurn— 3/2 = S(a — x,)S(y —z—,0 2,r ((a — z)2 + (y— y)2 + z)and thereforeurnrg(x,y)=I8(— —y,). (2.32)z,—’O 60Insertion of this charge density into equation (2.26) yields, in the absence of other inhomogeneities,the total potential= -•. (2.33)2iro r—This potential is double the whole-space primary potential and shows that, for a current sourceat the surface of a flat earth, the surface charge density provides an additional potential whichis of equal strength as the whole-space primary potential. In problems formulated explicitly ashalf-space problems, it is convenient to regard (2.33) as the primary potential.The second part in (2.30), denoted by rg, is given by= ff r’)— x’)2 + (y— y’)2 + Z12)3/2(2.34)Integrating (2.34) over the surface to obtain the total charge, Q3, yields,= 2ir j7 ri(r)(JJ ((z — i)2 + (ui— yi)2 + Z12)3/2The integral in the bracket is recognized as the solid angle subtended by the earth’s surface ata point r’ underground and equal to 2ir. Substituting in and carrying out the integrals over Tyields,Qa = Q, (2.35)where Q is the net secondary charge on th inhomogeneity. Therefore, the sum of the secondarycharge on the surface due to the buried inhomogeneities is equal to the sum of the net secondaryCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 26charges underground. A similar conclusion was arrived by Alfano (1959) through a differentapproach.Using Coulomb’s law to compute the secondary potential produced by r) yields,= JJ r(f”)G(, r”)dz”dy”.Exchanging the integral over To and T2 and applying the integral identity (Stevenson 1934)0000I’ I’ zI III dxdyJJ ((x1 — z)2 + (y— y)2 +z12)3”((z2 — x)2 + (Y2—y)2 +2ir1/2’ forz1 >0 (2.36)((x1 — 2,2)2 + (y1 — y2)2 + (z1+ I z 1)2)yields= I ff r’) ((x — i)2 + y — y’)2 + 1z1 +z,)2)h,d8. (2.37)Above the earth’s surface where z <0, this is equal to the potential produced by the subsurfacecharges. Beneath the earth’s surface where z > 0, A4 is equal to the potential produced bythe subsurface charges residing at their image positions above the surface. It follows that at anypoint in space the effect of this part of the surface charge is completely represented by a set ofimage bodies coinciding with the buried inhomogeneities or at their mirror positions about theearth’s surface, depending upon whether the field point in consideration is above or beneath thesurface. This theoretically justifies the common practice in which the actual secondary potentialon a flat earth’s surface is obtained by doubling the secondary potential computed only from thesubsurface charge densities.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIWTY EXPERIMENTS 27Next, writing equation (2.29) to isolate ro(f), and specifying a flat earth’s surface, yieldsrj() = 1k2.VG(,i)÷ 4JJri()..vG(l)dS+e0 2ir 27r._F’(2.38),ffTi(7;’)ñ..VG(;*i;.I)dS +- JJro(’),• VG(,’)dx’dy’.—00—00Let /3 be the last term in the above equation. Substituting (2.30) into /3 and applying identity(2.36) yields=- ff Tn.vG(,i1t)dS, (2.39)where and r” are, respectively, the images of ii., and i’ about the surface plane. The first termin (2.39) is recognized as the effect of a primary source at the image location i and the secondterm is the effect of image bodies above the earth’s surface, as has been demonstrated. The finalexpression, obtained by combining (2.38) and (2.39), isr3() 1k3 / 1 1 \= + ,6o 2ro9 r — r —nIIr(rL_____—, ,,7zj.V/ ds+2ii-’---’JJ Ir—r’I:=1Wi ‘k- ft r-(’) 1 (2.40)-L ds+2irJJ c Ir—r’IIc, ffr(i’),.__ _n,.V . ds.27r. JJ Co1.This is the same equation as that of Snyder (1976), which is derived through formulating thepotential problem in a half-space and using the half-space Green’s function.The whole-space and half-space Green’s function approaches to calculating the charge densities are essentially equivalent. They differ in that (2.29) is slightly more general and can be usedCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 28to treat topographic problems whereas (2.40) is specifically suited to plane surface problems andis numerically more efficient in such cases.2.5 An Analytic Example of Charge AccumulationGenerally the solution of (2.40) demands numerical techniques but there are a few simplegeometries in which analytic solutions for the charge density are available, i.e., a point currentsource buried in a half-space, and a point source in a layer which overlies a half-space. Theformer is discussed in the preceding section within the context of the effect of surface charges. Ipresent the analysis for the latter example here. These examples provide enhanced insight into themanner in which charges are distributed on various surfaces and they also provide a foundation forunderstanding the relationship between the solution to potential problems via charge distributionand via the image method.I evaluate the charge distributions when a point source is located at depth d inside a layer ofthickness h overlying a half-space (see Fig. 2.4). I adopt the same coordinate system as before.To be consistent with the previous section, I take the layer in which the current source is locatedas the background medium and choose normals as shown in the Fig. 2.4. At the boundaries T0and F1,k0 =1,i 2T1 +if0I.(x,yd)a1 +1riFigure 2.4 A point source in a layer over a half-space. The normal vector of the twoboundaries are defined as pointing into background medium o.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 29To find the charge densities r0 andr1 on these surfaces, I employ (2.29) in an iterative manner.Equation (2.29) is first written as1k — , i = I when j = 0= n.VG(r,r)+e0 2iro(l—1).’ i = 0 when j = 1 (2.41)T121rjj o 1=1Tiwhere 1 denotes the iteration number. When 1 = 0, the second term on the right hand side vanishesand the estimated charge densities are those due to the primary potential, i.e.,de0 2ircri (12 +d2)3/’— 1k d (2.42)e0 2Toi (,2 + (h — d)2)31’2where , = ((x — x)2 + (y — y)2)2 is the radial distance from the source point.The presence of the primary charge distribution on F’0 affects the charge on F’1 (and viceversa); the next iteration attempts to account for this. Employing the integral identity (2.36)again, I obtainT(?1) r°() + 1k 2h — d2’ (2.43)60 o 2iro1 (,i2 + (2h— d)2) /r”(,i) ‘-°(‘i) + 1k h + d (2 44)60 6o 2iro1 (,2 + (h + d)2)3”2Proceeding with this iterative process, I obtain the final solutions for charge density as seriessummations:TO(11) I d +eo 2ir (2 + j2)3/2I °° I 2nh—d 2nh+dk’ I + I, (2.45)2iro1 (2 + (2nh — d)2))3 2 (2 + (2nh + d)2)3 2Ti(?)Ikj h—d +60 2TUj (,,2+(h_d)2)3/CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIWTY EXPERIMENTS 30I ,,( 2nh—(h—d) 2nh+h—dk1 +k1. (2.46)271t7 + (2nh — (Ii — d))2)3/2 (2 + (2nh + h — d)2)31’2The charge densities are observed to have maximum magnitude at = 0, i.e., just above and belowthe current electrode. Away from the central point the charge density decreases in magnitude butthe rate of decrease depends upon the magnitude of lci. As k1 I becomes larger the charges arespread further over the interface. It is clear from the two formulae that the interaction betweenthe charges on the upper and lower boundaries is primarily controlled by the magnitude of k1 andthe thickness of the layer.As the current source approaches the surface, i.e., as d —* 0, a delta-like charge distributionwill appear on the surface coinciding with the current source. In addition, there is a charge density00I 2nh= > k (2 ÷ (2nh))3/’which is due to the interaction with the bottom of the layer. This charge distribution will bepositive if the layer is more conductive than the underlying half-space and it will be negative,otherwise.Examining equations (2.45) and (2.46), we notice that every term has the form of equation(2.31); i.e., each term is like the charge density on a single plane interface induced by a pointsource. Furthermore, applying equation (2.36) to calculate the potential due to each of theseterms results in a potential which is the same as that produced by a point charge away from theinterface. This leads to an understanding of the physical basis for image method solution ofpotential problems.2.6 Image Method Explained by Charge AccumulationConsider the example of a point source buried in the lower half-space of conductivity cTJoverlain by an upper half-space of conductivity ao. Let the source be located at (x8,y,, d) in thesame coordinate system as before. The charge density along the interface, obtained by evaluatingCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 31equation (2.29), isr(x,y) 1k d (247)2iro1 ((a — x,)2 + (y — y.)2 + d2)3”2where k = (u1 — u0) / (01 + oO). As o —* 0 we get to the special case of a half-space whosecharge density is given by equation (2.31).The potential anywhere in the whole-space is the sum of the primary potential from thecurrent source and the secondary potential due to the charges on the interface. Evaluating (2.25)and using identity (2.36) yieldsIl 1k 1çf(x,y,z)= —+ 1/24ira1R 4iro ((z_xs)2+(y_ys)2+(d+IzI)2)where R = ((x — a)2 + (y—+ (z — d)2)1”2.Thus the potential in the lower medium wherez >0, is=1 (1+-), (2.48)4iru1 R Rwhere R’ = ((x — x3)2 + (y—y)2 + (z +d)2)”. This potential is equivalent to that arising from acurrent of strength I located at R and another current of strength kI located at R’. In the uppermedium where z <0, the potential may be written asI l+k(2.49)47ru0Rwith R” = ((x — x)2 + (y — y)2 + (z — d)2)W and t = 2o/ (oi + u0). The second form of thepotential in equation (2.49) is equivalent to that produced by a current of strength tI at a distance1?” in a whole-space of conductivity o. Equations (2.48) and (2.49) are identical to those derivedby the image method in Keller and Frischknecht (1966) and Telford et at. (1976) where ananalogy is drawn between electric current flow and geometrical optics based on the fact that theintensity of both current density and light emanating from a point source varies as the inversesquare of the distance. The plane interface between two media is viewed as a semi-transparentCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 32mirror. Based on this analogy, t and k defined above are referred as transmission and reflectioncoefficients, respectively. Although the correct equations are derived using the geometrical opticsformula, the physical understanding is not present. The derivation here has explicitly shown thatthe image method works because it was possible to find a fictitious point charge which producedthe same potential as the true distribution of charges on the boundary.Although the whole-space single interface example establishes the basic principles of theimage method, a slightly more complicated model is required so that internal interactions can bemodelled. Correspondingly, I consider the structure in Fig. 2.4 with the current source moved tothe surface. From (2.32), (2.45) and (2.46), the charge densities on the upper and lower surfacesare=--sx x8)8(y— !I)+ i 2nh3/2’ (2.50)CO 0•1n=i(,2 + (2nh)2)and_____(2n+l)h. (2.51)60 71O ((p2 + ((2n + 1)h)2)3/2where k = (u1— 02)/(01 + c72). By evaluating (2.29), the potential at any point in the space isI__+27ro1RI ic’327ru 2 ((2+(2nh+ I z 1)2)1/2 (2.52)I____________________________2iroi ((j2 + ((2n + 1)h+ z — h I)2)12where = ((z— x3)2 + (y — y)2)h/2. Equation (2.52) is identical to that derived from imagemethods. For that approach one needs to compute potentials due to two families of image sources.(see Keller and Frischknecht 1966, pp.108-11)To summarize the statements concerning the relationship between the charge density and theimage solution I return to equation (2.51). r represents the charge that is accumulated on theCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 33interface between the layer and the underlying half-space. The potential due to this charge isobtained by integrating with the whole-space Green’s function, i.e.,=——- If r’) 1 ds (2.53)4ir€ojj Ir—r’ITIat any position in the medium. The image method reproduces the same potential through aseries solution (the second summation in equation (2.52)). Each term in the series is recognizedas the potential due to a point charge located in the half-space and beneath the primary sourcelocation. Thus whereas Aç in (2.53) is obtained by the integration of a real charge distributionon the boundary, the series summation produces the same result by employing a set of fictitiousimage sources. Therefore the essence of the image method solution to the potential problemis to derive a set of fictitious sources which produce the same potential as does the true chargedistribution. It follows that the image method is viable only when the conductivity structure issuch that the effect of the accumulated charge can be represented by a set of point images.The image method is not applicable to all problems. In general, potentials from all Id conductivity functions can be found using images (Kunetz 1972, Szaraniec 1976, Levy et a!. 1988).Applicability of the image method to 2d or 3d problems is more restricted. Alfano (1959) considered conductivity structures composed of rectangular prisms of constant conductivity bounded bytwo orthogonal plane interfaces. Ifo1, 02, O3 U4 denote the conductivities in successive quadrantsdefining the prisms, Alfano proved that the image solution exists if and only if(2.54)Physically, this condition ensures that the charge density on each interface is continuous at theintersection point. When continuity of physical charge density exists then the potentials can berepresented as arising from a set of fictitious image charges. However when (2.54) is not satisfied,the charge density is discontinuous at the intersection point, no representation of potentials bypoint charges is available, and the image method is not applicable.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 34This understanding, plus the previous derivations, allows me to make general summarystatements regarding the image solution. Steady state current flow in the earth demands thephysical existence of charges which reside on the surfaces of buried bodies and on the earth’ssurface. The secondary potential at any point in the medium is obtained by integrating thesecharges with the whole space Green’s function. The secondary potential zqS arising from a chargedensity r on a boundary F is given by=—-— if r(r’)G(r,f’)ds. (2.55)4irqrSince r is a surface function it can be expanded in terms of an arbitrary set of basis asr() = ajj(i), (2.56)where a1’s are expansion coefficients. Substituting (2.56) into (2.55) yields=a, if1(’)G(, jds. (2.57)If the basis functions are chosen so that the integrals can be computed analytically, then (2.57)may be efficiently computed. The image method as presented in the literature is equivalent tochoosing the basis functions as a charge distribution on a plane interface arising from a pointcurrent source (see equation (2.47)). Suppose that the boundary F is the plane z = z and thatthe centre of mass of r(z, y, z) occurs at the point (, y, zr). Then a possible expansion of thecharge density is00 1r(x, y, z) = a,:= ((a—a) +(YYc) +(zj—z))where specification of z’s determines the basis function and the a’s are the strengths. Thepotential produced from this charge density is100=24ire0 ((a — Xc)2+(Y y)2+(I z — + z —z 1)2) 1CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 35The difficulty with this method is in choosing the basis functions, i.e., specifying z’s, andin evaluating their strengths, a. Ideally, the basis functions should form a complete set so thatevery surface distribution of charge can be represented. It is noted that if the basis functions arecontinuous and not complete it is impossible to construct an exact representation of a discontinuouscharge distribution. Therefore, an image method solution exists only for those structures withcontinuous charge distributions. It follows that it is not possible to have an image solution for aconductivity structure that violates condition (2.54).2.7 Numerical Examples of Charge DistributionFor an arbitrarily shaped conductivity structure, the charge density on all interfaces is foundby solving the integral equation (2.29). There are many ways to solve this Fredhoim integralequation but most methods center around first dividing all interfaces into a set of small arealsegments known as elements. A representation for the charge density on each element is thenassumed. Different choices are possible: the charge density may be assumed to be concentratedas a point charge located at the centre of gravity of the element; it may be assumed constant orlinearly varying over the element; or it may be represented by some higher order polynomial.Substitution of the representation for the charge density into the integral equation yields a systemof linear equation to be solved (e.g., Alfano 1959, Harrington 1968, Pratt 1972, Snyder 1976,Eskola 1979, Brebbia and Walker 1980, Okabe 1984, Oppliger 1984, Shulz 1985, Das andParasnis 1987).As an example, consider a pure topographic problem in which a single interface F overliesa homogeneous medium of conductivity 00. The surface F is divided into M elements and aconstant i- is assigned to the th element. Substitution into (2.29) yields= 1k. ôG(r,,r) + JJ ôG(i’)dS i = 1,... (2.58)CO 2iru0 on1 2ir . oniri34where is the center of gravity of the ith element, and r = r(r,) is the charge density of thiselement. The integral over ith element vanishes after the singularity is removed.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 36By reordering, the equations in (2.58) may be written in matrix form:(2.59)where A is an M x M coefficient matrix with components:=1,A1... JJ ôG(’ii,r’)d (2.60)is the M x I vector of the unknown charge densities and 0) has components:— 1k. 8G() (2 61)— 2ircro 8The linear system of equations can be solved directly by decomposing the matrix A, or it can besolved through iterative techniques.In (2.59)1 have denoted the right hand side by 0)• The reason for this is that 0) is a first orderapproximation to the charge density that I am seeking. It is the charge density that would exist ifthe potential at the boundary were equal to the primary potential. This approximation sometimesyields fairly good estimates of potential anomaly and delineates first order characteristics (Kellerand Frischknecht, 1966). The computational ease with which can be evaluated makes thisapproximation very desirable.To illustrate quantitatively the charge accumulation in more complicated geo-electric structures, I consider two examples. In the first, I compute the charge distribution on a surfacetopography overlying a uniform earth. The topographic model and the charge distribution arisingwhen a current electrode is located in the middle of the valley are shown in Fig. 2.5(a) andFig. 2.5(b), respectively. A negative charge is observed on the valley walls but everywhere elsethe charge is positive. The maximum negative charge densities occur on the wall immediatelybeside the current source. Along the strike, the charge density decays with the distance away fromCHAPTER 2. THEORETiCAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 37air (o=O)I -/ AV 50m—225. —200. —175. —150. —125. —100. —75.0 —50.0 —25.060.l2.180.240.-300.300. 500.1.02 1.04 1.06 1.08 1.10 1.12 1.14 1.16 1.18Figure 2.5 The charge density on a 2d topographic surface. The geometry of the cross-section is shown in (a). The charge density on the 2d topographic surface is shown in (b)(coulombs/m2scaled by eo). The y coordinate delineates offsets in the strike direction.Because of symmetry only charge densities for positive value of y are shown. Thesurface anomalous potential (mV) is in (c) and the surface apparent conductivity (mS/rn)is displayed in (d).-150-9 0 9150X(m)Oi=O.001 S/rn a—7.00 —6.00 —5.00 —4.00 —3.00 —2.00 —1.00 0. 1.00-500. -300. -100. 180.X (m)300. 500.S3--100. 100,X (m)-100. 100.X (m)CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 38the source. The charge distribution at the bottom of the valley is primarily a result of interactioneffects with the negative charge on the valley walls.This physical charge residing on the earth’s surface produces an anomalous potential whichis displayed in Fig. 2.5(c). When that potential is added to the primary potential and convertedto apparent conductivity using the formula cr0(r) = I/2irrc(r), I obtain the distorted apparentconductivity map shown in Fig. 2.5(d). In that figure the apparent conductivities are plotted atthe locations where the potential is observed. It is noted that the apparent conductivity in andnear the topographic depression is greater than the intrinsic half-space conductivity.This example shows that the physical charge caused by topographic variations gives rise to aperturbation potential that will distort the signal of buried bodies. However, the insight affordedby this example can also be used in reverse. A correction for topographic distortions can be madeby first estimating the accumulated surface charge and then subtracting the resultant potentialfrom the field observations. This approach has been used with success by Oppliger (1984).The second example is selected to illustrate the charge accumulation that exists on a conductiveprism buried in a uniform half-space. The geometry of the structure, location of the current source,the charge densities on the prism, the anomalous surface potential, and the surface apparentconductivity are shown in Fig. 2.6. Negative charges exist on the top and left faces of the cubeand this is in accordance with current flowing into the prism there. Positive charges are seen onall other faces. Near the upper left corner of the front face (and, by symmetry, on the back face),there is a small area with negative charge indicating that current is flowing into the prism. Thisis a vivid demonstration of current channelling. In terms of charge accumulation, that negativecharge is a consequence of interaction with the strong negative charges on the top and left facesof the cube.The secondary potential produced by the charge distribution on all faces of the cube is shownin Fig. 2.7(a). It is similar to a dipole field with a negative region centred above the prismedge near the current source and a positive region at the other side of the prism. The apparentconductivity map is shown in Fig. 2.7(b). As in Fig. 2.5, apparent conductivities are plotted at theCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 390.—15.0$ 480m•,—30.0_____25 m_____Ti_____ _ ____ _______100 ml 0.01 S/rn ‘ Iwo: 5.0A 1___I [i. —60.0LI____Oi.0OlS/mz y—105.a —120.17.5 17.516.0 16.014.5 14.513.0 13.011.5 11.510.0 10.08.50 8.507.00 7.005.50 5.50S-50.—30.—10.10.30.50.—50. —30. —10. 10. 30. 52.X (m)SN—50,-30.-12,25.45.65.85.105.125,—50.—30. —10. 10.Y (m)10.30.30. 50.50.25.45.65.25.—50. —30. —10. 10. 30. 50.X (m)85.105.125,—50. —30. —10. 12. 30. 50.Y (m)0. 15.0—15.0 45. 12.0—30.0 9.00—45.0 65. 6.00—60.0 ... 3.00N—75.0 85. 0.—90.0—3.00—105. 105.—6.00—120.—9.00125,—50. -30, —10. 10. 30. 50.X (m)Figure 2.6 The charge density on a buried prism. The geometry and conductivities aregiven in (a). (b)-(f) respectively shows the charge densities on the top, right, bottom, left,and front faces of the prism. Units are in coulombs/m2scaled by e0.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 4050.0 1.250. 1.20—50.0 1.15—100. 1.10—150. 1.05—200. 1.00—250 0.950240. —300. 0.900—350. 0.850400.-400. —240. —80. 80. 240. 400.X (m)Figure 2.7 (a) shows the secondary potential (mV) at the earth’s surface computed fromthe charge densities shown in Fig. 2.6. (b) is the corresponding apparent conductivitymap (mS/rn).locations of the potential electrodes. A conductivity high is centered near the current source butthere is a conductivity low at large source-observer distances to the right. It is perhaps contraryto intuition that the conductive prism has produced a resistive anomaly. This is understandable,however, after a somewhat more critical evaluation of the charge distribution. When the potentialelectrode is to the right of the prism it is close to the positive charge on those faces. Even thoughthe strength of the positive charge is less than the negative charge on the other side, the inversedistance decay in the Coulomb potential allows the secondary potential to be a positive quantityas shown in Fig. 2.7(a). Therefore, a conductive body can produce either a conductive or resistiveanomaly depending upon the relative position of potential measurement with respect to chargedistribution.2.8 Born Equation for d.c. ResponsesAs mentioned in the preceding section, the potential calculated using only the charge accumulation from the primary field may yield a reasonable approximation to the true secondarypotential. This is in fact the Born approximation. In this section, I explicitly outline the Bornequation for the d.c. potential in a 2d and a 3d environment. This is to establish the basis for theimaging and inversion algorithm which will be presented in the following chapters.—400.—240.—80.80.—400.—240.-80.CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIWTY EXPERIMENTS 41Let o(r) be the conductivity structure in a lower half-space V whose upper surface is flat.Following sections 2.4 and 2.5, the potential on the surface resulting from a surface point sourceof strength I is given by= + If/fVU (1’) V4 (“) dv, (2.62)27rr1.—rjo, 2r o(r) Ir—riVwhere i, and are the source and observation point respectively, and o, = cr(r,). I choosea Cartesian coordinate system with origin at the surface and z positive down. The integral in(2.62) represents the secondary potential which contains all the information about conductivityanomalies and is denoted by= 1//fVcT‘ c”.? dv. (2.63)2ir o(r) Ir —riVLet the conductivity be represented by cr(i) = o-op(). Here o-o is the conductivity of a uniformbackground and p() is a dimensionless function of spatial position i. Substituting into (2.63)yields= I f//Vdv. (2.64)The potential 4(r) is not known, but if the deviation of the conductivity from the backgroundis small over the entire model and if the surface conductivity is equal to 00 (thus o = 0-0), I canapply the Born approximation by replacing 4(f) in (2.64) with the half-space primary potentialI— -.2ir r— Ito obtainçf’,(r) = 4 /fJ(V1n(1) V1.., ) .,1dv. (2.65)This is a Fredholm equation of the first kind with the conductivity perturbation Inp as theunknown function. Physically, (2.65) approximates the actual secondary potential by consideringCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 42only the charge accumulation arising from the primary field. The internal interaction of theaccumulated charges is neglected. Because of the linearity between the conductivity perturbationand the secondary potential, (2.65) provides a basis for formulating linear inverse algorithms toimage or recover the subsurface conductivity using surface potential data. I shall use it in the 2dimaging algorithms and 3d inversion algorithms to be presented in the next two chapters.In the special case of 2d structures, the geo-electrical property remains constant in the strikedirection. Assume the strike of the 2d structure is in y-direction, then p() = (x, z). Substitutingthis into (2.64) and recognizing that terms involving j) are independent of variable y, I cancarry out the integration with respect to y and obtain=— 4ir2uo I j [8 In it(z’, z’) (‘ — + 8 in ,LV z’) z] g (i; ,) dxWz’,(2.66)where = (x, y, 0), = (x,, y, 0) and the kernel function g is—1 . 7 iga(r;rrs= j2 2 1/2J [(x — x’) + (y—y’) + z12jI 312dy’.[(au — a9)2 + (Y—y’)2 + z12]Performing the integration for the special case y = y yieldsg (7’;f,i8)=( 2(,12— ) [ii2E(s) — 7lK(s)] , ,1 > ‘i.;12 [K(t) — E(t)], 7;’hr(2.67)CHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 43where=x) + z2,=z,)2 + z2,= /2_ 2=and E and K are the complete elliptic integrals of the first and the second kind, respectively. Itshould be noted that the integral for a kernel in the general case of y ‘ y has also been carriedout so that off-line potentials can be evaluated.Since the Born equations are basic to all the algorithms in this thesis work, it is importantto examine their validity. I attempt to quantify this for a simple geologic situation in which aconductive or resistive prism, either 2d or 3d is buried in a half-space. The cross-section of thegeologic model and electrode geometry is shown in Fig 2.8(a). The true secondary potentialsand those evaluated from (2.65) are given in Fig. 2.8(b) for a 2d prism and in Fig. 2.8(c) for thecube. Generally the shape of the Born potential is very similar to that of the true potential but theamplitude can be significantly different. The Born potentials are greater than the true potentialsfor a resistive prism and smaller than the true potential for a conductive prism. Amplitudediscrepancies of 30% are noted. There are also sign differences. For a 2d conductive prism thetrue secondary potential is entirely negative but the Born solution shows a positive side lobe.The true secondary potential requires the interaction between charges; This reduces the positivecharge on the far side of the prism and makes the potential negative. Overall, however, theagreement between the true and Born potentials is quite good. In fact, the correspondence for the3d conductive prism is excellent. This provides optimism that analyses carried out using (2.65)or (2.66) will yield useful results.2.9 DiscussionThe basic theory of the charge accumulation in d.c. resistivity experiment is presented. Ihave clarified that permittivity plays no part in d.c. problems other than supply a portion ofCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTSxs= -150m0h = 50mlOOm1OOm4-a440.30.20. 1‘•— 0.0-0.2—0.30.150.100.05‘— 0.00—0.05—0.10-0.15bC—400 —200 0 200 400x (m)I ICI I—400 —200 0 200 400x (m)Figure 2.8 The comparison of the true and the Born approximation potentials over 2dand 3d prisms are shown in panel (b) and (c) respectively. The two prisms have theidentical cross-sections, as shown in panel (a), and are embedded in a uniform half-spaceof 0.01 S/rn. The conductivity contrast between the prism and the half-space is 10. Thesolid lines are the true secondary potentials. The dashed lines represent the potentialsfrom the Born approximation. Curves labeled C and R correspond to conductive andresistive prisms respectively.Cthe charges accumulated in the regions of conductivity variation. I establish the relationshipbetween the image method of calculating potentials and the effects of true charge accumulation,which states that the essence of the image method is to find a set of fictitious charges whichreproduces the potential from the physically accumulated charges. I also show that the chargeaccumulation is the direct cause for refraction of current at a boundary separating media withCHAPTER 2. THEORETICAL BASIS OF d.c. RESISTIVITY EXPERIMENTS 45different conductivities. I quantify the charge distributions through the use of integral equationswith half-space and whole-space Green’s functions, which provides a greater understanding aboutthe d.c. electrical problem.However, the most important benefit of this investigation of charge accumulations has beenin understanding the observed d.c. potentials. When the geologic structures are simple, it isoften possible to sketch the magnitude and sign of charges on buried bodies and on the earth’ssurface that result from an arbitrary specification of the current source location. Thus secondarypotentials and apparent conductivity perturbations can be estimated so as to delineate essentialfeatures of the conductivity anomaly. This helps to establish an intuitive understanding ofthe relationship between apparent conductivity anomalies and geo-electrical structures which isimmediately useful in making preliminary interpretations of observed potentials or in making firstorder judgements about the validity of results produced from numerical simulation.It is also shown that the Born approximation provides a reasonable representation of the truepotential. The approximation is especially good in 3d cases. This indicates that it may be viableto image or invert for subsurface structures using observed potential data based upon such andapproximation.CHAPTER 3IMAGING SUBSURFACE STRUCTURESE-SCAN® data are distinguished from the conventional d.c. data in that the potential field onthe earth’s surface is recorded at many sites for each location of the current source. The numberof data is further compounded by the large number of source locations which are used. This isin sharp contrast with the conventional format of d.c. resistivity experiments where a fixed arrayconfiguration is employed and all active electrodes move from one site to another in a group andonly the data corresponding to the chosen array configuration are acquired at those sites. Of allthe data sets from different surface d.c. resistivity experiments, an E-SCAN® data set providesprobably the most complete coverage of the survey area with a great amount of redundancy. Ateach source location, the input current field energizes the earth differently and hence each setof surface potential measurements associated with a current source provides us with differentinformation about the underground geo-electric structure. Such completeness in the data enablesone to extract information by many different means.For preliminary interpretation, only qualitative information need be extracted directly fromthe data set. This is a necessary step before inversion methods are applied. The informationobtained in this stage can sometimes answer the fundamental questions regarding the existenceof prospective anomalies or regional structures of interest. It can also provide a starting pointfor further quantitative work. Qualitative information is presented as images. The simplestis to plot the data as apparent conductivities. Features seen in the apparent conductivity maybe interpretable directly. The next step in sophistication would be to use the potential data toconstruct a structural image of the actual conductivity. The secondary potential may be invertedto find locations of accumulated charges which are associated with conductivity anomalies. Thisserves to image boundaries of the anomaly. A crude estimation of the conductivity may also beconstructed from the surface potentials. In this chapter, I shall present these imaging methods inthe order of increased processing involved.46CHAPTER 3. IMAGING SUBSURFACE STR UCTURES 473.1 Direct Imaging via Apparent ConductivityThe apparent conductivity is a convenient and intuitive format to display and to examine d.c.resistivity data. It is denoted by o. and defined byOa=O, (3.1)where Aq is the measured total potential difference between the two measuring electrodes, I is thestrength of the current source, and G is the geometrical factor determined by the inter-electrodedistances of a given array configuration. C is given in its general form by1/1 1 1G=————————+— J , (3.2)2ir \ RAM RAN RBM RBN Iwhere RAM, RAN, RBM, and RBN are the inter-electrode distances and (A, B) and (M, N) representthe current and potential eletrode pairs respectively. When an electrode is placed at infinity, theassociated terms vanish. For instance, the geometrical factor for a pole-pole array is given byC = (2IrRAM)’.The apparent conductivity data are usually presented in the form of plan view maps or pseudo-sections. They are indicative of the subsurface conductivity variation. The multiplicity of thedata allows apparent conductivities to be synthesized for arrays with different separations and indifferent directions. This allows apparent conductivities to be formed in 3d space. Such a datavolume provides a blurred 3d image of the actual conductivity variation and the plan-view andpseudo-section maps are easily extracted.Of all commonly used configurations, the pole-pole array has the largest signal amplitudeand effective depth. Apparent conductivities from all other arrays can be obtained by takinglinear combinations of pole-pole apparent conductivities. Therefore, I choose to form the datavolume from pole-pole arrays. A pole-pole array is distinctly defined by specifying its separationbetween the current and potential electrodes, and the orientation in which the two electrodes arealigned. The data for a pole-pole array are commonly registered at the mid-point of the array.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 48This is justified by the fact that measurements obtained by interchanging the current and potentialelectrodes are the same due to the reciprocity of the d.c. field.SinceE-SCAN® data are usually acquired over a regular grid, many different pole-pole arrayscan be defined. The orientations are constrained by the grid since the electrodes must be on gridpoints. The electrode separations are multiples of the basic spacing between grid points alongthe specified direction. For each chosen pole-pole array, many data points can be synthesizedthroughout the survey grid so that a data map can be formed. It is apparent that for a givenelectrode separation, there is a number of arrays with different directions.A set of pole-pole arrays is first chosen and the corresponding apparent conductivity data aregathered. The apparent conductivities of each array are then interpolated to form a data map overtheE-SCAN® grid. Since the positions of the array electrodes are limited by the grid, the apparentconductivity is defined only over an smaller region within the B-SCAN® grid. To avoid arbitraryextrapolation near the edge, I assign the mean value of the apparent conductivities gathered fromthe E-SCAN® data to the edges that have no defined apparent conductivities. The interpolationis then applied to the data consisting of the directly gathered data and the assigned points onthe edge. The maps with the same separation but different directions are then averaged to forma directionally averaged map. Each averaged map is assigned at a pseudo-depth proportionalto the separation. The composite of all such map forms a 3d prism. Further interpolationin the vertical direction is performed when necessary and this generates the desired apparentconductivity volume.The two horizontal dimensions of the data volume are consistent with the actual spatialcoordinates but the vertical dimension does not correspond to a real depth. Thus the data mayindicate the correct horizontal location of an anomaly but provide only a relative depth withrespect to other anomalies. With an appropriate choice of the mapping to assign a pseudo-depthZe from a given electrode separation, the image from the apparent conductivity data can beindicative of the depth of burial for simple anomalous structures. Such mappings have beenstudied by a number of authors (e.g., Roy and Apparao 1971, Edwards 1977). Edwards (1977)CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 49provided a rather complete list of the mapping factors for commonly used arrays. The mappingfactor for pole-pole arrays is 0.86, i.e., Ze =O.86l?AM.Based upon the identity derived by Roy (1978) for the surface potential over a half-space V,=- JJ V4’) . V-dv, (3.3)where,1 is the distance between the current and potential electrodes, Edwards defines a signalcontribution function of the depth zfe(’i, z)= _JJ Vq5(’). V!dxy. (3.4)so that,1) Jfe(,Z)dZ.He then defines the pseudo-depth to be the depth above which the conductivity is responsible forhalf of the potential signal. That is,Zmj= J f(,z)dz. (3.5)This definition is essentially based upon the electric field distribution in the lower half-space.If the background conductivity deviates greatly from a uniform half-space, the above mapping factor will no longer be valid. I have carried out a study of the depth mapping in the caseof layered Id background following the approach of Edwards (1977). I first extend (3.5) to thevarying id background. The signal contribution function f(i, z) is formed using the numericalsolution for the potential assuming the Id earth is composed of a sequence of layers with constantconductivities. The mapping factor given by the ratio of pseudo-depth over the electrode separation defined by Zd, is no longer a constant. Instead, it is a function varying with the separation.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 50I also show that an identity similar to (3.3) exists, i.e.,= 2iro JJ o(z)V4(i’) V!dv. (3.6)This is based upon the distribution of electric current field in the lower half-space. A signalcontribution function can be defined accordingly asfOi, z)= 2lrcTaJJcT(z)VqS(i’). V-dxy. (3.7)Results from numerical simulations indicate that a signal contribution function formed fromweighted sum of feO’l, z) and f(i, z) is necessary in order for the definition of the mappingfactor work for background conductivities with conductive and resistive layers. The criterion isthat the apparent conductivity curve plotted using the varying depth mapping should position theanomalies approximately at the depths of the corresponding conductivity layers. One combinationisf(,z) = fe(?1,Z) + pf(i,z) (3.8)l+ILa 1+t&where Pa = 0a(77)/0 and u0 is the surface conductivity. Simulations suggest that an embeddedlayer reduces the mapping depth for an array which has strong responses to that layer. Fig. 3.1shows the mapping depth for a three-layered background. However, the use of such varying depthmapping relies upon the knowledge of the background conductivity. The study of its effectivenessis inconclusive and so is its possible application to field data.The apparent conductivity data can also be displayed in the form of relative anomalies withrespect to the apparent conductivity of the pure id structure. Numerical experiments indicatethat when the relative anomaly data are plotted with the half-space mapping factor, they tend toindicate the depth of burial for simple anomalies.The data image can be enhanced by applying low-pass and directional filtering at each depthin the data volume. Low-pass filtering is used primarily to accentuate the coherent features inCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 51HE______100 I 11111111 I I 1111111 I I 1111111 I II100 101 102 io10-1I I lull I I 111111 I—I—t—•I—i I I101 102 ioZ, ZFigure 3.1 Varying depth mapping for a three-layered background. Panel (a) is thepseudo-depth as a function of the pole-pole array separation. The triangles are the idbackground mapping and the solid line is the half-space mapping. Panel (b) shows the Idapparent conductivity (solid line) plotted with pseudo-depth in (a). The triangle indicatethe actual data points. The dashed line is the true conductivity.a data map by eliminating the noise components. The coherent features are mainly in the lowwavenumber band. There are two major sources of data error. The first is the random noiseassociated with processes such as potential measurements, electrode positioning, telluric currentand cultural noise. The second is the electric signal from shallow or near-surface conductivityanomalies. Small scale conductivity anomalies close to electrodes can distort the electricalmeasurement even at large current-potential separations. These measurements are assignedto a greater pseudo-depth proportional to the separation. Thus the distortion of the near surfaceconductivity variations can appear as deep anomalies in the apparent conductivity data. Therefore,the response of an array to near surface anomalies may constitute valid signal at small arrayseparations but it becomes effective noise for the apparent conductivities measured at largeCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 52electrode separations. These distortions usually concentrate in the high wavenumber band. Whenapparent conductivities are displayed as a means of imaging, they should be smoothed to suppressthis type of noise. The degree of smoothing to be applied should reflect the decreasing resolutionof an array as the electrode separation of the array increases. Low-pass filtering in the wavenumberdomain is an effective way to carry out the smoothing. Geological features may exhibit stronglinearity in certain directions. Directional filtering is an effective way to either enhance or tosuppress such linear features.The pole-pole array is advantageous for its large effective depth and simple anomaly pattern.However, another good candidate is the square array. It was introduced by Habberjam (1979).A square array involves four active electrodes on the surface. The electrodes A, B, N, M arearranged so that each occupies a corner of a square in that order. As the array expands all theinter-electrode distances expand simultaneously. Although acquiring square array data is moredifficult in the field, synthesizing them from an E-SCAN® data set can be readily achieved. Inaddition, since E-SCAN® data are acquired over a surface with planar arrays (in contrast to lineararrays), it is reasonable to display the data using a planar array such as the square array. Thepseudo-depth of a square array is given by Ze = O.45a, where a is the length of one side of thearray.Despite its relatively shallow effective depth compared with a pole-pole array, there aretwo advantages with the square array. Firstly, the data synthesized for square arrays fromE-SCAN® data sets will not be affected by the placement error of the infinite electrodes B and N.Secondly, a square array measures a very compact anomaly over a body of finite geometry andthe response of shallow structures at large off-set is very small. This renders the array insensitiveto the near-surface variation in conductivity. Therefore, a square array may represent a goodtradeoff in the quest for a large effective depth and reduced sensitivities to near surface variations.To conclude the discussion, I present a simple example in Fig. 3.2. A conductivity modelis formed which has a conductive cube buried beneath a surface layer of random conductivityvariations in an otherwise uniform half-space. The cube is 200 metres on side and locatedFigure 3.2 Responses of the pole-pole array and the square array to a conductive cubeburied beneath a surface layer with correlated random conductivity variations. Theapparent conductivity from pole-pole array is shown in (a) for pseudo-depth Ze = lOOm.A low-pass filtered version of (a) is shown in panel (b), which has greatly suppressed theeffect of surface variation. The data from the square array at the same pseudo-depth areshown in panel (c). It is clear that the square array is less sensitive to surface variations.The grey scale is in log10 Ua.0.200.400.600.800.1000.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 53—1.83—1.86—1.89—1.92—1 .—1.98—2.01—2.04—2.0?0. 200. 400. 600. 800. 1000.X (m)0.200.400.600.800.1000.0. 200. 400. 600. 800. 1000.X (m)0.200.400.600.800.0, 200. 400. 680. 800. 1000.X (m)CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 5450 metres below the surface. The surface layer consists of correlated random perturbations tothe uniform background with a correlation length of 50 metres. The response to this randomlayer dominates at all separations in the data from pole-pole arrays. The slice at pseudo-depthZe = lOOm, as shown in Fig. 3.2(a), is characteristic of this “geological” noise. The same slicefrom the low-pass filtered version is shown in Fig. 3.2(b); the anomaly caused by the buriedconductive cube is clearly displayed. As a comparison, a slice of data at the same pseudo-depthfrom square arrays is shown in Fig. 3.2(c). No low-pass filtering was applied, yet a concentratedanomaly for the buried cube is well defined over the relatively weak response to the near-surfacevariations.3.2 Imaging via Equivalent Source from Secondary PotentialThe secondary potential arises from the accumulated charges associated with conductivityanomalies. An image of the charge constructed from the secondary potential, therefore, shouldprovide useful information. In this section, I seek to construct such an image by a sequence ofequivalent source layers.The potential field is related to the subsurface charge density by (see chapter 2)-,I iii p(r)-.., dv, (3.9)2ireoJjJ Ir—r IVwhere i and f’ are the observation and source points respectively, p is the charge density, and Vis the lower half-space. Adopting a right handed Cartesian coordinate system with origin at thesurface and the z-axis pointing downward, the above equation can be rewritten asç(z,y)=— fp(x,y,z’)øø dz’, (3.10)2lrfoJ /2+y2+zt20where the symbol 00 denotes a 2d convolution in the x-y domain. It follows from (3.10) thatthere exists a layer of equivalent source distribution at any depth which will reproduce the givenCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 55potential on the surface. Thus, I can proceed to construct a set of equivalent source layers with aunit thickness at depths z (i = 1,• , n) from1 (3.11)2ir0When these source layers are assembled together according to their depths, they form a 3dequivalent source image which may indicate the concentration of accumulated charges. A similartechnique is proposed by Karous and Hjelt (1983) for interpreting VLF-EM measurements.The solution of (3.11) is a linear deconvolution problem in a 2d domain and can be solved byspectral division. Define the 2d Fourier transform pair as (Bracewell 1978)f(p, q) =JJf(x, y)e”’thcdy,(3.12)f(x, y) =iJJ !(P, q)e”dpdq.where (p, q) are the wavenumbers in the x- and y-directions respectively. Applying the Fouriertransform to (3.11) yields/3(p,q,zj)=,(p,q)/p2+q2eV1 i= 1,•”,n (3.13)where and /3 are the Fourier transform of qS and p respectively. ,Q, q) can be computed byFourier transforming the secondary potential on the surface. Thus, given a set of depths z, /3 canbe calculated and the equivalent source layers in the spatial domain are obtained by applying aninverse Fourier transform.The equivalent source construction is essentially a two-step process. The first step downwardcontinues the surface potential (by the operator eV) and the second calculates the verticalderivative of the continued potential (by operator iJp2 +q2). The downward continuation is adivergent process. In order to avoid unphysically large amplitudes in the constructed equivalentCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 56source, I apply a low-pass filter to ,. The cut-off wavenumber decreases as the depth z increasesand is determined by(3.14)where S is a specified percentage. 1/S effectively limits the magnification allowed in the downward continuation and, therefore, eliminates higher wavenumber components from the equivalentsource at greater depth.Fig. 3.3 illustrates the imaging method using the secondary potential over a single prismburied in a uniform half-space. The secondary potential map generated by a single point currentsource over a conductive prism is shown in Fig. 3.3(a). The application of the above imagingmethod generates a 3d image of the equivalent source whose two sections are shown in Fig. 3.3(b) and (c). The boundary of the buried prism is outlined in each section. It can be seen that theminimum of the equivalent source corresponds to the top portion of the prism. From the viewpointof imaging, this simple processing of the secondary potential has provided useful information.In general, the extrema in such a 3d image serve to indicate the depth of burial of an anomalousconductivity, since the top portion has a greater density of accumulated charges.With multiple sources, the secondary potential maps can be stacked to simulate the potentialfrom near vertical or horizontal primary fields. The equivalent source obtained using these stackedsecondary potential maps provides better images as the effects from the strong asymmetry of theprimary field relative to the conductivity anomaly is reduced.The method is purely based upon the mathematical properties of the potential field. Assuch, its usefulness is very limited. The image is rudimentary but the approach can be furtherdeveloped. One can utilize the properties of the actual d.c. field to invert the potential for thecharge density, and hence, image the conductivity anomaly. This leads to the next level of theprocessing, namely, charge density imaging.8147.6-220.293.367.Figure 3.3 Imaging conductivity anomaly using equivalent sources. Panel (a) showsthe secondary potential (mV) generated by a point current source (marked by the cross)over a conductive prism in a uniform half-space. Panel (b) is the section at x = 500mfrom the constructed 3d equivalent source image. The magnitude given by grey scalerepresents a scaled charge density and only relative strength of the source. Panel (c) isanother section at y = 500m. The outline of the prism is shown by the white square inall panels.0.200.400.600.800.1000.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 570.—20.0—30.0—40.0-60.0—60 .0—70.0—80 .00.-6.00—10.0—15.0—20.0—25.0—30.0—35 .0—40.073.0. 200. 400. 600.Y (m)b800. 1000.86-0.1.10.73.147.220.293.367.0. 208. 400. 608. 800. 1000.X (m)CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 583.3 Charge Density Imaging of 2d StructuresIn practical applications, if the anomaly due to a geologic structure has a preferred strikeand if it is sufficiently elongated, the 2d approximation may be valid. Even though the Bornapproximation may be somewhat more inaccurate for the 2d models compared to a 3d model,from a computational viewpoint there are benefits to using a 2d formalism if it is possible. Forthis reason, I develop a 2d imaging algorithm using charge densities based upon (2.66). The terminvolving z) in (2.66) can be written asI Vlnp(x,z).,5c(x,z)=— 247ro0 ‘iswhere i = ( — x, z) is the projection onto the z-z plane of the vector from the source to a givensubsurface point. (2.66) can then be written as, (x)=J Jc(x,z)9 (x’, z’;z5,z) dx’dz’, (3.15)where g(x’, z’; z,, x) = sga is the new kernel function. The function c(x, z) represents the totaleffect of the accumulated charge and behaves like the accumulated charge density: it is non-zerowhenever there is a component of 4 parallel to the conductivity gradient V In It vanishesidentically away from any boundary of the conductivity structure.Equation (3.15) linearly relates the charge density c(x, z) to the secondary potential. Assuch, secondary potentials recorded at a number of locations x, and arising from a fixed source atx, can be inverted to recover an estimate of c(z, z). Given that physical charge accumulates atthe boundaries between blocks of different conductivities, it follows that an algorithm which iseffective in delineating the locations of the charge, can image the subsurface structure.This imaging is accomplished by setting up a linear inverse problem. The cross-sectionperpendicular to the strike is first discretized into m x m rectangular cells according to thepartitioning(zi,... ,,...and(z1,... ,Zmx). Ineach cell c(x,z)isassumedtobeCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 59constant. Equation (3.15) becomesme mx= (3.16)j=1 k=1where ç is the datum at location x, n is the number of the data available, and is the integralof th kernel function over !Cth cell AJk,Ak= JJ g (x’, z’; x, x) dx’dz’.‘jkThe number of data are generally far fewer than the number of cells and hence (3.16) isan underdetermined system. Equivalently, the solution is non-unique. To obtain a particularsolution, I seek a “simple” charge density which can explain the observed secondary potentialdata. The measure of simplicity is defined by the weighted l norm of the charge density. I donot fit the data exactly but rather introduce a misfit variable e for each data equation and permita total misfit bounded by a prescribed value. The linear inverse problem becomes the followingminimization problem,me m;mm. IEWjk C,k I,j=1 k=1me mx A:subjectto I, ,n (317)31 k=I ‘Ie:Iwhere 8, is the estimated standard deviation of the error associated with th datum, /3 is a fitparameter which controls the fit to the data is the expected value of the total misfitassuming a Gaussian distribution of the noise. w is a set of weighting coefficients whichultimately determine the type of model which is obtained. We are free to choose any weightingfunction. I have specifiedWJJg = (zak— ZO)312?la,CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 60where Zck is the central depth of the cell. The first part is introduced to overcome the naturaldecay of the kernel function with the depth and the second part to compensate the decay of theprimary field with the distance.The above minimization problem is solved using a linear programming technique (e.g., Murty1983). This choice suits the need for finding a sparse representation of the charge density, sincethe number of non-zero elements in the LP solution is bounded by the number of constraints.In order to accommodate the absolute value of variables in both the objective function and themisfit constraint, each bipolar variable needs to be expressed as the difference of two positivevariables. In this way, when the absolute value is replaced by the sum of the two correspondingpositive variables, the solution of the new minimization problem yields the solution for the originalproblem.When a current is put into the ground the physical charge density depends upon the sourceposition relative to the conductivity anomaly. The charge density cjk found from the aboveinversion reflects this relationship. The result only provides a partial image of the boundary.However, the field from a different current source position illuminates the boundary in a differentway. A relatively complete image of the boundary can be obtained by inverting for c(x, z) frommany different source positions using (3.15).Taking the advantage of the E-SCAN® data set, I can form many sets of common sourcepotential data. Each set can be inverted to yield a part of the boundary image. When all ofthe partial images are combined, they give a rather complete picture. Fig. 3.4 shows the resultobtained by inverting synthetic data generated from a simple prism buried in a uniform halfspace. Seven sets of secondary potential data corresponding to different source locations, shownin Fig. 3.4(a), are used in the inversion. Fig. 3.4(b) is the combined image. The top and thebottom boundary of the prism are correctly imaged. The charge density also has the correctsign. It is negative at the top where the current flows into the conductive body and is positiveat the bottom where the current flows out of the prism. Fig. 3.4(c)-(f) show some results fromCHAPTER 3. IMAGING SUBSURFACE STRUCTURESFigure 3.4 Imaging subsurface structure using charge density. Panel (a) shows thesecondary potentials corresponding to the labeled source locations over a 2d conductiveprism buried in a uniform half-space. Panel (b) is the composite of the individuallyrecovered charge density images corresponding to seven sources located from -300 to300 meters in equal spacing. The light region represents the negative charge density andthe dark region the positive charge density. The solid line in the image indicates theboundary of the 2d prism. Panels (c)-(f) are images corresponding to sources at 0, 100,200, and 300 meters, respectively.61—400 —200 0 200 400X (m)JII I I IC_n C_.C t —o o 0 0 0o o 0 0 0xI I I I I C C_nC_n (.C C__n — 0 0 0 0 00 0 0 0 0 0 0 0 0 0o 0 0 0 0x600-I C CC0 _n C.C 0 — 0 0 0 00 0 0 0 0 0 0 0 0 Co o 0 0 0x6000 — _nC_fl C__n C__n 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0x600! I I I I 0 C_n c_n CC_CC CC C — 0 0 0 0 C0 0 0 0 0 0 0 0 0 C0 0 0 0 0xCHAPTER 3. IMAGING SUBSURFACE STRUCTURES 62different individual inversions. They clearly show the different parts of the boundary imaged bythe potentials from the different source locations.From the viewpoint of physics, this approach is very appealing. The charge accumulationon the boundary is a physical reality in the d.c. experiment and any attempt to recover it isvaluable. The above example illustrates that the charge density imaging can work satisfactorilyfor simple structures. However, for a closed body in a uniform half-space, the sum of the chargeaccumulated on the boundary is zero. So the secondary potential is the field of dipolar andhigher-order sources. The recovery of monopole distributions from such fields is difficult and thesolution can be highly non-unique.Because the subsurface charge density depends upon the source position, there is no singledistribution of charges which is compatible with all the data. Therefore, in the recovery of thesecondary charge density, each inversion uses a single current field. Results from several suchindividual inversions form a composite image as the final result. A more stable approach wouldbe to invert for a single model using all the data simultaneously. This leads us to the directformulation to recover a conductivity image.3.4 Conductivity Imaging of 2d StructuresEquation (2.66) presents a linear relationship between the gradients of ln p(i) and the secondary potential. By forming a linear inverse problem to recover ln fL() directly, one cansimultaneously invert data from many different current locations and thereby greatly reduce thenon-uniqueness. I adopt the same type of discretization as in the charge density imaging. Theconductivity section is divided into m x m cells each of which is specified by a constant valueof In . The model for the inverse problem then becomesmk lfl(j.L3k). j = 1,...,m; k=CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 63Under these assumptions, the gradient becomes discrete. is non-zero only on the verticalinterfaces between cells and is non-zero only on the horizontal interfaces:ôln=(mk—m_1,k)S(z—xJ), (x,z)EF7k,ôx (3.18)ôlnttbz =(mJk—mj,k_1)S(z--zk), (x,z)EFk.Where 1 and denote the vertical and horizontal interfaces respectively. When the index jequals 0 or m + 1, or index k equals 0 or m + 1, the model parameter mJk is taken to be thebackground value, which is zero.Substituting the gradient functions (3.18) into (2.66) yields a discrete data equation of theformm mr= mjk7k (3.19)j=1 k=Iwherei — i i zi zi7k 7jk 7j÷1,k 7J1C 7Jk+1=( — ) Jga (x’,z’;x3, )dz’,(3.20)= —2zkJg’ (x’,z’;x8,2)dx’.I choose to minimize an objective function composed of a combination of the l norm of themodel and the 11 norm of the model gradient.=cr I mJk I +(1— a) [w7k I — I +Wk I — mJ,k_I i]j,k j,kThe weighting coefficients Wjk, Wik and 1.14k are introduced to control the relative importance ofa particular cell in the objective function based on the cell size and depth. The parameter a takeson a value in [0, 1] and determines the relative importance between the model and the gradient.Mathematically, the second part is a measure of the variation in the model in the sense of l norm.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 64When a approaches unity, the inversion produces a smallest model without much restriction onthe variation. As a approaches zero, the inversion produces a minimum-variation model (Dossoand Oldenburg 1989).The minimization problem can be solved using standard methods. Misfit parameters areintroduced, as they were in the charge density imaging, so that the final data are fit only towithin a global tolerance. Variables which can acquire both positive and negative values areagain replaced by the difference of two positive variables and the corresponding absolute valueis replaced by the sum of the two new variables. In addition, an upper and lower bound can alsobe imposed on the model elements as constraints. The kernel are computed numerically. Itcan be shown that the kernel is symmetric for interchanged source and receiver locations.Therefore, for every pair of electrodes in the E-SCAN® data, only one datum entry is required asthe d.c. potential satisfies the reciprocity.To illustrate the above algorithms, a set of secondary potential data from a synthetic modelis inverted. The model consists of a conductive prism and a resistive prism buried in a uniformhalf-space (see Fig. 3.5(b)). The half-space conductivity is lmS/m. The conductivities of thetwo prisms are 4mSIm and 0.25mSIm, respectively. A set of E-SCAN® data are generated overa traverse of 21 grid points. Fig. 3.5(a) shows the corresponding apparent conductivity pseudosection. 200 secondary potentials corresponding to this pseudo-section are used in the inversion.The model section is divided into 300 cells with thickness increasing with the depth. An upperbound of 10 for conductivity contrast is imposed on the solution. Fig. 3.5(c) shows the resultsof a minimum variation model with a = 0.9. The conductive prism is defined reasonably well.So is the top of the resistive prism. The poorer recovery of the bottom of the resistive prism isdue to the fact that the Born approximation differs more from the true potential in the case ofresistive prisms. However, for the purpose of imaging the structure, the algorithm has workedsatisfactorily.CHAPTER 3. IMAGING SUBSURFACE STRUCTURES—2.95_a.o’iJ—3.0122—3.31—3.543@.X(m)Figure 3.5 Imaging subsurface conductivity. The synthetic model shown in panel (b)consists of a conductive (4mS/m) and a resistive 2d prism (0.25mSIm) in a uniform half-space of 1 mS/rn. 200 secondary potential data on the surface are used in the inversion.Panel (a) is the corresponding apparent conductivity pseudo-section. Panel (c) is therecovered conductivity model. The upper grey scale is for panel (a) and the lower greyscale is for panels (b) and (c). The scale is in log10 u.65) N0 0 0 0o o 0 0X (m)—3.01— N CZo o 0 0o o 0 050016001I I I C — N C . ø0 N — 0 0 0 0 0o o 0 0 0 0 0 0 0 0o 0 0 0 0X (m)—2.51—2.74—3.01—3.2 CI I I I I0 C. N0 0 0 0 00 0 0 0 0— N C) 00 0 0 0 00 0 0 0 0CHAPTER 3. IMAGING SUBSURFACE STRUCTURES 663.5 DiscussionThe display of apparent conductivity for direct imaging is discussed. It is demonstrated thatwavenumber filtering can be used to suppress geological noise such as near-surface variationsin the conductivity. My study also indicates that a square array is superior in dealing with suchnoise. An algorithm is also presented to construct a 3d equivalent source image from commonsource secondary potentials on the surface. The image can indicate the depth of burial for simpleconductivity anomalies.Two imaging algorithms are then developed for 2d models using the Born equation. The firstapproach is to image the boundary of subsurface structures using the accumulated charge density.The secondary potential is produced by the accumulated charges on the subsurface boundariesand consequently the recovery of the charge density effectively images these boundaries. Linearprogramming is employed to invert the secondary potentials for a sparse representation of thecharge density due to each of the current sources. The composite of these images yields a relativelycomplete image for simple subsurface structures. The algorithm has worked satisfactorily forsimple models. However, the difficulty lies in combining a final result from several independentinversions. This is generally less stable than producing a final image using all the available dataat once.To overcome the difficulties encountered in the charge density imaging, a direct formulationis developed which inverts approximately for simple conductivity structures. The algorithm alsoemploys linear programming for the inverse solution. An approximate conductivity model isconstructed which satisfies all the available data. The method works well for relatively simplestructures. For the smallest l model construction, the algorithm is very efficient. As thecomponent of variation increases, the algorithm becomes slower due to the linear programmingsolution. Because of the approximations involved, the semi-quantitative result mostly serves asan image representation of the model. In general, the approximation is less accurate in the 2dcases than it is in 3d cases. Consequently, the application of such an inversion can be limited inpractice.CHAPTER 4APPROXIMATE 3d INVERSIONThe d.c. resistivity problem is by its nature a 3d problem. Firstly, the source in d.c. resistivityexperiments is always a point source so the electric field is 3d. Secondly, geo-electrical structuresare generally 3d. Although the study of id and 2d problems can provide much insight intotheoretical work and data interpretation, in practice their application is rather limited. Only whenthe 3d problem is tackled, can we use the d.c. resistivity method to its fullest advantage. Inprinciple, the techniques used in 2d analyses can be extended to 3d conductivities. In practice,however, the large numbers of model elements and data result in a large matrix to be invertedand the approach is impractical. I choose to develop a different methodology motivated by themultiplicity in E-SCAN® data.An integral equation is developed in which the secondary potential anomalies from a pole-pole array is expressed as a depth integral of the logarithmic conductivity anomaly convolvedwith a kernel function in the horizontal plane. The kernel takes the separation and orientationof specific current and potential electrode pairs as parameters. Fourier transforming the dataequation decouples wavenumber components so that the surface response of a given pole-poleconfiguration at any wavenumber is obtained by carrying out a single integral with respect todepth. This decomposition allows a full 3d inversion to be carried out by solving a id linearinverse problem at each wavenumber. In the following I shall first develop the integral equationand then proceed to the inversion.4.1 Data EquationI begin with the Born equation (2.65),67CHAPTER 4. APPROXIMATE 3d INVERSION 68Applying the vector identity V . (ba) = çbV +V ã for an arbitrary scalar b and an arbitraryvector a0, to the integrand in (2.65) yields(i)= JJJintt(’)V. ( i )civ4r Ob Ir —rob.I I’ —r,!+ JJJv. (in (p’) _, V _, dv.(4.1)4ir 0o Ir —r,&,I Ir —r.IjVWith the aid of Gauss’ theorem the last integral in the above equation can be expressed asI ff I I3=4ir2aJf lnJL(r) .nds,where S represents the entire boundary of the lower half-space with outward normal vector i.The contribution to /3 from the flat upper surface of the half-space is zero because- •%=Othere. Along the remaining boundaries,1 1 1—, V ., • n cc —i as Ir I —. oo.Ir —rI Ir_7’.I 1r1and hence the boundary integral vanishes. Thus /3 is identically zero.Applying the same vector identity to the first integral in (4.1) yieldsçS = — JJJin (i.1) V....,• V dv4ir 00 Ir —r0z, lr’—r,I+ 4 flJin (p’) ,1V2 1 1dv.(4.2)The second term in (4.2) vanishes identically sinceV2 -. 46(’j— r5ICHAPTER 4. APPROXIMATE 3d INVERSION 6921______0_______yFigure 4.1 Geometry of the surface pole-pole array is specified by a vector 2rpointingfrom the source location , to observation o is the midpoint of the array, where thepole-pole potential is recorded.and ln jL(f) = 0 because it has been assumed that the surface conductivity is crc. Thus I obtainthe following expression for the secondary potential on the surface,JJJlnhL(f)V ., I—, dv. (4.3)4ir 00—r—VEquation (4.3)is a Fredholm equation of the first kind relating the secondary potential to theperturbation of logarithmic conductivity. The equation is similar to that presented by Boernerand West (1989) which was derived using an integral equation for the potential on the boundaryof an isolated inhomogeneity with constant conductivity.Since E-SCAN® data are generally acquired on a regular grid, it is convenient to explicitlyconsider data generated from pole-pole arrays with fixed separations and orientations. As inFig. 4.1, let io= (x, y, 0) be the midpoint of the array, and let 2lbe the vector pointing fromsource electrode to potential electrode. Then-07o6a = ‘O +1.CHAPTER 4. APPROXIMATE 3d INVERSION 70Substituting into (4.3) and converting into relative anomalies 84 = cb3/ç yields= _!-JJJ 1n(L(f )V V dv, (4.4)I1—’I Io—l’Iwhere 1= Il’. Since io = (x, y, 0) and 1 = (li, li,, 0), it is recognized that (4.4) involves aconvolution operation in the x-y domain. Thus the surface pole-pole potential anomaly can beexpressed as00Sqf; 1)= J ln tt() ® Øg(i; f)dz, (4.5)0where g(; 1) is the kernel functiong(;1)=—--V I V 1 -I+lI Ir—l!1 [(a, + l) (a, — l) + (y + l) (y — i) + z2] (4.6)2 3/2 2 3/2[(x+l)2+(y÷l) +z2] [(X_l)2+(y_l) +z2jEquation (4.5) relates the surface pole-pole potential anomaly data from a pole-pole arrayto the subsurface conductivity structure. It is noticed that the anomaly is a depth integral of theconductivity perturbation convolved with a kernel in the z- and y-directions.Define the Fourier transform F over the infinite plane and its inverse transform asgiven by (3.12)f(p, q) =JJ f (z, y) e_dxdy,f(a,, y) = fp, q) edpdq.Taking the 2d Fourier transform of (4.5) and applying the convolution theorem yieldsj=I,.•.,ni (4.7)CHAPTER 4. APPROXIMATE 3d INVERSION 71where (p, q) are transform variables representing the wavenumber in z- and y-directions andj (p, q) = .:l:: [sc (i; ii)], j = i,... , flj,ih(p,q,z)=.:F {1nJL(],j=1,.”,nj.The index j identifies the th pole-pole array. In an E-SCAN® data set, many different pole-polearrays exist; each is specified by a unique separation and orientation of the current and potentialelectrodes. I denote the number of distinct pole-pole arrays by n.Equation (4.7) is a Fredholm integral equation of the first kind. The data ê (p, q) andmodel th(p, q, z) are complex because they are the Fourier transforms of arbitrary real functions.However, the kernels are real. This results because of the symmetry of the kernel function(4.6) about the origin. At each wavenumber, (4.7) shows that there are n constraints upon themodel th (y, q, z). Linear inverse theory can therefore be used to recover th(p, q, z) as a functionof the depth z (Backus and Gilbert 1967, Oldenburg 1984). Completing this inversion at allwavenumbers yields the model thQ,, q, z). The 3d spatial distribution of the conductivity isobtained by performing an inverse 2d Fourier transform in the horizontal directions at sampleddepths.The implementation of the above inversion method requires: (1) that a background conductivity be estimated in order to compute the relative potential anomaly; (2) that pole-pole datamaps corresponding to different separations and directions are gathered and interpolated to acommon grid so they can be Fourier transformed; (3) that kernel functions are evaluated; and(4) that id inverse problems are solved. The id inverse problem can be formulated to recovereither a continuous representation or a discrete form of th(p, q, z) as a function of depth z. Bothapproaches have merits in different aspects and will be discussed. These aspects are addressed inthe following sections.CHAPTER 4. APPROXIMATE 3d INVERSION 724.2 DataThe data required for the Id inversions are generated by applying a 2d spatial Fouriertransform to the relative potential anomaly 84i for each electrode configuration. The choice of therelative anomaly eliminates the effect of electrode separation and the data therefore reflect onlythe geo-electrical variation. The relative anomaly is also a quantity consistent with the model oflogarithmic conductivity perturbation.Since theE-SCAN® experiment measures the total potentials 4, the relative potential anomalyS4 is calculated bySç; 1 = cS(; 1) — q(i; 1)4,(; 1)This requires a reliable estimate of the background conductivity o so that the primary potentialcan be calculated accurately. The simplest approach to estimate o-o is to use the mean value ofthe apparent conductivity data. That is, I take the best fitting half-space model as the background.Experiments with both synthetic and field data sets show that this works well, especially forstructures with very localized inhomogeneities or where the conductivity value is distributedabout a regional value.Alternatively, one can formulate the id inversion at the zero wavenumber so that the background conductivity o is solved for as a unknown parameter together with the conductivityperturbation. This requires only the mean value of each potential data map. I shall discuss thisapproach following the Id inversion.The process of data preparation also requires that the potential data for each pole-pole arraybe gathered from initial E-SCAN® field data and be interpolated and extrapolated. The procedureis best explained by using Fig. 4.2. An E-SCAN® survey is generally carried out over a regulargrid. A basic grid is shown by the solid lines in Fig. 4.2. Each node in that grid would havebeen occupied by a current or potential electrode as the survey data are collected. Suppose thata data map is desired for i = (tx/2, zy). This pole-pole configuration is given by the arrows inFig. 4.2. The potential datum qf for each electrode pair having this orientation can be obtainedCHAPTER 4. APPROXIMATE 3d INVERSION 73extendedregion4, assignedE-SCANsurvey gridFigure 4.2 Specification of the data map. The intersections of solid lines indicate theE-SCAN® grid points. The potential data are acquired over all grid points but the pole-pole potential is defined only in a smaller area. One possible array configuration is shownby the arrows in the diagram, which point from the source to the potential location. Thepotential data for this array can be formed only at locations indicated by solid dots. Thisactivated area is smaller than the survey area and decreases as l increases. To generate adata map suitable for inversion, a mathematical area is defined by extending the originalB-SCAN® survey region half of the maximum array separation in x- and y-directions.and assigned to the midpoint between these two electrodes. The solid dots indicate the locationsat which potential data can be gathered for this particular electrode pair. Notice that this areaactivated by the solid dots is smaller than the initial survey area. The size of the activated areadecreases as increases. The gathered data in this area need to be interpolated and extrapolatedto a common grid over the domain on which the inversion is to be carried out.I have selected the domain to be a rectangle which is larger than the initial survey area andcontains the initial survey at its center. Typically, as illustrated in Fig. 4.2, the distance betweenavailable data pointsCHAPTER 4. APPROXIMATE 3d INVERSION 74the survey grid and the boundaries of the inversion domain have been set equal to the maximumvalues of la,, and l, in the x- and y-directions respectively.The remaining step is to define a rectangular grid on the surface of the inversion modeldomain. Because some data maps (for example that shown in Fig. 4.2) yield data at locationsmidway between the nodes of the survey grid, it is reasonable to select a spatial discretizationfor the Fourier transform grid which is half of the field survey grid spacing. With this choice ofgridding there will be many points on the inversion grid that are not defined directly by valuesfrom the survey grid. Both interpolation and extrapolation are required. Interpolation inside anactivated area (e.g. within the rectangle generated by the solid dots in Fig. 4.2) is not difficultbut extrapolation is always dangerous. In order to avoid the introduction of extraneous artifactsand also ensure that the data to be Fourier transformed are zero around the edges, I specifythe potential data around the boundary of the inversion grid to be c6, = I/4rlcro, the potentialcorresponding to the background conductivity. The evaluation of potentials at all nodes on theinversion grid can now be obtained by interpolating a data set consisting of the interior activatedpoints and the outside boundary points. I apply an interpolator consisting of a combination ofLaplacian and cubic spline components to ln(l’t/I). After interpolation, the relative anomaly 84is calculated and Fourier transformed to generate the data actually used in the inversion.4.3 Kernel FunctionsEquation (4.6) gives the kernel function in the spatial domain. For illustration, Fig. 4.3 showsthree slices through the spatial kernel function, g(; 15, for a given 1. The distance between theelectrodes is 50 m and the depth slices are at z = I , 20, and 40 m respectively. At shallowdepths the character of the kernel function is dominated by two large dipole-like features whichare directly beneath the current and potential electrodes. These features become singularities asz —f 0. At greater depths these peaks become smoother and spread out horizontally. It shouldbe noted, however, that for shallow depths, (e.g. Fig. 4.3(b), (c) ) the spatial kernel function hasinner regions of opposite sign. These observations are all in agreement with our understandingCHAPTER 4. APPROXIMATE 3d INVERSION—2x]—4x:75—2x107,— —4x1O—6x1O—8x10’ -0 0 0 00 0 0 0X (in)Figure 4.3 Spatial domain kernel function corresponding to (1,,, li,) = (25,0) at depthsz = 1,20,40 m respectively.X (m)X (in)CHAPTER 4. APPROXIMATE 3d INVERSION 76of the pole-pole dc resistivity experiment. For a given separation, the array is sensitive to smallerchanges in the conductivity structure from shallow regions but it responds only to a broaderaverage of conductivity structure at deeper regions. The positive region agrees with the fact thatthe array has a negative response to relatively shallow anomalies. For instance, it measures aresistive anomaly over a shallow conductor (cf chapter 2).The kernels (p, q, z) needed in the id inversion require that (4.6) be spatially Fouriertransformed. I have not been able to obtain an analytic solution and therefore have calculatedthem numerically. The approach is to evaluate each 3d kernel at a specified set of depths(z1,. .. , z,) and then Fourier transform each depth slice. The results are combined to yielda discretized kernel (p, q, zk) (k = 1,.. . , n). This discrete representation is interpolated toprovide a continuous representation of each kernel, which can be used directly in the continuousform of inversion or be further integrated to form the parametric kernels for the discrete form ofinversion.The computations proceed in the following manner. The set of depths, (z1,... , z,), for thekernel computation is spaced at intervals increasing with the depth. At each depth, the spatialkernel for a given i is sampled very finely. The kernel must be sampled frequently enough toavoid aliasing. Ascertaining the appropriate digitization interval is facilitated by the form ofthe kernel functions in the wavenumber domain. Consider the spatial transform of the kernelin Fig. 4.3(a). The kernel has a very broad bandwidth whose amplitude decays smoothlyas the wavenumber increases. A cut-off wavenumber is chosen so that the contribution fromwavenumbers higher than this cutoff point is negligible. It is observed that the envelope ofalong the p-axis for l = (li, 0) can be approximated by the function e_21)/lJ as shown in Fig. 4.4.For a given depth z, the sampling interval is chosen so that the envelope has decayed to a certainfracton at Nyquist wavenumber compared with zero wavenumber, i.e.,rzln,‘CHAPTER 4. APPROXIMATE 3d INVERSION 770.040.03•_ 0.020.01, 0.00—o.o—0.03—0.04___p (ract/m)Figure 4.4 Amplitude change of the wavenumber domain kernel function. The solidcurve shows the wavenumber kernel along the p axis, which corresponds to array =(25, 0)m. The envelope for the kernel along the axis is approximated by the functione_zP/l as shown by the dashed curve.where A is the sampling interval and the specified fraction. is usually taken to be 0.01. Withthis criterion, an increasingly larger sampling interval can be used for larger depth values. Thisensures adequate sampling at each depth level to avoid the aliasing error and at the same timeachieves the optimum computational efficiency.However, A becomes very small for small z and the number of function evaluations becomeslarge. This difficulty is overcome by specially treating the singularities. The function,f - —l 2l(x — l3,)+2l (y— i) +z2g, T,l) =(4l +4l ÷z2)3/21 [ — + ( — +—2l(x —l)—2l (y — i) +z2[(x — l)2-i- (y — l)2+z2]3/closely approximates the singularities for small z and has a known analytic solution for its Fouriertransform. The subtraction ofg3(, 15 from g(, 15 removes the major variation in the amplitudeand hence the high wavenumber components. Therefore, I only sample and perform digitalFourier transform on the residual (g— gj. The analytic transform of g3 is then added back to thetransform of the residual to give the complete Fourier transform of g(, 1).— I I I I0 1 2 3 4 5 6I I ICHAPTER 4. APPROXIMATE 3d INVERSION 78The spatial Fourier transforms of the slices shown in Fig. 4.3 are given in Fig. 4.5. TheFourier transform at shallow depth is characterized by numerous oscillations and large bandwidth. For greater depths, the transforms die out more rapidly with increasing wavenumber. Atany wavenumber, these three plots provide three point evaluations of (p, q, z). The completekernels (y, q, z) corresponding to (p, q) = (0.012,0) are shown in Fig. 4.6. Kernel functionscorresponding to small ij have more structure and amplitude toward the surface than do thosefor larger offsets. Kernel functions for large source-receiver separation decay more slowly withdepth. This is consistent with the fact that an array with a larger separation has a greater depthof penetration. It is noted also that the kernels are smooth functions of depth. Since the arrayseparations are increasing linearly, the kernels for large offsets become similar. Effectively thekernel for l = 200 m provides little new information that is not already in the kernel for l, 175m. That is, the kernels associated with the two largest offsets are almost linearly dependent forthis particular choice of (j, q). The smoothness and lack of structure exhibited by the kernels alsoadds insight into the nature of the inversion results. One would not expect for example that thesekernels could provide any vertical resolution about the conductivity below depths of about 200m.For a fixed array separation, the dominating portion of the kernel function shifts towards thesurface as the wavenumber increases. This can be characterized by the median depth, Zmed, ofthe absolute areaa2(p, q, z) covered by the kernel function. Zd is given byZm©jaj(p,q,00)= / (p,q,z)dz.At zero wavenumber this median depth is in fact the effective depth of pole-pole array definedby Edwards (1977). The decrease of this depth with increasing wavenumber indicates thatEdwards’ effective depth sets an upper bound on the effective depth. It becomes smaller as thelateral resolution requirement increases. For the same reason, the recovered conductivity modelCHAPTER 4. APPROXIMATE 3d INVERSIONNNNFigure 4.5 Wavenumber domain kernel function corresponding for 1= (25,0) at the depthof 1 m, 20m, and 40m. These three slices are the 2d Fourier transforms of the spatialkernels in Fig. 4.3. They provide three point evaluations of the kernel corresponding tothe given 1 used in the all wavenumbers. The complete kernel is obtained by evaluatingall depth slices.79p (rad/ifi)p (rad/m)p (rad/m)CHAPTER 4. APPROXIMATE 3d INVERSION 800.0040.0020.000—0.002—0.004—0.006—0.008—0.0 10—0.0 12—0.014______________________________________500Figure 4.6 A set of kernels at wavenumber (p, q) = (0.012, 0)for arrays in the j4= (li, 0).The curves correspond to arrays with l,, ranging from 25 to 200 m in equal increments of25 metres.becomes smoother as depth increases. Fig. 4.7 shows a set of kernels at different wavenumbersand Z normalized by the value at zero wavenumber.4.4 InversionWith the data and kernel functions evaluated in the wavenumber domain, the solution ofthe inverse problem can now be tackled. According to equation (4.7), at each wavenumber(p, q) a Fredholm equation of the first kind needs to be inverted to recover a complex modeliYQ,, q, z). There exists an infinite number of models which adequately reproduce the flj complexand inaccurate data e3 (j = 1,... , ni). I will formulate the inversion to generate a particularmodel of “minimum” structure with respect to the background conductivity.The mathematical development of the approximate 3d inversion made use of the Born approximation. That approximation is valid when the conductivity differs only slightly from thehalf-space value. For consistency, it is therefore reasonable to compute a conductivity modelwhich deviates as little as possible from this background model. Such a 3d conductivity model200m25mI I I0 100 200 300 400Z(m)I I I0.0 150.0100.0050.000—0.005—0.010—0.0 15—0.020—0.0251.00.90.80.70.60.50.40.30.2Figure 4.7 Panel (a) displays a set of kernel functions for array E=(25, O)m at wavenumber (p, q) = (n14p, O).The curves from top to bottom correspond to ii,, = 0,. . . , 8. Panel(b) is the normalized median depth Zmed of the absolute area under the kernels.is constructed by finding perturbations to a base model defined by the estimated half-space conductivity o, which is used in the data reduction. Therefore a model is sought which minimizes1)=JJJw(z) [ln1u(r)]2dxdydz, (4.8)z, zwhere w(z) is an arbitrary positive weighting function. Using Parseval’s theorem, (4.8) can bewritten as0CHAPTER 4. APPROXiMATE 3d INVERSION 81No40 80 120 160Z(m)2000 1 2 3 4 5 6 7 8NIñi(p,q,z)2dz.CHAPTER 4. APPROXIMATE 3d INVERSION 82In this form, it is observed that the total objective function is a sum of positive quantities.The objective function is effectively decoupled and the total objective function is minimized byminimizing the quantity= JwCz th, q, z) 12 dz (4.9)at each wavenumber.At each wavenumber the data equations (4.7) have the formj=1,”.,njwhere and (y, q, z) are complex and j(y, q, z) are real. The Id inverse problem is then tofind a solution for th minimizing c1(th) in (4.9) subject to fitting data ë (j = 1, . , mj). Either acontinuous or blocky th can be sought as a function of depth z. Each form has merit and I shallpresent both in the following.4.4.1 Continuous FormulationThe most general approach to the above-stated inverse problem is to formulate the inversionso that a continuous function iT(y, q, z) is recovered. This has the advantage that the model doesnot depend upon any specific a priori parameterization. The model and kernels are assumedto reside in a Hubert space M of complex functions defined on the region [0, where Zmaxis some maximum depth of interest and is sufficiently great so that all kernel functions havenegligible amplitude by that depth. If and 7 are any two functions in M, then their innerproduct is given by(,i) = Jw(z)g(z)ii*(z) dz, (4.10)where the * denotes complex conjugate. The norm on M is defined by1/2Zmax=(J w(z) I(z)I2 dz) . (4.11)CHAPTER 4. APPROXIMATE 3d INVERSION 83The data provide only mi constraints upon the model and hence Ti(p, q, z) cannot be obtainedunambiguously from them. There are, in fact, infinitely many models which fit the data exactly(e.g. Backus and Gilbert 1967). Given that the data are also inaccurate, I choose to find thatparticular model which minimizes= IIII2+e1—(4.12)where 8 is a Lagrange multiplier and ° (z) = (z) 1w (z). The minimization is straightforward.Any function iT E M can be written as th = thII + th-’- where ffi. E MI1 = asp {7} and M’is the complement of MI’. It is clear that fr’ does not contribute to the data misfit but it willcontribute to the model norm. As such, a minimization of c1 (th) demands thatth(p,q,z)=a7(p,q,z) (4.13)where a are coefficients to be determined. Substituting into (4.12) and writing the objectivefunction in a matrix/vector notation yields= + 8’ (è—T1(eTa), (4.14)where H denotes complex conjugate transpose and F is a real symmetric positive definite innerproduct matrix with elementsr=and = (a1,. . . , a,)T g = (ëb . . , e,)T. Perturbing in (4.14) and carrying out the minimization using a variational method, subject to the fact that cI(iT) must always be real valued,yields9_lf’H(ë_ F)=O.CHAPTER 4. APPROXIMATE 3d INVERSION 84Since F is positive definite and symmetric this may be written as(T+8I)c=ë°. (4.15)The solution of (4.15) is most easily obtained by writing F = RART, where A is a diagonalmatrix containing the eigenvalues of F in descending order, A = diag. .,, and R is aunitary matrix composed of corresponding eigenvectors R = (i1,. . .,. The solution of (4.15)is then given by= R(A+18)’ RTë. (4.19)The major difficulty in carrying out the inversion is specifying the value of 8 for each idinversion. If the data errors on ë. (j = 1,• . ., ni) were Gaussian, independent with zero mean andknown standard deviation there would be no problem. The procedure would be first to normalizeeach data equation by the standard deviation of the datum and then adjust 8 until a desired misfitvalue is reached. A desired misfit might be E [x2] n.The errors in the above inversion do not comply with the simple Gaussian assumption. Inaddition to measurement errors on the initial potential data, there are several other major sourcesof errors. Errors are introduced into the data by the extrapolation of the known data to fill themathematical area in which the inversion is carried out; errors exist in the model representationassociated with the Id Fredholm equation which is inexact due to the Born approximation; errorsare introduced into both data and model by estimation error of the background conductivity oo.Therefore, the data are biased and non-Gaussian in realistic problems. Nevertheless, I haveattempted to incorporate a Gaussian type strategy as much as possible.The approach has been to first estimate an approximate error for each datum. These will beinterpreted either as approximate standard deviations or at least used as a relative weighting forthe different data to be inverted. The estimation of the errors are obtained in the following manner.The errors on initial potential data in each map are converted to the variance of the wavenumberdomain datum ë,(p, q). Since the initial data errors are assumed Gaussian with zero mean andCHAPTER 4. APPROXIMATE 3d INVERSION 85uncorrelated and the Fourier transform kernel has a unit amplitude, this variance is equal to thesum of the variance of the initial errors. This error is then scaled by the ratio of the total area ofthe inversion to activated area. In the estimation of the initial error, an upper and lower boundon the possible error are imposed. The lower bound is especially useful for those data with smallmagnitude.The above errors have been used in three ways. In the first these values have been acceptedas valid standard deviations, used to normalize the data equations, and 9 is chosen so that the chisquared misfit was equal to nj. In all the test examples, this method failed to work satisfactorily.The amplitude of the data j(p, q) decays rather rapidly with the increasing wavenumber whilethe estimated “standard deviation” stays the same. This results in the domination of the finalmodel by a very few low wavenumber components and the conductivity model is geophysicallyunacceptable.In the second method the estimated standard deviations are used only as relative weightingsfor the data equations and 9 is chosen in accordance with its effect on the model norm. The modelnorm increases monotonically as 9 decreases and is easily evaluated in the spectral domain. Usingthe decomposition of F, we can define a set of rotated basis functions (Parker, 1977),‘(z)= rjj7(p,q,z), (4.17)where r1 is the element of eigenvector . The t,b(z)’s form an orthogonal set, i.e.,(b,&1)= (4.18)where is the Kronecker delta. The regularized solution is given by(4.19)CHAPTER 4. APPROXIMATE 3d INVERSION 86where ê is the th element of the rotated data vector ê = RTè. It follows from the orthogonalityof the i,b(z)’s thatIIthII2=()2(4.20)The curve Ifiz.I 2(9) is usually smooth and often characterized by a near vertical slope forsmall 9 and near zero slope for large 9. The value of 9 corresponding to the onset of rapid normincrease can serve as estimates for 9 used in the inversion. This is implemented by normalizingboth 9 and IthII2 to interval [0, 1]. The value of 9 where the tangent to the normalized curve hasa given slope k (e.g. k = —1.0) yields the desired estimate. This approach has proven very stableand produced models very representative of the true models in synthetic examples. However,without referring to the data misfit, this approach is somewhat arbitrary in the choice of the slopek within this permissible small range. Therefore, when choosing 9 based upon a slope criterion,the value of the slope k is best chosen by incorporating a global misfit as done in conjunctionwith using a single value for9 as presented next.The third approach is an attempt to find the value of 8 according to the misfit to the originaldata in the spatial domain. A single value of 9 is used to regularize all Id inversions. A given 8will generate a model from which the Born data can be calculated. The total chi-squared misfit,evaluated only for the original field data is given by2=(4)2(4.21)where and are observed and predicted total potential values respectively, S is the estimatederror of each datum and t is the total number of the potential data used in the inversion. Thecalculation is efficiently carried out by first evaluating the Fourier transform of the relativeanomaly S4 by (4.7). Applying the inverse Fourier transform yields S4, from which the totalpotential 4?’ is obtained. x2 is the misfit for the approximate equation rather than a true misfit.Carrying out the inversion for a number of values of 8 yields a curve x2 (9) which can beinterrogated to estimate an optimum value or to find that 8 which corresponds to an expected x2.CHAPTER 4. APPROXIMATE 3d INVERSION 87This approach has the advantage that the regularization is controlled by misfit to the original data.It does require that the inversion be carried out a number of times but that poses no computationaldifficulties. The CPU time needed for all id inversions is only a small fraction of that used forcomputing the kernel functions. Once the kernels are computed, they can be stored and additionalinversions can be carried out efficiently.However, the eigenvalues of the inner product matrix F decreases with increasing wavenumber since the amplitude of the kernel function decreases. When a single value of 9 is used toregularize all id inversions, it tends to over-damp the solution at higher wavenumbers. As aresult, the recovered conductivity model has few details.Experiments with synthetic examples indicate that a better model is produced when theregulanzation is applied to each individual inversion based upon the model norm. Therefore, Ichoose to use second approach of regularization discussed in this section.4.4.2 Discrete FormulationIt is often desirable to formulate the inversion so that a blocky model is obtained directly.This is especially advantageous when a numerical modelling scheme such as a finite differencealgorithm is used for forward modelling. Introduce a depth partitioning (z0,... , z) where zo = 0and z, = Zm. A set of boxcar functions can be defined according to the partitioning as the basisfunction for the model th(p, q, Z),(1, z E [zk_j,Zk],bk(Z) = (4.22)1. 0, otherwise.The model can then be expanded asffi(y, q, z) = thkbk(z), (4.23)CHAPTER 4. APPROXIMATE 3d INVERSION 88where the coefficients ñZk are now constants to be determined. To be consistent with the modelbeing sought, a weighting function expandable in {bk(z)} would be reasonable:w(z)=Wkbk(Z),where wk are a set of positive real numbers. The basis {bk(z)} spans a portion of the model spaceM, which can be denoted by Mb = asp{bk}. Since bk(z) are not orthogonal to any of kernelfunctions, Mb intersects with M”. It is in this intersection where the model th(p, q, z) will beresiding.Substituting (4.23) into (4.12) yields the new objective function= + 8 (Gñ — ë)”(Gni — (4.24)where ñ. = (Ti,• . . , ,T’i,j, G = (,. . . , g )T and the discrete kernel = (,... , )T (j =1,.. , n,) is given by(p, q)= / j(p, q, z)dz.Zk_ IW is the weighting matrix and W = diag(wi(zi — zo),. ,w(z — z_1)). Minimizing (4.24)with respect to model parameters ñ. in the similar manner as minimizing (4.14) yields(GTG + 6W)ñ = GTë.Since the weighting matrix W is positive definite and diagonal, W exists and is diagonal.Therefore, the above equation may be written as(GGW + 9I) = (4.25)1where G = W 2G and m = W2m.Equation (4.25) is solved efficiently using the singular value decomposition (SVD) of the matrix G. Let G = USVT, where S is a diagonal matrix of singular values, S = diag(s1,. • , s)CHAPTER 4. APPROXIMATE 3d INVERSION 89and n is the effective rank of the matrix G. U and V are left and right singular vector matrix,respectively. The solution of (4.25) is then given by= VS(S2+ 9J)_I UT e, (4.26)and the final model is obtained by unweighting ñ:-. _i_•m=W 2m,,.As in the continuous formulation, the crucial step is to choose a proper value for 9. I haveemployed the same approaches discussed in the continuous formulation. As can be expected, itdoes not work to take the estimated errors as the valid standard deviations and to misfit the dataaccordingly. The second and the third approaches have worked in this discrete formulation asin the continuous formulation. In the second approach the model norm needs to be computedefficiently for a series of 9 values. The norm of the model expressed by (4.26) is given byflr /—2 p.21 S ‘Mmli =leI j,29) (4.27)where ê are the elements of the rotated data vector ê = UTe. With (4.27), an estimate for 9 canbe obtained by exactly the same procedure described in the continuous formulation.4.5 Calculation of Background Conductivity oWhen the background conductivity is uncertain, a desirable alternative approach is to estimateit through inversion. Note that the data used in Id inversions are the transform of the percentageanomaly Sçf= (ç5—q5,)/4,. Since the primary potential (for unit current strength)14ir1o0is constant for a given array configuration, 00 contributes only to the data at zero wavenumber andacts as a scale factor at all non-zero wavenumbers. Substituting Sçb into the spatial data equationCHAPTER 4. APPROXIMATE 3d INVERSION 90(4.5) and rearranging the terms yields= 1 + 10 ®g(r;i)dz.4irloo 47thT00The zero-wavenumber component after Fourier transforming is, therefore,= + .—!—-!- / ñz(z)(z)dz, (4.28)41rlj00 4irl3uo 0where==and= i JJg(r;i)dxdY.iiS denotes the mathematical area defined for the inversion (see Fig. 4.2). Equation (4.28) hasboth o and 7(z) as unknowns and provides the opportunity to recover both simultaneously.Equation (4.28) can be written= --- +-- / n(z)(z)dz, (4.29)O•0 O•O 0where ë1 = 4rl1ct. Now, (4.29) is a slight variation of the Fredholm equation of the first kindin which 00 appears as an extra unknown. As at any other wavenumber, I look for a weightedsmallest model (z) subject to (4.29) as constraints. For generality, I shall present the solution inthe continuous form in the following. Choose an objective function consistent with the previouslyused form,cI(, II,II2+61 !Lroëj— I — (7Yi,7(z))l2, (4.30)CHAPTER 4. APPROXIMATE 3d INVERSION 91where = g(z)/w(z). 9 is a Lagrange multiplier. Based upon the similar argument, theformal minimization of c1 with respect to ñz. and a0 yields(4.31)(T+90)=o0è— Q, (4.32)andeTa= 0, (4.33)where I’ is again the inner product matrix of kernel functions as defined before, è = (ë,•• . , )T,and Q = (1,.. ., )T• It should be noted, however, that I = (1/Si,... , l/S,1)” when the dataequations are normalized by their estimated errors 83.Equations (4.31) and (4.32) represent a similar solution as that given by (4.13) and (4.15).But (4.33) imposes a further constraint on the kind of solution permissible. Solving (4.32) and(4.33) for oo yieldsT(T÷ 9I)_1° T(r+9I)-1 (4.34)Since the value of 9 is to be determined based upon its effect on either the data misfit or the modelnorm, (4.32) and (4.34) need to be solved simultaneously. One approach would be to start withan initial 9 and solve the two equations iteratively. Once oo is obtained through the id inversionat zero-wavenumber, it can be used in the inversions at non-zero wavenumbers as the estimatedbackground conductivity.The oo obtained through this is the optimum value for the smallest 12 norm ñi(0, 0, z). Thisprocedure amounts to finding the background conductivity which represents the surface data asmuch as possible so that the conductivity perturbation is minimized. A similar technique hasbeen used by Aidridge et al. (1991) in the construction of the flattest models.The question remains as how to calculate the mean of each potential data map which isconsistent with the mathematical area of the inversion. One can use the mean of the data gatheredCHAPTER 4. APPROXIMATE 3d INVERSION 92from the E-SCAN® data set as the mean for the inversion region. Since the area activated by thesedata always dominates the inversion region, such extrapolation is not unreasonable. Alternatively,one can take a two-step approach. The mean value of all apparent conductivities gathered fromthe E-SCAN® data can be used as a first estimate of oo. This first estimate is used to interpolatethe data maps by augmenting them with primary potential on the boundary of the mathematicalarea (see section 4.2). The mean values of these interpolated maps can then be inverted to obtainthe final estimate for o.4.6 Synthetic ExampleThe algorithm has been tested on several synthetic data sets. As an illustration, I presentthe inversion results in the continuous form for one synthetic data set. The data set is generatedover a 21 x 21 grid of 50m spacing for a model consisting of five prisms buried in a uniformhalf-space. Table 4.1 gives the parameters of the model. The perspective view of the modellayout is shown in Fig. 4.8. Prisms Bi and B2 are the buried targets. Prisms Si, S2, and S3simulate near surface variations in the conductivity. 8 potential data maps are gathered in both- and y-directions. The electrode separations are equal to nA (n = 1,.. . , 8), where A is thegrid spacing and equal to 50m. 5 maps are also gathered in both zy- and yx-diagonal directionswith separations equal to n/A (n = 1,. . . , 5). These 26 data maps are used in the inversion(i.e., there are 26 complex data for each id inversion in the wavenumber domain). No a prioriinformation except the background conductivity is given. A weighting function w(z) = 1 /(z + ZO)is applied, where zo = 5m. With such a weighting function, the inversion is effectively carriedout on a logarithmic depth In(z). The regularization parameter 6 for each id inversion is chosenaccording to the norm of the model. A slope of tan(ir/3) is specified.Fig. 4.9 is the comparison of the true and recovered model in a section at x = 450m whichcuts through the four major prisms. Fig. 4.10 and Fig. 4.11 show the comparison at depth z = 30mand z = 1 50m respectively. Geometrically, the definition of the surface prisms are excellent.The image becomes less sharp as the depth increases. The magnitude of the recovered anomalyCHAPTER 4. APPROXIMATE 3d INVERSION 93Table 4.1 Parameters of the 5-prism model.. x-dimension y-dimension z-dimension ConductivityPrism (m) (m) (m) (mS/rn)Si 225-325 225-325 0-40 10S2 275 - 675 375 - 475 0-40 5S3 325-675 625-675 0-40 0.5Bi 275-625 275-375 50-250 0.5B2 375-475 475-675 75-275 10(0,0,0) X’AXISy-f15- -NFigure 4.8 The perspective view of the synthetic model consisting of five prisms embedded in a uniform half-space of ImS/m. The surface prisms SI (5mS/m), S2 (lOmS/m),and S3 (0.5mS/m) simulate the near surface variations in the conductivity. The Bi(0.5mS/rn) and B2 (lOmS/m) are the buried targets.CHAPTER 4. APPROXIMATE 3d INVERSION 940.100.200.300.400.500.0.100.200.F300.400.500.0. 200.Figure 4.9 The comparison of the true modelsection x=450 m (scale in log10 cr).600. 800. 1000.1000.—2.13—2.28—2.39—2.52—2.65—2.78—2.91—3.04—3.17(a) and the recovered model (b) in thedecreases with the depth. In general, the conductivity of the buried blocks is under-estimated.Because of the loss of resolution and magnitude at the depth, the bottoms of the buried prisms arepoorly defined. Instead, there is a smooth transition from the anomaly to the background. Theloss of resolution is due to the limited bandwidth of the data and the fact that the kernel functionin the wavenumber domain decays more rapidly with depth at higher wavenumbers. The loss ofamplitude is due to the nature of smallest model construction, where a ridge regression parameteris used to regularize the solution. Overall the recovered model appears as a depth varying filteredversion of the true model. Nevertheless, all anomalies are well resolved and despite the influenceof the surface conductivity variations the buried prisms are clearly defined.0. 200. 400.Y (m)400. 600. 800.X (m)CHAPTER 4. APPROXIMATE 3d INVERSION 950.200.400.600.800.1000.0.200.400.600.800.1 000.0. 200. 400. 600.X (m)800. 1000.—2.13—2.26—2.39—2.52—2.65—2.78—2.91—3.04—3.170. 200. 400. 600. 800. 1000.X (m)Figure 4.10 The comparison of the true model (a) and the recovered model (b) at depthof 30 m (scale in log10 o).0.200.400.600.800.1000.0.200.400.600.800.1000.0. 200. 400. 600. 800. 1000.X (m)Figure 4.11 The comparison of the true model (a) and the recovered model (b) at depthof 150 m (scale in log10 o).CHAPTER 4. APPROXIMATE 3d INVERSION 96—2.13—2.26—2.39—2.52—2.65—2.78—2.91—3.04—3. 170. 200. 400. 600. 800. 1000.X (m)CHAPTER 4. APPROXIMATE 3d INVERSION 97It is noteworthy that such good results have been obtained even though three conductivityblocks with high contrasts are placed at the surface. This is apparently in violation of theassumption that the surface conductivity be constant and the conductivity deviation is small. Thisshows that the algorithm is quite robust.4.7 DiscussionA 3d approximate inversion algorithm has been developed based upon the integral equationfor surface pole-pole potential anomaly. The algorithm is designed to work for general 3d modelscomposed of perturbations to a background conductivity. Synthetic examples show that it canhandle fairly large conductivity contrasts and, to a certain degree, surface variations. The approachI have taken has some general implications. The algorithm can be viewed from different aspects.By applying a linear operator (2d Fourier transform) to the basic data equation, the 3d inverseproblem is broken into a set of independent Id inverse problems. Each I d problem can be solvedeasily and efficiently. This results in an algorithm which is much more efficient than solving theoriginal problem directly.From the view point of inversion, this approach also has merits. The Fourier transformseparates the large scale features (low wavenumber components) and small scale features (highwavenumber components) in the data. Features of different scales are fit by finding the modelcomponents of corresponding scales. The horizontal variability in the conductivity is controlledby the largest wavenumber included in the inversion. Therefore, the method only attempts togenerate features in the model which are required by the data. Variations beyond the resolutionprovided by the data are excluded. This process makes the inversion stable and less prone toproducing spurious structures in the model.The parameter o only acts as a scaling factor for the data at non-zero wavenumbers, butprovides an additive term to the data at zero wavenumber. Thus an inaccurate choice of oOwill only scale the model components at non-zero wavenumbers and the structural image is notaltered in the horizontal directions. When the error on the estimated oO is large, it will produce anCHAPTER 4. APPROXIMATE 3d INVERSION 98incorrect mean value (from the inversion at zero wavenumber) of the model at each depth. Thiswill result in a series of false layers in the spatial domain. The final result would be the correctstructural image (with scaled amplitude) superimposed on these false layers. Because of this fact,a small variation in the estimate of o does not significantly alter the resultant conductivity model.It should also be noted that the application of any transform introduces an extra level ofcomplication in data processing. This can cause difficulties. For example, in order to apply theFourier transform, the observed potential data must be interpolated and extrapolated to a regulargrid. This introduces additional errors whose statistical properties are unknown. Consequently,the statistical properties of the transformed data are uncertain. If the original data are aliasedupon acquisition, the transformed data will also be incorrect.As the algorithm is based upon the Born approximation, its result is semi-quantitative. Therecovered conductivity model is not likely an acceptable model since it does not reproduce theobserved data via a full forward modelling algorithm. Therefore, the result of the approximate3d inversion is best regarded as an image. However, the success with the test examples providesa base for optimism for this approximate inversion. The next step will be to input the recoveredconductivity model into a full 3d forward modelling algorithm to quantify the misfit betweenobserved and predicted data. The approximate inverse mapping can then be incorporated into aniterative algorithm to yield a quantitative and rigorous solution to the 3d resistivity problem. Thisleads to the application of the AIM (Approximate Inverse Mapping) inversion to the 3d resistivityproblem using the 3d approximate inversion as an approximate inverse mapping.CHAPTER 5AIM INVERSIONThe interpretation of a data set often requires that an acceptable conductivity model beconstructed. The inversion of d.c. data is a non-linear problem. The most common approach forsuch problems has been to first linearize the problem (by Taylor expansion) and then to constructthe model iteratively using Newton’s method. This requires that, at each iteration, the sensitivitymatrix (or the Fréchet derivatives) be computed so that the change in the data can be mappedto a change in the model. In the cases of a Newton-type method, the matrix needs also to beinverted. These operations involve a great amount of computation and therefore pose difficultiesin practical applications.However, the 3d approximate inversion is a very efficient method which maps a set ofpotential data to a conductivity model without resorting to the generation and inversion of a3d sensitivity matrix. If an algorithm is constructed to utilize such an inverse mapping, it willhave a great computational advantage. The formalism of AIM (Approximate Inverse Mapping)inversion developed by Oldenburg and Ellis (1991) provides the framework for such an approach.In the AIM inversion, an approximate inverse mapping is used in conjunction with an exactforward mapping to update the model successively by finding either a perturbation to the modelor a modification for the data to be inverted. The algorithm iterates towards an acceptable modelwithout necessarily computing and inverting a large sensitivity matrix. In this chapter, I shall firstbriefly review the AIM formalism and then present its application to 3d d.c. resistivity problemusing the approximate 3d inversion.5.1 Review of the AIM FormalismSuppose a geophysical experiment is described by a forward operator F which maps anelement m in model space to an element ë in data space,F[m]=ë. (5.1)99CHAPTERS. AIM INVERSION 100For the inverse problem of this experiment, one would ideally like to have an inverse operatorso that it maps the data ëback to the model m,.F’[ëJ=m. (5.2)If the problem is non-linear, it is most commonly tackled by first linearizing it using mathematicaltools such as a Taylor series expansion. By such linearization, the Fréchet derivatives (orsensitivity matrix) are used to relate the change in the data to the change in the model. The modelis updated by applying the inverse of the matrix to the difference between the observed and thepredicted data. Because of the intensive computations involved, such approaches are marginallyapplicable on present day workstations to large scale problems in the 3d environment. In order toovercome this difficulty, Oldenburg and Ellis (1991) developed the AIM (Approximate InverseMapping) inversion as an alternative approach. The AIM inversion utilizes an approximate inversemapping and an exact forward mapping to update the model successively without resorting to thefull sensitivity matrix. Thus it can avoid the need for intensive computations.The philosophy of the AIM inversion is to break up a non-linear inverse problem accordingto the different levels of physical processes involved, instead of using mathematical tools. Thisapproach enables the physical understanding of the problem be introduced optimally into thederivation of the solution. It is recognized that most geophysical problems can have their physicsdivided into primary and secondary (or residual) processes. For example, the charge accumulationdue to the primary electric field would constitute the primary physical process (whose field isgoverned by the Born approximation). The interaction between the accumulated charges would bethe secondary process. Once the primary process is identified, an approximate inverse mapping,denoted by P’, can be defined accordingly. The division between the primary and secondaryphysics is problem dependent. They can be chosen so that the primary physics has either alower dimensionality or a reduced non-linearity and its inverse P’ can be realized with lesscomputational effort. Thus one can use the approximate inverse mapping to obtain that portionCHAPTER 5. AIM INVERSION 101of the model which can be seen “through” the primary process and find the remaining portionfrom the difference between the exact and approximate mappings.There are two ways of utilizingP to update the model. The first seeks a model perturbationat each iteration so that the updated model will reproduce the observed data according to F. Thesecond attempts to find a modification to the data so that, when the approximate inverse mappingis applied to the updated data, the resulting model will reproduce the observed data. Oldenburgand Ellis refer to these two approaches as AIM-MS and AIM-DS, respectively, for the reasonthat the perturbation is sought in either model space (AIM-MS) or data space (AIM-DS). Thereare several ways to derive the recursive formulae for the AIM algorithms. In the following, Ishall present a very brief derivation. More detailed discussions can be found in Oldenburg andEllis (1991).5.1.1 AIM-MSIn the AIM-MS, a sequence of model perturbations is sought to update the model successively.Let ë’ be the observed data. Let m be the model and ë°” be the predicted data at the thiteration,= (n) (5.3)Application of the approximate inverse mapping P to ë and è(1 yields ‘h°” and th(’). Thatis,P [ë0b5] = (5.4)[(fl)] = (n) (5.5)Since F’ is only approximate, its application to the data produced by a model m according to(5.1) will not recover the original model. The difference is quantified byA[m] = (Im — P’F)[m], (5.6)CHAPTER 5. AIM INVERSION 102where Tm is an identity mapping on model space. A[m] is called the mapping error on modelspace. Suppose m is the model which reproduces the observed data e°”, then by definition,m = th°1’ + zi[m].Since m is the unknown to be found, the above equation has to be solved iteratively by replacingthe m on the two sides of the equation by the models at two successive iterations,= + ii[m]. (5.7)Combining (5.6) and (5.7) and substituting in (5.3) and (5.5) yields,= + m — (5.8)Equation (5.8) is the final formula for the iterative solution of m. The iteration starts with aninitial model which can be supplied by iii°’.The equation (5.8) can be viewed in two ways. Firstly, it can be written asm+D = m + (thobs —and shows that the term (th°1 — th) constitutes a perturbation at (n + I )th iteration, which isthe difference between the models obtained by applying the approximate inverse mapping to theobserved data and the data predicted by the current model. Therefore, the model is built up witha series of perturbations. Secondly, as (5.7) shows, the final model is made of two parts. The firstpart, fi°1’, is the portion derivable directly using the approximate inverse mapping. The secondpart is a corrective term to account for the incompleteness of the approximate inverse mapping.In fact, the operator t = Tm — PF quantifies the portion of the model which is annihilated bythe approximate inverse mapping P.It is necessary for the algorithm to converge that the sequence of perturbations(W1’th) —.40 as the iteration progresses. Equivalently, the mapping error i[mJ should approach a constant,CHAPTERS. AIM INVERSION 103i.e., (A[m’]— [m”]) —* 0. This requires that the approximate inverse mapping be stable sothat a small change in the data space is mapped to a small change in the model space.5.1.2 AIM-DSIn the AIM-DS algorithm, a sequence of data perturbations is sought so that the applicationof the approximate inverse mapping to the updated data will generate a model which reproducesthe observed data. For clarity, I present the derivation in parallel to that of AIM-MS. Againdenote the observed data by é*0b1 and let m be the model and è be the modified data at thenth iteration, thus,m =F1[ë]. (5.9)Since P is only approximate, a model m obtained by applying P to data element e will notreproduce that data element through the forward mapping. The difference is quantified by= (Z— .FP)[e1, (5.10)where I is an identity mapping on data space. A[ë] is called the mapping error in the dataspace. Suppose m is a mode] derived by applying P’ to the modified data ë, and m adequatelyreproduces the observed data ë°1”. Then by definition,e= sobsSince we seek to form the modified data to recover the model m, this equation must be solvediteratively. Replacing on the two sides with the modified data in two successive iterations yields(n+I) = oba +[‘], (5.11)which can be further expanded by substituting (5.3), (5.9), and (5.10),(n+I) = + (n) — (5.12)CHAPTERS. AIM INVERSION 104This is the final formula for iteratively updating the data in AIM-DS algorithm. Applying theapproximate inverse mapping F yields the updated model in each iteration. The inversionagain can start with any reasonable initial data However, the choice = e°1 is appropriate.There is an alternative derivation for AIM-DS, which relies upon the explicit definition ofan approximate forward mapping and offers some intuitive understanding of the algorithm. Iformally introduce the approximate mapping P which describes the primary physical process andsatisfies P-’P = I. ThusP[m] = (n) (5.13)Suppose m is the model solution reproducing the observed data then the data equation canbe expressed asobi= P[rn] + (.1 — P)[m]. (5.14)That is, the data consist of the part predicted by the primary physics and a residual from thesecondary physics. Rearrange the terms in (5.14) intoF[m] = + F[m] —The unknown m is to be solved iteratively by replacing it with m’ and m” on the left andright hand sides respectively. Substituting in (5.13) yields the same equation as (5.12),(n+l) = (n) + oba —Thus the AIM-DS approach is equivalent to an iterative inversion where the data are successivelyupdated by the difference between the exact and approximate responses from the current model.Equation (5.12) can be viewed in either of the following two ways. The modified data areupdated by adding to the current data “ the perturbation term which is the difference betweenthe observed data and the predicted data, (e0l — e). Therefore, the data are built up by asequence of perturbations. Alternatively, the data are composed of the observed data ë°1 and aCHAPTERS. AIM INVERSION 105corrective term given by the mapping error in data space ( (n) — e). In order for the AIM-DSalgorithm to converge, it is necessary that the perturbation (e° — —* 0. Therefore it isrequired that the difference between the data mapping errors at successive iterations approachzero. This means that ë’ + z[ë (n)j must converge to a fixed point in the data space.As the iteration progresses, the models from both AIM-MS and AIM-DS will likely tohave increased structure. In the AIM-MS algorithm, the model is built up by a sequence ofperturbations. The final model cannot be regularized explicitly in general. In the AIM-DS, thedata are built up with a sequence of perturbations thus they will exhibit more and more structures.This increased structure will be translated to the resulting model. However, one advantageof the AIM-DS is that the final model is always obtained by a single application of the inversemapping. Therefore, it is possible that a norm of the model can be minimized so that the structuralcomplexity of the model be controlled to a certain extent. If both the approximate forward andinverse mappings are linear, it can be shown that AIM-MS and AIM-DS degenerate to the samealgorithm. This is an uncommon occurence since the regularization in the inversion renders themapping non-linear. However, there are situations in which the inverse mapping is nearly linear.Then the results from the two algorithms will not differ significantly.5.2 AIM Inversion of 3d d.c. Resistivity DataIn this section, I apply the AIM formalism to the 3d d.c. resistivity problem and developan iterative inversion algorithm to invert E-SCAN® data sets. The two essential components forforming an AIM inversion algorithm are an approximate inverse mapping and an accurate forwardmapping. The approximate 3d inversion presented in Chapter 4 is chosen as the approximateinverse mapping P. For the forward mapping, F, a 3d finite difference modelling algorithm isused.The finite difference forward modelling algorithm presented by Dey and Morrison (1979)is chosen. The conductivity model represented by a finite rectangular region is discretized intoprismatic cells by a 3d orthogonal mesh. Each cell is assigned a constant conductivity value. TheCHAPTERS. AIM INVERSION 106algorithm can work with general 3d conductivity models but the current source must be placedon a nodal point. The resulting potential is calculated at all other nodal points. The current sinkis assumed to be at infinity.In order to deal with the finite domain, a mixed boundary condition is adopted on the verticaland the bottom sides of the rectangular region, which simulate the infinitely distant planes of thereal model. The potential on those surfaces are assumed to have the general form of 4(e) = A/rwhere A is a constant, and r is the length of the position vector, assuming that the coordinateorigin is at the centre of the grid on the surface. By differentiating the potential with respect tothe outward normal of the bounding plane, the mixed boundary condition is derived to be+cos8 = 0, (5.15)on rwhere 9 is the angle between the normal vector and the position vector .The discrete equation for potentials is formed by integrating the governing differential equation (2.12),V (oV) = —I6(r— it),in a volume around each node. The volume extends halfway into the eight cells adjacent to thenode. The resultant coefficient matrix is banded and is solved using a conjugate gradient method.Several tests have been made to compare this finite difference code with other forward modellingcodes. For some models and some data I found discrepancies of 10%, but usually discrepancieswere less than 5%.The approximate 3d inversion as discussed in Chapter 4 can take either a continuous ordiscrete form in the vertical direction. In theory, they can both be used in the AIM application.A discrete form is chosen here because a single vertical grid can then be used for both theapproximate inverse mapping and the forward mapping. This avoids rediscretization of theresults from inverse mapping for input of the forward modelling.CHAPTERS. AIM INVERSION 107Two different meshes are used for the forward mapping and for the approximate inversemapping. For the forward problem, the mesh consists of a core portion and an extension region.The E-SCAN® survey grid is used as the horizontal grid of the mesh in the core. A verticalpartitioning is then chosen. It extends to a maximum depth so that the kernel functions for theapproximate 3d inversion all have negligible amplitude at that depth. Surrounding this core mesh,an extension mesh is added so that the boundary condition in the forward modelling algorithm canbe handled. This is usually a three—cell extension with increasing intervals towards the boundary.The mesh for the inverse mapping is chosen with the horizontal grid conforming to thatdescribed in Chapter 4. That is, the horizontal area consists of the original survey grid at thecentre and is extended by one half of the maximum array separation in each direction (Fig. 4.2).The grid spacing is chosen to be half of the original survey grid spacing. A vertical partitioningis chosen to coincide with that of the forward modelling mesh.The resultant conductivity model of inverse mapping attains a discrete point representation inthe horizontal directions by the nature of the fast Fourier transform (FF1’). Beneath each horizontalpoint, the model possesses a constant value along the depth within each vertical interval of themesh. The logarithmic conductivity ln(u) is assumed to be represented by piecewise bilinearinterpolations in the horizontal directions. A conductivity value is then derived from the integratedaverage of ln(o) within each cell of the forward mesh. This value is assigned to the correspondingcell as the recovered conductivity value.With the forward and inverse mapping and their respective meshes being defined, the nextstep is to define the data and the model of the AIM inversion. Although P takes the relativepotential anomaly as data and produces a conductivity model to which the forward mappingF is applied, the data and the model for the entire process of AIM inversion can be defineddifferently. For AIM-MS, it is better to define the model as the logarithm of the conductivityvalues. Therefore, the model for AIM-MS inversion is composed of the set of the logarithmicconductivity values ln(oJk) (i = 1, . . , M, j = 1,.. . , M, k = 1,.. . , Mi), where M, M, Mare respectively the number of cells in the x-, y-, and z-direction. Working with logarithmicCHAPTER 5. AIM INVERSION 108conductivities is advantageous because the variation in the earth’s conductivity commonly rangesover several orders of magnitude and also it ensures the positivity of the conductivity model. Inaddition, it is consistent with the logarithmic conductivity perturbation sought in the approximateinverse mapping. Thus, the iterative formula (5.8) takes the formln(o’’) = ln(&°1’)+ ln(u) — ln(&), (5.16)Given the total potential qj(x, y) associated with n pole-pole array configurations, the datafor the AIM-DS inversion is defined as the logarithm of the potential values, ln(ç$3(x, y)) (j =1,... ,n,). The formula (5.13) becomesIn(1+D)= 1n(’) + ln(q50b6) — ln(q’), (5.17)Adopting the logarithmic potential as the data eases the task of data interpolation and ensures thatthe updated potential is always positive.The AIM inversion of 3d d.c. resistivity data proceeds by the following steps:AIM-MS:1. Apply the approximate 3d inversion (P’) to the observed data to obtain &°. Let= u° and generate the predicted data 1) by applying the forward modelling.2. Apply P’ to to obtain model3. Update the current model by (5.16), ln(o’) = ln(&0b1) +ln(o) — ln(&’), to obtain a newmodel (n)4. Apply F to compute the predicted data5. Compute the misfit between the observed data qS and the predicted data q5. Exit thealgorithm if the misfit is acceptable. Otherwise, go to step-2 and begin a new iteration.AIM-DS:1. Take the observed dataq as the initial data and apply the approximate inverse mappingto produce model Compute the predicted data by applying forward modelling toCHAPTERS. AIM INVERSION 1092. Update the data by (5.17), ln() = ln() + ln(q5°6’)— ln(4), to obtain new data3. Compute the new model o by applying the inverse mapping to4. Compute the predicted data from the new model.5. Calculate the data misfit. Exit the algorithm if the misfit is acceptable. Otherwise, go tostep-2 and begin a new iteration.5.3 Examples of the AIM InversionTo illustrate the inversion of the E-SCAN® d.c. resistivity data using the AIM formalism,both AIM algorithms are applied to two sets of synthetic data. These examples demonstrate theconvergence of the algorithm, the resolution of the algorithm given the grid spacing, and thedepth of investigation given the grid spacing and the maximum array separation.The synthetic data are computed from two cellularized conductivity models using the finitedifference forward modelling algorithm. The distribution of the computed data simulates thatof the field data over a 21 x 21 grid with a 50m spacing in both the - and y-directions. Thefirst model is a 5-prism model similar to that used in Chapter 4 (Fig. 4.2). The locations of thefive conductivity anomalies are changed slightly so that their boundaries are aligned with the cellboundaries used in the finite difference program.The second model consists of correlated random perturbations to a uniform half-space.Two sets of random numbers with different spatial correlation lengths in the 3d space are firstgenerated. A weighted sum of the two sets are used as the conductivity perturbation 1n(t) toform the model. The component with longer correlation length is tapered to zero near the surfaceso that the model simulates the increased spatial variation in the conductivity near the surface.The resultant conductivity model is then converted into prismatic cells by integrations. Thecellularized conductivity constitutes the “true model” in this example and will be referred as thegeo-statistical model.CHAPTER 5. AIM INVERSION 1105.3.1 The 5-prism modelI now proceed to the results of AIM inversion applied to these examples. For the 5-prismmodel, I choose a vertical grid of 21 nodes in the core portion with intervals increasing with thedepth. The horizontal grid is the same as the data grid. A three-cell extension zone is appendedto the core mesh. The AIM-MS algorithm is first applied to the data set and twelve iterationsare performed. The regularization of the Id inversions in the wavenumber domain is based uponthe norm of the model in each individual inversion. (This method of choosing regularizationparameter is used for all the examples in this chapter.) Throughout the inversion, the half-spaceconductivity has been held at the correct value of lmS/m. Fig. 5.1 shows the comparison ofmodels at the section x = 460m, which cuts through the four major conductivity anomalies in themodel. The gray scales in this and all the following figures are in log10( ). The top panel is thetrue model. The middle panel is the model after first iteration (i.e., after the simple applicationof the approximate inverse mapping), and the bottom panel shows the recovered model fromthe 12th iteration. Fig. 5.2 shows the similar comparisons at two depth slices at z = 20m andz = 150m respectively. The first iteration recovers an image of the conductivity model. Thesurface anomalies and the buried conductive anomaly are clearly defined but the buried resistiveanomaly is not defined very well. The amplitudes of the anomalies are far less than that of the truemodel. However, further updates by the AIM algorithm greatly improves the result. The finalmodel as shown in the bottom panels of the two figures recovers all the conductivity anomaliesand has a dynamic range comparable with that in the true model. In terms of the model recovery,the inversion is quite successful.Fig. 5.3 attempts to illustrate the convergence property of this particular inversion. The toppanel is the RMS misfit of the total potentials at each iteration. The RMS misfit is adopted becausethe data are accurate. The half-space model (Oth iteration) and the model from the approximateinverse mapping have a misfit of 12% and 4.2%, respectively. As the iteration progresses, themisfit is steadily reduced. By the end of the l2 iteration, an RMS misfit of 0.5% is achieved. Itis noted that the misfit decreases monotonically with the iteration. The middle panel of Fig. 5.3CHAPTERS. AIM INVERSION 111Figure 5.1 Comparison of the 5-prism model and the AIM-MS inversion result in sectionx =460m. The upper panel is the true model, the middle panel is the model from the firstiteration, and the lower panel is the final model.0.100.200.300.400.500.0. 200. 400. 600. 800.Y (m)1000.E1-0..0.100.200.300.400.500.—2.13—2.26—2.39—2.52—2.65—2.78—2. 91—3.04—3.170. 200.0.100.400. 600. 800.Y (m)1000.200.300.400.500.0. 200. 400. 600. 800.Y (m)1000.CHAPTERS. AIM INVERSION 1120.200.400E>600.800.—2.131000.—2.260.—2.39200.—2.52400.—2.85600.—2.78800.—2.911000.—3.040.—3.17200.400.600.800.1000.z=20m z=150m0. 200. 400. 600. 800. 1000.X (m)0. 200. 400. 600. 800. 1000.X (m)0. 200. 400. 600. 800. 1000.X (m)0.200.400.S680.800.1000.0.200.400.S>-600.800.1000.0.200.400.S600.800.1000.0. 200. 400. 600. 800. 1000. 0. 200. 400. 600. 800. 1000.X(m) X(m)Figure 5.2 Comparison of the 5-prism model and the AIM-MS inversion result at depthsz = 20 and 150m. The true model and the models from the first and last iteration areshown from top to bottom.0. 200. 400. 600. 800. 1000.X (m)CHAPTER 5. AIM INVERSION— 120100= 800 2 4 6 8 10 12113Figure 5.3 The change of the RMS data misfit (a), the model norm (b) and the norm ofthe model perturbation (c) with the iteration in the AIM-MS inversion for 5-prism model.shows the norm of the model ln(())II, where p() = o(i)/o-o. The bottom panel shows thenorm of the model perturbation at each iteration. The model norm increases monotonically andapproaches an asymptotic value. The norm of the model perturbation decreases in a similarmanner and approaches a negligible value compared with that of the model norm (The ratio ofthe two quantities is 0.0 12 at the last iteration). Thus the necessary condition outlined in thederivation of AIM-MS algorithm that the model perturbation approach zero is met. Together withthe decrease of the misfit, this tends to suggest that this particular application of the AIM-MSinversion is stable and on a path towards convergence.1200 2 4 6 8 10 123020100 2 4 6 8 10 12ITERATION NUMBERCHAPTERS. AIM INVERSION 114Next, the AIM-DS inversion is applied. The mesh is the same as that used in the AIM-MSinversion. Again, the inversion has been run for twelve iterations. Fig. 5.4 is the final AIM-DSmodel at slices x = 460m, z = 20m, and z =1 50m. The recovery is good in terms of both locationsand amplitudes of the anomaly. It is noticed that this conductivity model is very similar to thatobtained by AIM-MS inversion. A detailed comparison finds that the difference between therecovered conductivity values from AIM-MS and AIM-DS can differ by up to 20%. The normof the AIM-DS model is also slightly smaller than that of the AIM-MS model.Fig. 5.5 is again to illustrate the convergence properties of the inversion. The top panel in thefigure shows the RMS misfit of the total potential. Similar to the result in the AIM-MS inversion,the reduction of the misfit is along a monotonic curve. The final misfit achieved is below 0.5%.It is worth noticing that the major reduction of the misfit is achieved with first few iterations inboth the AIM-MS and the AIM-DS inversions. Once the misfit is reduced below a level ( 1.0%in this case), further iteration does not improve it significantly.The middle and the bottom panel ofFig. 5.5 display respectively the norm of the modified dataand the norm of the data perturbation S at each iteration. Here the data are the relative anomaly,(c — 4o)/4’o, since that is the quantity used in the inverse mapping. The norm of the modified dataincreases monotonically towards an asymptotic value while the norm of the perturbation decreasesmonotonically and levels off at a negligible fraction of the data norm (2.8%). This is very similarto the behavior of model norm and the norm of the model perturbation in AIM-MS. For this idealcase of accurate data, the result conforms with the condition outlined in the derivation of theAIM-DS algorithm that the perturbation of the modified data approach zero.The above tests have used only the accurate data. A more realistic test, however, is touse data that are contaminated with noise. For the d.c. resistivity experiment there are twomajor types of noise, namely the geological noise and experimental errors. The former is mostlycaused by features such as near surface variations in the conductivity. This type of noise mostlyconcentrates in the high wavenumber band and can be suppressed in the inversion by limitingthe highest wavenumbers. The second type of noise is the error associated with all stochastic‘0001•008‘009•001’‘0001008(in)x‘009‘001’•00‘0‘0001woct=zpun‘vgg=z‘111o9v=soogsaqiAloAnoadsaJamwonoqodoiwoijsioundaqjiapowwspd-aqi.iojiinsaiuO!5JOAU!su-wiVp•çain2j(in)x‘0001‘008‘009‘001’‘00‘01’LI‘C—to‘C—16’3—81/3—38’3—6C’3—$93,3—CI‘3—tn‘001:z‘003$ ‘001citNOISI3ANIMLIVcKILLJVH3‘008‘009‘001’‘003(at)A‘0001‘008‘009‘001’‘003‘0‘00s‘001’CHAPTER 5. AIM INVERSION12— 8014T 12100 2 4 6 8 10 121164.0—3.02.01.00 2 4 6 8 10 12ITE]?/i TION NUMBERFigure 5.5 The change of the RMS data misfit (a), norm of the modified data (b) and thedata perturbation (c) with iteration in the AIM-DS inversion for the 5-prism model.processes involved in the data collection. Such noise is often uncorrelated and distributed overthe entire wavenumber band. It is often simulated by adding to the data independent randomnumbers with prescribed standard deviation and mean.For the data from the 5-prism model, I added independent Gaussian noise with zero mean andstandard deviations which amount to 5% of the corresponding total potential. A more realisticsimulation requires that there be a minimum for the standard deviation of the noise so that theerrors on the data from large array separations can be modelled. However, the separations of thearray in this numerical experiment only span the short range of 50 to 400 meters. Consequently0 2 4 6 8 10 12CHAPTERS. AIM INVERSION 117the measured potential varies with the separation only in a small range and adding a percentageis reasonable.The synthesized inaccurate data are input into the AIM-DS algorithm and eight iterations arerun. The x2 misfit of the total potential is reduced below the expected value after 4 iterations.Since the data now have genuine errors on them, the x2 misfit can be computed. Fig. 5.6 showsthe model at the 4th iteration. It is immediately clear from the figure that the noise in the datahas introduced some spurious structures into the resultant conductivity model. Such distortion isespecially strong near the surface. At depth, the shape of buried resistive prism, which is a veryweak anomaly, has been altered significantly. However, all anomalies in this particular model arerecovered and their definitions are clear.Because only the norm of the model is minimized in the inversion, there is no explicit controlover the structural complexity of the model. Consequently, the spurious structures introducedby the noise cannot be suppressed directly through the regularization in the inversion. However,these structures seem to be uncorrelated and become progressively weak as the depth increases.This is judged to be a direct result of the approximate inverse mapping. From Chapter 4 we knowthat the kernel functions in the wavenumber domain decays more rapidly with depth at higherwavenumbers. Meanwhile, the energy of the data is concentrated in the lower wavenumber bandwhereas the energy of the noise is spread throughout the entire band. Thus at lower wavenumber,the signal-to-noise ratio is much higher and the recovered model components have less spuriousfeatures. Structures at depth are dominated by these low wavenumber components. On theother hand, the noise component in the data becomes stronger as the wavenumber increases.The corresponding model component becomes more and more affected by the noise. Since thekernel functions are sensitive to shallower regions at higher wavenumbers, the constructed modelcomponents are non-zero in the progressively shallower region as the wavenumber increases andthe recovered conductivity model in the spatial domain has more noisy features near the surface.This also suggests that these undesirable features can be treated by applying a low-pass filteringto the model with the cut-off wavenumber decreasing with the depth. The type of filter and theCHAPTERS. AIM INVERSIONS—2.13Figure 5.6 AIM-DS inversion result of the noisy data from 5-prism model. The panelsfrom top to bottom are respectively the slices at x =460m, z =20m, and z = 150m.1180.100.200.300.400.500.0. 200. 400. 600.Y (m)800. 1000.0.200.400.600.800.1000.0. 200. 400. 600.X Cm)0.—2.28—2.39—2.52—2.85—2.78—2.91—3.04—3.17800. 1000.S200.400.600.800.1000.0. 200. 400. 600.X (m)800. 1000.CHAPTERS. AIM INVERSION 119cut-off wavenumber will no doubt depend upon the spectrum of data and the noise and, hence,will be problem dependent.Fig. 5.7 gives some details on the progress of the inversion as a function of iteration. The firstpanel shows the x2 misfit at each iteration. The dashed line indicates the expected value of x2which is equal to 8804. Again we observe a steady reduction of the misfit. After four iterations,the x2 misfit has dropped below the expected value. Further iterations do not result in significantimprovement. In fact from the 6th iteration there is a slight increase in x2• This is to be expected.It has been shown (Bertero eta!., 1988) that the degree of regularization is inversely proportionalto the number of iterations for an iterative inversion without explicit regularization. Such is thecase with the AIM inversion as applied here.The norm of the modified data (panel b) linearly increases with iteration number. The norm ofthe data perturbation (panel c) is first on a decreasing curve and then starts to increase from the 6iteration. This turning point seems to coincide with that of the x2 misfit curve and indicates that theinversion has begun to diverge. This suggests that, in practical applications, the inversion shouldbe terminated once the expected misfit is achieved. If the expected misfit cannot be reached,the onset of the increase in data perturbation should serve as the criterion for termination. Sincethis specific implementation of the AIM inversion does not explicitly minimize the model normsubject to a prescribed data misfit, the prolonged iterations increase the structural complexity ofthe model but do not necessarily improve the data misfit.It is interesting that the norm of the model increases smoothly with decreasing rate up to thelast iteration in spite of the near-divergent behavior of the modified data. Careful examination ofthe modified data reveals that the large increase of the data norm is mostly due to the outliers inthe data where the added noise is strong. These outliers seem to be enhanced by data updating.However, their energy is spread over the wavenumber band and the regularization based on themodel norm tends to suppress their effect so that the recovered model is relatively stable.CHAPTERS. AIM INVERSION30205.2Figure 5.7 Progress of the AIM-DS inversion applied to noisy data from 5-prism model.Panel (a) shows the change of the misfit with the iteration. The dashed line indicates theexpected x2 misfit. Panels (b) and (c) are the norm of the modified data and the norm ofthe data perturbation respectively. Panel (d) shows the change of the model norm.1200402 4 6 8100 26.05.64 6 84.80 2 4 61401201008080 2 4 6 8ITERATION NUMBERCHAPTER 5. AIM INVERSION 1215.3.2 The geo-statistical modelThe 5-prism model clearly illustrates both methods of the AIM inversion and demonstratestheir ability to derive a conductivity model which fits the data within the error tolerance. It alsoshows the performance and the complications in the presence of the noise in the data. Howeverthe model is simple and the conductivity variations are limited in a relatively shallow region. It ishoped that the geo-statistical model represents a geologically more realistic test which can providesome insight into the performance in more complicated situation. For this model, I choose aninversion mesh similar to the one used for the 5-prism model, but with the core portion beingpartitioned uniformly with depth at 25 m intervals.The ATM-MS inversion is first applied. The inversion is terminated after twenty-two iterations. The recovered model after the 1 6th iteration is shown in Fig. 5.8 and Fig. 5.9 in comparisonswith the true model in four slices at different depths. The three panels in each column of thefigures display, respectively from top to bottom, the true model, the models at the first and 16thiteration. The first column in Fig. 5.8 is the conductivity in the first layer (between depths 0 to25m). The second column is the third layer (between 50 to 75m). The conductivity in theselayers is recovered with great fidelity, although fewer details are recovered in the third layer.The dynamic range of the recovered conductivity is very close to that of the true model at thesedepths. The model from the simple application of the inverse mapping is plotted to demonstratethe improvement achieved by the iterative process.The two columns in Fig. 5.9 show respectively the comparisons in layer-6 (between depths125 to I 50m) and layer- 10 (225 to 250m). It is clear that both the detail and the amplitude of therecovered conductivity have decreased considerably. As depth increases, we only see larger scalefeatures with smaller dynamic range present in the recovered model. However, strike directionof the major features are still evident and this could be important in practical applications.In general it is again observed from the four slices, that the resolution of the model decreasesrapidly with the depth. The resultant model appears to be a smoothed version of the true modeland the smoothing is increased with the depth. A question of both theoretical and practicalCHAPTERS. AIM INVERSION 1221000.0.layer-i0. 200. 400. 600. 800. 1000.X (m)0. 200. 400. 600. 800. 1000.X (m)layer-30. 200. 400. 600. 800. 1000.X (m)Figure 5.8 Comparison of the geo-statistical model and the AIM-MS inversion results inlayer-i (z = (0, 25)m) and layer-3 (z = (50, 75)m). The panels from top to bottom in eachcolumn show respectively the true model and the models at the first and 16th iteration.0.200.400.600.800.1000.0.200.400.600.0.1000.200.400.0. 200. 400. 600. 800. 1000.X Cm)EBB0.200.400.600.800.600.—1.17—1.34—1.51—1.68—1.85—2.02—2.19—2.36—2.53800.1000.0.0. 200. 400. 600. 800. 1000.200.400.600.800.200.400.600.800.1000.1000.CHAPTERS. AIM INVERSION 123layer-6200. 400. 600. 800. 1080.X (m)0. 200. 400. 600.X (m)0.200.layer- 100. 200. 400. 600. 800. 1000.X (m)0. 200. 400. 600. 800. 1000.X Cm)0. 200. 400. 600. 800. 1000.X (m)—1.17—1.34—1.51—1.68—1 .85—2.19—2.36—2.53Figure 5.9 Comparison of the geo-statistical model and the AIM-MS inversion results inlayer-6 (z=(l25, 150)m) and layer-b (z=(225, 250)m). The panels from top to bottomin each column show respectively the true model and the models at the first and 16thiteration.200.400.E600.800.0.200.400.600.800.1000.1000.0.0. 200. 400. 600. 800. 1000.X (m)200.400.E200.600.800.400.600.800.1000.0.0.—2 . 02E200.400.600.800.1000.400.1000.800. 1000.CHAPTER 5. AIM INVERSION 124importance is: “What is the maximum depth above which the algorithm can recover someconductivity variations from the given data set?” This question is closely related to the questionabout the depth penetration of the surface d.c. resistivity experiment and its answer is far moreinvolved. Parker (1984) demonstrates that the depth penetration of id d.c. resistivity datais infinitesimal without additional constraints. This has direct bearings on the 3d experiment.However, if the value of the conductivity is limited to a finite range and arbitrarily rapid variationsare restricted, the experiment should still have a finite depth of penetration. The limited Fouriercomponents and the finite partitioning in the vertical direction used in the AIM inversion imposesuch constraints implicitly. Therefore, I attempt to provide a qualitative evaluation.Fig. 5.10 and Fig. 5.11 display the comparison in a cross-section of the true and recoveredmodel. The cross-section is taken at z =420m and is plotted at the scale of the true model. Overthe depth extension of the model, the inversion result seems to be a very poor representationof the true model. However, close examination shows that the recovery above the depth ofapproximately 200m is actually rather good. Below this depth, oniy the gross features causedby conductive anomalies are present. The recovered model exhibits few features below 300m. The Fig. 5.11 is a detailed comparison in the region above 200m. It is apparent that therecovered model and the true model agree well. Thus, the maximum depth of the model recoveryis estimated at approximately 200 m, one half of the maximum array separation. This result is inagreement with the observation at the end of Section 4.3 that the kernel function does not exhibitsignificant structure below 200m for an array with 400m separations. To further verify the abovemaximum depth, the conductivity below 200 m in the true model is muted to the background valueand the new responses are computed by the same forward modelling program. The RMS relativedifference between the new and original data is 4.5%. The contribution to the difference increaseswith the array separation. This suggests that the inability to recover conductivity variations belowthe maximum depth is partly inherent to the data set.Fig. 5.12 shows the change of misfit, norm of the model and model perturbation. The RMSmisfit exhibits the familiar steady reduction up to iteration- 16 (this is the iteration whose model is200.E300.200.300.Figure 5.10 Comparison ofresult (b) in section x =420m.depth of 200m.the geo-statistical model (a) and the AIM-MS inversionThe recovered model only represent the true model aboveCHAPTERS. AIM INVERSION 1250.100.——1. 17—1.34—1.51—1.68—1 .85—2.02—2.19—2.36—2.53400.OO.0.1 00.0. 200. 400. 600. 800.Y (m)1000.400.500.0. 200. 400. 600. 800.Y (m)1000.—. 0.%40.80.B— 120.fr.1 160.200.0.40.80.E-4 120.ç 160.200.0. 200. 400. 600. 800.Y (m)1000.0. 200. 400. 600. 800.Y (m)CHAPTERS. AIM INVERSION 126—1.40—1.50—1.60—1.70—1.80—1.90—2.00—2 101000.—2.20Figure 5.11 Detailed comparison of the geo-statistical model (a) and the AIM-MSinversion result (b) at x = 420m above z = 200mplotted in the preceding figures). An RMS misfit of 0.5% is achieved. However, the misfit startsto increase slightly beginning from iteration-17. So does the norm of the model perturbation. Thenorm of the model increase approximately linearly with the number of iteration after the initialrapid increase. This indicates that the inversion is starting to diverge.Next the AIM-DS inversion is applied. Twenty-two iterations are performed. Unlike theAIM-MS results, both the misfit and the norm of the data perturbation are still decreasing. For aconsistent comparison, I again plot the model after 16 iterations. Four slices at the same depths asbefore are shown in Fig. 5.13. They are to be compared with the true and final recovered modelin Fig. 5.8 and Fig. 5.9. The AIM-DS model is very similar to the AIM-MS model; differencesoccur only in the fine details. Calculation shows that the AIM-DS model has a slightly smallernorm. Fig. 5.14 shows the misfit reduction and the change of the norm of the data and the data20‘0CHAPTERS. AIM INVERSION 1271612cij800 5 10 15 20140130= 1101000 5 10 15 200 5 10 15 20ITERATION NUMBERFigure 5.12 The change of the data misfit (a), the model norm (b) and the norm of themodel perturbation (c) with the iteration in the AIM-MS inversion for the geo-statisticalmodel.perturbation. It is noticed that, at the later iterations, the norm of the modified data again increaseslinearly as does the norm of the model in the AIM-MS inversion.This behavior of the norm of model and data is different from that in the case of the 5-prismmodel, where we saw convergence towards an asymptotic value in both model and data. This isdue to the fact that the current model has more anomalous structure and the approximate inversemapping based upon the Born approximation deteriorates. However, this may not be a concernin practical applications, since the misfit is reduced to below 1.0% by the fourth iteration in bothcases. Further iteration has not reduced the misfit significantly but only increased the complexityof the model.CHAPTERS. AIM INVERSION 128layer-i layer-30.200.400.600.800.1000.0.200.400.600.800.1000.0.200.400.600.800.1000.0. 200. 400. 600. 800. 1000.X (m)layer-6—1.17—1.34—1.51—1.68—1.85—2.02—2. 19—2.36—2.530.200.400,E600.800.1000.0. 200. 400. 600. 800. 1000.X (m)0. 200. 400. 600. 800. 1000.X (m)Figure 5.13 AIM-DS inversion result for the geo-statistical model.CHAPTER 5. AIM INVERSION 1291612CI)80181614ITERATION NUMBERFigure 5.14 The change of the data misfit (a), norm of the modified data (b) and the dataperturbation (c) with iteration in the AIM-DS inversion for the geo-statistical model.To complete the example, I add 5% Gaussian independent noise to the data and apply bothAIM-MS and AIM-DS inversion. After four iterations, both algorithms produced models whichreduced the x2 misfit to the expected value. The two models are very similar, although theAIM-MS model has a slightly larger norm. Fig. 5.15 displays the AIM-DS model in four slicesat the same depths as before. It is evident that there are small scale features introduced by theadded noise. However, the first layer is not as noisy as that obtained by inverting the noisy datafrom 5-prism model. This difference is related to the difference in vertical partitioning. Theinversion mesh is uniform in the vertical direction for the geo-statistical model whereas it wasfinely partitioned near the surface for the 5-prism model. The very thin cells near the surface are0 5 10 15 200. 200. 400. 600. 800. 1000. 0. 200. 400. 600. 800. 1000.X(m) X(m)Figure 5.15 The result from AIM-DS inversion of the noisy data from geo-statisticalmodel.more prone to the effect of the noise in the data, since they allow more structure in the verticaldirection. This has direct bearings on the application to field data. Overly fine partitioning shouldbe avoided in the mesh. Fig. 5.16 shows the reduction of the x2 misfit and the change of themodified data, its perturbation, and the model in terms of their respective norms. Similar behavioris observed as in the case of 5-prism model. The norm of the data increases linearly but the modelnorm is more stable.CHAPTERS. AIM INVERSION 1300.layer-i layer-3200.400.E1-600.800.1000.200.400.600.800.1000.0. 200. 400. 600. 800. 1000.X (m)layer-60.200.0. 200. 400. 600. 800. 1000.X (m)layer- 10400.—1.17—1.34—1.51—1.68—1.85—2.02—2.19—2.36—2.53600.800.200.400.600.800.1000.1000.CHAPTERS. AIM INVERSIONITERATION NUMBER131Figure 5.16 Progress of the AIM-DS inversion applied to the noisy data from the geostatistical model. Panel (a) shows the change of the misfit with the iteration. The dashedline indicates the expected x2 misfit. Panels (b) and (c) are the norm of the modified dataand the norm of the data perturbation respectively. Panel (d) shows the change of theioio0 2 4 6 84030205.6= 5.24.8130120— 110//d0 2 4 6 8model norm.CHAPTERS. AIM INVERSION 1325.4 DiscussionThe formalism of the AIM inversion is successfully applied to the 3d d.c. resistivity inverseproblem and both AIM-MS and AIM-DS iterative inversion algorithms are constructed to invertE-SCAN® data for a 3d conductivity model. The model is assumed to have a factorized formof o(i) = o-oj(), where o is the background conductivity and i’() is a scale factor. Thealgorithm employs the approximate 3d inversion developed in Chapter 4 based upon the Bornapproximation, which generates the 3d inverse solution by decomposing the problem into asequence of id problems in the wavenumber domain. A finite difference forward modellingalgorithm is used as the forward mapping. Both accurate and noisy data from two syntheticmodels are used to test the AIM-MS and AIM-DS inversion algorithms. In all test cases thealgorithms have succeeded in producing acceptable models. These models adequately reproducethe data within the tolerance of the associated errors, and represent the true models with reasonablefidelity.In general, the AIM-DS is advantageous in comparison with AIM-MS inversion. AIM-MSgenerates the final model by a series of model perturbations and there is no explicit control overthe norm of the model. As a result, the model from AIM-MS inversion can acquire excessivestructure. The AIM-DS inversion, however, generates the final model from a single applicationof the approximate inverse mapping. A norm of the model is explicitly minimized to reduce thestructural complexity. Thus AIM-DS can potentially produce a better model and is preferable inpractical applications. This prompted me to invert the data sets from two synthetic model withboth AIM-MS and AIM-DS algorithm. I had expected that the models from the two algorithmswould be different. However, the results show that the difference is not large. It is only notedthat the norm of the model generated by AIM-DS is slightly smaller than that from AIM-MS andthe conductivities of individual cells between the two models can differ up to 20%. Beside thesedifferences, the two algorithms have generated almost identical results. The similarity betweenthe results of the AIM-MS and AIM-DS inversions suggests that the approximate 3d inversionmay be nearly linear when regularized according to the model norm.CHAPTERS. AIM INVERSION 133The algorithm has combined the advantages of the approximate 3d inversion and the AIMformalism and therefore is very fast. The kernel functions for the approximate 3d inversionare computed at the first iteration and stored for subsequent iterations. As a result, the inversemapping at each iteration is very rapid, since its solution in the wavenumber domain takes only asmall fraction of the time needed for computing kernel functions. In addition, the AIM formalismrequires one forward mapping at each iteration. This reduces the computation effort at eachiteration to a minimum. This is the advantage of the algorithm.The iterative application of the approximate inverse mapping has been able to constructmodels which fit the data. The reduction of the misfit to the expected level is achieved rapidlywithin the first few iterations. The corresponding model is greatly improved compared with thatfrom a simple application of the approximate inverse mapping. However, the algorithm has itslimitations. The misfit cannot be further reduced after the initial reduction in the first iterations.It is not guaranteed that these first few iterations will achieved the expected misfit. Prolongedexecution of the inversion only increases the structural complexity of the model. Therefore, themethod is best used as a means to generate a conductivity model which reduces the data misfitto a low level with minimum computational efforts. The resulting model can be used to drawgeologic conclusions or it may be used as an initial model for more rigorous inversion based uponlinearization.CHAPTER 6ANALYSIS OF FIELD DATAIn this chapter, a set of field data is processed and inverted using the AIM inversion algorithmdeveloped in the preceding chapter. These data are acquired by a survey designed to assistexploration for epithermal type deposits. I first describe the field procedure of data acquisitionand the data format pertaining to theE-SCAN experiments. Next I proceed to processing the datain preparation for the inversion. Then I will present the inversion results and make a comparisonwith the geology in the area based upon the information available.6.1 Geologic ObjectivesThe survey which acquired this data set is a part of an exploration effort for epithermaldeposits. They are a class of the hydrothermal mineral deposits. In general, epithermal depositsare considered as having been formed “by hot ascending waters of uncertain origin, but chargedwith igneous emanations. Deposition and concentration of ore minerals occur at slight depth”(Lindgren, 1933). Such deposits most commonly occur within the upper 600 m in areas withwell developed fracture and fault systems. Existing faults or fractures act as feeders of hot fluidwhich alters the host rock and deposits minerals. Consequently, ore is formed in the alterationzones, which often expand in the upper portion to form cone-like or mushroom-shaped featureswith deep “roots”. Precious metal deposits are frequently associated with silicification broughtabout by hydrothermal alteration.Because of the highly resistive nature of the veins and silicified zones, the ore deposits exhibit localized resistive anomalies as typical signatures in the d.c. resistivity experiment. TheE-SCAN® survey is therefore employed to map such silicification zones and possible fault structures. The area is covered with eluvium on the surface. As a result, there are few surfaceexpressions of the geology in this region but it is hoped that the E-SCAN® experiment willprovide information about the subsurface conductivity so that ore bodies can be located.134CHAPTER 6. ANALYSIS OF FIELD DATA 1356.2 Data AcquisitionE-SCAN® experiments, especially when applied to mineral exploration, are usually carriedout on a regular grid with constant grid spacings in each direction. The spacing can be the samefor both directions. Theoretically, each grid point is occupied in turn by a current source andthe surface potentials are measured within a prescribed radius around the source location. Themeasurement of the potential within that radius is over a specified pattern of the grid points. Sincethe potentials measured with slightly different source-receiver orientations at a large separationdo not provide very different information, omitting some of these data points will not reduce theinformation content of the data set significantly. Such omission is usually for economic reasons.The “infinite” current and potential electrodes B and N are placed at sufficiently long distancesaway.The field data acquisition is controlled by a digital computer. Each electrode at the grid pointshas an electronic switch box attached to it. Along each line these switch boxes are connectedin series and the end of each line is connected via a switch box to a main cable leading to thecentral control. Each box has a unique address within this network. The central control is thusable to instruct each individual box either to connect to the attached electrode (or the attachedline for those on the main cable), or to act as a relay connecting further down the network. Withthis setup, the controlling computer can measure the potential at any node within the grid. Theenergizing current is transmitted through a different cable and the current electrode is usuallymoved to each site manually. The waveform of the source current is the standard periodic squarewave with reversing polarities. The total potential is measured just before the source is cut off.The experiment also measures the IP (induced polarization) potential during the decay periodafter the current is cut off.Since the measurement radius is usually a fraction of the distance covered by the entiregrid, the wire network is laid out only over a group of adjacent lines in a portion of the grid.As the survey progresses, new lines are added in front of the network and lines along whichmeasurements have been taken are removed. To reduce the measurements needed in the field, theCHAPTER 6. ANALYSIS OF FIELD DATA 136reciprocity of the d.c. field is utilized so that potential measurements are made only in the forwarddirection. That is, if the current source is moving from left to right on the line sweeping throughthe grid from top to bottom, the potentials are measured only on the grid points to the right of thesource location and below the source line. This process is explained in Fig. 6.1.current source moving direction—00——- - - -———-/-- ----7:::::zzz__Figure 6.1 Field acquisition procedure of the E-SCAN® experiment. The solid gridrepresents the survey grid. The circles represent the grid points activated by measuringnetwork. For any current source location A, potential data are collected only within aprescribed radius (shown by the large circle) at the nodes in front of the source. Thesenodes are shown in the shaded area. The potential at the nodes behind the source, forinstance at P, are substituted by the reciprocal data, i.e, the potential measured at A whenthe current is at P. This group of activated nodes sweeps through the grid and the currentsource moves from one end to another so that all grid points are occupied by the sourcein turn.I—CHAPTER 6. ANALYSIS OF FIELD DATA 137The solid lines in Fig. 6.1 represent the E-SCAN® survey grid. The grid points marked bycircles denote activated electrodes (connected into the network). The current source is marked bythe solid dot. As the survey continues from the top of the grid to the bottom, the line below thisgroup will be activated and the line on which the current source presently lies will be deactivatedwhen the current source has occupied the last point on this line. So the group moves downwardas indicated by the arrow on the left side. The circle centred at the current source marks theboundary in which potential data are collected for this particular source location. The activatedgrid enclosed by the dashed line consists the points in the “forward” direction. Therefore, thefield measurements are made only within the intersection of the circle and the region bound bythe dashed line, which is highlighted by the shading. This area of measurements is often extendedso that potentials are acquired to the array extremities on the x- and y-axes and 45° diagonaldirections.The particular experiment in consideration was carried out on a 32 x 35 grid with two smallareas missing on the margin. The grid has a spacing of 91.44 metres and covers an area of2.5km x 3km. Fig. 6.2 shows the grid layout. The solid lines represent the E-SCAN® grid. Thegrid points marked by crosses constitute a primary grid with a grid spacing of 182.44 metres. Themajority of the data are collected on the primary grid following the procedure described above.Fig. 6.3(a) shows the pattern of points at which the potentials are measured for each currentlocation. The remaining grids are used to acquire supplemental measurements so that near offset data are collected. Only current sources are placed on these points and their potentials aremeasured on the primary grid. For these sources, the measurement pattern is specially designedand varies with the type of source locations in the grid. Only a small number of potentials aremeasured for these sources. The exact patterns are shown in Fig. 6.3 by panels (b), (c), and(d). For these special cases, all the potentials are measured in both the forward and backwarddirections.The lines AA’ and BB’ in Fig. 6.2 mark the positions of two sections whose apparentconductivity pseudo-sections and the conductivity section from the constructed model will beCHAPTER 6. ANALYSIS OF FIELD DATA-1500B1380 500 1000X (m)1500 2000 2500Figure 6.2 Grid layout of the field data set. The solid lines denote the E-SCAN® gridand grid points marked by cross form the primary grid. Potentials are measured only onthe primary grid. The coordinates are reproduced from the field positions. The lines AA’and BB’ indicate respectively the position of two sections to be discussed. They are tobe denoted by section-A and section-B, respectively.studied. These sections are chosen to coincide with two geological sections compiled fromdrillhole information. They help to provide some limited verification for the inversion result.Each datum is a stacked result of measurements taken with both polarities over a numberof current supply cycles. Usually the mean value is taken. Ideally, a standard deviation of theassociated error should be estimated from the repeated measurements. However, the data setfrom the field only has a numerical indicator of the noise for each datum. These indicators arequalitative summaries of the measurement repeatability, signal strength, and the background noisesuch as the strength of cultural interferences and telluric currents. These indicators are useful in0-500-1000•B’AE A’-2000-2500-3000BCHAPTER 6. ANALYSIS OF FIELD DATA 139H•1— — H H —p— -— -—-1—— HH———H -——H ———H——H H—(a) on primary grid—HL--H H——H--H(c) on vertical line of primary gridr. I-——--“H HH H(b) on horizontal line of primary gridL— —HHH(d) not on primary gridFigure 6.3 Four different patterns of data collection. The solid lines are the primary gridand the dashed lines are the remaining grid lines. The source is marked by the cross andits location is indicated by the labels. The circles are the nodes at which the potentialsare measured.identifying bad data points, but they cannot be used as standard deviations of the data error dueto the lack of quantitativeness.Each datum in an E-SCAN® data set consists of 19 quantities. These include, in the order ofthe file format:1. serial number of the datum2. grid type indicating the direction of the grid axes3-5. (x, y, z) coordinates of the current source A6-8. (x, y, z) coordinates of the potential site MCHAPTER 6. ANALYSIS OF FIELD DATA 1409-li. (x, y, z) coordinates of the current sink B (at “infinity”)12-14. (x, y, z) coordinates of the reference potential site N (at “infinity”)15. the transmitted current strength16. measured total potential17. measured IP potential (in percentage of the total potential)18. noise indicator for the total potential19. noise indicator for the IF potentialThese 19 quantities provide the complete information about each datum. These data recordscollectively define the survey grid and they also indicate the sequence in which the experiment iscarried out.6.3 Data ProcessingThe field data set in the above format needs to be processed before it can be used in theinversion to construct a conductivity model. The first problem is the topographic distortion.Since the inversion algorithm assumes a flat earth surface, any significant topographic relief inthe survey area will be translated by the inversion into spurious conductivity structures beneaththe surface. Fortunately, the topographic relief in the survey area is small. Any distortion dueto such small topography can be effectively treated as geological noise and be suppressed in theprocess of inversion by limiting the high wavenumber components. Thus, the data set is treatedas being collected from plane surface and the inversion is applied directly.The “infinite” electrodes, B and N, in this particular experiment are placed over 7km fromthe centre of the grid. This is sufficiently distant for measurements from short source-receiverseparations, but they are relatively close in comparison to large separations and have significantcontributions to the measured potential. When the pole-pole data are formed, this contribution tends to reduce the potential and to introduce spurious large scale apparent conductivityhighs. These spurious features will lead to false structures in the constructed conductivity model.CHAPTER 6. ANALYSIS OF FIELD DATA 141Therefore, the data must be corrected first for such biases before they can be used to recover aconductivity model. However, the difficulty is in that the correction relies upon the knowledgeof the conductivity structure in the region. Thus, I only attempt to make the first order correction.A best fitting half-space model oo of the region can be derived from the data set. At this stage thedata should be correctly treated as measurements from four-electrode configurations (the data sethas complete information on the electrode locations). A value of 0.04S/m is estimated for thisregion. The data set is then corrected for the “infinite” electrodes using this value by subtractingaway the contribution of the B and N electrodes over the half-space. The results show that thecorrections can be as large as 20% and they are mainly for the large off-set data. Once thiscorrection is done, the data are dealt with as true pole-pole potentials.For most data processing, it is necessary to have a complete “common current gather” forevery current site in the data set. Such a gather consists of all measurements of the potentialgenerated by a common current source. The direct measurement from the field experiment onlyprovides half of the gather as a result of the acquisition procedure. Therefore, the reciprocity ofthe d.c. potential field is invoked to complete the gather. As shown in Fig. 6.1, suppose A is thesource location and P is a point without a direct measurement. This datum is substituted by thepotential measured at A with a source at P. Such measurements are called indirect measurements.Fig. 6.4 illustrates the composition for one current site, where the circles and squares representthe direct and indirect measurements, respectively.With the completed common current gathers, the data can be easily sorted into differentformats for subsequent use. For example, data maps for pole-pole arrays can be synthesizedand the apparent conductivity volume can formed. The pole-pole data maps also constitute thedata for 3d inversion. 26 data maps with good areal distributions are gathered in the horizontal,vertical and 45° diagonal directions. The corresponding pole-pole separations range from 1 to 10basic grid spacings, which are sufficiently large for the interested depth. These 26 maps are to beused in the AIM inversion.CHAPTER 6. ANALYSIS OF FIELD DATA 142H -1 H—1_... --1-- -.----H--HH--HH.._..rt__...,-.----H--HH4-HH--F--- - - -. F- -- F - -- - - -..J_ .4.. .I...H—-pH+F+.--—-.- -. - -. F-1 -- F - -. -. - --— — —- H — — — - — — — — — — — — —. ..•r rt,_.__.-I- -I— — —- H — —— H — — — H—---—---H-- --HH — — — — — — — —H H —Figure 6.4 Composition of potential measurements in a common current gather. Thesource location is marked by the cross. The direct and reciprocal measurements areindicated by circles and squares respectively.There is a highly resistive anomaly in the northwest corner of the grid. This is caused byresistive volcanic tuff, but it has no bearing on the structures of interest in the rest of the region.Because its magnitude is much higher than that of other anomalies, it could have an adverse effecton the inversion by masking the signal from those anomalies. Therefore, the potential at the pointsassociated with this anomaly have been replaced by the mean value from the rest of data pointsin each map. This replacement should not cause problems as long as no significance is attachedto the conductivity in this region in interpreting the conductivity model from the inversion.It is also clear that some data maps are very noisy. In the extreme cases, the data maps havemany isolated spikes. As mentioned earlier, each datum does not have an estimated error soerroneous data cannot be isolated on that basis. The potential data are to be Fourier transformedand any erroneous data points in the spatial domain will affect all data in the wavenumber domain.To improve the quality of the data, I smooth each data map (with the resistive anomaly replaced)using spline interpolation. The spline algorithm constructs a smooth representation, f(x, y), ofCHAPTER 6. ANALYSIS OF FIELD DATA 143the data, which minimizes the functionalq) = _ (f(x, y) — ej2+ JJ ((L)2 + + (_)2) dxdy,where e is the quantity to be smoothed, which is usually the logarithm of the total potential. ). isthe tradeoff parameter determining the degree of the smoothing. Since the data error is uncertain,the generalized cross-validation (GCV) technique (Craven and Wahba 1979) is used to estimatethe optimum value of ). The smoothed data are then used in the inversion as the observed data.Fig. 6.5 shows two apparent conductivity maps formed directly from the raw data. Thecorresponding arrays are in x-direction and the electrode separations are 1 and 4 grid spacings.Fig. 6.6 displays the same two maps after the aforementioned processing steps. The resistiveregion at the upper left corner of the grid in Fig. 6.5 is the anomaly which has been replaced. It isevident that the GCV smoothing has improved the data quality by suppressing the outliers. Fig. 6.7shows two apparent conductivity pseudo-sections in section-A and section-B, respectively. Thepseudo-sections are formed from the arrays aligned with the orientation of the section and theprocessed data are used.6.4 InversionThe AIM-DS algorithm is applied to the processed data maps to construct a 3d conductivitymodel. Again, a dual mesh system is used. The mesh for the inversion is designed as described inChapter 5. The area of the inversion is obtained by expanding the survey area by half maximumarray separation (457.2 m) in horizontal directions. The horizontal mesh over this area has aspacing equal to half of the survey grid spacing. The vertical mesh has a uniform partitioningwith a spacing of 30 metres in the interval [0, 600]m. Three cells with increasing spacing are thenappended to extend the mesh down to 900 metres.For the forward modelling mesh, the horizontal grid has, in the core portion, the field surveygrid of 32 x 35 with 91 .44m spacing (the two missing areas are filled in). Three cells withincreasing spacings are appended to four sides to complete the grid. The vertical grid coincidesCHAPTER 6. ANALYSIS OF FIELD DATA 1440.—622.—1244.—1865.-2487.—3109.—274.0.—622.—274. 293. 859. 1426.X (m)—1244.E-1865.—2487.-3109.1993. 2560.—1.15—1.20—1.25—1.30—1.35—1.40—1.45—1.50—1.55Figure 6.5 Two apparent conductivity data maps from the raw field data. Panels (a) and(b) correspond to arrays aligned in the x-direction and with separations equal to 91.44mand 365.76m respectively.C.,859. 1426.X (m)CHAPTER 6. ANALYSIS OF FIELD DATA—3109—274.Figure 6.6 Two processed apparent conductivity data maps from the field data set. Panels(a) and (b) correspond to arrays aligned in the x-direction and with separations equal to91.44m and 365.76m respectively. These two maps are to be compared with those in1450.-622.—1244.2—1865.—2487.—3109.0.-622.—1244.2—1865.—2487.-274. 293. 859. 1426.X (m)1993. 2560.—1.15—1.20—1.25—1.30—1.35—1.40—1.45—1.50—1.55293. 859. 1426.X (m)1993. 2560.Fig. 6.5.79.220.362.L 503.645.786.220.362.503.645.786.—274. 293. 859. 1426. 1993.X (m)79.2560.—3109.CHAPTER 6. ANALYSIS OF FIELD DATA 146—1.14—1.18—1.22—1.26—1.30—1.34—1.38—1 42—1.46Figure 6.7 Two apparent conductivity pseudo-sections from the processed field data.Panels (a) and (b) show respectively the section-A and section-B. Each pseudo-sectionis formed from pole-pole arrays aligned with the section.with the partitioning in the inversion mesh. This results in a mesh with 38 x 41 x 24 nodal points.Thus the AIM inversion has a total of 34,040 unknowns.As described in section 6.2, in the original survey the current source is placed at every gridpoint in turn but the potential is measured only over the primary grid. Because of the reciprocityin the d.c. field, this data set is equivalent to one in which the potential is measured over the entiregrid while the current source is placed only at the primary grid locations. Thus I treat the originalpotential sites as the source locations and compute the potentials at the original source sites.This reduces the computation required for the forward modelling by half, since the operation isproportional to the number of the source locations.-2487. —1865. —1244. -622. 0.Y (m)CHAPTER 6. ANALYSIS OF FIELD DATA 147353o25ITERATION NUMBERFigure 6.8 Reduction of the data misfit in the AIM-DS inversion of the field data set.An AIM-DS inversion is carried out for 6 iterations. A weighting function w(z) = 1/(z + ZO)is applied, where zO = 20m. The half-space conductivity is fixed at O.0526S/m, which is thebest fitting half-space model from the 26 data maps. The half-space model has a RMS misfit of36%. This initial inversion achieved a lowest misfit of 12.3%. However, the misfit curve exhibitsstrong oscillations as shown in the Fig. 6.8. This indicates that the modified data ë are likelyover-corrected. A second attempt is made in which a relaxation parameter is introduced to limitthe stepsize of the data update.The quantity (ë°’ — e) in equation (5.13),(n+I) = (n) + (*ObS —can be considered to provide both the direction and the stepsize for the perturbation. We canperturb the data in this direction but control the stepsize. That is, the data can be updated by,=+ w(ë°—where w is the relaxation parameter. A choice of w less than unity should prevent an overcorrection of the data to a certain extent and make the iterative process more stable. A similarrelaxation parameter can be applied to the model perturbation in the AIM-MS algorithm.In the second attempt, a relaxation with w = 0.6 is applied to the data update. The stabilityis apparently improved. This is indicated by the smooth decrease of the misfit. Fig. 6.9 shows0 1 2 3 4 5 6ThCl)6040760720680640600Figure 6.9 Reduction of the data misfit in the AIM-DS inversion with a relaxationparameter. (a) is the RMS misfit of the predicted data. (b) and (c) are respectively thenorms of the modified data and the data perturbation. (d) is the norm of the model. Notethat both data misfit and data perturbation start to increase after 4 iterations.CHAPTER 6. ANALYSIS OF FIELD DATA 14830200 1 2 3 4 5 616.015.014.013.0i*0 1 2 3 4ITERATION NUMBER5 6CHAPTER 6. ANALYSIS OF FIELD DATA 149the reduction of the misfit, and the change of the model, data and data perturbation measuredin their respective norms as the inversion progresses. It should be noted that the achieved bestfit to the observed data is about the same for both inversions, although the stepsize of the dataupdate is different. It is also noted that the major reduction of the misfit is achieved in the firsttwo iterations. The remaining iterations have little effect on the misfit, but the model norm doesincrease considerably.The inversion is apparently on a path towards divergence. This is indicated by the increaseof the data misfit at the last iteration as well as the commencement of the rapid increase in dataperturbation after four iterations. Since after four iterations the data misfit does not improve, theinversion should be terminated at this point. If the desired misfit is not achieved or is unknown,the onset of the increase in data misfit, data perturbation (for AIM-DS inversion) or modelperturbation (for AIM-MS inversion) can serve as the criterion for termination. The prolongedexecution of the inversion here is to examine the convergence properties of the algorithm. Inkeeping with this argument, I take the model from the 4thi iteration as the final conductivity modelfrom this inversion.Fig. 6.10 shows two apparent conductivity maps predicted by this model. Fig. 6.11 shows twopseudo-sections of the predicted data. These are to be compared with the observed data shown inFig. 6.6 and Fig. 6.7. It is seen that the large scale features in the field data are reproduced quitewell. However, some fine details are different. Noticeably, the predicted data have smoothercontour lines and fewer isolated small peaks. These isolated points may be contributing a greatdeal to the data misfit. Thus the structural fit to the data might be better than that indicated by theRMS misfit.Fig. 6.12 shows the constructed conductivity model at section-A and section-B. The recovered conductivity model has many structural details at shallow depth but becomes increasinglysmooth at greater depth. This is similar to the results from synthetic examples presented inChapter 5. These conductivity sections can be compared with the apparent conductivity pseudosections displayed in Fig. 6.7. The benefits of carrying out the inversion can now be seen. TheFigure 6.10 Two apparent conductivity maps from the data predicted by the constructedconductivity model. To be compared with the maps in Fig. 6.6.CHAPTER 6. ANALYSIS OF FIELD DATA 1500.-622.—1244.—1.152—1.20—1865.—1.25—2487.—1.30—3109.-274. —1.35—1.40—1.45—1.50—1.55I I I1426.X (m)0.—622.—1244.—1865.—2487.—3109.-274. 293. 859. 1426.X (m)1993. 2560.79.220.362.503.645.786.220.362.503.645.786.—274.79.293. 859. 1426. 1993. 2560.X (m)—3109.CHAPTER 6. ANALYSIS OF FIELD DATA 151—1.14—1 . 18—1.22—1.26—1.30—1.34—1.38—1 42—1.48Figure 6.11 Two apparent conductivity pseudo-sections from the data predicted by theconstructed conductivity model. To be compared with the sections in Fig. 6.7.model obtained from the inversion displays many structures not visible or clear in the apparentconductivity pseudo-sections.Known geology in this area seems to support some features observed in the conductivitymodel. The layer of eluvium covering the area is resistive. This may have been correctlyrepresented in the model by the resistive layering near the surface. Borehole information in thesection-A indicates that there is a fault between x = I 000m and x = I 200m. The sharp conductivitycontrast at the similar location (with the steep contour lines) and the structural change across thisposition may reflect the existence of the fault. In the section-B, the sharp conductivity contrastneai y = —400m seems to indicate the position of a fault known from borehole information.—2487. —1865. —1244. —622. 0.Y (m)Figure 6.12 Two sections from the constructed conductivity model. Panels (a) and (b)show the conductivity in the section-A and section-B, respectively. These sections arechosen to coincide with the apparent conductivity pseudo-sections displayed in Fig. 6.7.Some geological information is available at these locations.The model exhibits clear definitions of various conductive and resistive units. Ideally, thesedefinitions should be related to different lithologic units using the information obtained fromdrilling. Although five major lithologic units have been identified in this region, their geo-electricproperties are still uncertain. This, together with the limited availability of drilling information,makes such verification difficult at this stage.From the viewpoint of the algorithm performance, the model may also prove to have itsvalidity. The inversion starts from a best fitting half-space model with a misfit of 36%, andconstructed a sequence of minimum norm models which steadily reduce the misfit to 12%. Sucha fit is not overly loose for a field data set. Due to the nature of the algorithm, there is little chance—274.CHAPTER 6. ANALYSIS OF FIELD DATA 152—.7400.—.880120.‘ 240.—1.02360.480.—1.16600.—1.300.—1.44- 120.240.—1.5836O.480.—1.72600.—1.88293. 860. 1426. 1993.X (m)2560.—3109. -2487. -1865. —1244. -622. 0.Y (m)CHAPTER 6. ANALYSIS OF FIELD DATA 153that large scale spurious structures are introduced. Therefore, features appearing in the model aremost likely required by the observed data and may reflect true large scale conductivity variations.6.5 DiscussionA set of field E-SCAN® data is analyzed using the techniques developed in this thesis. Thedata set was acquired as a part of an exploration effort for epithermal deposits. The field procedureof data acquisition is first reviewed. I then apply a first order correction for the placement errorof the “infinite” electrode using an estimated value for the regional background conductivity.The maximum correction is about 20%, which indicates that the correction is necessary for fielddata sets. The common source gathers of the potential measurements are then formed from thecorrected data by invoking the reciprocity of the d.c. potential field. Pole-pole data maps areformed and smoothed using the generalized cross-validation (GCV) smoothing technique.The AIM-DS iterative inversion is then applied to 26 data maps with a total of 6825 potentialdata. The constructed conductiiity model is represented by 34,040 parameters and achieved aRMS relative misfit of 12%. The model exhibits various structures and geo-electrical units whichare not visible in the apparent conductivity images. However, due to the lack of the informationabout the true geology in the area, the conductivity model is difficult to verify, but it seems that afew known structures are reflected in the model.Of course, the interpretation of a data set has at least two steps. The first is to construct aconductivity model which satisfies the data and reflects the true conductivity structures beneaththe surface. The second is then to correlate the recovered conductivity structure with variouslithologic units and to extract information about unknown geological structures. The first stagefalls in the realm of the inversion. Although one always like to see the verification of the modelby the ground truth in the second stage, such a procedure is beyond the scope of an inversionalgorithm.For the purpose of finding an acceptable conductivity model, the AIM algorithm has workedwell. It constructs a model which fits the field data reasonably. The reduction of the misfit isCHAPTER 6. ANALYSIS OF FIELD DATA 154achieved within 4 iterations. Each iteration takes less than 300 minutes CPU time on a 4/330SUN workstation with 32 MB memory. This shows that the application of the algorithm to fielddata is practical and efficient.CHAPTER 7CONCLUSIONSThe goal of the research in this thesis is to develop techniques for interpreting 3d d.c. resistivity data to extract information about subsurface conductivity structures. Several imaging andinversion algorithms have been developed based upon the understanding of charge accumulationin d.c. resistivity experiments and the Born approximation of the surface electrical potential. Allwork on imaging and inversion has concentrated on aspects of model construction.Various aspects of charge accumulation in the d.c. resistivity experiment are studied inChapter 2. When electric current is introduced into a conductive medium, charges accumulatein the region where there is a non-zero component of electric field parallel to the conductivitygradient. It is these accumulated charges which give rise to the observed potentials. The behaviorof a d.c. field can be explained with better understanding by examining the distribution of thesecharges. It is demonstrated that the refraction of current across a boundary separating two mediawith differing conductivity is a direct effect of the charges accumulated on the boundary. Forsimple geologic structures, it is often possible to sketch out the approximate potential responsebased upon the understanding of how the accumulated charges will distribute. This is useful inmaking preliminary interpretations of observed data, in understanding the result generated by aninversion algorithm, or in making first order judgements about the validity of the result from anumerical forward simulation. The traditional image method of solving potential problems isshown to be finding the set of fictitious point charges which produces the same potential as doesthe true accumulated charge. For systems consisting of isolated bodies of constant conductivity ina uniform background, the integral equation for the accumulated charges is generalized to includean undulating surface of the earth. This allows the solution of problems with surface topography.Both analytic and numerical examples are presented in which the charge accumulation is quantifiedusing an integral equation.155CHAPTER 7. CONCLUSIONS 156In Chapter 3, different imaging techniques are developed for extracting first order informationdirectly from the data. The display of apparent conductivity data volume with a pseudo-depth as ameans of imaging is examined. It is demonstrated that filtering techniques can be used to enhancethe data image. However, such data images only provide an indication of the lateral variationof the conductivity and possible relative depth of different anomalies. It provides neither theabsolute depth nor the actual structural information. A second method is developed to image thedepth of burial of simple structures by equivalent source images. The secondary electric potentialis produced by subsurface charges. Potential field theory allows calculation of an equivalentsource layer at any depth. Therefore, a 3d image can be obtained by constructing such equivalentsources at a sequence of depths. Numerical results show that the depth of burial for simpleanomalies can be estimated from such images.In the special case of 2d structures, the Born approximation of the potential field is used todevelop two algorithms for imaging the conductivity structures. Under the approximation, thesecondary potential on the surface is linearly related to a “pseudo-charge” which is a functionof the conductivity gradient and the position of the current source. This function represents theintegrated effect of accumulated charges distributed along the strike direction, hence the namepseudo-charge. Since the charge only accumulates in the region where the conductivity changes,the recovery of this function can yield a structural image for the conductivity. This is the basisof the first algorithm. A sparse representation of the pseudo-charge for each source location isconstructed by inverting the corresponding secondary potentials. Linear programming is used tosolve the inverse problem. The composite of these pseudo-charges from multiple source locationsforms the desired structural image. The drawbacks of the method are the severe non-uniquenessin recovering the charge distribution from its potentials, and the instability in combining severalindividual inversions.In the second approach, the same Born equation is used to formulate the inversion so thatan image of the conductivity is recovered from the secondary potentials generated by all currentsources. The pseudo-charge density can be expressed explicitly as a function of the conductivity.CHAPTER 7. CONCLUSIONS 157In this way, the surface secondary potential becomes a linear functional of the conductivity.Unlike estimation of the charge location, the inversion now has a single function to recover.Consequently, all surface data can be inverted at once. The algorithm constructs a piecewiseconstant conductivity model which minimizes a combination of the 11 norm of the model andthe 11 norm of the model gradient in a section parametrized into rectangular cells. Differentproportions of the two model components results in models having different flatness. Linearprogramming is again used to solve the problem. Due to the approximation involved, the solutionis to be interpreted as an image representation of the true conductivity.An efficient 3d inversion algorithm is developed in Chapter 4 based upon the Born approximation. The algorithm can be used to construct general 3d conductivity models from pole-poledata maps. The model is assumed to be composed of perturbations over a given uniform half-space, i.e., o) =0(i), where Ob is the background conductivity and ji() is the dimensionlessfunction representing the perturbation. With this model assumption and the Born approximation,an integral equation is derived for the potential anomalies of the surface pole-pole array configuration. The equation expresses the pole-pole potential anomaly as a depth integral of the 2dconvolution between the conductivity perturbation, ln(), and a kernel function in the horizontalplane. This convolutional representation then enables the data equation to be decomposed into asequence of independent components in the wavenumber domain by the Fourier transform. Thatis, each component of the potential data from any pole-pole array is related to the same componentof model perturbation by a single depth integral with a wavenumber domain kernel. It follows thatthe 3d inverse problem of recovering a multiplicative model perturbation is solved by a sequenceof id inversions in the wavenumber domain. Each available pole-pole data map in the spatialdomain provides one (complex) datum at every wavenumber via a 2d Fourier transform. At eachwavenumber, the id inversion is solved using linear inverse theory and a smallest model as afunction of depth is constructed for the conductivity perturbation. By Parseval’s theorem, thisgives a smallest model for the perturbation in the spatial domain. Thus the conductivity model isessentially a smallest deviatoric model with a uniform base model ln(o-o).CHAPTER 7. CONCLUSIONS 158This formulation of the 3d inverse problem has several advantages. It is computationallyefficient. With the convolutional representation of the surface data, the application of a Fouriertransform decomposes the 3d problem into a set of independent id problems. Each id problemis solved with minimal computational effort. The solution of the id problems collectively yieldthe solution for the original 3d problem through an inverse Fourier transform. This approachis much more rapid than solving the 3d problem directly. For a problem consisting of 20,490(32x32x20) cells and 8,800 potential data, the CPU time needed is about 100 minutes on a SUN4/330 workstation with 32 MB memory. Most of the CPU time is used in computing the kernelfunctions.The algorithm is capable of solving large problems. Since the problem is solved at eachwavenumber through id inversions, the computing resources and time required for the inversesolution is greatly reduced. Therefore, it can cope with large numbers of model parameters andlarge numbers of data. This is an important step toward practical applications.The algorithm is stable and less likely to produce large scale spurious structures. From analternative view, the Fourier transform decomposes a single large inverse problem into manysmaller problems each of which is sensitive to a different spatial scale length. The horizontalvariability is limited by the smallest scale (highest wavenumber) of component included. Inaddition, a model of smallest norm is sought at each component. The kernel function used formodel construction decays more rapidly with depth for higher wavenumbers. Therefore, thesolution eliminates from the model those features which cannot be resolved by the data. At thesame time, the algorithm naturally reflects in the model the decrease of resolution with the depth.In Chapter 5, the 3d approximate inversion and AIM formalism is used to formulate aniterative algorithm which constructs a 3d conductivity model reproducing the observed data. The3d approximate inversion is used as an approximate inverse mapping under the AIM formalism.Together with a finite difference forward modelling algorithm, it updates successively a initialmodel to produce a conductivity model which predicts the observed data within the toleranceof the data errors. The model obtained by simple application of the approximate 3d inversionCHAPTER 7. CONCLUSIONS 159has been used as an initial model in all test cases, however, any appropriate model can servethe purpose. Both AIM-MS and AIM-DS algorithms are constructed and tested with severalsynthetic data sets.The performance of the AIM algorithm depends to a great extent upon that of the approximateinverse mapping. Thus the constructed iterative inversion algorithm for 3d d.c. resistivity inheritsthe advantages of the approximate 3d inversion. It is rapid, stable, and capable of solving largeproblems. The kernel functions in the wavenumber domain are computed at the first iteration andstored for use in subsequent iterations. The computational effort for the solution of id systems isa very small fraction of that needed to calculate kernels. Thus the CPU time needed for applyinginverse mapping is almost negligible in each iteration. In addition, the forward mapping needsto be performed only once at each iteration. Therefore, the total computation is minimal — oneforward modelling is required in order to compute the data misfit in any iterative inversion. In allthe synthetic test cases, the algorithm has steadily reduced the data misfit to the expected levelfor both accurate and noisy data sets. The reduction of the misfit is achieved within the first fewiterations. The original conductivity models are reproduced with great fidelity in the depth rangeto which the given data set has sufficient responses.The AIM inversion algorithm is then applied in Chapter 6 to a field E-SCAN® data set toconstruct a 3d conductivity model. The data set is acquired in an epithermal deposit area. Thepurpose of the E-SCAN® survey is to locate the fault structures and silicified zones with whichore deposits may be associated. The inversion constructs a 3d model which fits the observed data.The model defines various conductive and resistive units. Due to the lack of available geologicalinformation, the validity of the model is difficult to verify. However, limited geological sectionsseem to support several features observed in the constructed model.The iterative inversion algorithm developed in this thesis offers a rapid approach to constructing 3d conductivity models by inverting large d.c. resistivity data sets fromE-SCAN® experiments,or any surface experiment which can provide multiple set of pole-pole data maps. The implementation is stable and efficient. The resultant conductivity model may fit the observed dataCHAPTER 7. CONCLUSIONS 160adequately and suffice to provide the qualitative and quantitative information needed in the datainterpretation. However, it is observed that the algorithm is likely to diverge after the initialreduction of the data misfit. The algorithm does not explicitly minimize the norm of the modelsubject to a desired misfit. As the number of iterations becomes large, the structural complexityof the model will increase and this will lead to the increase in the data misfit. As a result, thealgorithm does not necessarily construct a model which fits the data to an expected accuracy.In such cases, the resultant model from this algorithm can be used as a starting model for otherlinearized inversions. Such a scheme would be able to produce a model which approximately fitsthe data with minimal effort and then to carry out detailed analysis with a more rigorous inversionif necessary.REFERENCESAidridge, D. F., Dosso, S. E., Endres, A. L., and Oldenburg, D. W., 1991. New methods forconstructing flattest and smoothest models, Inverse Problems, 7, 499-513.Alpin, L. M., 1947. Source of the field in the theory of electroprospecting by direct current,Applied Geophysics, 3, 57-100.Alfano, L., 1959. Introduction to the interpretation of resistivity measurements for complicatedstructural conditions, Geophys. Prosp. 7, 311-368.Alfano, L., 1960. The influence of surface formations on the apparent resistivity values in theelectrical prospecting, part I, Geophys. Prosp., 8, 576-606.Alfano, L, 1961. The influence of surface formations on the apparent resistivity values in theelectrical prospecting, part II, Geophys. Prosp., 9, 213-241.Backus, G. E. and Gilbert, J. F., 1967. Numerical applications of a formalism for geophysicalinverse problems, Geophys. J. R. astr Soc. , 13, 247-276.Bertero, M., De Mol, C., and Pike, E. R., 1988. Linear inverse problems with discrete data: II.Stability and regularisation, Inverse Problems, 4, 573-594.Boerner, D. E. and West, G. F., 1989. Fréchet derivatives and single scattering theory, Geophys.J. mt., 98, 385-390.Bracewell, R. N., 1978. The Fourier transform and its applications, McGraw-Hall, Inc.Brebbia, C. A. and Walker, S., 1980. Boundary element techniques in engineering, Butterworth& Co. Ltd., London.Brown, B. H., 1986. Applied potential tomography, Phys. Bull., 37, 109-112.Coen, S. and Yu, M. W., 1981. The inverse problem of the direct current conductivity profile ofa layered earth, Geophysics, 46, 1702-17 13.Craven, P. and Wahba, G., 1979. Smoothing noisy data with spline function: Estimating thecorrect degree of smoothing by the method of generalized cross-validation, NumerischeMathematik, 31, 377-403.Das, U. C. and Parasnis, D. S., 1987. Resistivity and induced polarization responses of arbitrarilyshaped 3-D bodies in a two-layered earth, Geophys. Prosp., 35, 98-109.161REFERENCES 162Dey, A. and Morrison, H. F., 1979. Resistivity modelling for arbitrarily shaped three-dimensionalstructures, Geophysics, 44, 753-780.Dieter, K., Paterson, N. R., and Grant, F. S., 1969. IP and resistivity type curves for three-dimensional bodies, Geophysics, 34, 615-632.Dines, K. A. and Lytle, R. J., 1981. Analysis of electrical conductivity imaging, Geophysics , 46,1025-1036.Dosso, S. E. and Oldenburg, D. W., 1989. Linear and non-linear appraisal using extremal modelsof bounded variation, Geophys. J. mt. , 99,483-495.Edwards, L. S., 1977. A modified pseudosection for resistivity and IP, Geophysics , 42, 1020-1036.Eskola, L., 1979. Calculation of galvanic effects by means of the method of sub-areas, Geophys.Prosp., 27, 616-627.Goldberg, S., Loewenthal, D. and Rotstein, Y., 1982. An improved algorithm for magnetotelluricand direct current data interpretation, J. Geophys., 50, 151-158.Grant, F. S. and West, G. F., 1965. Interpretation theory in applied geophysics, McGraw-Hill,New York.Habberjam, G. M., 1979. Apparent resistiviiy observations and the use ofsquare array techniques,Borntraege, Berlin.Harrington, R. F., 1968. Field computation by momentum methods, MacMillan Co., New York.Inman, J. R., 1975. resistivity inversion with ridge regression, Geophysics , 40, 798-817.Jeans, Sir James, 1966. The mathematical theory of electricity and magnetism, CambridgeUniversity Press, Cambridge.Jupp, D. L. B. and Vozoff, K., 1975. Stable iterative methods for the inversion of geophysicaldata, Geophys. J. R. astr Soc., 42, 957-976.Karous, M. and Hjelt, S. E., 1983. Linear filtering of VLF dip-angle measurements, Geophys.Prosp., 31, 782-794.Kaufman, A. A., 1985. Distribution of alternating electrical charges in a conducting medium,Geophys. Prosp., 33, 171-184.REFERENCES 163Kaufman, A. A. and Keller, G. V., 1985. Inductive mining prospecting, Elsevier, New York.Keller, G. V., and Frischknecht, F. C., 1966. Electrical methods in geophysical prospecting,Pergamon Press Inc., New York.Koefoed, 0., 1976. Progress in the direct interpretation of resistivity soundings: An algorithm,Geophys. Prosp., 24, 233-240.Kunetz, G., 1972. Processing and interpretation of magnetotelluric soundings, Geophysics, 37,1005-1021.Lee, T., 1972. A general technique for the direct interpretation of resistivity data over two-dimensional structures, Geophy. Prosp., 20, 847-859.Levy, S., Oldenburg, D., and Wang, J., 1988. Subsurface imaging using magnetotelluric data,Geophysics, 53, 104-117.Lindgren, W., 1933. Mineral deposits, McGraw-Hill, New York.McGillivray, P. R. and Oldenburg, D. W., 1990. Methods for calculating Fr&het derivatives andsensitivities for the non-linear inverse problem: A comparative study, Geophys. Prosp.,38,499-524.Morse, P. M. and Feshbach, H., 1963. Metlwds of theoretical physics, McGraw-Hill, New York.Murty, K. G., 1983. Linear programming, John Wiley & Sons Inc., New York.Okabe, M., 1984. Boundary element method for the arbitrary inhomogeneities problem inelectrical prospecting, Geophys. Prosp., 29, 39-59.Oldenburg, D. W., 1978. The interpretation of direct current resistivity measurements, Geophysics, 43, 610-625.Oldenburg, D. W., 1984. An introduction to linear inverse theory, IEEE Trans. Geosci. RemoteSensing, GE-22, 665-674.Oldenburg, D. W. and Ellis, R. G., 1991. Inversion of geophysical data using an approximateinverse mapping, Geophys. J. Int., 105, 325-353.Oppliger, G. L., 1984. Three-dimensional terrain corrections for mise-á-la-masse and magnetometric resistivity surveys, Geophysics, 49, 17 18-1729.REFERENCES 164Park, S. K. and Van, G. P., 1991. Inversion of pole-pole data for 3-D resistivity structures beneatharrays of electrodes, Geophysics, 56, 951-960.Parker, R. L., 1977. Understanding inverse theory, Ann. Rev. Earth Planet. Sci., 5, 35-64.Parker, R. L., 1984. The inverse problem of resistivity sounding, Geophysics, 49, 2143-2158.Pelton, W. H., Rijo, L., and Swift, C. M., Jr., 1978. Inversion of two-dimensional resistivity andinduced-polarization data, Geophysics, 43, 788-803.Petrick, W. R., Jr., Sill, W. R., and Ward, S. H., 1981. Three-dimensional resistivity inversionusing alpha centers, Geophysics, 46, 788-803.Pratt, D. A., 1972. The surface integral approach to the solution of the 3-D resistivity problems,ASEG Bulletin, 3, 33-50.Roy, A, 1978. A theorem for direct current regimes and some of its consequences, Geophys.Prosp., 26, 442-463.Roy, A. and Apparao, A., 1971. Depth of investigation in direct current methods, Geophysics,36, 329-340.Scollar, I., Tabbagh, A., Hesse, A., and Herzog, I., 1990. Archeological Prospecting and RemoteSensing, Cambridge University Press, Cambridge.Shulz, R., 1985. The method of integral equation in the direct current resistivity method and itsaccuracy, Journal of Geophysics, 56, 192-200.Smith, N. C. and Vozoff, K., 1984. Two-dimensional DC resistivity inversion for dipole-dipoledata: IEEE Trans. Geosci. Remote Sensing, GE-22, 2 1-28.Snyder, D. D., 1976. A method for modelling the resistivity and IP responses of two-dimensionalbodies, Geophysics, 41, 997-1015.Stefanescu, S. and Stefanescu, D., 1974. Mathematical models of conducting ore bodies for directcurrent electrical prospecting, Geophys. Prosp., 22, 246-260.Stevenson, A. F., 1934. On the theoretical determination of earth resistance from surface potentialmeasurements, Physics, 5, 114-124.Stratton, J. A., 1941. Electromagnetic theory, McGraw-Hill, New York.REFERENCES 165Szaraniec, E., 1976. Fundamental functions for horizontally stratified earth, Geophys. Prosp.,24, 528-548.Szaraniec, E., 1980. Direct resistivity interpretation by accumulation of layers, Geophys. Prosp.,28, 257-268.Telford, W M., Geldart, L. P., Sheriff, R. E., and Keys, D. A., 1976. Applied Geophysics,Cambridge University Press, Cambridge.Tripp, A. C., Hohmann, G. W., and Swift, C. M., Jr., 1984. Two-dimensional resistivity inversion,Geophysics, 49, 1708-1717.Van Ziji, J. S. V. and Joubert, S. J., 1975. A crustal geoelectrical model for south AfricanPrecambian granitic terrains based on deep Schiumberger soundings, Geophysics, 40,657-663.Vozoff, K., 1960. Numerical resistivity interpretation: General inhomogeneity, Geophysics, 25,1184-1194.Ward, S. H. and Fraser, D. C., 1971. in Mining geophysics - theory, Chapter 2, Part 2, Editor D.A. Hansen et al., Society of Exploration Geophysicists, Tulsa.Ward, S. H. (Ed.), 1990, Geotechnical and environmental geophysics, Society of ExplorationGeophysicists, Tulsa.Ward, S. H. and Hohmann, G. W., 1988. in Electromagnetic method in applied geophysics -theory, Chapter 4, Editor M. N. Nabighian, Society of Exploration Geophysicists, Tulsa.Wright, P. M., Ward, S. H., Ross, H. P., and West, R. C., 1985. State-of-the-art geophysicalexploration for geothermal resources, Geophysics, 50, 2666-2699.Yorkey, T. J. and Webster, J. G., 1987. A comparison of impedance tomography reconstructionalgorithms, Clin. Phys. Physiol. Meas. Suppl. A, 8, 55-62.Zohdy, A. A. R., 1974. Use of Dar Zarrouk curves in the interpretation of vertical electricalsounding data, U.S. Geol. Sun’. Bull., 1313-D, 1-41.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Inversion of three-dimensional direct current resistivity...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Inversion of three-dimensional direct current resistivity data Li, Yaoguo 1992
pdf
Page Metadata
Item Metadata
Title | Inversion of three-dimensional direct current resistivity data |
Creator |
Li, Yaoguo |
Date Issued | 1992 |
Description | A direct current (d.c.) resistivity experiment investigates subsurface geo-electrical structures by measuring the electric field set up by introducing current into the earth. Information about geo-electrical structures is extracted by inverting the observed data to generate an image of the conductivity or to construct a conductivity model. The goal of this thesis is to develop efficient inversion techniques for the interpretation of three-dimensional (3d) d.c. resistivity data. The study assumes data consisting of pole-pole potentials measured over a regular grid on the surface for many current locations. The Born approximation is employed to linearize the inverse problem. The source of the electric field measured in the d.c. resistivity is the accumulated electric charges. Different aspects of the charge accumulation are reviewed, enlarged with new insights and presented in a unified notation. This provides the basis for understanding the fundamentals of d.c. resistivity experiments. Two algorithms are developed to image simple 2d conductivities. The first constructs a structural image by combining the charge density images obtained by inverting multiple sets of common current potentials. The second constructs a conductivity image directly. Processing and displaying the apparent conductivity, and constructing equivalent sources from secondary potentials are studied as the means of imaging. Assuming a multiplicative perturbation to a uniform half-space, the potential anomaly of pole-pole arrays is expressed as a depth integral of the logarithmic perturbation convolved with a kernel function in the horizontal directions. Applying the Fourier transform decomposes the data equation for a 3d problem into a set of id equations. A rapid approximate 3d inversion is developed based upon this decomposition by solving a sequence of id inversions in the wavenumber domain. The approximate 3d inversion is used to construct iterative inversion algorithms using the AIM (Approximate Inverse Mapping) formalism. The approximate inversion and an exact forward mapping are used to update the model successively so that the final result reproduces the observed data. The AIM inversion is applied to analyse a set of field data. |
Extent | 7178648 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085637 |
URI | http://hdl.handle.net/2429/3142 |
Degree |
Doctor of Philosophy - PhD |
Program |
Astronomy |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1992-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1992_spring_li_yaoguo.pdf [ 6.85MB ]
- Metadata
- JSON: 831-1.0085637.json
- JSON-LD: 831-1.0085637-ld.json
- RDF/XML (Pretty): 831-1.0085637-rdf.xml
- RDF/JSON: 831-1.0085637-rdf.json
- Turtle: 831-1.0085637-turtle.txt
- N-Triples: 831-1.0085637-rdf-ntriples.txt
- Original Record: 831-1.0085637-source.json
- Full Text
- 831-1.0085637-fulltext.txt
- Citation
- 831-1.0085637.ris
Full Text
Cite
Citation Scheme:
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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085637/manifest