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, 2012Abstract 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 eld dataset from a RESOLVE survey and compare the inversion results to those obtained from a well-established 1D implementation. iiTable 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 di erential operators . . . . . . . 19 2.4.3 Discretization of the material averaging . . . . . . . . . 25 iii2.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 eld data . . . . . . . . . . . . . . . . . . . . . . . 50 5 Conclusion and future work . . . . . . . . . . . . . . . . . . . 59 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 ivList 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 di erent cell sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 vList 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 eld from the transmitter induces eddy currents in the conductor and a secondary eld is produced as a result. The induced eld 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 ight direction and depth. Each conductor has a conductivity of i and its extent in the strike direction is in nite. . . . . . . . . . 7 Figure 2.1 Discretization of the electromagnetic eld and the conduc- tivity model on a 2D staggered grid. Black dots denote the location of the nodes, boxes are the x-edges (de ned as edges whose normal points in the x-direction), diamonds are the z-edges and crosses represent the cell-centers. . . . 21 viFigure 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 mis t is underrepresented. The opposite happens when is small. When the noise level is known, the target mis t, 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 viiFigure 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: mis t of predicted data normalized by noise standard de- viation. Middle and right: The response curves are o set 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 mis t for the inversions on the lines shown in Fig. 4.5: a)-b) Line 10300, c)-d) Line 20590. The curves are o set 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 viiiAcknowledgments 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 di culties 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 eld 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 eld of applied mathematics and computer science, as well as in geophysics and geology. ixChapter 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 elds 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 in nitely 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 1reasonable 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- ts to the observed data and satis es 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- 2mitter loop to produce a magnetic dipole eld. If a conductive body is buried in the ground, the primary eld will induce eddy currents in the conductor and a secondary magnetic eld will be generated as a result. This secondary eld changes the magnetic ux through the receiver loop and its magnitude can be quanti ed 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 elds. In typical frequency-domain AEM systems, each conductivity depth sound- ing (mapping of conductivity as a function of depth) consists of magnetic eld measurements at ve 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 di usion length of the magnetic eld, which is given by the skin depth = r 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 pro le the conductivity distribution as a function of depth. To map the horizontal direction, soundings are carried out at regularly spaced intervals along a ight 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 own over the area. The line separation can range from tens to hundreds of meters, depending on the size of the geological features under study. 3bird 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 eld from the transmitter induces eddy currents in the conduc- tor and a secondary eld is produced as a result. The induced eld measured at the receiver loop is then used to determine the conductivity of the body. 41.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 di erent source locations. Compared to ground survey methods, this is about two orders of magnitude more sources. The rst 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 nite 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 eld measurement does not extend notice- ably into the adjacent lines and so each ight line can e ectively 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 ight direction and depth. Each conductor has a conductivity of i and its extent in the strike direction is in nite. 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 - nite element (FEM) or nite di erence (FD) approaches. But none used the nite 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 nite volume approach to implement an algo- rithm for inversion of frequency-domain AEM data. The algorithm, which is developed in Matlab, was rst tested on synthetic data generated by a 2D model and later on a 3D AEM eld 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. 71.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 eld data. A conclusion and discussion of future work can be found in the last chapter. 8Chapter 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 elds are rewritten as the sum of 2D elds 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 nite volume discretization of the elec- tromagnetic eld 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 eld. 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 eld and the primary eld produced by the transmitter in free space. For horizontal coplanar loops, the response is given 9by Resp = ~Hz(rtr) ~Hp;z(rtr) j ~Hp;z(rtr)j ; ~Hp;z(r) = m 4 r3 (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 eld and ~Hz(rtr) is the total eld measured at the receiver. To model the magnetic response due to harmonically changing electro- magnetic elds, ~F (~x; t) = ~F (~x)e i!t, we solve Maxwell’s equations in the frequency domain: r ~E i! 0 ~H = ~Ms (2.2a) r ~H ~E = ~Js on (2.2b) n^ ~H = 0 on @ (2.2c) where we have imposed a vanishing magnetic eld 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 elds, ! is the angular frequency of the elds, 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 simpli ed to the real-valued conductivity . From now on, we will use for the conductivity. To solve Maxwell’s equations, we decompose the total elds ~E; ~H into primary and secondary parts: ~E = ~E0 + ~E1; ~H = ~H0 + ~H1; (2.3) so (r ~E0 i! 0 ~H0) + (r ~E1 i! 0 ~H1) = 0 (2.4a) (r ~H0 ~E0) + (r ~H1 ~E1) = ~Js on (2.4b) 10where we have set ~Ms = 0, since only electric sources are used in our system. This decomposition separates Maxwell’s equations into a primary problem, r ~E0 i! 0 ~H0 = 0 (2.5a) r ~H0 = ~Js (2.5b) and a secondary problem, r ~E1 i! 0 ~H1 = 0 (2.6a) r ~H1 ~E1 = ~E0 (2.6b) Setting ~Js to the current in a HCP, the primary problem describes the eld 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 eld. 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 eld dominates the induced eld 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 in nite extent in the strike direction y: = (x; z). Despite the 2D con- ductivity, each eld component remains 3-dimensional as the source eld 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 eld components are either symmetric or antisymmetric in y. Therefore, we can either use the sine or cosine integral 11for the transformation: ~F (x; ky; z) = 8 < : i R1 0 ~F (x; y; z) sin(kyy)dy if ~F is anti-symmetric in y R1 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: ~rky ~E i! 0 ~H = 0 (2.8a) ~rky ~H ~E = ~E0 (2.8b) where ~rky is the 2.5D curl operator in the wavenumber domain and it acts on the eld ~F as follows: ~rky ~F = 2 6 4 iky ~Fz @z ~Fy @z ~Fx @x ~Fz @x ~Fy + iky ~Fx 3 7 5 (2.9) where ~F = ( ~Fx; ~Fy; ~Fz). The Fourier transform of the source eld, ~E0, 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) = 8 < : 2i R1 0 ~F (x; ky; z) sin(kyy)dky if ~F is anti-symmetric in ky 2 R1 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 12transforms: ~F (x; y; z) = 8 < : 2i PN j=1 ~F (x; ky;j; z) sin(ky;jy)gj if ~F is anti-symmetric in ky 2 PN 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 elds ~E0; ~H0 are the free-space electromagnetic elds due to the current density in the source loop, ~Js: r ~E0 i! 0 ~H0 = 0 (2.12a) r ~H0 = ~Js (2.12b) Instead of solving (2.12) directly, the primary eld ~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 eld 13is given by the vector potential of the dipole, ~A [11]: ~E0(x; y; z) = @t ~A = i! ~A = i!m 0 4 r3 2 6 4 y x xs 0 3 7 5 (2.13) where r = p (x xs)2 + y2 + (z hs)2. For every source on the same survey line, the Fourier transform of ~E0 is given by ~E0(x; ky; z) = i! 0m 4 2 6 4 ikyK0(jkyjs) jkyjK1(jkyjs)=s 0 3 7 5 (2.14) where s = p (x xs)2 + (z hs)2 and Kn is the n-th order modi ed 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 xed 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 nding the ky;j’s and gj’s that minimize the error in approximating a vertical dipole eld 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 eld due to a homogeneous half- space. For the range of conductivities and frequencies typically found in 14AEM surveys, the low-frequency approximation applies, i.e., jkrj 1, where k = p i ! 0 and r is the distance from the center of the source. The quasi-static eld 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 4 r3 1 k2r2 4 (2.15) where r = p 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 4 jkyj x K1(jkyjx) k2 4 K0(jkyjx) (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 X 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 ~H lz(xi; 0; 0) X 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 k k22 = 1 2 k1 Mgk22 (2.19) where Mij = 2 ~H lz(xi; ky;j; 0) ~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. 15To 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 k1 M(MTM) 1MT1k22 = 1 2 k1 v(ky)k22 (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 k1 v(k(0)y ) @v @ky kyk22 (2.24) gives the search direction ky: ky = (JTJ) 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 di erentiation: J(i) = @v @ky (i) = h v k(i)y;1 : : : v k(i)y;N i ; 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 nite volume approach. The linear system to be solved is set up by discretizing the eld 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 elds ~Ex and ~Ez on the edges of each cell, as shown in Fig. (2.1). This con guration places constraints on the possible locations for the rest of the eld quantities. To see this, we expand (2.8) to obtain the system of di erential equations for the eld components: 17 iky ~Ez @z ~Ey i! 0 ~Hx = 0 (2.28a) @z ~Ex @x ~Ez i! 0 ~Hy = 0 (2.28b) @x ~Ey + iky ~Ex i! 0 ~Hz = 0 (2.28c) iky ~Hz @z ~Hy ~Ex = ~E0;x (2.28d) @z ~Hx @x ~Hz ~Ey = ~E0;y (2.28e) @x ~Hy + iky ~Hx ~Ez = ~E0;z (2.28f) Firstly, Eqns. (2.28a) and (2.28c) imply that the elds ~Hx and ~Hz should be discretized at the same location as ~Ez and ~Ex, respectively. Secondly, Eqn. (2.28c) places the x-directional derivative of ~Ey on the z-edges (edges whose normal points in z-direction). And since the two-point formula is used for the derivative, this con guration puts ~Ey 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 ~Hy should be discretized at the cell-centers. Unlike the secondary elds, which can only be approximated numerically, the source eld ~E0 is continuous and its value is known everywhere. However, although ~E0 is an electric eld, we will not discretize it like ~E. The reason is that ~E0 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 ~Hz can be computed at the receiver center without averaging, we should avoid discretizing ~E0 on z-edges. As a result, we will evaluate ~E0 at the cell-centers instead, and then use averaging to bring the current density, ~E0, to the same mesh as the electric eld. 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 18discontinuous 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 di erential operators To construct the discrete version of the 2.5D curl operator (2.9), CURL, and the gradient operator, GRAD, we discretize the di erence 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 di erence 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 di erence which maps nodes to cell-centers. And similarly, the edge!node and cell-center!edge mappings follow the structure of the 1D cell-center di erences. Moreover, the boundary conditions (BC) given in (2.2c) require a Dirichlet type BC for the magnetic eld quantities, while the electric elds follow BC of Neumann type. From this information, we can see that the 1D cell-center di erences with Dirichlet BC and the 1D node di erences with Neumann BC are needed to build CURL and GRAD. In this section, we will demonstrate how to discretize the 1D di erences 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 19Table 2.1: Overview of the discretization of the physical quantities and the operators used in the Maxwell’s equations Quantity Notation Discretization electric eld 2 4 ~Ex ~Ey ~Ez 3 5 2 4 z-edge node x-edge 3 5 magnetic eld 2 4 ~Hx ~Hy ~Hz 3 5 2 4 x-edge cell-center z-edge 3 5 source electric eld 2 4 ~E0;x ~E0;y ~E0;z 3 5 2 4 cell-center cell-center cell-center 3 5 conductivity cell-center Operator Notation Mapping curl of electric eld CURL ~E 2 4 z-edge node x-edge 3 5! 2 4 x-edge cell-center z-edge 3 5 curl of magnetic eld CURL ~H 2 4 x-edge cell-center z-edge 3 5! 2 4 z-edge node x-edge 3 5 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 ~Ex; ~Ez edge ! cell-center ~Ey node ! edge ~Hx; ~Hz edge ! node ~Hy cell-center ! edge 20x z 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 eld and the conduc- tivity model on a 2D staggered grid. Black dots denote the loca- tion of the nodes, boxes are the x-edges (de ned as edges whose normal points in the x-direction), diamonds are the z-edges and crosses represent the cell-centers. 21column 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 eld quantities are stored as: ~F = 2 6 4 ~Fx ~Fy ~Fz 3 7 5 (2.30) where each 2D eld 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 di erence operator: D1Dx;node!cell-center = 2 6 6 6 6 6 4 1 hx;1 1 hx;1 1 hx;2 1 hx;2 . . . . . . 1 hx;m 1 hx;m 3 7 7 7 7 7 5 ; (2.31) and the Kronecker product, A B := 2 6 6 6 6 6 4 a11B a12B : : : a1;nB a21B a22B . . . ... ... . . . . . . an 1;nB an;1B : : : an;n 1B annB 3 7 7 7 7 7 5 (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) 22where In is the n n identity matrix and we have used the second-order two-point formula d dx u(xi+ 12 ; 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 di erence 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 = 2 6 6 6 6 6 6 6 4 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 3 7 7 7 7 7 7 7 5 ; (2.36) to build the 2D di erence 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 23condition @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 eld ~E: CURL ~E(ky) = 2 6 4 0 Dxend ikyI Dccze 0 D cc xe ikyI Dzend 0 3 7 5 (2.38) and the curl for the magnetic eld ~H: CURL ~H(ky) = 2 6 4 0 Dzecc ikyI Dndxe 0 D nd ze ikyI Dxecc 0 3 7 5 : (2.39) For clarity, we have simpli ed 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 = (D 1D x;cc!nd) T , and the two curl operators become adjoints of each other: CURL ~H = CURL T ~E: (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) 242.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 rst-order 1D cell-center averaging: i;j+ 12 = i 12 ;j+ 12hi + i+ 12 ;j+ 12hi 1 hi 1 + hi +O(h): (2.42) However, when the conductivity has jumps at the cell boundaries where its derivative is not de ned, 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 di erence. And so to obtain the 2D average, we use the sparsity pattern of Eqn. (2.36) to get A1Dx;cell-center!node = 2 6 6 6 6 6 6 6 4 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 3 7 7 7 7 7 7 7 5 ; (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 A 1D x;cell-center!node (2.44c) 252.4.4 Matrix representation of Maxwell’s equations For each combination in the set of parameters fxs; !; kyg, the following dis- crete 2.5D Maxwell’s equations have to be solved: CURL ~E ~E i! 0 ~H = 0 (2.45a) CURL ~H ~H diag(A1 )~E = A2( )~E0 (2.45b) where A1 := 2 6 4 Azecc Andcc Axecc 3 7 5 (2.46) and A2( ) = 2 6 4 Azeccdiag( ) 0 0 0 Andcc diag( ) 0 0 0 Axeccdiag( ) 3 7 5 : (2.47) To solve the system, we can isolate for ~E in (2.45) and rearrange the equations to obtain: [CURL ~HCURL ~E + cdiag(A1 )] ~E = cA2( )~E0 (2.48a) ~H = CURL ~E~E= c (2.48b) where c = i! 0. Once the matrix A = CURL ~HCURL ~E + cdiag(A1 ) is set up, the system can conveniently be solved with the Matlab backslash: ~E = An( cA2( )~E0). Using (2:48) to solve for all eld 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 elds. This approach reduces the size of the linear system to be solved. Once a solution is found, the other eld components can be computed straightfor- wardly using numerical di erentiation. To demonstrate, we start by writing 26out explicitly the matrices acting on each eld component: iky~Ez Dxend ~Ey i! 0 ~Hx = 0 (2.49a) Dccze ~Ex D cc xe ~Ez i! 0 ~Hy = 0 (2.49b) Dzend ~Ey + iky~Ex i! 0 ~Hz = 0 (2.49c) iky ~Hz Dzecc ~Hy diag(A ze cc )~Ex = A ze ccdiag( )~E0;x (2.49d) Dndxe ~Hx D nd ze ~Hz diag(Andcc )~Ey = A nd cc diag( )~E0;y (2.49e) Dxecc ~Hy + iky ~Hx diag(A xe cc )~Ez = A xe ccdiag( )~E0;z (2.49f) Now, after some algebra, we can eliminate all eld components except ~Ey and ~Hy from the equations and express the linear system for these remaining unknowns as " ~Ey ~Hy # = " S11 S12 S21 S22 # 1 " ikyDndze 1 ze I ikyD nd xe 1 xe Dccze 1 ze 0 D cc xe 1 xe # A2( )~E0 (2.50) Here, the blocks of S are de ned as S11 = Dndxe xeD xe nd + D nd ze zeD ze nd diag Andcc (2.51a) S12 = ikyDndxe 1 xe D xe cc ikyD nd ze 1 ze D ze cc (2.51b) S21 = iky c Dccxe 1 xe D xe nd + iky c Dccze 1 ze D ze nd (2.51c) S22 = Dccxe 1 xe D xe cc + D cc ze 1 ze D ze cc Im n (2.51d) where xe = diag (Axecc ), 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 27Eqn. (2.50) is solved, the remaining eld components can be recovered from " ~Ex ~Ez # = " 1ze 0 0 1xe #" ikyDzend cD ze cc ikyDxend cD xe cc #" ~Ey ~Hy # + c " 1ze 0 0 0 0 1xe # A2( )~E0; (2.52) " ~Hx ~Hz # = " xeDxend iky 1 xe D xe cc zeDzend iky 1 ze D ze cc #" ~Ey ~Hy # + iky " 0 0 1xe 1ze 0 0 # A2( )~E0: (2.53) Compared to solving for all components simultaneously, this method only requires solving for two eld components. Cutting the number of variables by a third signi cantly 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 rst reduced and then discretized, the order of accuracy of the remaining eld components is not guaranteed to be the same as that of ~Ey and ~Hy. 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 e ective 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 elds 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) 28This 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~ECURL T ~H cdiag(A1 ) = 2 6 4 0 Dxend iky Dccze 0 D cc xe iky Dzend 0 3 7 5 T 2 6 4 0 Dzecc iky Dndxe 0 D nd ze iky Dxecc 0 3 7 5 T cdiag(A1 ) = 2 6 4 0 (Dccze) T iky (Dxend) T 0 (Dzend) T iky (Dccxe) T 0 3 7 5 2 6 4 0 (Dndxe ) T iky (Dzecc) T 0 (Dxecc) T iky (Dndze ) T 0 3 7 5 cdiag(A1 ) (2.55) And now by noticing that the blocks of CURLT~E and CURL ~H have the same sparsity pattern, as well as those of CURLT~H and CURL ~E, 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 = (D A B) 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 identi ed 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 29Table 2.3: Substitution of operators for the setup of the adjoint prob- lem A ! AT Dzecc ! (D cc ze) T Dndxe ! (D xe nd) T Dndze ! (D ze nd) T Dxecc ! (D cc xe) T Dxend ! (D nd xe ) T Dccze ! (D ze cc) T Dccxe ! (D xe cc) T Dzend ! (D nd ze ) 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 Trhs) 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 e ort 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 di erent 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- 300 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 e ects, are not shown in the gure. 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 di erence of a meter in the altitude of the source can a ect 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. 31Chapter 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 de ned 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 di erent application dependent requirements on the solution can be achieved by adding the appropriate regularizations to the data- tting via mathematical constraints. 323.1 Objective function Given the magnetic eld response, dobs, observed above a survey line, the goal of an EM inversion is to recover a smooth conductivity model that minimizes the data mis t. Similar to other EM inverse problems, the log of the conductivity m = log10( ) is used and the data tting 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 mis t: min m (m) = min m d(m) + m(m) (3.1) where d(m) = 1 2 X i dpredi (m) d obs i stdi 2 2 (3.2) is the data mis t and i is an index for all combinations of ! and xs. m(m) = 1 2 skm mrefk 2 2 + xk@x(m mref)k 2 2 + zk@z(m mref)k 2 2 (3.3) is the model mis t and is the regularization parameter. The predicted data, dpred, is computed by extracting the secondary eld due to = em at the receivers using the measurement matrices Pi: dpredi = Pi ~H(xs; !) = Pi NkyX j gj ~H(xs; !; ky;j) (3.4) The data mis t term is the squared sum of the mis t 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 33terms, respectively and minimizing them has the e ect 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 kQ 1 dpred(m) dobs k22 + 2 skm mrefk22 + xkDx(m mref)k 2 2 + zkDz(m mref)k 2 2 : (3.5) where Q = diag( std), and Dx = Dxecc and Dz = D ze cc as shown in Chapter 2.4.2. Next, we di erentiate the objective function with matrix calculus to nd the gradient: r (m) = @dpred @m T Q 2 dpred dobs + s(m mref) + xDTxDx(m mref) + zD T z Dz(m mref) = X i 0 @ NkyX j gjJ T i;jP T i Q 2 i (d pred i d obs i ) 1 A + s(m mref) + xDTxDx(m mref) + zD T z Dz(m mref) (3.6) In this equation, the sensitivity matrix Ji;j is given by the derivative of the magnetic eld with respect to the model: Ji;j = @ ~H @m = @ ~H @ @ @m = CURL ~E c @~E @ @ @m (3.7) 34where ~E = ~E(xs; !; ky) is given by the forward solution: ~E = A 1( cA2( )~E0) (3.8) To nd its derivative w.r.t to the conductivity model, we apply implicit di erentiation to (2.48): @ @ [CURL ~HCURL ~E + cdiag(A1 )] ~E = @ @ cA2( )~E0 (3.9) A @~E @ + c @ @ diag(A1 )~E = @ @ cA2( )~E0 (3.10) Now by applying the property diag(a)b = diag(b)a twice, once on each side of (3.10), we obtain: A @~E @ + c @ @ diag(~E)A1 = c @ @ 0 B @ 2 6 4 Azeccdiag(~E0;x) Andcc diag(~E0;y) Axeccdiag(~E0;z) 3 7 5 1 C A (3.11) A @~E @ + cdiag(~E)A1 = cA3(~E0) (3.12) @~E @ = cA 1 A3(~E0) + diag(~E)A1 (3.13) where we have de ned A3(~E0) := 2 6 4 Azeccdiag(~E0;x) Andcc diag(~E0;y) Axeccdiag(~E0;z) 3 7 5 : (3.14) Now, by going back to Eqn. (3.7), we nd that: Ji;j = @ ~H @m = CURL ~EA 1 A3(~E0) + diag(~E)A1 diag(ln(10)10m) (3.15) 35At last, we take the transpose of the sensitivity to obtain JTi;j = @ ~H @m !T = diag(ln(10)10m) AT3 (~E0) + A T 1 diag(~E) (AT ) 1 CURLT~E (3.16) where ~E represents the complex conjugate of ~E. 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 eld ~E 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~E P T i dpredi d obs i 2i;std (3.18) Together with the evaluation of the mis t 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 mis t 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 nds 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) 36where 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 r (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 o ers 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 + xDTxDx + zD T z Dz; (3.22) and the model and gradient updates sl = ml+1 ml (3.23) yl = r (ml+1) r (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]. 373.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 elds 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 elds, 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 rst 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 elds. The second dimension is used to di erentiate the sources, the third dimension separates the frequencies and the fourth dimension the di erent wavenumbers. To extract the eld ~Ex computed for source number s, frequency number w and k-th wavenumber, we type the command: Exfw; kg(:; 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 con ned 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 di erent types. 38Most 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 con nes 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) ; = fm 2 Rn : mmin m mmaxg (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-de ned lower and upper bound model, respectively. By bounding each cell individually, one can de ne the conduc- tivity in those areas that one has con dence 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 de ne spatial regions in the model that can be up- dated in each iteration. This option is implemented in our algorithm by rst de ning an active region and then modifying the computed gradient as 39follows: r (mn)modk = 8 < : r (mn)k if k $ (i; j) and cell ij is in the active area 0 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 e ect of placing a target in di erent locations on the data tting. 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 di erent . In geophysical inversions, the most commonly used strategies are the discrepancy principle, the L-curve and generalized cross-validation (GCV). The rst 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 rst 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 mis t, 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 mis t in this way can be justi ed by the assumption of Gaussian noise and the fact that d(m ) 40follows 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 ts the noise, i.e., over- tting. On the other end, if is too large the optimization renders an under- tted and featureless model as only the model mis t is minimized and the recovery does not make much use of the observed information. The e ect of the regularization parameter on the two mis ts at convergence is illustrated by the well-known Tikhonov curve (Fig. 3.1). To nd 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 mis t is reached within a de ned 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 bene cial 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 mis t 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 ts to a larger (smaller) error. Also, the same strategy can be used to pick a model that under ts or over ts 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 mis t is underrepresented. The opposite happens when is small. When the noise level is known, the target mis t, d( ;m ), can be used to determine the regularization parameter. 42log(φ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 ts 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. 433.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; j (mn+1) (mn)j=j (mn)j tol the norm of the gradient is small enough; jjr (mn+1)jj=jjr (m0)jj tol 44Chapter 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 eld 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 45Tx 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 di erent 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 su ciently small for the purpose of inversion. Now that the forward problem can be accurately solved, we can proceed to test the inversion algorithm. 46102 103 104 105 106 0 500 1000 1500 2000 2500 Frequency (Hz) Inphase and Quad (ppm ) 1D Inphase 1D Quad 2D Inphase 2D Quad 102 103 104 105 106 0 1 2 3 4 5 Frequency (Hz) Relative error (% ) 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 di erent 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 474.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 rst. Using synthetic models that replicate complex geological structures and topology can help in nding 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 oor 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 mis t 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 mis t 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 rst source which can be explained by the 48smaller 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 nding 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 e ect 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 di er- 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 con dent 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 eld 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) Depth (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 eld data To check the e ectiveness of our algorithm in inverting eld 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 own 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) Depth (m ) −3 −2.5 −2 −1.5 −1 −0.5 0 0 10 20 0 2000 4000 6000 8000 10000 12000 source # Inphase (PPM ) 400Hz 1800 8200 40k 140k 0 10 20 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 source # Quad (PPM ) 400Hz 1800 8200 40k 140k 0 10 20 0 5 10 15 20 25 30 35 source # Data misfit normalized by stn dev +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: mis t of predicted data normalized by noise standard deviation. Mid- dle and right: The response curves are o set by 1000 ppm for clarity. Crosses represent the synthetic inphase and quadrature responses and solid curves are the predicted data. 51covered 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 ve 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 out tted with two 6-core Intel Xeon X5660 CPUs at 2.8 GHz and 64 GB of memory. 5210300 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]. 537549000 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 Inph as e (p pm ) Inph as e (p pm ) Quad (p pm ) Quad (p pm ) d) b) 8000 10000 400Hz 1800 8200 40k 140k Figure 4.6: Predicted data mis t for the inversions on the lines shown in Fig. 4.5: a)-b) Line 10300, c)-d) Line 20590. The curves are o set by 1000 ppm for clarity. Solid lines: predicted response. Crosses: observed data. Observed data are from [21]. 540 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. 55No rt hi n 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 rt hi n 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 Figur e 4.8 : 2D (a),c) ) an d 1D (b), d) ) in verte d lo g conductivi ty m odel s fo r th e tw o sur ve y line s high - lig hte d in Fig . 4.5 . a)-b ) Lin e 10300 , c)-d ) Lin e 20590 . 56In 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 uctuations 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 uctuations 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 t to a large portion of the observations and they reproduced the uctuations caused by the ight 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 e ect 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- ow 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 57a more conclusive answer. To summarize the test results, our implementation of the 2.5D EM for- ward modeling and inversion scheme were successful when veri ed against more established 1D counterparts. The forward modeling was able to com- pute the magnetic response to a level of accuracy ( 1%) that is su ciently 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. 58Chapter 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 eld 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 de nite 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 re ned 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 e ective 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. 59Therefore 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. 60Bibliography [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 e. 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 nite volume methods for distributed non-smooth parameter identi cation. 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 nite 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 elds 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. Green eld. 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 nite 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
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Inversion of airborne electromagnetic data in 2.5D
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Inversion of airborne electromagnetic data in 2.5D Yu, Wing Wa 2012-10-16
pdf
Page Metadata
Item Metadata
Title | Inversion of airborne electromagnetic data in 2.5D |
Creator |
Yu, Wing Wa |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | In this work, we implement an inversion algorithm for airborne electromagnetic (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. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-10-16 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0073275 |
URI | http://hdl.handle.net/2429/43416 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2012_fall_yu_wing-wa.pdf [ 4.06MB ]
- Metadata
- JSON: 24-1.0073275.json
- JSON-LD: 24-1.0073275-ld.json
- RDF/XML (Pretty): 24-1.0073275-rdf.xml
- RDF/JSON: 24-1.0073275-rdf.json
- Turtle: 24-1.0073275-turtle.txt
- N-Triples: 24-1.0073275-rdf-ntriples.txt
- Original Record: 24-1.0073275-source.json
- Full Text
- 24-1.0073275-fulltext.txt
- Citation
- 24-1.0073275.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.24.1-0073275/manifest