Inversion of Airborne Electromagnetic data in 2.5D by Wing Wa Yu B.Sc., University of Iceland, 2006 M.Sc., University of British Columbia, 2009 a thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in the faculty of graduate studies (Mathematics) The University Of British Columbia (Vancouver) October 2012 c© Wing Wa Yu, 2012 Abstract In this work, we implement an inversion algorithm for airborne electromag- netic (AEM) data in the frequency domain by using 2D conductivity models. First, we discretize the 2.5D Maxwell’s equations on a staggered grid and test the numerical accuracy of the forward solution. The inverse problem is then solved by regularized minimization approach using the limited memory BFGS variant of the quasi-Newton method. Next, EM responses from a synthetic 2D conductivity model are inverted to validate the algorithm. Finally, we use the algorithm on an AEM field dataset from a RESOLVE survey and compare the inversion results to those obtained from a well-established 1D implementation. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Inverse problems . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 AEM system . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Motivation and Goal . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . 8 2 Forward Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1 Maxwell’s equations in 2.5D . . . . . . . . . . . . . . . . . . . 9 2.2 Fourier transform of the source term . . . . . . . . . . . . . . 13 2.3 Optimization of ky and g . . . . . . . . . . . . . . . . . . . . . 14 2.4 2D discretization of Maxwell’s equations . . . . . . . . . . . . 17 2.4.1 Discretization of the physical quantities . . . . . . . . . 17 2.4.2 Discretization of the differential operators . . . . . . . 19 2.4.3 Discretization of the material averaging . . . . . . . . . 25 iii 2.4.4 Matrix representation of Maxwell’s equations . . . . . . 26 2.5 Solving the adjoint problem . . . . . . . . . . . . . . . . . . . 28 2.6 Solving for multiple RHS . . . . . . . . . . . . . . . . . . . . . 29 3 Inversion framework . . . . . . . . . . . . . . . . . . . . . . . 32 3.1 Objective function . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Sensitivity matrix . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.3 Method of optimization . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Implementation of the optimization . . . . . . . . . . . . . . . 38 3.4.1 Data structures . . . . . . . . . . . . . . . . . . . . . . 38 3.4.2 Constrained optimization . . . . . . . . . . . . . . . . . 38 3.4.3 Regularization parameter selection . . . . . . . . . . . 40 3.4.4 Stopping criteria . . . . . . . . . . . . . . . . . . . . . 44 4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.1 Testing the forward solution . . . . . . . . . . . . . . . . . . . 45 4.2 Inversion of synthetic data . . . . . . . . . . . . . . . . . . . . 48 4.3 Inversion of field data . . . . . . . . . . . . . . . . . . . . . . . 50 5 Conclusion and future work . . . . . . . . . . . . . . . . . . . 59 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 iv List of Tables Table 2.1 Overview of the discretization of the physical quantities and the operators used in the Maxwell’s equations . . . . . . . . 20 Table 2.2 Mapping of the four types of gradient operator . . . . . . . 20 Table 2.3 Substitution of operators for the setup of the adjoint problem 30 Table 4.1 Solution error of the 2.5D forward modeling for different cell sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 v List of Figures Figure 1.1 An AEM system and the magnetic response from a buried conductor. The transmitter (Tx) and receiver (Rx) of the system are housed in a bird which is towed behind and be- low a helicopter. In the presence of a conductive body, the primary field from the transmitter induces eddy currents in the conductor and a secondary field is produced as a result. The induced field measured at the receiver loop is then used to determine the conductivity of the body. . . . 4 Figure 1.2 A cross-sectional view of the layered earth model. The model is made of layers of earth, each with a thickness of hi and a conductivity of σi. . . . . . . . . . . . . . . . . . 6 Figure 1.3 A cross-sectional view of the 2D earth model. The model is a function of both position along the flight direction and depth. Each conductor has a conductivity of σi and its extent in the strike direction is infinite. . . . . . . . . . 7 Figure 2.1 Discretization of the electromagnetic field and the conduc- tivity model on a 2D staggered grid. Black dots denote the location of the nodes, boxes are the x-edges (defined as edges whose normal points in the x-direction), diamonds are the z-edges and crosses represent the cell-centers. . . . 21 vi Figure 2.2 An example of a mesh generated for AEM inversion. The color shows the initial guess of the log of the conductivity model and the crosses represent the source positions. . . . 31 Figure 3.1 Tikhonov curve. For large β, the model norm φm domi- nates the objective and the data misfit is underrepresented. The opposite happens when β is small. When the noise level is known, the target misfit, φd(β ∗,m∗), can be used to determine the regularization parameter. . . . . . . . . . 42 Figure 3.2 The L-curve. It can be used to select the regularization parameter β when the noise level is unknown. In this strategy, β is chosen at the point of maximum curvature of the L-curve. . . . . . . . . . . . . . . . . . . . . . . . . 43 Figure 4.1 Schematic of the mesh and the conductivity model used to test the 2D forward modeling. . . . . . . . . . . . . . . 46 Figure 4.2 Comparison of 2D response with 1D analytical solution. a) Computed responses from the synthetic model shown in Fig. 4.1. Circles are the 1D analytical responses and the crosses are the 2D modeled responses. b) Relative error of the 2D responses. . . . . . . . . . . . . . . . . . . 47 Figure 4.3 Synthetic model used in testing the 2D inversion algo- rithm. The model is 2D and it consists of a conductive prism and an overburden that is conductive on one end and resistive on the other. Responses were recorded 10 m away from the sources and 40 m above the ground. Loca- tions of the sources are indicated by the red crosses. . . . . 50 vii Figure 4.4 a) Log of the recovered 2D conductivity model for the syn- thetic data generated by the model shown in 4.3. b) Left: misfit of predicted data normalized by noise standard de- viation. Middle and right: The response curves are offset by 1000 ppm for clarity. Crosses represent the synthetic inphase and quadrature responses and solid curves are the predicted data. . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 4.5 Satellite image of the Red Dog mine area in Alaska. Red line: outline of the survey area. White line: outline of the waste area. Blue line: outline of the tailing pond. Black lines: ore deposits. Green lines: survey lines inverted. Image from [21]. . . . . . . . . . . . . . . . . . . . . . . . 53 Figure 4.6 Predicted data misfit for the inversions on the lines shown in Fig. 4.5: a)-b) Line 10300, c)-d) Line 20590. The curves are offset by 1000 ppm for clarity. Solid lines: predicted response. Crosses: observed data. Observed data are from [21]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Figure 4.7 Normalized bird altitude and a measured response curve. The two observables are strongly correlated. . . . . . . . 55 Figure 4.8 2D (a),c)) and 1D (b),d)) inverted log conductivity models for the two survey lines highlighted in Fig. 4.5. a)-b) Line 10300, c)-d) Line 20590. . . . . . . . . . . . . . . . . . . . 56 viii Acknowledgments From the beginning to the end of this Master’s program, I have received a lot of supports from many good people and I am very thankful to each one of them. My thesis advisor, Eldad Haber, has always been extremely helpful when- ever I encounter difficulties with my research. He is always very supportive of the directions his students would like to take and he always leaves them with plenty of freedom to create and develop their own ideas. For that I am very grateful. I would like to express my gratitude to Doug Oldenburg, Philip Loewen for being in my thesis committee and for providing invaluable helps in making a successful completion of my degree possible. Many thanks also go to Joel Jansen and Boris Lum of Teck Resources for sharing the Red Dog field dataset as well as giving me the chance to gain industrial experience through an internship with Teck. Last but not least, I am very thankful to all my fellow students and colleagues in UBC and in Teck for broadening my understanding in the field of applied mathematics and computer science, as well as in geophysics and geology. ix Chapter 1 Introduction 1.1 Inverse problems The need to reconstruct a property of a system by using observations made on the system arises in many fields of science and engineering: for exam- ple geophysics, medical imaging, machine learning, remote sensing, particle physics, astronomy and many more. In many cases, despite the relative ease to conduct experiments on a sys- tem we are interested in, acquiring enough information such that an property can perfectly be predicted is usually not possible. This problem is most ap- parent in physical systems, where usually only a fraction of the total number of degrees of freedom is easily accessible, for example in a 3D system where measurements can only be made on the boundaries. This lack of measurable parameters makes the system under-determined: the number of variables in the property that we wish to reconstruct is greater than the number of data points available. Parameter estimation for an under-determined system is essentially an ill-posed problem, meaning that its solution space is infinitely large and a unique solution does not exist. This type of problem is solved by reducing the solutions to only those that are both good predictors of the observations and 1 reasonable for the system according to some criteria. Information previously learned about the system, or a priori information, is used to devise the criteria to make educated guesses on the solutions. Studies of these problems, known as inverse problems, have given us a general framework to convert any type of observations on a system to infor- mation about any type of properties of the system. But before any data can be inverted, the forward solution has to be known. That is, given a prop- erty model, we need to be able to predict the outcome of observations on the model. Once the forward problem is solved, the inverse framework can subsequently be employed to reconstruct the model. The goal of an inverse algorithm is to rebuild the system property by choosing a model that both best-fits to the observed data and satisfies the constraints formulated by a priori information on the system. Solving inverse problems therefore involves both solving systems of equations as well as optimizing objective functions with constraints. 1.2 AEM system In recent geophysical exploration projects, airborne electromagnetic (AEM) surveys have become the method of choice for many mining companies to determine the conductivity model under the area of investigation. The pop- ularity of AEM systems comes from their ability to cover a large area in a relatively short time while requiring only a small number of human opera- tors. They are also ideal when the survey area has complicated terrain or is not accessible by ground transportation. In most AEM systems, horizontal coplanar loops (HCP) are used as trans- mitters (Tx) and receivers (Rx) to detect conductivity contrasts in the earth. In frequency-domain AEM systems, the loops are typically separated by 10 meters and they are housed in a bird that is towed behind and beneath an aircraft at a height of 30-70 meters above the earth’s surface. To determine the conductivity contrast, a time-varying current is passed through the trans- 2 mitter loop to produce a magnetic dipole field. If a conductive body is buried in the ground, the primary field will induce eddy currents in the conductor and a secondary magnetic field will be generated as a result. This secondary field changes the magnetic flux through the receiver loop and its magnitude can be quantified by measuring the emf induced in the receiver (see Fig 1.1). In this way, the AEM system probes the conductivity distribution next to the loops by measuring the earth’s response to changing magnetic fields. In typical frequency-domain AEM systems, each conductivity depth sound- ing (mapping of conductivity as a function of depth) consists of magnetic field measurements at five to six frequencies. The frequencies usually span over three orders of magnitude and are spread logarithmically from hundreds of Hz to hundreds of kHz. In this frequency range, the diffusion length of the magnetic field, which is given by the skin depth δ = √ 2 σωµ0 , (1.1) is on the order of 1 to 100 meters. Thus, current can be induced from a few meters up to a few hundreds of meters deep into the ground, and the depth of penetration drops as the inverse square root of the frequency. The relation between frequency and penetration depth allows AEM systems to profile the conductivity distribution as a function of depth. To map the horizontal direction, soundings are carried out at regularly spaced intervals along a flight line. Typically, data are sampled every few tenths of a second which is equivalent to a sounding every few meters along the line. To create a pseudo-3D map of the conductivity of the survey area, multiple lines are flown over the area. The line separation can range from tens to hundreds of meters, depending on the size of the geological features under study. 3 bird Tx Rx Conductor Secondary eld Primary eld Figure 1.1: An AEM system and the magnetic response from a buried conductor. The transmitter (Tx) and receiver (Rx) of the sys- tem are housed in a bird which is towed behind and below a helicopter. In the presence of a conductive body, the primary field from the transmitter induces eddy currents in the conduc- tor and a secondary field is produced as a result. The induced field measured at the receiver loop is then used to determine the conductivity of the body. 4 1.3 Motivation and Goal Inverse algorithms that are both reliable and fast enough to process an entire AEM survey dataset are still challenges to many researchers at the moment. The constant movement of the source over a large area of investigation means a vast number of sources has to be solved in each forward modeling. It is not unusual for an airborne survey to have more than hundreds of thousands of different source locations. Compared to ground survey methods, this is about two orders of magnitude more sources. The first generation of inverse algorithms was built on a 1D conductivity model that assumes the earth is made of layers of conductors stacked on top of each other (see Fig. 1.2). An analytical solution to the 1D problem can be obtained for HCP sources [18] and numerical implementations [4, 23] of the solution have been extensively studied and they are well-established. This simple layered earth model has been used in many surveys with good results and on today’s processors the 1D algorithm is extremely fast to run. However, the main shortcoming of the model is its simplicity. There are cases where the 1D implementation has troubles resolving correctly 2D or 3D geological structures that exhibit rapid lateral change in the immediate region under the probes (e.g., [5, 14, 25, 27]). To obtain a more accurate description of the earth, inversion techniques based on the full 3D conductivity model are now widely being studied and tested. At the moment, 3D inversion algorithms based on the finite volume (FV) method [9] and the integral equation approach [3] are available. Although the inversion of an entire AEM dataset is now possible with these algorithms, they still face numerous problems. A full inversion with the FV method can take hours to complete, and the amount of computational resources required is still huge and expensive. Although the integral equation technique is faster, it may not correctly resolve models with high conductivity contrasts. To try to resolve some of the problems that existing 1D and 3D inversions have, we tested and developed an inversion algorithm using a 2D conductivity 5 σ2 σ3 σΝ σ1 h2 . . . h1 h3 hΝ Figure 1.2: A cross-sectional view of the layered earth model. The model is made of layers of earth, each with a thickness of hi and a conductivity of σi. earth model(see Fig. 1.3). In AEM surveys, the spatial separation between survey lines is often large, and the survey is done perpendicular to geological strikes. As a result, each magnetic field measurement does not extend notice- ably into the adjacent lines and so each flight line can effectively be inverted independently without too much loss in resolution. This motivates the use of 2D earth models to invert AEM datasets. Compared to 1D algorithms, 2D inversions have the advantage of enabling more precise interpretation of conductive anomalies; and compared to 3D strategies, they require much less computing resources and are relatively straightforward to parallelize for multi-core CPUs. Over the past decades, the well-known 2.5D formulation of Maxwell’s equations has been used to design numerous 2D EM algorithms ([6, 20, 6 σ2 σ1 Strike direction σ3 Figure 1.3: A cross-sectional view of the 2D earth model. The model is a function of both position along the flight direction and depth. Each conductor has a conductivity of σi and its extent in the strike direction is infinite. 22]). These algorithms, which have been applied to various airborne and controlled-source EM problems (e.g., [1, 13, 15, 24]), are either based on fi- nite element (FEM) or finite difference (FD) approaches. But none used the finite volume method, which we have preferred over FEM or FD, due to its simplicity in generating the matrix for the linear system, as well as its ability to handle high conductivity contrast. In this work, we used the finite volume approach to implement an algo- rithm for inversion of frequency-domain AEM data. The algorithm, which is developed in Matlab, was first tested on synthetic data generated by a 2D model and later on a 3D AEM field dataset from a Fugro RESOLVE sur- vey. Our 2D inversions on the dataset showed strong agreement with models recovered by a well-established 1D algorithm. 7 1.4 Outline of the thesis The remainder of this thesis is divided into four chapters. In chapter two, we formulate and solve the forward problem by discretizing the 2.5D Maxwell’s equations in the frequency domain on a 2D orthogonal mesh. In chapter three, we present the inversion framework by discussing and constructing the objective function and the sensitivity matrix, followed by a look into the technique used in optimizing the objective function. Finally, in chapter four, we show the accuracy of the forward model and then present the results of our 2D inversions on both synthetic and field data. A conclusion and discussion of future work can be found in the last chapter. 8 Chapter 2 Forward Modeling In this chapter, we start by presenting the well-studied 2.5D formulation of the Maxwell’s equations. By assuming a 2D conductivity model that is in- variant in the strike direction, 3D electromagnetic fields are rewritten as the sum of 2D fields in the wavenumber domain by applying Fourier transfor- mations w.r.t. the strike direction. Once the 2.5D Maxwell’s equations are solved for an optimized set of wavenumbers, the inverse transformation is used to obtain the 3D solution. Subsequently, we will present a finite volume discretization of the elec- tromagnetic field on a 2D staggered grid, followed by a discretization for the CURL operator and the material averaging. Finally, the discrete system for the 2.5D Maxwell’s equations is set up to compute the electromagnetic field. 2.1 Maxwell’s equations in 2.5D Conventionally, measurement data in AEM surveys are presented as complex numbers: the real component of the magnetic response, which is called in- phase, and the imaginary part, called quadrature. The response is the ratio between the secondary magnetic field and the primary field produced by the transmitter in free space. For horizontal coplanar loops, the response is given 9 by Resp = ~Hz(rtr)− ~Hp,z(rtr) | ~Hp,z(rtr)| , ~Hp,z(r) = − m 4pir3 (2.1) where rtr is the transmitter to receiver separation, m is the magnetic moment of the transmitter, ~Hp,z is the z-component of the primary field and ~Hz(rtr) is the total field measured at the receiver. To model the magnetic response due to harmonically changing electro- magnetic fields, ~F (~x, t) = ~F (~x)e−iωt, we solve Maxwell’s equations in the frequency domain: ∇× ~E − iωµ0 ~H = ~Ms (2.2a) ∇× ~H − γ ~E = ~Js on Ω (2.2b) n̂× ~H = 0 on ∂Ω (2.2c) where we have imposed a vanishing magnetic field on the boundary of the spatial domain Ω, which consists of both the earth and the air above it. In these equations, ~E and ~H are the electric and magnetic fields, ω is the angular frequency of the fields, µ0 the magnetic permeability of free space and ~Js and ~Ms denote the external electric and magnetic sources. The conductivity γ = σ− iω is complex valued, but for most applications in geophysics, where the quasi-static condition ω σ usually applies [23], it can be simplified to the real-valued conductivity σ. From now on, we will use σ for the conductivity. To solve Maxwell’s equations, we decompose the total fields ~E, ~H into primary and secondary parts: ~E = ~E0 + ~E1, ~H = ~H0 + ~H1, (2.3) so (∇× ~E0 − iωµ0 ~H0) + (∇× ~E1 − iωµ0 ~H1) = 0 (2.4a) (∇× ~H0 − σ ~E0) + (∇× ~H1 − σ ~E1) = ~Js on Ω (2.4b) 10 where we have set ~Ms = 0, since only electric sources are used in our system. This decomposition separates Maxwell’s equations into a primary problem, ∇× ~E0 − iωµ0 ~H0 = 0 (2.5a) ∇× ~H0 = ~Js (2.5b) and a secondary problem, ∇× ~E1 − iωµ0 ~H1 = 0 (2.6a) ∇× ~H1 − σ ~E1 = σ ~E0 (2.6b) Setting ~Js to the current in a HCP, the primary problem describes the field due to the current loop in free space. With the source now replaced by the induced current density evaluated on the entire conductivity model σ, the secondary problem is the Maxwell’s equations for the induced magnetic field. This secondary form enables the magnetic response to be computed accurately which is otherwise very hard as the discretization error in com- puting the total field dominates the induced field strength. From now on, we will work with the secondary form and drop the secondary notation: ~E1, ~H1 → ~E, ~H. In this work, we assume an isotropic 2D conductivity model that has infinite extent in the strike direction y: σ = σ(x, z). Despite the 2D con- ductivity, each field component remains 3-dimensional as the source field is a three dimensional function. The advantage of using a 2D model is that now the 3D Maxwell’s equations can be reduced to 2D equations by Fourier transforming Eqn. (2.6) with respect to y. Since our system is invariant un- der translation in the y-direction, all field components are either symmetric or antisymmetric in y. Therefore, we can either use the sine or cosine integral 11 for the transformation: F̃ (x, ky, z) = −i ∫∞ 0 ~F (x, y, z) sin(kyy)dy if ~F is anti-symmetric in y∫∞ 0 ~F (x, y, z) cos(kyy)dy if ~F is symmetric in y (2.7) where F̃ denotes ~F in the Fourier space which is parameterized by the wavenumber ky. For antisymmetric functions, the sine transform is used, while the cosine transform is applied to symmetric functions. In the wavenumber domain, we have the following 2D Maxwell’s equations for each ky: ∇̃ky × Ẽ − iωµ0H̃ = 0 (2.8a) ∇̃ky × H̃ − σẼ = σẼ0 (2.8b) where ∇̃ky× is the 2.5D curl operator in the wavenumber domain and it acts on the field F̃ as follows: ∇̃ky × F̃ = −ikyF̃z − ∂zF̃y∂zF̃x − ∂xF̃z ∂xF̃y + ikyF̃x (2.9) where F̃ = (F̃x, F̃y, F̃z). The Fourier transform of the source field, Ẽ0, can be approximated by assuming a vertical dipole source and the formulation is shown in the next section. Once a solution to (2.8) is obtained for each ky, the 3D solution can be retrieved with the inverse Fourier transform: ~F (x, y, z) = 2ipi ∫∞ 0 F̃ (x, ky, z) sin(kyy)dky if F̃ is anti-symmetric in ky 2 pi ∫∞ 0 F̃ (x, ky, z) cos(kyy)dky if F̃ is symmetric in ky (2.10) To evaluate these integrals, we approximate them using the discrete Fourier 12 transforms: ~F (x, y, z) = 2ipi ∑N j=1 F̃ (x, ky,j, z) sin(ky,jy)gj if F̃ is anti-symmetric in ky 2 pi ∑N j=1 F̃ (x, ky,j, z) cos(ky,jy)gj if F̃ is symmetric in ky (2.11) where gj is a weight associated with the j-th wavenumber ky,j. Later in the chapter, we will show how the choice of the wavenumbers and weights can be optimized for our problem in order to minimize the computational cost. Now, we have seen how Maxwell’s equations in 2.5D are formulated. Not only is the computing time reduced because of the smaller size of the problem, the 2.5D implementation also makes the forward calculation easily paralleliz- able as each 2D equation in (2.8) can be solved independently. Once the 2D conductivity cross-section directly underneath each survey line is recovered, interpolation between adjacent lines can be used to create a pseudo-3D model of the survey area. 2.2 Fourier transform of the source term The primary fields ~E0, ~H0 are the free-space electromagnetic fields due to the current density in the source loop, ~Js: ∇× ~E0 − iωµ0 ~H0 = 0 (2.12a) ∇× ~H0 = ~Js (2.12b) Instead of solving (2.12) directly, the primary field ~E0 can be obtained by modeling the transmitter loop as a vertical magnetic dipole source. For a source centered at (xs, 0, hs) and oscillating at frequency ω, the electric field 13 is given by the vector potential of the dipole, ~A [11]: ~E0(x, y, z) = −∂t ~A = iω ~A = iωmµ0 4pir3 −yx− xs 0 (2.13) where r = √ (x− xs)2 + y2 + (z − hs)2. For every source on the same survey line, the Fourier transform of ~E0 is given by Ẽ0(x, ky, z) = iωµ0m 4pi ikyK0(|ky|s)|ky|K1(|ky|s)/s 0 (2.14) where s = √ (x− xs)2 + (z − hs)2 and Kn is the n-th order modified Bessel function of the second kind. It is worth mentioning that with this formula- tion, vertical loops can easily be modeled as well by simply switching the x and z coordinates. 2.3 Optimization of ky and g Given a fixed number of wavenumbers to be used to solve (2.8), we seek the optimum ky,j’s and the associated weights such that the error of approximat- ing the Fourier transform (Eqn. (2.11)) is minimized. In this work, this is done by finding the ky,j’s and gj’s that minimize the error in approximating a vertical dipole field on the surface of a homogeneous half-space. The same optimized values are then used to compute responses due to inhomogeneous half-spaces. Since the receivers in most AEM systems are located in-line with the sources (we will use yRx = yTx = 0), we only need to be concerned with minimizing the error for the locations on the survey line. For this, we will follow the optimization technique proposed in [17, 26]. We start by Fourier transforming the field due to a homogeneous half- space. For the range of conductivities and frequencies typically found in 14 AEM surveys, the low-frequency approximation applies, i.e., |kr| 1, where k = √−iσωµ0 and r is the distance from the center of the source. The quasi-static field of a vertical dipole source on the surface of a half-space with conductivity σ in this regime is given by [19]: ~H lx(x, y, 0) = ~H l y(x, y, 0) = 0, ~H l z(x, y, 0) = −m 4pir3 ( 1− k 2r2 4 ) (2.15) where r = √ x2 + y2 for a source centered at the origin. Substituting (2.15) into the Fourier transform for symmetric functions (2.7) gives H̃ lz(x, ky, 0) = −m 4pi ( |ky| x K1(|ky|x)− k 2 4 K0(|ky|x) ) (2.16) Our objective is to minimize the relative error in approximating (2.15) with the discrete Fourier cosine transform (2.11): ~H lz(x, 0, 0) ≈ 2 pi ∑ j H̃ lz(x, ky,j, 0)gj (2.17) for a range of x. Given a location xi, the relative error i is: i = 2 pi ~H lz(xi, 0, 0) ∑ j H̃ lz(xi, ky,j, 0)gj − 1 (2.18) and so our objective function can be written in the matrix form: Φ = 1 2 ‖‖22 = 1 2 ‖1−Mg‖22 (2.19) where Mij = 2H̃ lz(xi, ky,j, 0) pi ~H lz(xi, 0, 0) , (2.20) and g is a vector whose entries are gj and 1 is a vector of ones. Notice that the objective function (2.19) is a function of both ky and g. 15 To reduce it to a function that depends only on k, we optimize for g using the normal equation of least squares, g = (MTM)−1MT1 (2.21) and substitute it back into (2.19) to obtain Φ = 1 2 ‖1−M(MTM)−1MT1‖22 = 1 2 ‖1− v(ky)‖22 (2.22) Since the function v(ky) is nonlinear in ky, we linearize it by expanding it in Taylor series around a point k (0) y v(ky) = v(k (0) y ) + ∂v ∂ky δky, δky = ky − k(0)y . (2.23) Now the least squares solution for the objective function Φ = 1 2 ‖1− v(k(0)y )− ∂v ∂ky δky‖22 (2.24) gives the search direction δky: δky = (J TJ)−1JT (1− v(k(0)y )) (2.25) where J = ∂v ∂ky . At the i-th iteration, a new k (i+1) y is thus given by k(i+1)y = k (i) y + αδk (i) y (2.26) where α is a line search parameter. This iterative procedure reduces v(ky) at each iteration and is repeated until a certain error tolerance is reached. At this point, the corresponding ky is used to compute g using (2.20) and (2.21). It should be noted that since we do not have an analytical expression for 16 ∂v ∂ky , it has to be approximated by numerical differentiation: J(i) = ( ∂v ∂ky )(i) = [ ∆v ∆k (i) y,1 . . . ∆v ∆k (i) y,N ] , v(k(i)y + ∆k (i) y )− v(k(i)y ) (2.27) where we have set ∆k (i) y = 0.05k (i) y . 2.4 2D discretization of Maxwell’s equations In this section, we discretize the 2.5D Maxwell’s equations on a 2D staggered grid using a finite volume approach. The linear system to be solved is set up by discretizing the field quantities on either the edges, nodes or cell-centers of the staggered grid and then forming the discrete version of the curl operator and the material averaging for these grid variables. The discretization of the gradient, which appears in the regularization of the inverse problem, will also be presented here. 2.4.1 Discretization of the physical quantities To discretize the physical quantities in Maxwell’s equations, we consider an m×n staggered grid. Each cell is labeled by the indices i and j and has a size of hx,i by hz,j. To start, we put the electric fields Ẽx and Ẽz on the edges of each cell, as shown in Fig. (2.1). This configuration places constraints on the possible locations for the rest of the field quantities. To see this, we expand (2.8) to obtain the system of differential equations for the field components: 17 −ikyẼz − ∂zẼy − iωµ0H̃x = 0 (2.28a) ∂zẼx − ∂xẼz − iωµ0H̃y = 0 (2.28b) ∂xẼy + ikyẼx − iωµ0H̃z = 0 (2.28c) −ikyH̃z − ∂zH̃y − σẼx = σẼ0,x (2.28d) ∂zH̃x − ∂xH̃z − σẼy = σẼ0,y (2.28e) ∂xH̃y + ikyH̃x − σẼz = σẼ0,z (2.28f) Firstly, Eqns. (2.28a) and (2.28c) imply that the fields H̃x and H̃z should be discretized at the same location as Ẽz and Ẽx, respectively. Secondly, Eqn. (2.28c) places the x-directional derivative of Ẽy on the z-edges (edges whose normal points in z-direction). And since the two-point formula is used for the derivative, this configuration puts Ẽy midway between the z-edges in the x-direction, namely on the nodes. Following a similar reasoning, but now using the constraint set by Eqn. (2.28f), we can see that H̃y should be discretized at the cell-centers. Unlike the secondary fields, which can only be approximated numerically, the source field Ẽ0 is continuous and its value is known everywhere. However, although Ẽ0 is an electric field, we will not discretize it like Ẽ. The reason is that Ẽ0 is singular at the center of a source, (x, z) = (xs, hs). And since we would like to place the center on an z-edge, so that H̃z can be computed at the receiver center without averaging, we should avoid discretizing Ẽ0 on z-edges. As a result, we will evaluate Ẽ0 at the cell-centers instead, and then use averaging to bring the current density, σẼ0, to the same mesh as the electric field. In our modeling, the conductivity is assumed to be piecewise constant and therefore it should be discretized at the cell-centers. This discretization assumes that conductivity is constant throughout a cell but jumps may occur at the boundaries. This way, we can allow for highly varying and highly 18 discontinuous conductivity which is usually the case in geophysical settings. As required by the Maxwell’s equations, the conductivity has to be evaluated at the other grid locations as well. For this, the material averaging matrix, which is discussed later in this chapter, is used to transfer the cell-centered conductivity to the desired locations. An overview of the discretization of the physical quantities and the dis- crete operators that act on them is given in Table 2.1. 2.4.2 Discretization of the differential operators To construct the discrete version of the 2.5D curl operator (2.9), CURL, and the gradient operator, GRAD, we discretize the difference operator ∂x (discretization of ∂z is then found by symmetry). We begin by realizing that for our 2D staggered grid, there are four types of difference operators. Each of them can be characterized by one of the mappings shown in Table (2.2). Next, we notice that the gradient operators that map edge→cell- center or node→edge use the same 1D nodal difference which maps nodes to cell-centers. And similarly, the edge→node and cell-center→edge mappings follow the structure of the 1D cell-center differences. Moreover, the boundary conditions (BC) given in (2.2c) require a Dirichlet type BC for the magnetic field quantities, while the electric fields follow BC of Neumann type. From this information, we can see that the 1D cell-center differences with Dirichlet BC and the 1D node differences with Neumann BC are needed to build CURL and GRAD. In this section, we will demonstrate how to discretize the 1D differences and then how to use them to construct the discrete version of the 2D operators ∂x and ∂z. To build the matrices for the discrete operators, we will follow the column vector convention. Let U be a M×N matrix whose entries, uij, are the values of the function u(x, z) associated with the cell at (xi, zj), 1 ≤ i ≤ M, 1 ≤ j ≤ N . Instead of storing U as a 2D array, we vectorize it and store it in the 19 Table 2.1: Overview of the discretization of the physical quantities and the operators used in the Maxwell’s equations Quantity Notation Discretization electric field ẼxẼy Ẽz z-edgenode x-edge magnetic field H̃xH̃y H̃z x-edgecell-center z-edge source electric field Ẽ0,xẼ0,y Ẽ0,z cell-centercell-center cell-center conductivity σ cell-center Operator Notation Mapping curl of electric field CURLẼ z-edgenode x-edge → x-edgecell-center z-edge curl of magnetic field CURLH̃ x-edgecell-center z-edge → z-edgenode x-edge x-edge average Axecc cell-center → x-edge nodal average Andcc cell-center → node z-edge average Azecc cell-center → z-edge Table 2.2: Mapping of the four types of gradient operator Acts on Mapping Ẽx, Ẽz edge → cell-center Ẽy node → edge H̃x, H̃z edge → node H̃y cell-center → edge 20 xz EyEx, Hz Hy, σ Ez, Hx ~ ~ ~ ~~ ~ (xi, zj) (xi, zj+1) (xi+1, zj+1) (xi+1, zj) (xi+1/2, zj+1/2) (xi+1/2, zj) (xi, zj+1/2) cell ij Figure 2.1: Discretization of the electromagnetic field and the conduc- tivity model on a 2D staggered grid. Black dots denote the loca- tion of the nodes, boxes are the x-edges (defined as edges whose normal points in the x-direction), diamonds are the z-edges and crosses represent the cell-centers. 21 column format such that: U(k) = uij, k = (N − 1)× i+ j (2.29) where U(k) is the k-th component of the column form of U . In this way, the field quantities are stored as: F̃ = F̃xF̃y F̃z (2.30) where each 2D field component is in the vector column form. With this convention, the 2D operator that maps edges to cell-centers and nodes to edges can easily be constructed from the 1D nodal difference operator: D1Dx,node→cell-center = −1 hx,1 1 hx,1 −1 hx,2 1 hx,2 . . . . . . −1 hx,m 1 hx,m , (2.31) and the Kronecker product, A⊗B := a11B a12B . . . a1,nB a21B a22B . . . ... ... . . . . . . an−1,nB an,1B . . . an,n−1B annB (2.32) in the following way [10]: Dx,x-edge→cell-center = In ⊗D1Dx,node→cell-center (2.33a) Dx,node→z-edge = In+1 ⊗D1Dx,node→cell-center (2.33b) 22 where In is the n × n identity matrix and we have used the second-order two-point formula d dx u(xi+ 1 2 , zj) = u(xi+1, zj)− u(xi, zj) hx,i +O(h2x,i). (2.34) for the 1D derivative. It is easy to verify that the 2D difference matrices for the z-direction are: Dz,z-edge→cell-center = D1Dz,node→cell-center ⊗ Im (2.35a) Dz,node→x-edge = D1Dz,node→cell-center ⊗ Im+1 (2.35b) where D1Dz,node→cell-center is constructed by replacing the hx’s in (2.31) with hz’s. In a similar way, we use the 1D cell-center derivative, D1Dx,cell-center→node = 1 hx,1 −2 (hx,1+hx,2) 2 (hx,1+hx,2) . . . . . . −2 (hx,m−1+hx,m) 2 (hx,m−1+hx,m) −1 hx,m , (2.36) to build the 2D difference matrices for the mappings edge→node and cell- center→edge: Dx,z-edge→node = In+1 ⊗D1Dx,cell-center→node (2.37a) Dx,cell-center→x-edge = In ⊗D1Dx,cell-center→node (2.37b) Dz,x-edge→node = D1Dz,cell-center→node ⊗ Im+1 (2.37c) Dz,cell-center→z-edge = D1Dz,cell-center→node ⊗ Im (2.37d) Note that the use of the matrices (2.31) and (2.36) ensures that the Neumann 23 condition ∂n̂u = 0, and the Dirichlet condition, u(x, z) = 0, are imposed automatically at the boundaries and so the boundary condition (2.2c) is already taken care of. Now, with all four types of derivatives in our arsenal, it is straightforward to construct the curl operator that acts on the electric field Ẽ: CURLẼ(ky) = 0 −D xe nd −ikyI Dccze 0 −Dccxe ikyI D ze nd 0 (2.38) and the curl for the magnetic field H̃: CURLH̃(ky) = 0 −D ze cc −ikyI Dndxe 0 −Dndze ikyI D xe cc 0 . (2.39) For clarity, we have simplified DA→B to DBA and abbreviated x-edge, z-edge, node and cell-center with xe, ze, nd and cc, respectively. It is worth mentioning that if all the cells are of the same size, then D1Dx,nd→cc = −(D1Dx,cc→nd)T , and the two curl operators become adjoints of each other: CURLH̃ = CURL T Ẽ . (2.40) And as a consequence, the resulting matrix of the system is also symmetric. As we will see later, the regularization of the inverse problem requires the gradient of the conductivity to be computed. To achieve this, we put the cell-center→edge derivatives (2.37b) and (2.37d) together to obtain GRADcell-center→edge = [ Dxecc Dzecc ] . (2.41) 24 2.4.3 Discretization of the material averaging In this section, we construct the material averaging matrix, which averages the cell-centered conductivity and maps it to the edges and nodes where it is evaluated. For this, we use the nearest neighbor extrapolation to obtain the first-order 1D cell-center averaging: σi,j+ 1 2 = σi− 1 2 ,j+ 1 2 hi + σi+ 1 2 ,j+ 1 2 hi−1 hi−1 + hi +O(h). (2.42) However, when the conductivity has jumps at the cell boundaries where its derivative is not defined, this approximation is only O(1). To build the averaging matrix, we note that the 1D averaging has the same sparsity pattern as the 1D cell-center difference. And so to obtain the 2D average, we use the sparsity pattern of Eqn. (2.36) to get A1Dx,cell-center→node = 1 2 hx,2 (hx,1+hx,2) hx,1 (hx,1+hx,2) . . . . . . hx,m (hx,m−1+hx,m) hx,m−1 (hx,m−1+hx,m) 1 2 , (2.43) and then use the Kronecker product like before to construct the 2D averaging operator: Acell-center→x-edge = In ⊗A1Dx,cell-center→node (2.44a) Acell-center→z-edge = A1Dz,cell-center→node ⊗ Im (2.44b) Acell-center→node = A1Dz,cell-center→node ⊗A1Dx,cell-center→node (2.44c) 25 2.4.4 Matrix representation of Maxwell’s equations For each combination in the set of parameters {xs, ω, ky}, the following dis- crete 2.5D Maxwell’s equations have to be solved: CURLẼ Ẽ− iωµ0H̃ = 0 (2.45a) CURLH̃ H̃− diag(A1σ)Ẽ = A2(σ)Ẽ0 (2.45b) where A1 := A ze cc Andcc Axecc (2.46) and A2(σ) = A ze ccdiag(σ) 0 0 0 Andcc diag(σ) 0 0 0 Axeccdiag(σ) . (2.47) To solve the system, we can isolate for Ẽ in (2.45) and rearrange the equations to obtain: [CURLH̃CURLẼ + cdiag(A1σ)] Ẽ = −cA2(σ)Ẽ0 (2.48a) H̃ = CURLẼẼ/−c (2.48b) where c = −iωµ0. Once the matrix A = CURLH̃CURLẼ + cdiag(A1σ) is set up, the system can conveniently be solved with the Matlab backslash: Ẽ = A\(−cA2(σ)Ẽ0). Using (2.48) to solve for all field components at once can be computation- ally expensive when the mesh is dense. Another, simpler, way to a solution is to express (2.45) in terms of the y-component of the electric and magnetic fields. This approach reduces the size of the linear system to be solved. Once a solution is found, the other field components can be computed straightfor- wardly using numerical differentiation. To demonstrate, we start by writing 26 out explicitly the matrices acting on each field component: −ikyẼz −DxendẼy − iωµ0H̃x = 0 (2.49a) DcczeẼx −DccxeẼz − iωµ0H̃y = 0 (2.49b) DzendẼy + ikyẼx − iωµ0H̃z = 0 (2.49c) −ikyH̃z −DzeccH̃y − diag(Azeccσ)Ẽx = Azeccdiag(σ)Ẽ0,x (2.49d) Dndxe H̃x −Dndze H̃z − diag(Andcc σ)Ẽy = Andcc diag(σ)Ẽ0,y (2.49e) DxeccH̃y + ikyH̃x − diag(Axeccσ)Ẽz = Axeccdiag(σ)Ẽ0,z (2.49f) Now, after some algebra, we can eliminate all field components except Ẽy and H̃y from the equations and express the linear system for these remaining unknowns as[ Ẽy H̃y ] = [ S11 S12 S21 S22 ]−1 [ ikyD nd ze Γ −1 ze I ikyD nd xe Γ −1 xe −DcczeΓ−1ze 0 DccxeΓ−1xe ] A2(σ)Ẽ0 (2.50) Here, the blocks of S are defined as S11 = D nd xe ΣxeD xe nd + D nd ze ΣzeD ze nd − diag ( Andcc σ ) (2.51a) S12 = ikyD nd xe Γ −1 xe D xe cc − ikyDndze Γ−1ze Dzecc (2.51b) S21 = −iky c DccxeΓ −1 xe D xe nd + iky c DcczeΓ −1 ze D ze nd (2.51c) S22 = D cc xeΓ −1 xe D xe cc + D cc zeΓ −1 ze D ze cc − Im×n (2.51d) where Γxe = diag (A xe ccγ), Γze = diag (A ze ccγ), Σxe = diag (A xe ccσ) Γ −1 xe , Σze = diag (Azeccσ) Γ −1 ze and γ = k 2 y1 − cσ where 1 is a vector of ones. Once the 27 Eqn. (2.50) is solved, the remaining field components can be recovered from[ Ẽx Ẽz ] = [ Γ−1ze 0 0 Γ−1xe ][ ikyD ze nd cD ze cc ikyD xe nd −cDxecc ][ Ẽy H̃y ] + c [ Γ−1ze 0 0 0 0 Γ−1xe ] A2(σ)Ẽ0, (2.52)[ H̃x H̃z ] = [ ΣxeD xe nd ikyΓ −1 xe D xe cc −ΣzeDzend ikyΓ−1ze Dzecc ][ Ẽy H̃y ] + iky [ 0 0 −Γ−1xe Γ−1ze 0 0 ] A2(σ)Ẽ0. (2.53) Compared to solving for all components simultaneously, this method only requires solving for two field components. Cutting the number of variables by a third significantly improves the computation time as the number of matrix entries is now shrunk to half of that in the full system. Additionally, performing the elimination after the system is discretized preserves the order of accuracy of the original approximation, as rearranging the matrices and variables with algebra does not change the order of accuracy. On the contrary, if the system is first reduced and then discretized, the order of accuracy of the remaining field components is not guaranteed to be the same as that of Ẽy and H̃y. In this work, we have used Matlab’s backslash direct solver for the solution of the matrix system. For many right hand sides like we have here, it is more effective to use direct solving methods. However, iterative methods like conjugate gradient (CG) and multigrid become necessary when dealing with larger systems, for example when the fields are 3D functions. 2.5 Solving the adjoint problem As we will see later, solving the inverse problem requires solving the adjoint problem, which has the form: ATx = rhs. (2.54) 28 This problem is relatively easy if matrix A is small and Matlab’s backslash is used to obtain the numerical solution, as the adjoint matrix AT can be formed easily. However, if the reduced system is used for solution, the setup of the problem becomes a bit more involved. To set up the adjoint problem for the reduced system, we start by writing out AT explicitly, AT = CURLT Ẽ CURLT H̃ − cdiag(A1σ) = 0 −D xe nd −iky Dccze 0 −Dccxe iky D ze nd 0 T 0 −D ze cc −iky Dndxe 0 −Dndze iky D xe cc 0 T − cdiag(A1σ) = 0 (D cc ze) T −iky −(Dxend)T 0 (Dzend)T iky −(Dccxe)T 0 0 (D nd xe ) T −iky −(Dzecc)T 0 (Dxecc)T iky −(Dndze )T 0 − cdiag(A1σ) (2.55) And now by noticing that the blocks of CURLT Ẽ and CURLH̃ have the same sparsity pattern, as well as those of CURLT H̃ and CURLẼ, we can now solve the adjoint system by making the substitutions indicated in Table (2.3) to the operators in Eqns. (2.50)-(2.53). If the grid is uniform in size, then DBA = −(DAB)T for any type of mapping, and only the last substitution is needed. 2.6 Solving for multiple RHS Compared to other survey methods, the number of data points acquired for each source is much smaller in AEM systems. Each source in a frequency domain AEM survey is identified by its location and frequency. And for each source, only one EM sounding measurement is made. If an iterative method is used, this one to one ratio between the number of sources and 29 Table 2.3: Substitution of operators for the setup of the adjoint prob- lem A → AT Dzecc → −(Dccze)T Dndxe → −(Dxend)T Dndze → −(Dzend)T Dxecc → −(Dccxe)T Dxend → −(Dndxe )T Dccze → −(Dzecc)T Dccxe → −(Dxecc)T Dzend → −(Dndze )T c → −c the number of receivers can make the computation quite costly, as the total number of forward solves (x = A−1rhs) or transposed solves (x = A−T rhs) scales with the number of data points. However, using direct solver, mul- tiple sources can be solved at about the same time as solving for a single source. In this case, the majority of the computational effort goes into LU factorizing the forward matrix and the cost of solving for each source with forward/backward-substitution is only a fraction of the cost of the factoriza- tion. However, to take advantage of this, the sources have to share the same linear system. Namely, the mesh needs to be set up such that it contains all source locations. With such a mesh, the number of forward solves required does not scale linearly with the number of sources, but only with the number of sounding frequencies Nω and the number of optimum wavenumbers Nky . In Matlab, multiple right hand sides can be solved by concatenating the different RHS column vectors into a RHS matrix. The matrix is passed to backslash and the outputs are stored in a matrix whose columns are the solutions of the corresponding columns in the RHS matrix. A typical mesh generated for multiple source locations can be seen in Fig- 30 0 500 1000 1500 2000 2500 −200 −150 −100 −50 0 50 100 x (m) z (m ) −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 0 1 Figure 2.2: An example of a mesh generated for AEM inversion. The color shows the initial guess of the log of the conductivity model and the crosses represent the source positions. ure (2.2). The coarser padding cells, which are used to extend the core mesh boundary so to reduce boundary effects, are not shown in the figure. In the mesh, we can see that in order to incorporate the source altitudes accurately, cells around the sources have a hz less than 1 meter. This high tolerance is necessary, as a difference of a meter in the altitude of the source can affect the computed response by a couple of percent, creating discretization errors that could dominate other sources of error in the data. If the discretization error is not carefully controlled, serious error underestimation can occur in the inversion process. This can lead to unsuccessful inversions. Now that we know how to predict the magnetic response from any given conductivity distribution, we can proceed to apply the standard inversion framework to build an algorithm to invert AEM data to conductivity models. 31 Chapter 3 Inversion framework The objective of an inversion algorithm is to convert observed data into a model of a system property such that the data can be explained by the model. But since most inverse problems are ill-posed, meaning that they are not well- posed in the sense defined by Jacques Hadamard, it is necessary to restrict the possible recovered solutions to the most probable and most reasonable ones. What counts as a good recovered model depends on the application. For example, in geophysical settings where physical properties of the earth are recovered, a good model should show few features on a length scale longer than the expected size of the anomalies, while the recovered anomalies should have a size typical for structures created by similar geological formations. By contrast, in image deblurring, a good reconstruction adds sharper contrasts to the features of the image and makes it richer in details. In inverse prob- lems, the different application dependent requirements on the solution can be achieved by adding the appropriate regularizations to the data-fitting via mathematical constraints. 32 3.1 Objective function Given the magnetic field response, dobs, observed above a survey line, the goal of an EM inversion is to recover a smooth conductivity model that minimizes the data misfit. Similar to other EM inverse problems, the log of the conductivity m = log10(σ) is used and the data fitting is regularized by incorporating a priori information into the objective function. Here, we regularize with the most commonly used constraints in geophysical inversions: we require that the model is smooth and close to a reference model. To accomplish this, we add the Tikhonov regularization to the data misfit: min m φ(m) = min m φd(m) + βφm(m) (3.1) where φd(m) = 1 2 ∑ i ∥∥∥∥∥dpredi (m)− dobsiσstdi ∥∥∥∥∥ 2 2 (3.2) is the data misfit and i is an index for all combinations of ω and xs. φm(m) = 1 2 ( αs‖m−mref‖22 + αx‖∂x(m−mref)‖22 + αz‖∂z(m−mref)‖22 ) (3.3) is the model misfit and β is the regularization parameter. The predicted data, dpred, is computed by extracting the secondary field due to σ = em at the receivers using the measurement matrices Pi: dpredi = Pi ~H(xs, ω) = Pi Nky∑ j gjH̃(xs, ω, ky,j) (3.4) The data misfit term is the squared sum of the misfit of the data points normalized by the noise standard deviation σstdi . The second term is the smallness regularization and the corresponding parameter, αs, controls the length scale of the details in the inverted model. The smoothness in x and z-directions is controlled by the parameters αx and αz in the third and forth 33 terms, respectively and minimizing them has the effect of reducing sharp jumps in the solution so that a smooth model is produced. 3.2 Sensitivity matrix As we will discuss later, the objective function is minimized by a quasi- Newton optimization technique. In this section, we derive the gradient of the objective function needed to compute the search step for the optimization. To start, we discretize the objective function and obtain φ(m) = 1 2 ‖Q−1(dpred(m)− dobs)‖22 + β 2 ( αs‖m−mref‖22 + αx‖Dx(m−mref)‖22 + αz‖Dz(m−mref)‖22 ) . (3.5) where Q = diag(σstd), and Dx = D xe cc and Dz = D ze cc as shown in Chapter 2.4.2. Next, we differentiate the objective function with matrix calculus to find the gradient: ∇φ(m) = ( ∂dpred ∂m )T ( Q−2 ( dpred − dobs)) + β ( αs(m−mref) + αxDTxDx(m−mref) + αzDTz Dz(m−mref) ) = ∑ i Nky∑ j gjJ T i,jP T i Q −2 i (d pred i − dobsi ) + β ( αs(m−mref) + αxDTxDx(m−mref) + αzDTz Dz(m−mref) ) (3.6) In this equation, the sensitivity matrix Ji,j is given by the derivative of the magnetic field with respect to the model: Ji,j = ∂H̃ ∂m = ∂H̃ ∂σ ∂σ ∂m = CURLẼ −c ∂Ẽ ∂σ ∂σ ∂m (3.7) 34 where Ẽ = Ẽ(xs, ω, ky) is given by the forward solution: Ẽ = A−1(−cA2(σ)Ẽ0) (3.8) To find its derivative w.r.t to the conductivity model, we apply implicit differentiation to (2.48): ∂ ∂σ ( [CURLH̃CURLẼ + cdiag(A1σ)] Ẽ ) = ∂ ∂σ ( −cA2(σ)Ẽ0 ) (3.9) A ∂Ẽ ∂σ + c ∂ ∂σ ( diag(A1σ)Ẽ ) = ∂ ∂σ ( −cA2(σ)Ẽ0 ) (3.10) Now by applying the property diag(a)b = diag(b)a twice, once on each side of (3.10), we obtain: A ∂Ẽ ∂σ + c ∂ ∂σ ( diag(Ẽ)A1σ ) = −c ∂ ∂σ A ze ccdiag(Ẽ0,x) Andcc diag(Ẽ0,y) Axeccdiag(Ẽ0,z) σ (3.11) A ∂Ẽ ∂σ + cdiag(Ẽ)A1 = −cA3(Ẽ0) (3.12) ∂Ẽ ∂σ = −cA−1 ( A3(Ẽ0) + diag(Ẽ)A1 ) (3.13) where we have defined A3(Ẽ0) := A ze ccdiag(Ẽ0,x) Andcc diag(Ẽ0,y) Axeccdiag(Ẽ0,z) . (3.14) Now, by going back to Eqn. (3.7), we find that: Ji,j = ∂H̃ ∂m = CURLẼA −1 ( A3(Ẽ0) + diag(Ẽ)A1 ) diag(ln(10)10m) (3.15) 35 At last, we take the transpose of the sensitivity to obtain JTi,j = ( ∂H̃ ∂m )T = diag(ln(10)10m) ( AT3 (Ẽ0) + A T 1 diag(Ẽ) ) (AT ) −1 CURLT Ẽ (3.16) where Ẽ represents the complex conjugate of Ẽ. We plug this result back into (3.6) to get the gradient of our optimization objective. From Eqns. (3.6) and (3.16), we can see that each time the gradient is evaluated, two forward solves are needed: one computes the secondary electric field Ẽ by solving the regular forward problem given by Eqn. (3.8), and the other solves the adjoint problem given by: x = (AT ) −1 rhs (3.17) rhs = CURLT Ẽ PTi dpredi − dobsi σ2i,std (3.18) Together with the evaluation of the misfit function in each line search iter- ation, which we will discuss in the next section, the total number of solves involving A and AT is 3NωNky if one assumes that the misfit function only needs to be checked once for each new search direction. 3.3 Method of optimization In this work, the objective function is minimized using the quasi-Newton iterative method L-BFGS ([16]). This method finds the search direction by approximating the inverse Hessian using updates of the argument and the gradient. As in all iterative method, the optimization is initialized with a starting guess m0. Then, in each subsequent iteration, a new model is computed using mn+1 = mn + α`δmn (3.19) 36 where the line search parameter α` decreases exponentially with the line search iteration number `: α` = const ` , ` = 0, 1, ..., and δmn is the search direction given by the quasi-Newton approximate Hessian Hn and the gradi- ent: δmn = −H−1n ∇φ(mn). (3.20) Using the L-BFGS algorithm to approximate the inverse Hessian has the advantage that the sensitivity does not have to be computed explicitly but instead only matrix-vector products involving the Jacobian: JTi,jx, (3.21) are needed [8]. In situations where the full Hessian is too costly to be eval- uated, L-BFGS offers a cheap way to compute a good approximation to the inverse Hessian. At each L-BFGS iteration, the search step δmn, is com- puted from an initial guess of the Hessian, for which we use the Hessian of the regularization term Hm = αsI + αxD T xDx + αzD T z Dz, (3.22) and the model and gradient updates sl = ml+1 −ml (3.23) yl = ∇φ(ml+1)−∇φ(ml) (3.24) collected from the previous k iterations n − 1, ..., n − k. A history of only k = 10 recent update vectors usually does a good job in obtaining a descent search direction. For those interested, details on implementing the L-BFGS algorithm can be found in [12]. 37 3.4 Implementation of the optimization 3.4.1 Data structures In our inversion scheme, we seek to recover the 2D conductivity distribution of the earth directly underneath a survey line by simultaneously inverting soundings made along the line. Good data structures for storing the ob- served data and the electromagnetic fields are important to making the code readable and easy to run and maintain. The observed data are characterized by the source position, xs and the frequency ω. As an example, the data point measured with a source located at position number s at w-th frequency can be retrieved in Matlab with: dataObs(:, s, w). The fields, which are stored as column vectors, are functions of the wavenumbers ky as well. The most convenient way is to adapt a four di- mensional array structure for the storage. The first dimension will have the length of the stored variable: 1 for the observed data and the total number of edges/nodes/cell-centers for the fields. The second dimension is used to differentiate the sources, the third dimension separates the frequencies and the fourth dimension the different wavenumbers. To extract the field Ẽx computed for source number s, frequency number w and k-th wavenumber, we type the command: Ex{w, k}(:, s) (the cell structure is used to allow for storage of sparse matrices as Matlab does not support sparse matrices with more than two dimensions). 3.4.2 Constrained optimization By incorporating additional constraints into the optimization, the solution space can be further confined and the recovered solution will resemble more to a realistic real-world model. As mentioned earlier, constraints are de- rived from a priori information and they can come in many different types. 38 Most commonly, geophysical constraints come from borehole measurements, physical properties of the survey area extracted from previous geological or geophysical studies, and inversion results from other surveys. In geological settings, the range of observed conductivity/resistivity can span over eight orders of magnitude. Putting restrictions on the value of the conductivity during the optimization is therefore very helpful in producing a realistic model. In the algorithm, the restriction is implemented by the projected gradient method [2]. The method confines the conductivity in each cell of the recovered model to a given interval by applying a projection function to the updated model at each iteration: mn+1 = PΘ(mn + α`δmn) ,Θ = {m ∈ Rn : mmin ≤m ≤mmax} (3.25) where PΘ(m(xi, zj)) = { mmax(xi, zj) if m(xi, zj) > mmax(xi, zj) mmin(xi, zj) if m(xi, zj) < mmin(xi, zj) (3.26) with mmin and mmax being the user-defined lower and upper bound model, respectively. By bounding each cell individually, one can define the conduc- tivity in those areas that one has confidence in. For example, conductivity data from a drillhole sample can be used to narrow down the conductivity distribution in the region immediately around the borehole core, soil/rock samples on the surface can be used to make better inversion of the overbur- den and regional geophysical surveys can be used to place a global bound on the entire survey area. In addition to the ability to restrict the conductivity to be within a range of values, we can also define spatial regions in the model that can be up- dated in each iteration. This option is implemented in our algorithm by first defining an active region and then modifying the computed gradient as 39 follows: ∇φ(mn)modk = ∇φ(mn)k if k ↔ (i, j) and cell ij is in the active area0 otherwise. (3.27) By doing this, parameters outside the active region are not updated during the optimization and are left the same as in the initial guess. The use of active areas allows one to, for example, control the location in which a potential target can be found. In addition, it also provides an easy way to test the effect of placing a target in different locations on the data fitting. 3.4.3 Regularization parameter selection Choosing the appropriate regularization parameter, β, is essential for a good solution to the regularized minimization problem (3.1). Each strategy of reaching a solution to the problem produces a different β. In geophysical inversions, the most commonly used strategies are the discrepancy principle, the L-curve and generalized cross-validation (GCV). The first approach is usually used when the noise is Gaussian and its standard deviation is known. When the noise level cannot easily be determined, as is often the case, one can turn to the L-curve or GCV. In our algorithm, only the first two approaches are implemented. (The GCV approach requires knowing the full Hessian which is too computationally expensive to obtain for the size of our problem.) The discrepancy principle recovers a model by selecting the β that con- verges to a target data misfit, tmisfit, equal to tmisfit := φd(β ∗,m∗) = Ndata (3.28) where Ndata is the number of measurement points and the star ∗ denotes the optimum value of the variable. Constraining the data misfit in this way can be justified by the assumption of Gaussian noise and the fact that φd(m ∗) 40 follows a χ2 distribution and its expected value is E(φd(m ∗)) = Ndata. (3.29) If β is too small, the minimization will try to reduce φd further, resulting in a solution that fits the noise, i.e., over-fitting. On the other end, if β is too large the optimization renders an under-fitted and featureless model as only the model misfit is minimized and the recovery does not make much use of the observed information. The effect of the regularization parameter on the two misfits at convergence is illustrated by the well-known Tikhonov curve (Fig. 3.1). To find the β that gives tmisfit, it is advised to start with a large guess for β and then decrease its value in subsequent optimization trials until the desired target misfit is reached within a defined tolerance. This is because more iterations are required to reach convergence when φd is the dominant term. Also, since φm changes logarithmically with β, it is beneficial to generate a series of solutions for logarithmically spaced β’s. The search for β can be automated and the time it takes reduced by using a cooling strategy [7]. However, in the current version of the algorithm, β can only be chosen manually. In geophysical inversions, the errors come from uncertainties in the mea- surements as well as error in the forward solution. Estimating the error precisely can therefore be a tough task. In the cases where the noise level is uncertain, one can account for the uncertainty in the estimation by multi- plying the target misfit by a factor, chifact, such that the tmisfit is now tmisfit = chifact×Ndata. (3.30) Setting the chifact to above (below) 1 allows us to choose another model that fits to a larger (smaller) error. Also, the same strategy can be used to pick a model that underfits or overfits the current noise level. In the occasions when the noise level is completely unknown or the un- 41 β → ∞ β → 0 φm φd φd(β∗,m∗) Figure 3.1: Tikhonov curve. For large β, the model norm φm dom- inates the objective and the data misfit is underrepresented. The opposite happens when β is small. When the noise level is known, the target misfit, φd(β ∗,m∗), can be used to determine the regularization parameter. 42 log(φm) lo g (φ d ) β∗ β → ∞ β → 0 Figure 3.2: The L-curve. It can be used to select the regularization parameter β when the noise level is unknown. In this strategy, β is chosen at the point of maximum curvature of the L-curve. certainty in determining it is too high, β is generally chosen at the point of maximum curvature of the L-curve (see Fig. 3.2). The idea behind this is that in the absence of a good estimate of the noise, we prefer a model that both fits the data well and has a small model norm. In the current implemen- tation of the inversion, this method can only be used by manually generating a series of solutions for a range of β like in the discrepancy strategy. 43 3.4.4 Stopping criteria Once the quasi-Newton optimization is started, the inverted model is up- dated at each iteration until a stopping criteria or a convergence condition is reached. In our algorithm, the iteration is terminated when one of the following conditions is met: • the maximum number of quasi-Newton iterations, maxit, is exceeded • the size of the search step is small enough; max(abs(α`δmn)) ≤ tol step • the change in the objective function is small enough; |φ(mn+1)−φ(mn)|/|φ(mn)| ≤ tolφ • the norm of the gradient is small enough; ||∇φ(mn+1)||/||∇φ(m0)|| ≤ tol 44 Chapter 4 Results To evaluate the numerical accuracy of our forward modeling, a known 1D analytical solution was used to benchmark the 2.5D implementation. To assess the reliability of our inversion, we tested the algorithm on synthetic data and field data from a RESOLVE AEM survey. The inverted models were compared to the synthetic 2D models and also to models from a well- established 1D inversion supplied by Fugro. 4.1 Testing the forward solution The numerical accuracy of our algorithm was tested by computing the re- sponse from a synthetic layered earth model and then comparing the results against the 1D analytical solution [4, 23]. The model is shown in Fig. 4.1; it consists of a 20 m thick conductive layer of 10−1 S/m buried 30 m under the surface of a 10−2 S/m half-space. Figure 4.2 shows the 2D inphase and quadrature response. The frequencies used in the computation spanned a range that is typically used in AEM systems and the response was computed 10 m away from the transmitter and 40 m above the ground. Our results show that with only six optimized wavenumbers and a cell size of 10 m × 5 m, the error in computing the 2D response can be reduced 45 Tx Rx σ = 0.01 S/m σ = 0.1 S/m σ = 0.01 S/m 10 m 40 m 30 m 20 m 5 m 10 m 350 m 1200 m Figure 4.1: Schematic of the mesh and the conductivity model used to test the 2D forward modeling. to within 1% of the 1D analytical value. The response was also computed for different cell sizes and the discretization errors are listed in Table 4.1. For cell sizes under 10 m×10 m, the relative error averaged over all inphase and quadrature responses was found to be less than 1%. Since the measurement uncertainties in AEM systems are usually well above this level, these solution errors are sufficiently small for the purpose of inversion. Now that the forward problem can be accurately solved, we can proceed to test the inversion algorithm. 46 102 103 104 105 106 0 500 1000 1500 2000 2500 Frequency (Hz) In ph as e an d Q ua d (p pm ) 1D Inphase 1D Quad 2D Inphase 2D Quad 102 103 104 105 106 0 1 2 3 4 5 Frequency (Hz) R el at iv e er ro r ( % ) Inphase Quad a) b) Figure 4.2: Comparison of 2D response with 1D analytical solution. a) Computed responses from the synthetic model shown in Fig. 4.1. Circles are the 1D analytical responses and the crosses are the 2D modeled responses. b) Relative error of the 2D responses. Table 4.1: Solution error of the 2.5D forward modeling for different cell sizes. Cell size Max. Inphase error(PPM) Max. Quad error(PPM) Mean rel. error(%) 10×10 38.35(2.89%) 20.71(3.10%) 1.00 10×5 12.28(0.92%) 8.88(1.00%) 0.37 5×5 10.76(0.81%) 6.68(0.76%) 0.31 2.5×2.5 3.75(0.28%) 1.61(0.18%) 0.11 47 4.2 Inversion of synthetic data We next test our inversion code on synthetic data. Before inversions are performed on actual geophysical datasets, it is always a good idea to test the method on known solutions first. Using synthetic models that replicate complex geological structures and topology can help in finding limitations and shortcomings of the algorithm. Here, we tested our algorithm on a 200 m wide and 50 m deep conductive 2D prism that is buried inside a 0.01 S/m half space (see Fig. 4.3). The prism has a conductivity of 0.1 S/m and is buried 30 m under a 20 m thick overburden that is resistive (0.003 S/m) on one half and conductive (0.03 S/m) on the other. Two small inclines were added to the topology to mimic variations in the height of earth’s surface. In order to avoid committing the “inversion crime”, the synthetic response was generated using a twice denser mesh with cell size of 5 m×5 m. The response 10 m away from the transmitters was computed for 20 horizontal coplanar sources and 5 logarithmically spaced frequencies. Each source-receiver pair is separated by 20 m and is placed 40 m above the zero depth level. Gaussian noise with a standard deviation of 5% of the response value plus a floor of 10−10 A/m was used to contaminate the data before the inversion. In the inversion, the data were assigned a standard deviation of 8%+10−10 A/m and the mesh cell size of 10 m×10 m generates a forward system of about 8000 variables. The other inversion parameters were Nky = 6, αs = 10−4, αx = αz = 4 × 10−2. The reference and initial model were set to a 0.01 S/m half-space. The inversion was run using the Parallel Computing Toolbox in Matlab with 8 parallel processes and a laptop equipped with a 2.30GHz Intel i7-2820QM CPU. For a β of 1000 and a chifact of 1, the target misfit of 100 was reached in 23 iterations and 20 minutes. The results of inversion on the synthetic data are shown in Fig. 4.4. The data misfit showed that the optimization converged to a model that pre- dicted the majority of the observed response data to within the noise level, except for several points near the first source which can be explained by the 48 smaller data per parameter density on the edges. The recovered model re- vealed three distinct regions of conductivity contrast in accordance with the synthetic model. The relative position of the regions are consistent with the true model, but the sizes and conductivities displayed some discrepancies. A slightly better recovery was observed for the resistive half of the overbur- den than the conductive half, while in total only about 60-70% of the area of the true overburden was reconstructed with the correct property value. This is actually a decent result considering that only about a third of the length of the overburden was covered by the measurements. This finding is evidence that the range of detection for conductors of similar conductivity does not extend beyond 100 m at the frequencies tested, or in other words the footprint of the AEM system setup is at most 100 m. As for the inverted target, its conductivity plateaued at a value that is 6% less than the true value. From the peak value, the conductivity dropped continuously until it reached the reference value about two cells away from the edge of the true target, which puts the uncertainty on the position of the edge at about 20 m. This uncertainty is the effect of regularizing with the mean-squared error (Eqn. 3.3) which produces models with smooth edges by heavily penalizing sharp jumps in the model. Although it is possible to reduce it with a differ- ent regularizing function (such as the total variation), the uncertainty does not pose major problems for mineral targeting purposes, as the dimension of many economically viable ore deposits is often an order of magnitude larger. With these results, we are now more confident that the algorithm is func- tioning like a proper inversion is supposed to. However, before we can label it as a working code, we need to test the inversion once more, but now on a 3D field dataset. Details and results of this test are described in the next section. 49 −200 −100 0 100 200 300 400 500 600 −200 −150 −100 −50 0 50 100 150 200 x (m) D ep th (m ) −3 −2.5 −2 −1.5 −1 −0.5 0 Figure 4.3: Synthetic model used in testing the 2D inversion algo- rithm. The model is 2D and it consists of a conductive prism and an overburden that is conductive on one end and resistive on the other. Responses were recorded 10 m away from the sources and 40 m above the ground. Locations of the sources are indicated by the red crosses. 4.3 Inversion of field data To check the effectiveness of our algorithm in inverting field data, data from a Fugro RESOLVE survey were inverted. Since a unique solution does not exist, we can only verify the algorithm by comparing to results obtained with other well tested inversions. For this purpose, we compared our results to models from an established 1D inversion tool supplied by Fugro. In the survey, a total of 684 line kilometers were flown over the Red Dog mine site in Alaska. The survey lines were spaced 50 m apart and they 50 −200 −100 0 100 200 300 400 500 600 −200 −150 −100 −50 0 50 100 150 200 x (m) D ep th (m ) −3 −2.5 −2 −1.5 −1 −0.5 0 0 10 20 0 2000 4000 6000 8000 10000 12000 source # In ph as e (P P M ) 400Hz 1800 8200 40k 140k 0 10 20 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 source # Q ua d (P P M ) 400Hz 1800 8200 40k 140k 0 10 20 0 5 10 15 20 25 30 35 source # D at a m is fit n or m al iz ed b y st n de v +1 −1 a) b) Figure 4.4: a) Log of the recovered 2D conductivity model for the syn- thetic data generated by the model shown in 4.3. b) Left: misfit of predicted data normalized by noise standard deviation. Mid- dle and right: The response curves are offset by 1000 ppm for clarity. Crosses represent the synthetic inphase and quadrature responses and solid curves are the predicted data. 51 covered the site in two orthogonal directions: 13◦ and 103◦ azimuth. To map the conductivity contrasts of the survey area, the AEM system measured the magnetic response of the area in five frequencies at heights of 40-60 m above the ground and 7.9 m away from the sources. The objective of the survey was to use the conductivity detected in and around the mine’s tailing pond to identify any possible leakage of tailings from the waste area located to the east into the pond. Here, instead of running a full inversion using every data point of each line, the responses were resampled at 40 m intervals. Equivalently, the data were decimated to 7% of the original sampling rate. The area of inversion had a 200 m depth extent and it was discretized with a mesh of 8 m× 8 m cells. The paddings added 300 m to both ends of the line and 300 m to the depth of the mesh. The data were assigned a noise standard deviation of 10%×dobs+10−9 A/m and a half-space of 10−3 S/m was used as the reference model. By using a β determined from L-curve strategy and a convergence criteria of tol step = 0.02, convergence was usually reached within 50 iterations of L- BFGS with a history of 20 recent updates. Each L-BFGS iteration took about 3 mins on a cluster node outfitted with two 6-core Intel Xeon X5660 CPUs at 2.8 GHz and 64 GB of memory. 52 10300 Tailing pond waste area Figure 4.5: Satellite image of the Red Dog mine area in Alaska. Red line: outline of the survey area. White line: outline of the waste area. Blue line: outline of the tailing pond. Black lines: ore deposits. Green lines: survey lines inverted. Image from [21]. 53 7549000 7550000 7551000 7552000 0 2000 4000 6000 Northing, UTM 3N a) 7549000 7550000 7551000 7552000 0 1000 2000 3000 4000 5000 6000 Northing, UTM 3N 400Hz 1800 8200 40k 140k 588000 589000 590000 591000 0 2000 4000 6000 8000 10000 Easting, UTM 3N c) 588000 589000 590000 591000 0 1000 2000 3000 4000 5000 Easting, UTM 3N In ph as e (pp m ) In ph as e (pp m ) Qu ad (p pm ) Qu ad (p pm ) d) b) 8000 10000 400Hz 1800 8200 40k 140k Figure 4.6: Predicted data misfit for the inversions on the lines shown in Fig. 4.5: a)-b) Line 10300, c)-d) Line 20590. The curves are offset by 1000 ppm for clarity. Solid lines: predicted response. Crosses: observed data. Observed data are from [21]. 54 0 500 1000 1500 2000 2500 3000 3500 0 0.2 0.4 0.6 0.8 1 Distance from 1st source (m) A .U . Normalized bird altitude Normalized Inphase at 140kHz Figure 4.7: Normalized bird altitude and a measured response curve. The two observables are strongly correlated. 55 No rth in g (m ) Orthometric height (m) a) b) c) d) 75 49 00 0 75 49 50 0 75 50 00 0 75 50 50 0 75 51 00 0 75 51 50 0 75 52 00 0 10 0 20 0 30 0 - 3 - 2 - 10 - 3 - 2 - 10 Ea st in g (m ) 58 80 00 58 85 00 58 90 00 58 95 00 59 00 00 59 05 00 59 10 00 20 0 40 0 - 3 - 2 - 10 -3-2-10 Ea st in g (m ) 58 80 00 58 85 00 58 90 00 58 95 00 59 00 00 59 05 00 59 10 00 20 0 40 0 No rth in g (m ) 75 49 00 0 75 49 50 0 75 50 00 0 75 50 50 0 75 51 00 0 75 51 50 0 75 52 00 0 10 0 20 0 30 0 F ig u re 4 .8 : 2D (a ), c) ) an d 1D (b ), d )) in ve rt ed lo g co n d u ct iv it y m o d el s fo r th e tw o su rv ey li n es h ig h - li gh te d in F ig . 4. 5. a) -b ) L in e 10 30 0, c) -d ) L in e 20 59 0. 56 In the test, we inverted several 3 km long survey lines from both orthogo- nal directions. As seen in Fig. 4.6, the observed responses are highly varying between measurement points. A quick comparison of the length scale of these changes and that of fluctuations observed in the height of the bird (see Fig. 4.7) suggests that the variations are not due to instrument noise, but are rather caused by the changing bird altitude. As a consequence, one would expect the predicted data to exhibit the same fluctuations as in the observed data if the forward modeling takes the varying altitude into account. This was indeed the case: the predicted data showed a close fit to a large portion of the observations and they reproduced the fluctuations caused by the flight altitude. Inversions for the two lines highlighted in Fig. 4.5 are shown in Fig. 4.8. When compared to 1D inversions, our inverted models displayed the same ability to resolve areas of high conductivity contrasts. The conductive re- gions were found to be in the pond and on top of the waste dump area, which respectively correspond to the water in the pond and the conductive waste pile-up on the dump site. Features in the 1D models that are likely contributed by conductive geological structures were all recovered in the 2D models. However, most of these structures appeared smoother in our 2D inversions. We relate this to our choice of large smoothness parameters αx and αz, which have the effect of penalizing sharp features in the solution. Despite the size of the anomalies showed some discrepancies, their conduc- tivity values found in the 2D models were in good agreement with the 1D results as well as their location and depth extent. As for locating possible channels of tailing drainage, no clear evidence of their presence was found in any of the lines inverted. But like every good in- version work-flow will suggest, inversion should always be repeated multiple times such that decisions are made based on a collection of models rather than on a single model. Another run of inversion with better adjusted reg- ularization parameters and smaller cell sizes is therefore necessary to obtain 57 a more conclusive answer. To summarize the test results, our implementation of the 2.5D EM for- ward modeling and inversion scheme were successful when verified against more established 1D counterparts. The forward modeling was able to com- pute the magnetic response to a level of accuracy (±1%) that is sufficiently high for inversion purposes. In addition, there is strong evidence that the inversion algorithm is correctly resolving conductive anomalies at least as accurately as a standard 1D inversion is doing. 58 Chapter 5 Conclusion and future work In this project, we developed an inversion algorithm to recover a 2D conduc- tivity model of the earth from frequency domain AEM data generated by horizontally coplanar vertical dipoles. We demonstrated that our implemen- tation of the 2.5D EM problem can accurately predict the magnetic response of a three layer model and the inversion scheme can successfully and reliably be applied to field data to correctly locate and identify conductive anomalies. However, the computation speed still poses a big problem and improve- ments are necessary to give our algorithm a more definite edge over other 1D or 3D inversion tools that are now available. This could be done by using adaptive meshes [10] which can greatly reduce the size of the linear system. The technique saves computing time by starting out with a coarse grid, which is then successively refined for only those regions in the grid that are rich in details. Another solution to the problem is by speeding up the matrix decomposition in the direct solve step. This can be achieved by using parallelized sparse direct solver algorithms such as MUMPS and SuperLU. Moreover, whether the inversion is effective in correctly resolving 3D structures is still a question to be answered. As one might have quickly noticed by inspecting the aerial map of the mine site in Chapter 4, the conductive structures under the two inverted lines are rather 2D in nature. 59 Therefore in the future, in order to make more correct interpretations of higher dimensional geologies, our algorithm should also be tested on data produced by strongly 3D structures, for example by a conductive block with a strike extent comparable to the size of the AEM footprint. For future work, a direct continuation of this project would be to imple- ment the same 2.5D formulation but for time domain EM data. For this, the same inversion scheme shown earlier can be applied once the forward problem is solved. Solving time domain problems is however a little more involved as it requires time stepping, which cannot be easily parallelized. Alternatively, the time domain problem can be reformulated and solved in the frequency domain by applying a Fourier transform. But this approach is only advantageous if the number of frequencies required for an accurate transformation can be kept small. 60 Bibliography [1] A. Abubakar, T. M. Habashy, V. L. Druskin, L. Knizhnerman, and D. Alumbaugh. 2.5D forward and inverse modeling for interpreting low-frequency electromagnetic measurements. Geophysics, 73: F165–F177, 2008. → pages 7 [2] P. H. Calamai and J. J. Moré. Projected gradient methods for linearly constrained problems. Mathematical Programming, 39:93–116, 1987. → pages 39 [3] L. H. Cox, G. A. Wilson, and M. S. Zhdanov. 3D inversion of airborne electromagnetic data using a moving footprint. Exploration Geophysics, 41:250–259, 2010. → pages 5 [4] A. Dey and S. H. Ward. Inductive sounding of a layered earth with a horizontal magnetic dipole. Geophysics, 35:660703, 1970. → pages 5, 45 [5] R. G. Ellis. Inversion of airborne electromagnetic data. 68th SEG meeting, Expanded Abstracts, pages 2016–2019, 1998. → pages 5 [6] M. Everett and R. Edwards. Transient marine electromagnetics: the 2.5-D forward problem. Geophys. J. Int., 113:545–561, 1992. → pages 6 [7] E. Haber. Numerical Strategies for the solution of inverse problems. PhD thesis, UBC, 1997. → pages 41 [8] E. Haber, U. Ascher, and D. Oldenburg. On optimization techniques for solving nonlinear inverse problems. Inverse Problems, 16: 1263–1280, 2000. → pages 37 61 [9] E. Haber, U. Ascher, and D. Oldenburg. Inversion of 3D electromagnetic data in frequency and time using an inexact all-at-once approach. Geophysics, 69:1216–1228, 2004. → pages 5 [10] E. Haber, S. Heldmann, and U. Ascher. Adaptive finite volume methods for distributed non-smooth parameter identification. Inverse Problems, 23:1659–1676, 2007. → pages 22, 59 [11] J. D. Jackson. Classical Electrodynamics. Wiley, 3 edition, 1999. → pages 14 [12] C. T. Kelley. Iterative Methods for Optimization. Siam, 1999. → pages 37 [13] K. Key and J. Ovall. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling. Geophys. J. Int., 186: 137–154, 2011. → pages 7 [14] A. Ley-Cooper, J. Macnae, and A. Viezzoli. Breaks in lithology: Interpretation problems when handling 2D structures with a 1d approximation. Geophysics, 75:WA179–WA188, 2010. → pages 5 [15] Y. Mitsuhata, T. Uchida, and H. Amano. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source. Geophysics, 67:1753–1768, 2002. → pages 7 [16] J. Nocedal and S. Wright. Numerical Optimization. Springer, New York, 1999. → pages 36 [17] A. Pidlisecky and R. Knight. FW2 5D: A MATLAB 2.5-D electrical resistivity modeling code. Computers & Geosciences, 34:1645, 2008. → pages 14 [18] J. Ryu, H. F. Morrison, and S. H. Ward. Electromagnetic fields about a loop source of current. Geophysics, 35:862–896, 1970. → pages 5 [19] B. R. Spies and F. C. Frischknecht. Electromagnetic sounding: in Nabighian, M. N., editor, Electromagnetic Methods in Applied Geophysics–Applications. Soc. of Expl. Geophys., 1 edition, 1991. → pages 15 62 [20] C. Stoyer and R. J. Greenfield. Numerical solutions of the response of a two-dimensional earth to an oscillating magnetic dipole source. Geophysics, 41:519–530, 1976. → pages 6 [21] Teck Resources Ltd. Technical report. → pages viii, 53, 54 [22] M. J. Unsworth, B. J. Travis, and A. D. Chave. Electromagnetic induction by a finite electric dipole source over a 2-d earth. Geophysics, 58:198–214, 1993. → pages 7 [23] S. H. Ward and G. W. Hohmann. Electromagnetic theory for geophysical applications: in Nabighian, M. N., editor, Electromagnetic Methods in Applied Geophysics-Applications. Soc. Expl. Geophys., 1988. → pages 5, 10, 45 [24] G. A. Wilson, A. P. Raiche, and F. Sugeng. 2.5D inversion of airborne electromagnetic data. Exploration Geophysics, 37:363–371, 2006. → pages 7 [25] X. Xie, J. Macnae, and J. Reid. The limitations of 1-D aem inversion for 2-D overburden structures. 68th SEG meeting, Expanded Abstracts, pages 760–763, 1998. → pages 5 [26] S. Xu, B. Duan, and D. Zhang. Selection of the wavenumbers k using an optimization method for the inverse fourier transform in 2.5D electrical modelling. Geophysical Prospecting, 48:789–796, 2000. → pages 14 [27] D. Yang and D. W. Oldenburg. Three-dimensional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit. Geophysics, 77:B23–B34, 2012. → pages 5 63