A Cooperative Magnetic Inversion Method With Lp-normRegularizationbyDominique FournierBSc Honours Geophysics, The University of British Columbia, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Geophysics)The University of British Columbia(Vancouver)October 2015c© Dominique Fournier, 2015AbstractIn this master’s thesis, I implement a Cooperative Magnetic Inversion (CMI) al-gorithm for the 3-D modeling of magnetic rocks at depth. While in most casesit is assumed that the magnetic response is purely induced, certain rocks have theability to retain a permanent magnetic moment in any orientations, also known asremanent magnetization. The effect of remanence has long been recognized as anobstacle for the geological interpretation and modeling of magnetic data. My ob-jective is to improve current magnetic inversion methods to recover simpler andbetter defined magnetization models.The CMI algorithm brings together three inversion techniques previously intro-duced in the literature. First, magnetic data are inverted for an equivalent-sourcelayer, which is used to derive magnetic amplitude data. Next, amplitude data areinverted for an effective susceptibility model, providing information about the ge-ometry and distribution of magnetized objects. Finally, the effective susceptibilitymodel is used to constrain the Magnetic Vector Inversion (MVI), recovering theorientation and magnitude of magnetization. All three algorithms are formulatedas regularized least-squares problems solved by the Gauss-Newton method.I further constrain the solution by imposing sparsity constraints on the modeland model gradients via an approximated lp-norm penalty function. I elaborate aScaled Iterative Re-weighted Least-Squares (S-IRLS) method, allowing for a stableand robust convergence of the algorithm while combining different lp-norms on therange 0 ≤ p ≤ 2. The goal is to reduce the complexity of magnetization modelswhile also imposing geometrical constraints on the solution.As a final test, I implement the CMI algorithm on an airborne magnetic surveyover the Ekaty Property, Northwest Territories. I formulate a tiled inversion schemeiiin order to reduce the computational cost and increase the level of parallelization.The final merged magnetization model provides insights into the distribution ofdyke swarms and kimberlite pipes. Following the regional inversion, I focus myanalysis on sixteen known kimberlite pipes. Magnetization vector inclinations arecompared to the expected polarity of rocks inferred from radiometric dating.iiiPrefaceAll the work in this thesis is my own. The algorithms and ideas used from othersources are duly cited.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Magnetic methods . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Magnetostatic Problem . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Magnetization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.1 Regularized inversion . . . . . . . . . . . . . . . . . . . . . . . . 163.1.1 Iterative solver . . . . . . . . . . . . . . . . . . . . . . . 203.1.2 Preconditioner . . . . . . . . . . . . . . . . . . . . . . . 223.1.3 Bound constraints . . . . . . . . . . . . . . . . . . . . . 243.2 Magnetic susceptibility inversion . . . . . . . . . . . . . . . . . . 253.2.1 Synthetic example . . . . . . . . . . . . . . . . . . . . . 26v3.2.2 Effect of remanent magnetization . . . . . . . . . . . . . 313.3 Magnetic vector inversion . . . . . . . . . . . . . . . . . . . . . . 343.4 Magnetic amplitude inversion . . . . . . . . . . . . . . . . . . . . 383.5 Equivalent source method . . . . . . . . . . . . . . . . . . . . . . 433.5.1 Comment for future research . . . . . . . . . . . . . . . . 474 Mixed Lp-norm Regularization . . . . . . . . . . . . . . . . . . . . . 494.1 Least-norm and regularized least-squares . . . . . . . . . . . . . 494.2 Lp-norm and iterative re-weighted least squares . . . . . . . . . . 534.3 IRLS solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.3.1 Basic IRLS algorithm . . . . . . . . . . . . . . . . . . . 584.3.2 1-D synthetic example . . . . . . . . . . . . . . . . . . . 604.3.3 Regularization scaling . . . . . . . . . . . . . . . . . . . 654.3.4 Threshold parameter ε . . . . . . . . . . . . . . . . . . . 684.4 Scaled-IRLS method (S-IRLS) . . . . . . . . . . . . . . . . . . . 744.4.1 Cell-based weights (wr, wm) . . . . . . . . . . . . . . . . 794.5 Mixed lp-norm regularization . . . . . . . . . . . . . . . . . . . . 824.5.1 Localized S-IRLS . . . . . . . . . . . . . . . . . . . . . . 824.5.2 2-D example . . . . . . . . . . . . . . . . . . . . . . . . 874.5.3 3-D example . . . . . . . . . . . . . . . . . . . . . . . . 924.6 Case study - Tli Kwi Cho kimberlite complex . . . . . . . . . . . 954.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985 Cooperative Magnetic Inversion (CMI) . . . . . . . . . . . . . . . . 1015.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.2 CMI with S-IRLS regularization . . . . . . . . . . . . . . . . . . 1075.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1076 Case Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1106.1 Ekati Property . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1116.1.1 Geology . . . . . . . . . . . . . . . . . . . . . . . . . . . 1116.1.2 Airborne magnetic data . . . . . . . . . . . . . . . . . . . 1136.2 Regional scale inversion . . . . . . . . . . . . . . . . . . . . . . 1146.3 Deposit scale inversion . . . . . . . . . . . . . . . . . . . . . . . 118vi6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1227 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1277.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130viiList of TablesTable 1.1 Magnetic susceptibility for various rock types. . . . . . . . . . 3Table 3.1 Conjugate Gradient algorithm . . . . . . . . . . . . . . . . . . 22Table 3.2 Line-search . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Table 3.3 Inversion parameters. . . . . . . . . . . . . . . . . . . . . . . 28Table 4.1 IRLS Solver 1: Fix parameters . . . . . . . . . . . . . . . . . 59Table 4.2 IRLS Solver 2: β -search . . . . . . . . . . . . . . . . . . . . . 65Table 4.3 IRLS Solver 3: Scaled regularization . . . . . . . . . . . . . . 68Table 4.4 IRLS Solver 4: ε-Cooling . . . . . . . . . . . . . . . . . . . . 70Table 6.1 Dyke swarms of the Lac de Gras region listed in increasing ageas published by Buchan et al. (2009) . . . . . . . . . . . . . . 113Table 6.2 Regional-scale inversion parameters. . . . . . . . . . . . . . . 116Table 6.3 Deposit-scale inversion parameters. . . . . . . . . . . . . . . . 121Table 6.4 Average magnetization amplitude and direction from the localinversion over 16 known kimberlite pipes. Uncertainties werecalculated from a three-cell cube standard deviation around theapproximate location of each kimberlite pipe. . . . . . . . . . 121Table 6.5 Summary of published results over selected kimberlite pipes . . 125Table 6.6 Radiometric age and inverted magnetization direction for 11 pipes.126viiiList of FiguresFigure 1.1 (a) Dipolar magnetic field of the Earth in its normal orientation.(b) Polarity of the field inferred from the geological record overthe past 90 Ma (Cande and Kent, 1995). . . . . . . . . . . . . 4Figure 1.2 Magnetic field data~b observed on a plane over a magnetizedsphere located at the origin. The orientation and magnitude ofmagnetization of the sphere is function of the magnetic sus-ceptibility κ , the inducing field ~B0 and remanent components~MNRM, neglecting self-demagnetization effects. . . . . . . . . 6Figure 3.1 (a) Synthetic susceptibility model consisting of a folded anomaly(κ = 0.075 SI) arching around a discrete block (κ = 0.05 SI) . 27Figure 3.2 (a)Data generated from the synthetic susceptibility model sub-ject to a vertical 50,000 nT inducing field. (b) Data are thencorrupted with (c) random Gaussian noise, 1 nT standard devi-ation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 3.3 Convergence curve showing the data misfit φ (k)d and modelnorm φ (k)m as a function of β iterations. The inversion achievestarget misfit after the 7th iteration, which in this case also corre-sponds to the point of maximum curvature on the misfit curve.Attempting to further lower the data residual comes at the riskof fitting some of the Gaussian noise. . . . . . . . . . . . . . 28ixFigure 3.4 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered susceptibility model for a purely induced response. Themodel is smooth but recovers the arc and block anomaly atroughly the right depth. . . . . . . . . . . . . . . . . . . . . . 30Figure 3.5 Comparison between (a) observed and (b) predicted data fromthe recovered susceptibility model. (c) The normalized dataresiduals appear to be correlated with the location of the mag-netic body. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.6 Perspective view and sections through the synthetic magneti-zation model. The arc-shaped anomaly is magnetized at 45◦from horizontal and with variable declinations directions be-tween [−45◦N ; 45◦N]. . . . . . . . . . . . . . . . . . . . . 32Figure 3.7 (a) Data generated from the synthetic magnetization model.(b) Observed data are corrupted with (c) random Gaussian noise,1 nT standard deviation. . . . . . . . . . . . . . . . . . . . . 32Figure 3.8 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered susceptibility model assuming no remanence. The arc-shaped anomaly is poorly recovered and magnetic susceptibil-ities are pushed at depth and outwards. . . . . . . . . . . . . . 33Figure 3.9 Comparison between (a) observed and (b) predicted data fromthe recovered susceptibility model assuming a purely inducedresponse. (c) The inversion has a harder time fitting the largenegative fields along the arc. . . . . . . . . . . . . . . . . . . 33Figure 3.10 (a) Iso-surface (κe =0.001) and (b) sections through the recov-ered magnetization model from the MVI method. The inver-sion recovers the true orientation of magnetization inside theblock, but the thin arc is poorly resolved. . . . . . . . . . . . 37Figure 3.11 Comparison between (a) observed and (b) predicted data fromthe recovered magnetization model from the MVI method. Themodel can replicate the data at the same level achieved by thepurely induced problem. . . . . . . . . . . . . . . . . . . . . 37xFigure 3.12 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered effective susceptibility model. The arc-shaped and blockanomalies are recovered at the right location, but smoothlystretched vertically. . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.13 Comparison between (a) observed and (b) predicted data fromthe recovered effective susceptibility model. The inversion canpredict most of the data within one standard deviation. . . . . 42Figure 3.14 (a) Synthetic model consisting of 200 unit cubes of suscepti-ble material in a non-susceptible background. Data are gener-ated on a plane one unit above the source location, assuming apurely vertical inducing field. Various components of the fieldsare shown in figure (b) to (f). . . . . . . . . . . . . . . . . . 45Figure 3.15 (a) Recovered equivalent source layer from TMI data using apositivity constraint. The residuals between observed and pre-dicted data are shown in figure (b) to (f) for various compo-nents of the field. Each component is well recovered withinthe noise level. . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 3.16 (a) Recovered equivalent source layer from TMI data after re-moving a portion of data over the corner of the magnetic anomaly.The residuals between observed and predicted data are shownin figure (b) to (f) for various components of the field. Notethe large correlated artifacts recovered on the bx, by and |b|components. . . . . . . . . . . . . . . . . . . . . . . . . . . 48xiFigure 4.1 Comparative contour maps for various objective functions overa range of model values [m1; m2]. (a) The mininum of the mis-fit function φd forms a line spanned byN (F) (red dash). Theleast−norm solution is marked as a solid dot. (Middle) Reg-ularization functions and gradient directions (blue arrows) forapproximated lp-norm measures of the model for (b) p = 2,(d) p= 1 and (f) p= 0. The gradient directions are shown fortwo different starting models (black arrows). (bottom) Contourmaps of the initial objective functions φ(m) = φd+φm for thesame two starting models: m(0)1 (solid) and m(0)2 (dash). (c) Inthe l2-norm case , the function has a global minimum regard-less of the starting model, while for non-linear functions for(e) p = 1 and (g) p = 0, the objective function changes withrespect to the starting model. . . . . . . . . . . . . . . . . . . 52Figure 4.2 Contour maps for various objective functions after convergenceof the IRLS algorithm. (a) Final model obtained with the l2-norm penalty on the model for two starting models at m(0)1 =[0.2;0.4] and m(0)2 = [0.6;0.2] for a fixed trade-off parameter(β = 1e− 4). In both cases, the solution converges to theglobal minimum, which is also the least − norm solution atmln = [0.2;0.4]. (b) Solution with the l1-norm penalty for thesame starting models and trade-off parameter, converging to aglobal minimum at m∗ = [0;0.5]. This solution is sparse andcan reproduce the data. (c) The same experiment is repeatedfor the l0-norm penalty, converging to two different solutionsdepending on the relative magnitude of the starting model pa-rameters. Both solutions are sparse and honor the data. . . . . 53xiiFigure 4.3 (a) Penalty function ρ(x) for different approximate lp-normmeasures, and enlarged view near the region of influence of ε ,making the lp-norm continuous at the origin. (b) IRLS weightsr(x) as a function of model function x(m) and p-values and fora fix threshold parameter (ε = 1e− 2). (c) Gradients g(x) ofthe model penalty function for various p-values. Note that thegradients are on a logarithmic scale due to the rapid increaseas xi→ ε for p< 1. . . . . . . . . . . . . . . . . . . . . . . . 57Figure 4.4 (a) Synthetic 1D model made up of a rectangular pulse and aGaussian function. (b) Kernel functions consisting of expo-nentially decaying cosin functions of the form f j(x) = e− j x ·cos(2pi jx). (c) Data generated from d= F m , with 5% ran-dom Gaussian noise added. . . . . . . . . . . . . . . . . . . 61Figure 4.5 (a) Recovered models from the smooth l2-norm regularization,and (bottom) measure of model misfit and model norm as afunction of iterations. Both the rectangular pulse and Gaus-sian function are recovered at the right location along the x-axis, but the solution is smooth and dispersed over the entiremodel domain. (b) Solution obtained from the IRLS with sp-asity constraints on the model gradients (q=0), using the l2-norm model (a) to initiate the IRLS steps. The algorithm usesa fixed threshold parameter (ε = 1e−8) and fixed trade-off pa-rameter β . The final solution is blocky, as expected from thenorm chosen, but fails to achieve the target data misfit. Clearlythe influence of the regularization function has overtaken theminimization process. . . . . . . . . . . . . . . . . . . . . . 64xiiiFigure 4.6 (Top) Recovered models from two different algorithms used toimplement the IRLS method for q = 0 and a fix threshold pa-rameter (ε = 1e− 8). (Bottom) Measure of model misfit andmodel norm as a function of iterations. (a) The first algorithmsearches for an optimal trade-off parameter β (k) between eachIRLS step, requiring a solution to multiple sub-inversions. (b)The second algorithm only adjusts β (k) once after each IRLSstep. A new scaling parameter γ(k) is added to smooth the tran-sition between the IRLS updates. The algorithm recovers asimilar blocky model but is computationally cheaper, as indi-cated by the total number of beta iterations and CG solves. . . 67Figure 4.7 (Top) Recovered models from two different algorithms usedto implement the IRLS method for q = 0 with cooling of thethreshold parameter ε . (Bottom) Measure of model misfit andmodel norm as a function of iterations. Both algorithms adjustβ and scaling parameter γ(k) between each IRLS iteration. (a)In the first case, the threshold parameter ε is monotonicallyreduced until reaching the convergence criteria. The solutionis sparser than the previous algorithm, even though ε is muchlarger than machine error. (b) In the second case, ε is decreaseduntil reaching the target threshold ε∗. The solution is blocky,penalizing small model gradients. . . . . . . . . . . . . . . . 71Figure 4.8 (a) Distribution of model parameters and (b) model gradientsrecovered from the smooth l2-norm regularization. Both curvesshow a sharp corner around which the model functions varyrapidly. Similarly, a measure of (c) model norm sMS and (d)model gradient norm sMGS can be computed over a range of εvalues, yielding a similar L− curve. The point of maximumcurvature can be used to determine the optimal e f f ective zeroparameter εp and εq. . . . . . . . . . . . . . . . . . . . . . . 73xivFigure 4.9 (a) (Top) Recovered model using a mixed-norm regularizationfor p= 0, q= 2, ε = 1e−3. (Bottom) Measure of model mis-fit and model norm as a function of iterations. The inversionconverges to a sparse solution without apparent penalty on thegradient, indicative of an imbalance in the regularization. (b)Recovered model and convergence curves after rescaling of theregularization, yielding a model that is both sparse and smoothas expected from the applied mixed-norm . . . . . . . . . . . 75Figure 4.10 Contour maps for (a) the misfit function φd , (b) the model normφs and (c) the norm of model gradients φx. (d) The total ob-jective function φ(m) has a global minimum located at m =(0.5,1.0) for a given small trade-off parameter (β = 1e− 3).The direction of update is shown for two starting models m(0)(black and white dot). (e) Section through the objective func-tion along the minimum of φd . The global minimum occurswhere the partial gradients of ∂φs∂ m1 and∂φx∂ m1 have equal and op-posite signs. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Figure 4.11 (a) Partial gradients of approximated lp-norm penalties for afix stabilizing parameter ε = 1e− 2. Gradients for p < 1 areconsistently larger on the interval [0 < xi <√1− ε2], makingit hard to combine multiple norm penalties within the same ob-jective function. (b) Function gradients after applying a scaleof ε(1−p/2), forcing each lp-norm to intersect at m=√ε . . . . 78Figure 4.12 Recovered models for two different depth weighting formula-tions: (red) weighted sensitivity φ(mˆ), (black) weighted regu-larization φ(m). (a) True and recovered models using the φ(mˆ)and φ(m) formulations for penalty applied on the model gradi-ents for q= 0 and (b) for p= 0, q= 2. The weighted sensitiv-ity formulation φ(mˆ) increases the influence the regularizationfunction with distance along the x-axis, skewing the model to-wards the right. . . . . . . . . . . . . . . . . . . . . . . . . . 81xvFigure 4.13 (a) Model error ‖m−m∗‖1 and (b) misfit function for the 441inverted models using a range of regularization with mixed-norm penalty on the model for 0≤ p≤ 2 and on model gradi-ents for 0 ≤ q ≤ 2. (c) The largest model error (‖δm‖1) wasobtained with the mixed-norm for p = 0, q = 2, compared to(d) the optimal solution found with p= 1.5 and q= 0.4. . . . 83Figure 4.14 (a) Nine of the 441 inverted models for a range of mixed-normpenalties on the model and its gradient for 0 ≤ p ≤ 2 and 0 ≤q≤ 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 4.15 (Left) Improved solution for the 1-D problem after applying alocalized mixed-norm penalty, where the regularization is di-vided into two regions with independent lp-norm regulariza-tion: (left) p= q= 0, (right) p= 1,q= 2. (Right) Convergencecurves for the mixed-norm S-IRLS inversion. . . . . . . . . . 86Figure 4.16 (a) Synthetic 2-D model made up of a square block and asmooth Gaussian function. (b) Example of a kernel functionfor e−ω r · cos(2piωr) and (c) data generated from d= F m.Five percent random Gaussian noise is added. . . . . . . . . 88Figure 4.17 (a) Recovered model for p = qx = qz = 1 penalizing finitedifference gradients in orthogonal directions, yielding right-angled anomalies. (b) Recovered model for the same normsbut penalizing the absolute gradient of the model (|∇m|) re-covering round edges. . . . . . . . . . . . . . . . . . . . . . 90Figure 4.18 Distribution of lp-norm on the model and model gradients overthe 2-D model domain. The original boundary of each regionwas smoothed in order to get a slow transition and reduce vis-ible artifacts. Regions were chosen to cover a larger area thanthe anomalies to simulate a blind-inversion. . . . . . . . . . . 90xviFigure 4.19 (a) Smooth l2-norm solution used to initiate the IRLS itera-tions. (b) Recovered model using the mixed-norm regular-ization after seven S-IRLS iterations. The contour line (red)marks the value of εp, judged to be the effective zero valueof the model (mi ≤ 5e-2). Both models (a) and (b) fit thedata within 2% of the target misfit φ ∗d . (c) Dual plots show-ing the distribution of model parameters and the gradient ofthe l0-norm penalty function gˆp(m) as a function of S-IRLSiterations. High penalties are applied to progressively smallermodel values. The final model nicely recovers both the blockyand smooth Gaussian anomaly. . . . . . . . . . . . . . . . . . 91Figure 4.20 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered susceptibility model after five IRLS iterations for (p =0, q = 2). The final model is substantially closer to the truesolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Figure 4.21 Comparison between (a) observed and (b) predicted data fromthe recovered susceptibility model using compact norms for(p = 0, q = 2). (c) Normalized data residuals are within twostandard deviations. . . . . . . . . . . . . . . . . . . . . . . . 93Figure 4.22 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered susceptibility model after nine IRLS iterations for (p =0, q= 1) . Sparsity constraints on the model and model gradi-ents yield a simple and blocky model. . . . . . . . . . . . . . 94Figure 4.23 Comparison between (a) observed and (b) predicted data fromthe recovered susceptibility model using compact norms for(p = 0, q = 1). (c) Normalized data residuals are within twostandard deviations. . . . . . . . . . . . . . . . . . . . . . . . 94xviiFigure 4.24 (a) (Left) Horizontal section through the recovered susceptibil-ity model at 25 m depth below topography from the smooth l2-norm regularization. (Right) Iso-surface of susceptibility val-ues around 0.002 SI. (b) Recovered model using the mixed-norm S-IRLS algorithm. Magnetic dykes are better recovered,imaged as continuous plates and extending vertically at depth.Susceptibility values for DO-27 and DO-18 have increased,showing as compact vertical pipes. . . . . . . . . . . . . . . . 97Figure 4.25 Horizontal section through the mixed-norm models applied tofour sub-regions with smooth transition across zones. . . . . . 99Figure 4.26 (a) Observed and predicted data over the TKC kimberlite com-plex. (b) Residuals between observed and predicted data nor-malized by the estimated uncertainties (10 nT). Both the smoothand mixed-norm inversions reproduce the data within four stan-dard deviations. . . . . . . . . . . . . . . . . . . . . . . . . . 100Figure 5.1 Schematic representation of the Cooperative Magnetic Inver-sion (CMI) algorithm. Input and output parameters are indi-cated by a dash arrow. . . . . . . . . . . . . . . . . . . . . . 103Figure 5.2 (a) Inverted equivalent-source layer and (b-f) predicted xˆ, yˆ, zˆ-component, TMI and magnetic amplitude data for the syntheticmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104Figure 5.3 (a) Iso-surface (0.002 SI) and (b) sections through the recov-ered effective susceptibility model. This effective susceptibil-ity model is used to construct a weighting matrix to constrainthe MVI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105Figure 5.4 Comparison between (a) observed and (b) predicted data fromthe recovered effective susceptibility model. The inversion canpredict most of the data within one standard deviation. . . . . 105Figure 5.5 (a) Iso-surface (0.005 SI) and (b) sections through the recov-ered magnetization model from the CMI algorithm (p= 2, q=2) . The inversion recovers both the arc and the block anomalyas distinct objects. . . . . . . . . . . . . . . . . . . . . . . . 106xviiiFigure 5.6 Comparison between (a) observed and (b) predicted data fromthe recovered susceptibility model. (c) Normalized data resid-uals are within two standard deviations. . . . . . . . . . . . . 106Figure 5.7 (a) Iso-surface (0.01 SI) and (b) sections through the recov-ered effective susceptibility model from the amplitude inver-sion with sparsity constraint applied (p = 0, q = 1). The lp-norm constraint considerably reduces the complexity of themodel, although the model is still stretched vertically. . . . . . 108Figure 5.8 Comparison between (a) observed and (b) predicted amplitudedata from the recovered compact effective susceptibility model(p = 0, q = 1). (c) Normalized data residuals are within twostandard deviations. . . . . . . . . . . . . . . . . . . . . . . . 108Figure 5.9 (a) Iso-surface (0.01 SI) and (b) sections through the recov-ered magnetization model from the CMI algorithm. Compactnorms (p= 0, q= 2) were applied during the amplitude inver-sion. The lp-norm constraint considerably reduces the com-plexity of the model. . . . . . . . . . . . . . . . . . . . . . . 109Figure 5.10 Comparison between (a) observed and (b) predicted data fromthe recovered magnetization CMI model. (c) Normalized dataresiduals are within two standard deviations. Correlated resid-uals are no longer seen. . . . . . . . . . . . . . . . . . . . . . 109Figure 6.1 Topography and known kimberlites over the Lac de Gras re-gion. The location of the main operations for the Ekati andDiavik mines are shown (red dot). . . . . . . . . . . . . . . . 112Figure 6.2 Topography, aeromagnetic survey and known kimberlites overthe Lac de Gras region. The location of the main operationsfor the Ekati and Diavik mines is marked with a red dot. . . . 114Figure 6.3 Regional data (a) before and (c) after removal of a regionaltrend. The median value within selected regions (boxes) wereused to compute (b) a first-order polynomial trend. . . . . . . 115Figure 6.4 TMA data after removal of the regional signal and tiling con-figuration used for the regional inversion. . . . . . . . . . . . 117xixFigure 6.5 Interpreted dykes (line) from the property-scale magnetizationinversion and known kimberlite pipes location (circle). The 16pipes chosen for a deposit-scale inversion are labeled. . . . . . 118Figure 6.6 Comparison between (a) observed and (b) merged predicteddata from the recovered magnetization model over the EkatiProperty. (c) Normalized residual data show horizontal stria-tions due to variations in magnetic data between adjacent lines,which could not be accounted for by the inversion. Most of thepredicted data from the individual tiles can fit the observeddata within one standard deviation. (d) Forward modeled datafrom the merged magnetization model and (c) normalized dataresidual. It appears that a large portion of the the low frequencycontent has been lost during the merging step. More researchis required in order to preserve the long wavelength information.119Figure 6.7 (a) Local data collected over the Misery pipe showing a localwestern trend. (b) Regional field data are computed from a lo-cal inversion. Most of the signal comes from a dyke runningnorth-south along the western edge of the local tile. (c) Resid-ual data after regional field removal, showing a clear reverselymagnetized anomalie corresponding with the location of theMisery pipe. . . . . . . . . . . . . . . . . . . . . . . . . . . 122Figure 6.8 Horizontal sections through the recovered (a) induced and (b)remanent magnetization model. Several pipes with strong re-manence can easily be identified as discrete circular anomalies. 124Figure 6.9 Age of 11 pipes from the Ekati region with respect to Earth’spolarity reversal. . . . . . . . . . . . . . . . . . . . . . . . . 126xxAcknowledgmentsFirst of all I need to thank my supervisor Dr. Doug Oldenburg for his mentorshipand his contagious passion for research. It has been an honor and privilege to learnfrom a pillar in Earth sciences. Special acknowledgment to Dr. Kristofer Davis forhis support, beginning four years ago with my undergraduate honors thesis.Thanks to all the UBC-GIF members for the interesting conversations andfriendship. Sincere gratitude for the collaboration and mentorship that I have re-ceived over the years from several senior geoscientists: Peter Kowalczyk, NigelPhillips and all the staff at Mira Geoscience, Dr. Peter Fullagar at Fullagar Geo-physics, Dr. Randy Enkin at the GSC and Ste´phane Lorrain at EnvironnementIllimite´. Special thanks to Jon Carlson at Dominion Diamond Ekati Corporation,for granting me access to the rock physical property database. Also to Jennifer Pelland Brooke Clements at Peregrine Diamonds Ltd for providing the TKC dataset.Last but not least, extra special thanks goes to my editor-in-chief, part-timelover and full-time best friend Chantelle Bellrichard.xxiChapter 1IntroductionMagnetic methods have a long and rich history in Earth sciences (Nabighian et al.,2005). As early as the 19th century it was understood that particular rocks werethe source of local magnetic anomalies measurable at the surface. Since the devel-opment of the fluxgate magnetometer during World War II, magnetic surveys havecontributed greatly to what we know about the earth. At a global scale, marinesurveys over magnetically polarized oceanic plates have shaped our understandingof continental drift and plate tectonics (Vine and Matthews, 1963). By the 1990s,Global Positioning System (GPS) and increasingly sensitive instruments led to aproliferation of airborne surveys that gave us the ability to map geological featuresremotely. Later advancements in numerical computing contributed to the adop-tion of 3-D magnetic inversion as a tool to image magnetic bodies at depth. Thismaster’s thesis focuses on improving current inversion methods using a particularstrategy. I explore how various inversion codes can be used cooperatively, yieldinga more robust and flexible imaging tool in exploration geophysics.But first, a review of the magnetic properties of rocks. Assuming no free cur-rents and steady state, Maxwell’s equations can be written as:∇ ·~B= 0 (1.1)∇× ~H = 0 (1.2)~B= µ~H , (1.3)1where ~B is the magnetic flux density in Tesla (T) and ~H is the magnetic field inunits of amperes per meter (A/m). In matter, the magnetic permeability µ relatesthe magnetic flux density ~B to the field ~H such thatµ = µ0(1+κ) , (1.4)where κ is the magnetic susceptibility, a dimensionless positive number describingthe ability of certain material to become magnetized under an applied field. In freespace, the term magnetic field is often used interchangeably to describe ~B and ~Hsince linearly related by the magnetic permeability of free space µ0 (4pi×10−7 T·m/A).Factors responsible for rock magnetism can be divided into an induced and aremanent component such that:~M = κ~H+ ~MNRM (1.5)~H =1µ0(~B0+~BA) , (1.6)where ~M (A/m) is the magnetization per unit volume of a rock, the quantity ofinterest in mineral exploration. The inducing field can be further decomposed in theprimary geomagnetic flux ~B0 and anomalous local flux ~BA. The Natural RemanentMagnetization (~MNRM) describes the ability of matter to retain a net magnetizationcomponent in the absence of an inducing field.Generated within the Earth’s core (Campbell, 1997), the inducing field ~B0varies between 30,000 nT at the equator, to over 50,000 nT near the poles as il-lustrated in Figure 1.1. From the geochronological record, we also know that theEarth’s field went through at least 92 magnetic field reversals over the last 100 mil-lion years. The current orientation of the field, pointing down towards the Northpole, is defined as the normal polarity. In most cases, the induced componentfrom the primary field ~B0 is the main driver for the magnetic response of rocks.Secondary fields ~BA are generally much weaker than the primary field and arisefrom the interaction of neighbouring magnetic objects. Secondary fields may beresponsible for sel f − demagnetization effects as studied by Lelie`vre (2003). Inthe majority of cases, the strength of ~B0 largely dominates any other secondary2Table 1.1: Magnetic susceptibility for various rock types.Rock Type Magnetic Susceptibility κ×10−6SIGranite 20 - 40,000Slates 0 - 1,200Basalt 500 - 80,000Oceanic basalts 300 - 36,000Limestone (with magnetite) 10 - 25,000Gneiss 0 - 3,000Sandstone 35 - 950Hematite (ore) 420 - 10,000Magnetite (ore) 7×104 - 1.4×107Magnetite (crystal) 1.5×108fields, forcing the magnetization to be aligned parallel to the geomagnetic field. Asa rule of thumb, secondary fields only become important for rocks with magneticsusceptibility κ > 1. Table 1.1 summarizes the magnetic susceptibility of majorrock types, varying over several orders of magnitude. Magnetite is by far the mostabundant magnetic mineral, followed by pyrrhotite and ilmenite (Clark, 1991). Allthose minerals are commonly associated with iron oxides, and potentially goodindicators for mineral deposits.Natural Remanent Magnetization ~MNRM is a permanent magnetization momentpreserved in the absence of an inducing field. The reader is encouraged to refer toBlakely (1996) and Clark (1991) for a more in-depth explanation of chemical, ther-mal and biological processes responsible for the remanent component of variousrock types. While the induced component is linearly related to the inducing field,nothing can be assumed about the NRM orientation and it remains challenging toestimate. If aligned with the inducing field direction, the NRM is indistinguishablefrom a purely induced response, resulting in an over-estimation of the magneticsusceptibility of rocks. For large NRM components, perpendicular or anti-parallelto the inducing field, the magnetic response of compact objects may get distorted,potentially resulting in false interpretation about the distribution and geometry ofrock units. Direct geological interpretation of magnetic data that ignores the effect3Figure 1.1: (a) Dipolar magnetic field of the Earth in its normal orientation.(b) Polarity of the field inferred from the geological record over the past 90 Ma(Cande and Kent, 1995).of remanence has long been recognized as problem in mineral exploration. Thisbrings up an important question: How can we better recover the location and geom-etry of magnetized objects without knowledge of the total magnetization direction~M?From Gauss’s law, the relation between the observed magnetic field and rockmagnetization is expressed as:~b(r) =µ04pi∫V∇∇1r· ~M dV , (1.7)where~b is the magnetic field (T) as measured at some distance r from a magneticanomaly with magnetization per unit volume ~M (A/m). Since~b is usually small, itis commonly measured in units of nano-Tesla (nT). The majority of magnetic data4consist of Total Magnetic Intensity (TMI) measurements which can be written as:bTMI = |~B0+~BA| , (1.8)where we measure the magnitude of the field rather than the individual compo-nents. In most cases we are only interested in the anomalous local fields. Underthe assumption that |~BA| |~B0|, the anomalous field is approximated to be paral-lel with the direction of the inducing field Bˆ0. The Total Magnetic field Anomaly(TMA) is given by:bTMA = |~B0+~BA|− |~B0|' ~BA · Bˆ0 .(1.9)Figure 1.2 gives an example of TMA data measured on a plane above a magne-tized sphere. The main challenge in interpreting magnetic data is to characterizethe magnetic sources when nothing is known about the location and magnetic prop-erties (κ, ~M) of buried objects. Accurately imaging the location and geometry ofmagnetic objects is a core problem in exploration geophysics.Early geophysical studies relied primarily on filtering techniques to infer ge-ological structures and identify potential targets (Cowan et al., 2000). In order toreduce the complexity of the problem, the magnetization direction is in most casesassumed to be purely induced along ~B0, neglecting self-demagnetization and rema-nence effects. Inversion codes, such as the UBC-MAG3D from Li and Oldenburg(1996), also rely on this assumption. While the induced component may domi-nate in most cases, recent petrophysical studies seem to indicate that the effect ofremanent magnetization may be more important than previously thought (Enkin,2014). This is especially true for specific types of mineral deposits such as BandedIron-Formations (BIF) and diamondiferous kimberlite (Dransfield et al., 2003; Liet al., 2010) Strong remanent magnetic components can hinder the interpretationof magnetic data thus leading to false drilling targets.5Figure 1.2: Magnetic field data~b observed on a plane over a magnetized spherelocated at the origin. The orientation and magnitude of magnetization of thesphere is function of the magnetic susceptibility κ , the inducing field ~B0 andremanent components ~MNRM, neglecting self-demagnetization effects.1.1 Magnetic methodsSeveral studies have been dedicated to the inversion of magnetic data in the pres-ence of remanence. As summarized by Li (2012), the proposed inverse methodscan be divided into three categories. The first category estimates the magnetizationdirection from the data in a pre-processing step. Strategies such as the Helbig’smethod (Phillips, 2003) and the cross-correlation (Dannemiller and Li, 2006) areused upstream of standard inversion codes. Adapted from the 2-D analytical signalof Nabighian (1972), Roest et al. (1992) estimate the orientation of remanence in3-D. In most cases, data are transformed to the wavenumber domain in order to ex-tract magnetic field components. These methods are simple and robust for simpleand isolated anomalies, but become impractical when the data are acquired overrough terrain or complicated geology.The second category deals with magnetic data that are weakly dependent onthe orientation of magnetization. It has been shown that magnetic amplitude data6are weakly dependent on the magnetization direction in 3-D (Nabighian, 1972).Magnetic amplitude data are simply calculated as:|~b|= (b2x+b2x+b2x)1/2 , (1.10)where bx, by and bz are the three components of the anomalous magnetic field.Inverting magnetic amplitude data can provide a robust estimate for the locationof magnetized bodies as shown by Shearer (2005). Magnetic amplitude data aresuccessfully inverted over a kimberlite deposit, recovering the true location ofkimberlite pipes and associated magnetic dykes. Because most magnetic surveysonly record TMI data, the three components of the field must be derived in post-processing. Two strategies have been proposed in the literature to extract magneticfield components. The first relies on the fact that magnetic data are potential fielddata, hence knowledge of any component of the field on an infinite plane abovethe source is sufficient to calculate the remaining components. Data are interpo-lated onto a uniform grid, which is then used to calculate the components of thefield using Fourier transform techniques. This may be a problem however for air-borne surveys over steep topography or acquired at low latitude. Alternatively, theequivalent-source method has been proposed by Li and Oldenburg (2010) to alle-viate those constraints. Taking advantage of the non-uniqueness of potential fielddata, TMA data are inverted for a layer of magnetic sources then used to forwardmodel three-component magnetic data.The third method directly solves for the orientation of magnetization withoutmaking any assumptions about the location or geometry of causative bodies. In theMagnetic Vector Inversion (MVI) proposed by Lelie`vre and Oldenburg (2009), themagnetization vector is decomposed in its induced and orthogonal components.The MVI is closely related to the method proposed by Kubota and Uchiyama(2005), and later borrowed by Ellis et al. (2012). This method results in a largeunder-determined inverse problem, with roughly three times the number of free-variables over conventional susceptibility inversion codes. Recovered magnetiza-tion models are generally smooth and become overly complicated for direct geo-logical interpretation. As pointed out by Lelie`vre (2009), any a priori informationfrom surface or borehole measurements may greatly reduce the non-uniqueness of7the problem. Unfortunately this kind of information is rarely available in greenfieldsettings or is only available in sparse samples.Building upon the work of Shearer (2005), Liu et al. (2015) propose a hybridthree-step process to invert for the location, orientation and magnetic susceptibilitydistribution in 2-D. In the first step, the location of magnetic material is invertedfrom amplitude data. Next, the orientation of magnetization is found via a corre-lation method. The orientation of magnetization is then used in a standard inver-sion code to recover a susceptibility model. The combined method greatly reducesambiguity about the location and orientation of isolated targets, but becomes chal-lenging over multiple anomalies with variable magnetization directions. The samesynthetic model was applied to the MVI method, resulting in a poor recovery forthe location and intensity of magnetization. This result is surprising consideringthe satisfactory solution obtained in previous studies (Ellis et al., 2012; Lelie`vreand Oldenburg, 2009). The combined approach of Liu et al. (2015) is interestinghowever, as it links together complementary algorithms.1.2 Thesis outlineFollowing the same line of thought as Liu et al. (2015), I propose a CooperativeMagnetic Inversion (CMI) algorithm that directly incorporates the amplitude inver-sion of Shearer (2005) into the MVI algorithm of Lelie`vre and Oldenburg (2009).Magnetic field amplitude data are computed from the Equivalent Source techniqueproposed by Li and Oldenburg (2010), removing the need to grid the data. More-over, I introduce a Scaled Iterative Re-weighted Least Squares (S-IRLS) regular-ization method to recover blocky and sparse solutions. My thesis is divided intothe following chapters:In Chapter 2 and 3, I provide further details about rock magnetism and reviewthe theory related to inverse problems in the context of exploration geophysics. I re-visit three inversion codes from the literature: the magnetic susceptibility inversionfrom (Li and Oldenburg, 1996), the magnetic amplitude inversion from Shearer(2005), and the Magnetic Vector Inversion (MVI) from (Lelie`vre, 2009). Eachof these codes are tested on a synthetic example with complicated magnetizationdistribution. I also review the equivalent-source method from Li and Oldenburg8(2010).In Chapter 4, I introduce a mixed lp-norm regularization function for the recov-ery of compact and sparse solutions, increasing the flexibility of current inversioncodes. The method uses a Scaled Iterative Re-weighted Least Squares (S-IRLS)formulation to approximate any lp-norm penalties on the interval 0≤ p≤ 2. Sparseconstraints are applied on both the model and model gradients independently. Themixed-norm regularization is implemented on the magnetic problem, and tested onan airborne magnetic dataset from the Tli Kwi Cho kimberlite complex, NorthwestTerritories.In Chapter 5, I formulate a Cooperative Magnetic Inversion (CMI) method,combining the inversion codes presented in Chapter 3 and 4. The algorithm istested on the same synthetic example for comparison.Finally in Chapter 6, I apply the CMI algorithm to a large aeromagnetic dataset over the Ekati Property, Northwest Territories. I design a tiled inversion schemeto automate the inversion process and to reduce the computational cost. Individualmodels are then merged back onto a global mesh for analysis. A bulk estimateof magnetization direction is compared to values published in the literature. Thepolarities of magnetization over 11 pipes are compared to the estimated age ofemplacement. The analysis is extended to various dyke swarms, providing regionalgeologic information.The work presented in this thesis is significant for two reasons. First, themixed-norm regularization function introduced in Chapter 4 can considerably im-prove the flexibility of inversion algorithms, not limited to potential field problems.The S-IRLS method allows for a combination of sparse norms on the model andmodel gradients independently, granting access to a wider range of solutions thanpreviously offered by globally convex functions.Secondly, it is the first time that both the amplitude inversion and MVI algo-rithms are combined in a cooperative inversion method in 3-D. Structural infor-mation gained by the amplitude inversion is used directly to constrain the MVIalgorithm, reducing the non-uniqueness of the solution. From a practical stand-point, it is an automated process that streamlines the inversion workflow. This, inturn, can reduce the overall time required to process magnetic data and promotethe expansion of inversion methods in mineral exploration.9Chapter 2Magnetostatic ProblemThe general theory describing the magnetic static problem can be derived fromMaxwell’s equations. Assuming no free currents and no time varying electric field,the equations describing the electromagnetic field simplify to:∇ ·~B= 0 (2.1)∇× ~H = 0 (2.2)~B= µ~H , (2.3)where ~B is the magnetic flux density in Tesla (T) and ~H is the magnetic field inunits of amperes per meter (A/m). In matter, the magnetic permeability µ relatesthe magnetic flux density ~B to the field ~H such that:µ = µ0(1+κ) , (2.4)where µ0 is the magnetic permeability of free space. The magnetic susceptibilityκ is a dimensionless positive number describing the ability of certain material tobecome magnetized under an applied field. From 2.2, the magnetic field can beformulated via a potential field formulation such that:~H = ∇φ , (2.5)10Because Maxwell’s equations forbid magnetic monopoles from 2.1, we approxi-mate the potential field φ by a magnetic dipole. The potential from a dipole mo-ment ~m located at some location rQ as observed at location rP is given by:φ(r) =14pi~m ·∇(1r)(2.6)r = |rQ− rP|=√(xP− xQ)2+(yP− yQ)2+(zP− zQ)2~m=[mx, my, mz].Going from a discrete dipole moment to continuous magnetization allows us towrite:φ(r) =14pi∫V~M ·∇(1r)dV (2.7)~M =MxMyMz ,where ~M is the magnetization vector per unit volume in A/m. Taking the gradientof 2.7, the magnetic flux density can be expressed as:~b(P) =µ04pi∫V~M ·∇∇(1r)dV (2.8)~b(P) =bxbybz ,where~b(P) is the magnetic vector measurement at location P.In a geophysical context, we are interested in identifying discrete volumes ofmagnetic material. As presented in Sharma (1966), the integral equation can beevaluated analytically for the magnetic field of a rectangular prism. The integral 2.811then becomes:~b(P) = T · ~M , (2.9)where the tensor matrix T takes the form:T=Txx Txy TxzTyx Tyy TyzTzx Tzy Tzz . (2.10)Note that only five of the nine elements forming T are independent and need to becomputed. Expanding 2.9 in terms of components of the field yields:bx = Txx Mx+Txy My+Txz Mzby = Tyx Mx+Tyy My+Tyz Mzbz = Tzx Mx+Tzy My+Tzz Mz .(2.11)From 2.8, I divide the earth into nc discrete prisms, each of which has a constantmagnetization such that:~b(P) =nc∑j=1T j · ~M j . (2.12)From the superposition principal, the xˆ-component of magnetic flux density~b(P)corresponds to the cumulative contribution of all nc cells such that:bx(P) =[T 1xx . . . Tncxx T1xy . . . Tncxy T1xz . . . Tncxz]MxMyMzwhere, Mx =M1x...Mncx My =M1y...Mncy Mz =M1z...Mncz ,(2.13)and likewise for the by(P) and bz(P) components. In compact matrix vector nota-tion, 2.13 can be written as:~b(P) = T ~M , (2.14)12where the augmented tensor matrix T ∈ R3×(3nc) multiplies a vector of magnetiza-tion direction ~M ∈ R3nc.Because magnetic flux measurements can be recorded at multiple locations,2.14 is further augmented by a factor N such that:~b=bxbybzwhere, bx =bx(P1)...bx(PN) by =by(P1)...by(PN) bz =bz(P1)...bz(PN) .(2.15)This ordering defines the rows of the final matrix used for the forward calculationsof magnetic data:~b= T ~Mwhere, ~M=MxMyMz , T=Txx Txy TxzTyx Tyy TyzTzx Tzy Tzz~b ∈ R(3N) , T ∈ R(3N)×(3nc) , ~M ∈ R(3nc) .(2.16)Most geophysical surveys do not collect the components of the magnetic field~b, but rather its amplitude, or Total Magnetic Intensity (TMI) such that:bTMI = |~B0+~BA| , (2.17)where ~B0 is the primary geomagnetic field and ~BA are anomalous local fields. Inmineral exploration, we are only interested in the anomalous field arising frommagnetized rocks. The quantity of interest is refered to as Total Magnetic Anomaly,or bTMA :bTMA = |~B0−~BA| . (2.18)13Assuming that |~BA||~B0| 1. the anomalous field is approximated as:bTMA ' ~BA · Bˆ0= |~B0+~BA|− |~B0| .(2.19)The anomalous data bTMA can be obtained from 2.16 by a projection matrix actingon T such that:bTMA = PT ~M (2.20)P=1|~B0|[B0xI B0yI B0zI]P ∈ RN×(3N) ,(2.21)where the projection matrix P computes the inner-product between each magneticfield measurements and the inducing field, and I is an N×N identity matrix.2.1 MagnetizationWe have so far remained general and have omitted any details regarding the mag-netization vectors ~M. In matter, the total magnetization per unit volume can beexpressed as:~M = κ~H+ ~MNRM , (2.22)where ~MNRM is known as the Natural Remanent Magnetization and the magneticsusceptibility κ is the intrinsic physical property of rocks describing their ability tobecome magnetized under an applied field ~H. I will here assume that the inducedmagnetization is isotropic. The inducing field ~H can be further divided in twoparts:~H = ~H0+ ~HA , (2.23)where ~H0 is the Earth’s geomagnetic field and ~HA are the anomalous local fields.The geomagnetic field ~H0 is generally dominant and believed to be generatedwithin the core Campbell (1997). Smaller diurnal variations can be observed dueto the movement of charged particles within the upper atmosphere.Under the assumption of no free currents, secondary fields ~HA arise from the in-14teraction of magnetic objects, or self-demagnetization effects. For equidimensionalobjects, secondary induced fields oppose the geomagnetic field and reduce theoverall induced magnetization. For an elongated magnetic body, the total magneti-zation direction may be deflected towards the principal axis. Self-demagnetizationeffects usually become important as κ gets large (κ ≥ 0.1). For the remainder ofthis research project I will assume that self-demagnetization effects are negligible.Lastly, the Natural Remanent Magnetization, or ~MNRM component, is a perma-nent dipole moment preserved in the absence of an inducing field. Certain rocks,known as ferromagnetic material, can retain a net magnetization vector, which insome cases will reflect the orientation of the Earth’s field during formation. Theorientation and magnitude of the NRM component is generally unknown and dif-ficult to distinguish from the induced component. From laboratory measurements,the strength of remanence is often expressed as a ratio with respect to the inducedcomponent:Q=|~MNRM|κ|~H0|, (2.24)where Q stands for the Koenigsberger ratio. The NRM component has been foundto play an important role in many geological settings. Due to frequent changesin the Earth’s polarity and various geological processes, rock magnetization direc-tions can vary greatly within a given region. Accurately determining the orientationof magnetization is the principal focus of this research project.While most magnetic methods proposed in the literature assume a purely in-duced response, a more general approach must account for variability in magneti-zation direction. From a geophysical standpoint, we would like to infer the distribu-tion of magnetic material from the observed field data. Modeling magnetic objectsthough the inverse problem without knowledge about the magnetization directionhas proven to be difficult and remains a field of active research in geophysics.15Chapter 3Inverse Problem3.1 Regularized inversionIn chapter 2, I introduced the linear forward calculation for the anomalous magneticdata bTMA generated by volumes of magnetization. In a more general case, thediscretized system of equations relating the model to the geophysical data can bewritten as:F[m] = d , (3.1)where F is a generic operator relating the geophysical data d ∈ RN to discretemodel parameters m ∈ Rnc. In some cases, F is linear, but can be non-linear as isthe case with amplitude data. We are usually interested in finding a solution to theinverse problem:m= F−1d , (3.2)namely to recover some model parameters m responsible for a set of observationsd. At this point, the magnetization model m is general. In this thesis, m will takedifferent forms depending on assumptions made, coordinate system chosen or datatransformation.It is difficult to solve the inverse problem for two reasons. First, the inverseproblem is often ill− posed as the number of unknown parameters largely exceedsthe number of observations. The solution is highly non-unique. Secondly, fielddata are generally corrupted by random noise e such that the true linear system16should be written as:F[m] = dobsdobs = d+ e .(3.3)Even if the problem is linear and the matrix F is full ranked and invertible, theproblem is said to be ill-conditioned. Any solution that satisfies 3.3 is unstable,as small changes in the noise can induce large changes in model values. More-over, most geophysical methods have a strong spatial dependency related to thedistance between the source and observation location. Consequently, few com-ponents describing the system F can be significantly larger, which can adverselyimpact numerical solvers.The inverse problem described by 3.3 may admit an infinite number of solu-tions that are geologically unrealistic, and the solution may not be very stable withrespect to the observed data. As first introduced by Tikhonov and Arsenin (1977),the inverse problem can be formulated as a regularized least-squares problem ofthe form:minmφ(m)φ(m) = φd + βφmφd = ‖Wd(F[m]−dobs)‖22φm = R ,(3.4)where the optimal solution is found at the minimum of the objective function φ(m).The misfit function φd measures the residual between predicted and observed datadobs normalized by the estimated uncertainties Wd :Wd =1/σ1 0 . . . 00 1/σ2 0....... . . 00 . . . 0 1/σN , (3.5)where σi are assigned standard deviations. Assuming the noise to be Gaussian17and uncorrelated, the misfit function follows a chi-squared distribution with anexpected value of N. Under this assumption, the expected data misfit is also equalto the number of data N.The model objective function φm is added to stabilize and constrain the solu-tion. The trade-off parameter β balances the relative influence between the misfitfunction and any a priori information prescribed by a chosen regularization R.There has been much research done on designing robust and effective regulariza-tion functions. In the original magnetic inversion work of Li and Oldenburg (1996),the regularization function involves a measure of model smallness and smoothness.The general objective function takes the form:φ(m) = φd+β[αs∫Vws(r)∣∣m(r)−mre f ∣∣2 dV + ∑i=x,y,zαi∫Vwi(r)∣∣∣∣∂m(r)∂xi∣∣∣∣2 dV] ,(3.6)where αs, αx, αy and αz are adjustable constants balancing the relative contributionbetween the various components of the model objective function. The first integralmeasures the deviation from a reference model mre f . The three following integralsmeasure the spatial gradients of the model m(r) in Cartesian coordinates. Theweighting functions ws(r), wx(r),wy(r) and wz(r) are cell-based penalties addedto the system to reflect specific characteristics expended from the solution. Forpotential field problems, it is common to resort to a sensitivity-based weighting inorder to compensate for the natural decay of the kernel functions, such that:ws(r) = wr(r) w˜s(r)wi(r) = wr(r) w˜i(r) ,(3.7)where wr(r) is a general distance weighting, and w˜s, w˜i are customizable weightsthat reflect any available a priori information. Note that the distance weightingwr(r) is applied through the regularization function rather than directly to the sen-sitivity matrix as prescribed by Li and Oldenburg (1996). Reasons for this changewill be discussed in Section 4.4.1 while experimenting with different lp-norm reg-ularization functions.18I discretize 3.6 on a tensor mesh made of rectangular prims such that:φ(m) = φd+β[‖Ws (m−mref)‖22+ ∑i=x,y,z‖Wi Gi m‖22], (3.8)where the diagonal matrices Ws and Wi ∈ Rnc×nc contain dimensional scales thatarise from the discretization and cell-based weights (ws(r), wi(r)). The globalscaling constants αs and αx are also absorbed as cell weights. The measure ofspatial gradients are calculated by a generic forward difference scheme such that:Gi =−1 1 0 . . . 00. . . . . . . . ........ . . 0 −1 10 . . . 0 1 −1 , (3.9)where the gradient operator Gi ∈Rnc×nc calculates the horizontal changes in modelparameters m. I voluntarily removed all spatial dimensions out of the gradientoperators for reasons that will be discussed in Section 4.4.1. The last row is usinga backward difference in order to get a square matrix. The banding structure of thegradient operators is different for Gx, Gy and Gz to account for the ordering of thecells in the tensor mesh.For ease of notation, I write the model objective function in compact form as:φm = ‖Wm (∆m)‖22 , (3.10)such that:WTmWm =WTs Ws+ ∑i=x,y,zGTi WTi WiGi ,as well as:∆m=m−mre f ,to define the deviation between the model parameters m and the reference modelmre f .193.1.1 Iterative solverAs previously stated, my goal is to minimize the objective function described by3.8. The minimum solution is found where the partial gradients of the function arevanishing such that:g(m) =∂φ(m)∂m= 0 . (3.11)Taking the partial derivatives of 3.8 with respect to the model parameters yield:JTWTdWd[F[m]−dobs]+βWTmWm (∆m) = 0 , (3.12)where J is the Jacobian of the forward operator:J=∂F[m]∂m. (3.13)More generally, we are interested in solving the nonlinear system described by3.12, which can be done with second-order methods. Newton’s method computesa series of model updates δm by solving the system:H δm=−g(m( j)) , (3.14)where H is the Hessian, or second-order derivatives of the objective function. Theoptimization problem is solved iteratively such that:m( j+1) =m( j)+αδm , (3.15)where the superscript ( j) denotes the solver iterations for a model perturbation δmscaled by a step length α .Computing the true Hessian can be computationally expensive. For least-squares problems, the Gauss-Newton method can be used to approximate the Hes-sian such that:∂ 2φ∂m2≈ H˜= JTWTdWdJ+βWTmWm , (3.16)where the assumption is made that the change in Jacobian is small between eachiteration:J(m( j)+δm)≈ J(m)( j) , (3.17)20and that the model update is approximated to be along the gradient direction suchthat:F[m+δm]≈ F[m]+Jδm . (3.18)Combining 3.12 and 3.16, and ignoring the superscript ( j) yields:H˜ δm=−g(m))(JTWTdWdJ+βWTmWm)δm=−JTWTdWd[F[m]−dobs]−βWTmWm (∆m) .(3.19)Solving 3.19 gives a step direction. Each Gauss-Newton step requires a solu-tion to linear system of the form A x = b. In our case, the left-hand-side of 3.19can be assumed to be a symmetric positive definite matrix. Hence it is possiblefind a unique solution by directly computing the inverse of the pseudo-Hessian.From a practical standpoint however, the computational cost of this operation on adense matrix is O(nc3), which can rapidly become prohibitive for a large systemof equations.Krylov Space methods have been proposed as an alternative to direct solvers.Iterative solvers, such as the Conjugate Gradient (CG) method, apply a series oforthogonal steps towards the minimum of a function. Each step is said to be A-orthogonal to the previous ones, contributing to the fast convergence of CG oversimpler gradient descent methods. The reader is encouraged to read Shewchuk(1994) for a comprehensive review of the method. The algorithm is presentedin Table 3.1. It can be shown that for a symmetric positive definite matrix A,the CG method converges to a unique solution in at most O(nc) operations. It israrely required to solve the system exactly however. My algorithm uses a stoppingcriteria for the minimum CG update (δ > 1e− 4), after which a Gauss-Newtonstep direction δm is returned. The step length α is calculated by a line-searchmethod as presented in Table 3.2. The Gauss-Newton steps are repeated untilthe model updates falls below some threshold value, |αδm| < γ , where γ is somesmall value. In this thesis, I fix γ = 1e−4, which experimentally has proven to bea good compromise between computation cost and accuracy.In order to reduce the amount of memory required to store the dense matrix H˜21Table 3.1: Conjugate Gradient algorithmInitialize: d(0) = r(0) = b−A x(0)while : ‖r‖> δα(i) =rT(i)r(i)dT(i)A d(i)x(i+1) = x(i)+α(i)d(i)r(i+1) = r(i)−α(i)Ad(i)β(i+1) =rT(i+1)r(i+1)rT(i)r(i)d(i+1) = r(i+1)+β(i+1)d(i)Table 3.2: Line-searchInitialize: α = 1 , mˆ=m(j)while : φ(mˆ)≥ φ(m(j))α = α/2mˆ=m(j)+αδmm(j+1) = mˆin 3.19, the Gauss-Newton steps can be formulated as an overdetermined problemof the form: Wd J√βWs√βWx Gx√βWy Gy√βWz Gzδm=−Wd(F[m]−dobs)Ws (∆m)Wx Gx (∆m)Wy Gy (∆m)Wz Gz (∆m) . (3.20)The CG steps presented in Table 3.1 are altered slightly in order to calculate theresidual r(i) and descent direction d(i). Memory savings are important as I avoidforming a dense system H˜∈Rnc×nc, but only require sparse matrix-vector products.3.1.2 PreconditionerAs explained in Section 3.1.1, I want to solve iteratively a normal equation ofthe form Ax= b. The performance of iterative solvers depends strongly on the22condition number κ of the left-hand side matrix A such that:κ =λmaxλmin, (3.21)where λmax and λmin are the largest and and smallest eigenvalues of A. It canbe shown that an upper bound on the convergence rate of CG is function of thiscondition number such that:ω ≤ κ−1κ+1, (3.22)where the convergence rate ω gets worse as the condition number increases. Inorder to accelerate the convergence, it is common to resort to a preconditioner toimprove the spectral properties of A (van der Vorst, 2003). The basic idea is topre-multiply the linear system such that:P Ax= P b , (3.23)where the resulting condition number of P A is smaller than the original matrixA. Assuming that A is invertible, the perfect preconditioner is its inverse A−1, inwhich case A−1A is the identity matrix with a condition number of 1. Computingthis preconditioner would obviously be redundant, since if I know A−1, then Iwould have already solved the problem.I can decompose the matrix A in its lower and diagonal blocks:A= L+D+LT , (3.24)where in our case the matrix D is the diagonal of the pseudo Hessian describedby Equation 3.19. The Jacobi preconditioner is the simplest choice where I setP= D−1. In this research project, the matrix A is a positive definite matrix, henceD is assumed to have diagonal elements that are strictly non-zero positive numbers.The preconditioning matrix P is updated before each Gauss-Newton step and pre-multiplied ahead of the Conjugate Gradient solver presented in Table 3.1.233.1.3 Bound constraintsIn many cases, the solution m is expected to be bounded on a specific interval. Inparticular for the magnetic problem, the magnetic susceptibility κ can only be apositive real number on the interval κi ∈ [0,∞). Three strategies have been pro-posed in the literature. In Li and Oldenburg (2003), positive constraints are appliedvia a log-barrier penalty term such that:φ(λ ) = φd+βφm−2λnc∑i=1ln(κi) , (3.25)where the log barrier parameter λ controls the influence of the logarithmic penaltyfunction. The method solves a series of non-linear problems by monotonicallyreducing λ until convergence to a stable solution has been reached. As κi → 0,the log(κi) tends to be large and negative, hence strongly penalized during theminimization process. The same approach can be used to impose and upper boundon the model parameters.In a second approach, the model values are parameterized into a quantity thatcan only be positive. Lelie`vre (2003) proposes transforming the model parameterssuch that:mi = κ(1/2)i , (3.26)in which case if mi ∈ (−∞,∞), then κi ∈ [0,∞), hence a strictly positive value. Thismethod has the advantage of reducing the size of the system to be solved, but canonly be used as a lower bound.The third strategy, which is used in this work, is the Projected Gradient methodas presented in Vogel (2002). Within the line-search algorithm presented in Ta-ble 3.2, model parameters outside the bounds are replaced and set to be inactive forthe following Gauss-Newton iteration.All three strategies make the inverse problem non-linear with respect to m. TheProjected Gauss-Newton method is the simplest to implement however within aniterative solver framework. It does not require additional parameters to be adjustedand can be used to impose both an upper and a lower bound. All the algorithms pre-sented in this thesis were implemented with the Projected Gauss-Newton method.243.2 Magnetic susceptibility inversionAs formulated in 2.20, the observed Total Magnetic Anomaly data bTMA are re-lated to discrete magnetization parameters ~M by a large under-determined systemof equations P T ∈ RN×(3∗nc). A long standing strategy, as adopted by Li andOldenburg (1996), Pilkington (1997) and others, is to neglect the effect of self-demagnetization and remanence and invert directly for the magnetic susceptibilityof rocks. Assuming a uniform magnetization direction parallel to the Earth’s field~H0, the linear system 2.11 can be written as:bx =[Txx H0x +Txy H0y +Txz H0z]κby =[Txx H0x +Txy H0y +Txz H0z]κbz =[Txx H0x +Txy H0y +Txz H0z]κ ,(3.27)where the magnetic field response is linearly related to the effective susceptibilityparameter κ . For a larger system involving N observations over a model spacediscretized into nc cells, 3.27 is written as:bTMA =P T H κH=H0xIH0yIH0zI , κ =κ1κ2...κMH ∈ R(3nc)×nc , κ ∈ Rnc ,(3.28)where I is the identity matrix and H is a block diagonal matrix containing the orien-tation of the inducing field ~H0 in Cartesian coordinates. Since the inducing field ~H0is assumed to be constant, the matrix product P TH can be directly incorporatedinto a linear system of equations such that:bTMA = F κ , (3.29)25where F ∈ R(N)×nc is the linear forward operator relating the discrete magneticsusceptibility values κ to the magnetic field observations bTMA.3.2.1 Synthetic exampleHaving defined the linear system relating the anomalous field data bTMA to a dis-crete magnetic susceptibility model κ , I illustrate the forward calculation with asynthetic example shown in Figure 3.1. The model consists of a folded anomalywith induced magnetization of 3 A/m, arching around a discrete block with mag-netization of 2 A/m. The arc-shaped anomaly dips 20◦ towards the south. Fromthe relation between induced magnetization and ambient field:~Mind = κ~B0/µ0 , (3.30)the susceptibility of both anomalies are 0.075 and 0.05 SI respectively, subject to a50,000 nT inducing field ~H0. I create this model in order to simulate complex ge-ology with varying shapes and depth extent. The model is discretized on a uniform20 meters cell size mesh. A grid of 342 observation stations is placed on a plane 20meters above the top of the mesh. I generate bTMA data using 3.29 as shown in Fig-ure 3.2. Data are corrupted with random Gaussian noise, 1 nT standard deviation,to simulate field data conditions.I then follow with the inverse problem, namely to recover a distribution ofsusceptibility values from the bTMA data as presented in Li and Oldenburg (1996).Similar to the general system 3.4, the objective function to be minimized is writtenas:φ(κ) = ‖Wd( Fκ−dobs )‖22+β‖Wmκ‖22 . (3.31)In this case, Wd is simply the identity matrix since the noise is known to have astandard deviation of exactly one. I also set the reference mode mre f to zero, henceforcing the solution to be small and smooth. Table 3.3 summarizes the inversionparameters used.I minimize 3.31 iteratively following the procedure presented in Section 3.1.1.Figure3.3 shows the convergence curve for the data misfit φ (k)d and model norm φ(k)mas a function of β iterations. The inversion achieved the target misfit after seven26Figure 3.1: (a) Synthetic susceptibility model consisting of a folded anomaly(κ = 0.075 SI) arching around a discrete block (κ = 0.05 SI) .Figure 3.2: (a) Data generated from the synthetic susceptibility model subjectto a vertical 50,000 nT inducing field. (b) Data are then corrupted with (c) ran-dom Gaussian noise, 1 nT standard deviation.27Table 3.3: Inversion parameters.Core cell size 20×20×20 mnc cells 82,000N data 342αs, αx, αy, αz 2.5e-3, 1, 1, 1ws, wx, wy, wz 1, 1, 1, 1Uncertainties 1 nTiterations, after which the misfit function levels off. In cases where the true noiselevel is unknown, it is common practice to chose the point of maximum curvatureon the misfit curve as the optimal model. Letting the inversion progress furtherdown the misfit curve comes at the risk of fitting some of the noise, potentiallyintroducing artifacts in the model.Figure 3.3: Convergence curve showing the data misfit φ (k)d and model normφ (k)m as a function of β iterations. The inversion achieves target misfit after the 7thiteration, which in this case also corresponds to the point of maximum curvatureon the misfit curve. Attempting to further lower the data residual comes at therisk of fitting some of the Gaussian noise.Sections of the recovered susceptibility model from the 7th iteration are shown28in Figure 3.4. Both anomalies are recovered at roughly the right location and atdepth. As expected from the l2-norm regularization, the model is smoothly varyingspatially. Peak susceptibility values are underestimated by roughly 20% due to thesmooth constraint. The model can predict the data well within one standard devi-ation on average (Fig. 3.5). It is worth nothing that there are correlated residualsdirectly above the magnetic anomaly. Even though the global misfit criterion issatisfied, there is clearly information in the data that is not being captured. Oneoption would be to lower the target misfit, but at the risk of fitting some of therandom noise. I tested that strategy on this example, which returned speckly mod-els with isolated susceptibility anomalies near the surface. A second option wouldbe to manually reduce the uncertainties only over the correlated residuals. Thiskind of processing is time consuming and somewhat arbitrary. I will argue that thesmooth regularization prevents the inversion from fitting the high frequency con-tent. In Chapter 4 I introduce a compact norm regularization in order to test thishypothesis.29Figure 3.4: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredsusceptibility model for a purely induced response. The model is smooth butrecovers the arc and block anomaly at roughly the right depth.Figure 3.5: Comparison between (a) observed and (b) predicted data from therecovered susceptibility model. (c) The normalized data residuals appear to becorrelated with the location of the magnetic body.303.2.2 Effect of remanent magnetizationRecent studies have shown that that the assumption of a purely induced responsemay not be adequate in many geological settings (Buchan et al., 2009; Enkin,2014). As presented in Lelie`vre (2009), incorrect assumptions regarding the ori-entation of magnetization can make the inversion process unreliable. To simulatethis problem, I alter the synthetic model presented in Figure 3.1. I change the ori-entation of magnetization along the arc-shaped anomaly as shown in Figure 3.6 tosimulate a geological folding of a ferromagnetic object. The magnetization incli-nation is orientated at 45◦ from the horizontal, with variable declinations between[−45◦N ; 45◦N]. Data are generated from the new magnetization model and cor-rupted with random Gaussian noise, 1 nT standard deviation ( Fig 3.7 ). Comparedto the purely induced response, the positive anomaly from the arc shifts towardsthe center while creating a large negative anomaly on the outside boundary of thearc.Following the same procedure, I invert this new data set while still assuminga purely vertical magnetization direction. As presented in Figure 3.8, the inver-sion poorly recovers the location of the arc-shaped body. A broad susceptibilityanomaly is created on the periphery of the data set, pushing susceptibilities out-ward in order to account for the large negative data.. The signal from the arc alsoaffects the center block anomaly, that is now recovered deeper than the true loca-tion. As seen in Figure 3.9, the inversion algorithm has difficulty fitting the strongnegative data. In a mineral exploration context, the model presented in 3.8 couldresult in false drilling targets —costly both in time, resources and confidence ingeophysical methods. This is a well known problem in mineral exploration, thus aneed for inversion methods that can account for remanence.31Figure 3.6: Perspective view and sections through the synthetic magnetizationmodel. The arc-shaped anomaly is magnetized at 45◦ from horizontal and withvariable declinations directions between [−45◦N ; 45◦N].Figure 3.7: (a) Data generated from the synthetic magnetization model. (b)Observed data are corrupted with (c) random Gaussian noise, 1 nT standard de-viation.32Figure 3.8: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredsusceptibility model assuming no remanence. The arc-shaped anomaly is poorlyrecovered and magnetic susceptibilities are pushed at depth and outwards.Figure 3.9: Comparison between (a) observed and (b) predicted data from therecovered susceptibility model assuming a purely induced response. (c) The in-version has a harder time fitting the large negative fields along the arc.333.3 Magnetic vector inversionAs shown in the previous example, the assumption of a purely induced magneti-zation direction has its limitations. In his PhD thesis, Lelie`vre (2009) introducesa Magnetization Vector Inversion (MVI) procedure in order to recover both thelocation and orientation of magnetization. From the general expression for themagnetic field of a prism in 2.11, the forward problem is formulated as:~b= T · ~Mp+T · ~Ms+T · ~Mt , (3.32)where I define three orthogonal components of magnetization. The primary pˆ com-ponent is set to be aligned with the inducing field, the second sˆ component has thesame azimuth as the primary but is rotated 90◦ in inclination. The third componenttˆ lies on the xy-plane and perpendicular to pˆ and sˆ. I define the transformation fromCartesian to {p,s,t} coordinate system by a double rotation such that:[pˆ sˆ tˆ]= Rz(θ) Rx(φ)0 0 11 0 00 1 0 , (3.33)where θ and φ are the declination and inclination of the inducing field as definedin a Cartesian coordinate, positive counter-clockwise. For an inducing field ~H0oriented 0◦I, 0◦D, the components of magnetization pˆ, sˆ and tˆ would be alignedwith the Cartesian components yˆ, zˆ and xˆ respectively.For the bx component of the magnetic field due to a single prism magnetizedalong the {p,s,t} direction can be written as:bx =[Txx Txy Txz] HpˆxHpˆyHpˆzκp+HsˆxHsˆyHsˆzκs+HtˆxHtˆyHtˆzκt , (3.34)and likewise for the by and bz components of the magnetic field. The magnetic data~b is now function on three orthogonal components of effective susceptibility κp, κsand κt and their respective direction of magnetization ~Hp, ~Hs and ~Ht defined by the34{p,s, t} coordinate system such that:~Hp =|~H0| pˆ~Hs =|~H0| sˆ~Ht =|~H0| tˆ .(3.35)Just as I did for the magnetic susceptibility problem, I augment the linear sys-tem for the magnetic response of nc magnetized prisms as observed at N data loca-tions. Equation 3.32 becomes:bTMA = Fκˆ= P T[Hp Hs Ht]κpκsκt . (3.36)where Fˆ ∈ RN×(3nc) has now three times the number of variables as the suscepti-bility inversion described in 3.29. The block diagonal matrices Hp,Hs and Ht ∈R3nc×nc relate the three effective susceptibility parameters κp,κs and κt ∈ Rnc tothe system of equation T in Cartesian coordinates.The objective function to be minimized becomes:φ(κˆ) = ‖Wd( Fκˆ−dobs )‖22+β‖Wˆmκˆ‖22 . (3.37)where the weighting matrices and gradient operators comprised by Wˆm are blockdiagonal matrices such that:Wˆ =W . . . 0... W...0 . . . WGˆ =G . . . 0... G...0 . . . G .35Each operator is three times the size as those used for the susceptibility inversion.The inversion process follows the same method as the susceptibility inversionand uses the data presented in Figure 3.7. The only difference is that effectivesusceptibilities are allowed to go negative. Figure 3.10 presents sections of therecovered magnetization model, with color scale representing the total effectivesusceptibility κe where:κe = (κ2p+κ2s +κ2t )1/2. (3.38)The inverted model recovers the approximate orientation of magnetization inthe block anomaly and part of the arc. The result is smooth however, making it hardto distinguish the extremities of the arc. Effective susceptibilities associated withthe center block extend deeper than the true model. Large effective susceptibilityvalues are found near the top of the mesh, leaving holes in the iso-surface.The current inversion strategy may work well in simple cases, or when sub-stantial a priori information is provided. For more general cases, the smoothmodel presented in Figure 3.10 could lead to false interpretation due to the lackof boundaries. There is still need to develop a robust algorithm that could reducethe non-uniqueness of the MVI formulation with minimal input required from theuser.36Figure 3.10: (a) Iso-surface (κe =0.001) and (b) sections through the recoveredmagnetization model from the MVI method. The inversion recovers the trueorientation of magnetization inside the block, but the thin arc is poorly resolved.Figure 3.11: Comparison between (a) observed and (b) predicted data from therecovered magnetization model from the MVI method. The model can replicatethe data at the same level achieved by the purely induced problem.373.4 Magnetic amplitude inversionThe inversion of magnetic quantities, weakly dependent on the magnetization di-rection, have been proposed as an alternative to the full magnetic vector inversionapproach (Li et al., 2010). In this section I review the work done by Shearer (2005)in inverting magnetic amplitude data. Amplitude data have the interesting propertyof being weakly dependent on the orientation of magnetization (Nabighian, 1972),which can be used to get an estimate for the location of magnetized material. Iillustrate the idea with a single vertical dipole located at the origin. The magneticfield~b measured at some distance~r can be written as:~b=3µ0m4pi|r|3(sinφ cosφ [cosθ xˆ+ sinθ yˆ]+ (cos2 φ − 13) zˆ), (3.39)where θ is the angle on the xy-plane, φ is the angle with respect to the z-axis andm is the dipole moment. Computing the magnitude of the field yields:|~b|=√b2x+b2y+bz2=µ0m4pi|r|3√3 cos2 (φ)+1 .(3.40)The amplitude of the field |~b| is strongly dependent on the radial distance |r|, butweakly dependent on φ , the inclination of the dipole. If we were to collect magneticfield data at a constant height above the dipole, the location of maximum amplitudewould roughly occur above the location of the dipole, and the amplitude of the fieldwould vary by at most a factor two. Hence amplitude data may reduce some of theambiguity related to the horizontal distribution of magnetized material.Recall from Chapter 2, I have defined the component of the magnetic field atsome location P due to a prism with uniform magnetization ~M as:bx = [Txx Txy Txz] ~Mby = [Tyx Tyy Tyz] ~Mbz = [Tzx Tzy Tzz] ~M ,where the forward operator T maps the components of the magnetic field due to a38distribution of magnetized prisms ~M. Just as in the MVI method, the magnetizationdirection is unknown in most cases. The advantage of working with amplitude datais that the magnetization direction does not have to be known exactly in order toget a good estimate of |~b|. Under the assumption that the magnetization vectoris in most part parallel to the geomagnetic field direction, equation 2.22 can besimplified to:~M = κ(~H0+ ~Hs)+ ~Mrem≈ ~H κe ,where the effective susceptibility κe is a unitless quantity representing the totalmagnetization of each cell. Assuming a uniform inducing field ~H, the anomalousmagnetic field measurement due to a magnetized cell is approximated by:~b(P) =bxbybz=Txx Txy TxzTyx Tyy TyzTzx Tzy TzzκeHxκeHyκeHz=FxFyFz κe .(3.41)Note the parallel with the κe used in the MVI method to describe the magnetizationcomponents pˆ, sˆ and tˆ. Because the geomagnetic field is generally much larger thanany secondary fields, I approximate the magnetization direction to be parallel tothe inducing field (~H ' ~H0). From the linear relation described above, I write theforward calculation of magnetic amplitude data as:|~b|=[bx2+by2+bz2]1/2=[(Fx κe)2 +(Fy κe)2 +(Fz κe)2]1/2.(3.42)39Taking the partial derivative in terms of model parameter κe yields:∂ |~b|∂κe=∂∂κe[bx2+by2+bz2]1/2=1|~b|[bx∂bx∂κe+by∂by∂κe+bz∂bz∂κe]=1|~b| [bxFx+byFy+bzFz] .(3.43)The inverse problem is clearly non-linear with respect to the model parameter κe.A solution is found iteratively as outlined in section 3.1.1. At the kth iteration, thesensitivity relating the ith amplitude data due to the jth prism is given by:Ji j =~b(k)|b(k)| ·FxFyFz , (3.44)where the magnetic data ~b(k) is computed from an effective susceptibility foundat a previous iteration. In order to have a Jacobian defined at the first iteration, Ichoose a small number larger than zero for the starting model (κ(0) = 1e−4).In mineral exploration, amplitude data must be derived directly from TMAdata, which will be covered in Section 3.5. For this synthetic example, amplitudedata are generated directly from 3.42 as I already know the true magnetizationmodel. Data are corrupted with the same Gaussian noise, 1 nT standard deviation.Once again, positivity constraints are enforced on values of effective susceptibility.Figure 3.12 presents the recovered effective susceptibility model after reachingthe target misfit. Compared to the result presented in Figure 3.8, high κe valuesare recovered along the arc, and at the right depth inside the center block. Theamplitude model closely resembles the model obtained from the purely inducedresponse obtained in Section 3.2.1. I note however that the amplitude inversion hasthe tendency to stretch the model vertically. This remains an open question.Comparing the observed and predicted data, the magnetic amplitude inversionfits the data generally well, within one standard deviation as shown in Figure 5.4.Once again, the highest residuals are correlated with the location of the anomaly,40which will be tackled in Chapter 4.While neither the MVI or amplitude inversion managed to recover the locationand orientation of magnetization exactly, both methods brought complementaryinformation. From the MVI method, I am recovering a better estimate of the mag-netization direction. With the amplitude inversion, I get a closer estimate of the truelocation of magnetic anomalies. The next chapter explores different regularizationfunctions in order to further reduce the non-uniqueness of the MVI method. Mygoal is to combine those methods into a cooperative inversion work-flow in orderto a recover a simpler and more accurate solution.41Figure 3.12: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredeffective susceptibility model. The arc-shaped and block anomalies are recov-ered at the right location, but smoothly stretched vertically.Figure 3.13: Comparison between (a) observed and (b) predicted data from therecovered effective susceptibility model. The inversion can predict most of thedata within one standard deviation.423.5 Equivalent source methodThe magnetic amplitude inversion method introduced in 3.4 requires magnetic fieldcomponents [bx, by, bz] in order to calculate amplitude data such that:|b|=[bx2+by2+bz2]1/2. (3.45)While technically possible to measure three-component magnetic data, the vastmajority of current and past magnetic surveys consist of TMI measurements. Thisis largely due to the difficulty in determining the location and orientation of three-component receivers. Two methods have dominated the literature in order to extractvector components directly from TMI data, either in the frequency domain or bythe equivalent source method.As demonstrated by Bhattacharyya (1964), TMI data can be expressed as a 2-DFourier series expansion of the form:b˜(x,y,z) =∞∑n=o∞∑m=oe−2pi z(m2Lx2+ n2Ly2)1/2(Am cos2pimxLx+Bm sin2pimxLx)(Cm cos2pinyLy+Dm sin2pinyLy),(3.46)where Lx and Ly are the fundamental wavelengths in the x and y directions. Thecoefficients Am, Bm, Cn and Dn are computed from the data for each wavenumberm and n. The harmonic representation of the data in (3.46) can be seen as a productof two functions: an exponential function controlling the amplitude of the signaland an harmonic function for the spatial distribution.Converting data from the spatial to the wavenumber domain requires two im-portant assumptions: that the data are located on a plane and distributed over auniform grid. In most cases however, magnetic surveys are carried along unevenlyspaced grids and over rugged topography. Even if acquired on a plane above a flattopography, the data need to be interpolated and smoothed. The transformationbecomes even more complicated when data are collected at different elevations, asis often the case for airborne surveys with overlapping flight lines.43The Equivalent Source method has been suggested as an alternative to horizon-tal griding methods (Dampney, 1969). The technique makes use of the inherentambiguity of potential fields in determining the source location. It can be shownthat any magnetic response can be explained by an arbitrary distribution of sources.In discrete form, an equivalent source model mes can be calculated from the least-squares problem:minm‖F mes−d‖2 . (3.47)As demonstrated by Dampney (1969) and later revised by Li and Oldenburg (2010),the distance between the data and the equivalent source is limited by the frequencycontent of the signal. Most recent work by Li et al. (2014) addresses striation ar-tifacts when applying a reduction to the pole at low-latitude. They advocate for apositivity constraint formulation and prove the existence of an all-positive equiva-lent source layer in 2-D.I replicate the synthetic example presented in Li et al. (2014) as shown in Fig-ure 3.14. The model consists of 200 unit cubes with susceptibility of 0.01 SI placedin a non-susceptible background. Data are generated on a plane exactly 1 unitabove the anomaly, assuming a purely vertical inducing field of 35,000 nT. Ran-dom Gaussian noise with a standard deviation of 1 nT is added to the field compo-nents, from which TMI and amplitude data are calculated. Using the formulationin Li et al. (2014), I invert for an equivalent source with a positivity constraint.The equivalent source layer is placed at a depth that is half the data spacing, in thiscase at half a unit below the data plane. Figure 3.15 presents the equivalent sourcelayer, as well as the residual between the observed and predicted data. All threecomponents of the field, and consequently |b| data, are well recovered within thenoise level. Note however that some of the high frequency content related to theedges of the anomaly is lost in the process.44Figure 3.14: (a) Synthetic model consisting of 200 unit cubes of susceptible ma-terial in a non-susceptible background. Data are generated on a plane one unitabove the source location, assuming a purely vertical inducing field. Variouscomponents of the fields are shown in figure (b) to (f).45Figure 3.15: (a) Recovered equivalent source layer from TMI data using a pos-itivity constraint. The residuals between observed and predicted data are shownin figure (b) to (f) for various components of the field. Each component is wellrecovered within the noise level.463.5.1 Comment for future researchWhile testing the equivalent source method, I experimented on different data distri-butions to test the stability of the algorithm. Issues arose in cases where data cov-erage would partially cover the magnetic anomaly as illustrated in Figure 3.16. Inthis example, a portion of data are removed over one of the corners of the magneticsource. An equivalent-source layer is then calculated with the same parametersthat were previously used for Figure 3.15. Large correlated artifacts are created onthe bx and by component of the field along the missing portion of the data. Conse-quently, a strong narrow anomaly is recovered on |b| with amplitude in the 40 nTrange, well above the noise level. The same experiment was repeated for variousinducing field orientations, transferring the correlated artifacts to other componentsof the field accordingly. This suggests some level of ambiguity in the componentsof the field if the anomalous response is not fully captured by the data. Such ar-tifacts, if ignored during the inversion process, can have notable consequences onthe effective susceptibility model. It is to my knowledge the first time that this is-sue is raised, and it will require further research. For the remainder of this researchproject, I continue using the equivalent source method of Li et al. (2014), but I takespecial care in removing edge data with large magnetic anomalies.47Figure 3.16: (a) Recovered equivalent source layer from TMI data after remov-ing a portion of data over the corner of the magnetic anomaly. The residualsbetween observed and predicted data are shown in figure (b) to (f) for variouscomponents of the field. Note the large correlated artifacts recovered on the bx,by and |b| components.48Chapter 4Mixed Lp-norm RegularizationAs presented in Chapter 3, magnetic inverse problems are generally formulated asa regularized least-squares problem. The regularization has the dual purpose ofimproving the conditioning of the linear systems as well as to impose constraintson the solution. In this chapter, I explore different norm measures in order to im-pose sparse and blocky constraints on the solution. I introduce a new regularizationalgorithm based on the Iterative Re-weighted Least-Squares (IRLS) method to im-prove the flexibility and robustness of current inversion methods. The algorithm istested in 1-D, 2-D and 3-D on various synthetic examples. Finally, I implement thealgorithm on an airborne magnetic survey over the Tli Kwi Cho kimberlite deposit.4.1 Least-norm and regularized least-squaresAs a general introduction to inverse problems, I consider the linear system of equa-tion of the form:F m= d, (4.1)where the discretized forward operator F ∈ RN×nc is a linear system of equationsrelating a set of model parameters m ∈ Rnc to the observed data d ∈ RN . From ageophysical standpoint, we are interested in solving the inverse problem:m = F−1 d, (4.2)49where a set of unknown model parameters can be recovered from the observed data.In most cases, the inverse of F cannot be computed directly as N nc, giving riseto an undetermined system of equations. There are an infinite number of possiblesolutions satisfying 4.1, which corresponding to the null-space of F, with (nc−N)degrees of freedom. The choice of a specific answer is subjective and depends oncharacteristics expected from the true model.One simple option would be to find the smallest possible model that also min-imizes the data residual. The least-norm solution mln can be found from:mln = FT (F FT)−1 d, (4.3)where mln marks the point of shortest distance between the null-space of F and theorigin such that:mln ⊥N (F). (4.4)An equivalent answer can be found via an optimization problem of the form:minmφ(m) (4.5)φ(m) = φd + βφmφd = ‖Fm−d‖22φm = ‖m‖22 ,where φ is our objective function. As first introduced by Tikhonov and Arsenin(1977), a trade-off parameter β balances the relative importance between the dataresidual φd and the Euclidean norm of the model m. The minimum of the objectivefunction is found at the location of vanishing gradients such that:∂φ(m)∂m= 0(FTF+β I)m= FT d ,(4.6)where I is the identity matrix.I illustrate those concepts with a simple 2-variable example where I attempt to50solve:m1+2 m2 = 1 . (4.7)Equation 4.7 can be formulated as an under-determined linear system of the form:F=[1 2], d= 1. (4.8)Figure 4.1(a) displays contours along the surface formed by the l2-norm measureof data residual φd = ‖F m−d‖22 over a range of model values m. The misfitfunction can be seen as a trough with the minimum laying along the null space ofF. Any solution along the minimum of φd can satisfy 4.1. The least-norm solutionis marked with a solid dot, which is the closest distance between the null space of Fand the origin. Similarly, Figure 4.1(b) depicts the l2-norm regularization functionφm as a function of model values m, forming a symmetric paraboloid centered atthe origin.The full objective function φ(m) shown in Figure 4.1(c) can be interpretedconceptually as the sum of two competing objectives. On one hand we want to findthe best model reproducing the data, which in this case is any solution along thenull-space of F. On the other hand, we impose some constraints on the magnitudeof m in a l2-norm sense, pulling the solution towards zero. The shift between theleast-norm solution and the regularized least-squares depends on the regularizationparameter β . It can be shown that the solution to 4.6 converges to the least-normsolutions exactly as β → 0 as shown in Figure4.2(a). For a small enough β , thesolution converges to the global minimum at mln = [0.2,0.4], regardless of thestarting model.It is also important to note that the model norm regularization can be shiftedaway from the origin if we were to minimize:φm = ‖m−mref‖22 . (4.9)In this case, the minimum gradient of the objective function would occur at mref.We would therefore look for a solution that both minimizes the data residual, whilealso approximating the expected model value. For simplicity, I consider the casewhere the reference model is at origin.51Figure 4.1: Comparative contour maps for various objective functions over arange of model values [m1; m2]. (a) The mininum of the misfit function φd formsa line spanned by N (F) (red dash). The least − norm solution is marked asa solid dot. (Middle) Regularization functions and gradient directions (blue ar-rows) for approximated lp-norm measures of the model for (b) p = 2, (d) p = 1and (f) p= 0. The gradient directions are shown for two different starting mod-els (black arrows). (bottom) Contour maps of the initial objective functionsφ(m) = φd+φm for the same two starting models: m(0)1 (solid) and m(0)2 (dash).(c) In the l2-norm case , the function has a global minimum regardless of thestarting model, while for non-linear functions for (e) p = 1 and (g) p = 0, theobjective function changes with respect to the starting model.52Figure 4.2: Contour maps for various objective functions after convergence ofthe IRLS algorithm. (a) Final model obtained with the l2-norm penalty on themodel for two starting models at m(0)1 = [0.2;0.4] and m(0)2 = [0.6;0.2] for afixed trade-off parameter (β = 1e− 4). In both cases, the solution converges tothe global minimum, which is also the least−norm solution at mln = [0.2;0.4].(b) Solution with the l1-norm penalty for the same starting models and trade-offparameter, converging to a global minimum at m∗ = [0;0.5]. This solution issparse and can reproduce the data. (c) The same experiment is repeated for thel0-norm penalty, converging to two different solutions depending on the relativemagnitude of the starting model parameters. Both solutions are sparse and honorthe data.4.2 Lp-norm and iterative re-weighted least squaresI have so far only considered the Euclidean norm of the model parameters, penal-izing the square of model values. Likewise, Li and Oldenburg (1996) include anl2-norm penalty on the model gradients, yielding smooth models. In some cases,the solution may be expected to be sparse and blocky, either in terms of model pa-rameters or spatial gradients. There have been several studies dedicated to the useof non-l2 measures to recover models with sharp edges. Most methods proposedin the literature made us of the globally convex l1-norm, with proven convergenceto a global minimizer (Daubechies et al., 2010; Farquharson and Oldenburg, 1998;Sun and Li, 2014). Others have proposed approximations to the non-convex l0-53norm, methods such as the compact regularization of Last and Kubik (1983) andthe minimum support functional of Portniaguine (1999), to name a few. My goalis to further generalize the proposed methods in order to explore a wide range ofsolutions for any combinations of lp-norm penalty applied on the model and modelgradients independently.I begin with a generalized expression for the model objective function, suchthat (4.5) becomes:φm =nc∑i=1ρ(xi) , (4.10)where ρ is some norm measure of a model function x(m) commonly chosen to bethe model itself or some measure of spatial gradients. The general lp-norm measurecan be written asφm =nc∑i=1|xi|p , (4.11)where for p = 2, I recover the standard regularized inversion presented in (4.5).Several approximations to the lp-norm have been proposed, such as the Ekblomnorm (Ekblom, 1973):φm =nc∑i=1(x2i + ε2)p/2, (4.12)where a small number ε is added to guarantee that the function is continuous anddifferentiable as x→ 0. Figure 4.3(a) presents various norms for a fix thresholdvalues (ε = 1e−2). The derivative of (4.12) is given by:∂φm∂m=nc∑i=1ρ ′(xi)∂xi∂m= pxi(xi2+ ε2)1−p/2∂xi∂m.(4.13)Expression (4.13) is clearly non-linear with respect to the model function xi. Asfirst introduced by Lawson (1961), the norm can be linearized by the IterativelyRe-weighted Least-Squares (IRLS) method such that:φ (k)m =12nc∑i=1ri x2i (4.14)54∂φ (k)m∂m= ri xi∂xi∂m, (4.15)where we added the superscript (k) to denote the IRLS iterations. The weightsr(x) are computed from model values obtained at a previous iteration such that:ri =((xi(k−1))2+ ε2)p/2−1, (4.16)where r(x) ∈ Rnc.The goal of the IRLS method is to approximate the lp-norm by solving a seriesof locally convex least-squares problems. The general objective function to beminimized takes the form:φ(m) = φd+βφ(k)m= ‖Fm−d‖22+β‖R x(m)‖22 ,(4.17)where the diagonal matrix R ∈ Rnc×nc holds the IRLS weights such that:Rii = ri1/2 . (4.18)At each kth iteration, we seek a minimum along the gradient of the the objectivefunction. Replacing the regularization function from equation 4.6 we get:FTFm+β g(x) = FT dg(x) = RTR x(m)∂x(m)∂m,(4.19)where I explicitly define the gradient of the approximated lp-norm regularizationfunction g(x). I voluntarily neglect the constant of differentiation p from equation(4.15) for two reasons. First, note that in the special case where p= 0, the regular-ization function would vanish and reduce 4.19 to a simple least-squares problem.Secondly, for any p 6= 0, the constant would simply get absorbed by the trade-off parameter β . Other parameters will be introduced in the following section tohandle scaling issues arising from mixed-norm regularization functions.Figure 4.3(b) and (c) presents the IRLS weights r(x) and gradient functionsg(x) for a range of p-value. For p= 2 and x(m) =m, we obtain the smooth regu-55larization function presented in Figure 4.1(a). It is important to note that for a smalllp-norm (i.e p< 1), the IRLS weights and gradients rapidly increase as xi→ ε . Thebehavior of the regularization function around the threshold parameter ε is impor-tant for reasons that will be addressed in Section 4.3.4.I illustrate the various lp-norms on the same 2-variable inverse problem pre-sented in Section 4.1. Figure 4.1(d) and (e) show the regularization and the ob-jective function for p = 1, over a range of model parameters. I point out that,compared to the globally convex l2-norm, the shape of the objective function nowdepends on the starting model m1(0) and m2(0).For p = 0, I recover the compact regularization function put forward by Lastand Kubik (1983), which has been borrowed by many researchers in geophysicsAjo-Franklin et al. (2007); Barbosa and Silva (1994); Stocco et al. (2009). Fig-ure 4.1(f) and (g) presents the approximated l0-norm and corresponding objectivefunction over the same range of model parameters. Similarly, the minimum of theobjective function depends on the initial model m(0) used to construct the regular-ization.Although difficult to prove analytically, experiments have shown that the l0-norm can yield a sparser solution than the convex l1-norm (Chartrand, 2007). I usethe same 2-variable example to illustrate this idea. For simplicity, I impose a fixthreshold parameter (ε = 1e−8) and a fix trade-off parameter (β = 1e−4). Here,I am only interested in the final solution depending on the choice of lp-norm, andas a function of starting model m(0). In the next section, I will provide strategies toefficiently determine those parameters.In the first experiment, I choose the starting model m1(0) to be the least-normsolution at mln= [0.2;0.4], marked as a solid dot (Fig. 4.2). In both cases the l1-norm and l0-norm converge to the same optimal solution at m∗= [0.0;0.5]. Themodel m∗ is interesting as it is sparse with a single non-zero parameter, while alsohaving the smallest norm possible.In the second experiment, the initial model is also chosen to satisfy the targetdata misfit, but this time with a relatively smaller value on the second variablesuch that m2(0)= [0.6;0.2], marked as a white dot. Because both the l1-norm andl0-norm penalize small model parameters, the regularization forces the solution tobe sparse along N (F). Note the clear difference between the globally convex l1-56Figure 4.3: (a) Penalty function ρ(x) for different approximate lp-norm mea-sures, and enlarged view near the region of influence of ε , making the lp-normcontinuous at the origin. (b) IRLS weights r(x) as a function of model functionx(m) and p-values and for a fix threshold parameter (ε = 1e−2). (c) Gradientsg(x) of the model penalty function for various p-values. Note that the gradientsare on a logarithmic scale due to the rapid increase as xi→ ε for p< 1.57norm, converging back to m∗, compared to the non-convex l0-norm, reaching them1-axis at m(k)= [1;0]. The l0-norm has therefore two possible solutions that donot depend on the overall magnitude of m(0), but rather on the relative magnitudeof its components [m1 ; m2]. The l0-norm acts as a binary selector, influenced bythe largest elements of a starting model.The choice of specific norm should therefore reflect the expected character ofthe solution, and the chosen algorithm should allow access to the full range ofnorms for 0 ≤ p ≤ 2. While it is simple to solve (4.17) for a 2-variable problem,finding a solution to large systems for non-convex norms (p< 1) has proven to bedifficult and remains a field of active research. The following sections review po-tential complications for the non-convex cases, and review strategies to efficientlysolve the IRLS method.4.3 IRLS solverThe IRLS method presented in 4.17 is an iterative process, which depends on twokey parameters: the trade-off parameter β and the threshold parameter ε . In thissection, I explore four strategies for the implementation of the IRLS.4.3.1 Basic IRLS algorithmAs illustrated with the previous 2-variable problem, the choice of a starting modelm(0) is especially important for non-convex norms p< 1. According to Chartrand(2007), it may be possible to compute a global minimizer from non-convex func-tions, as long as the initial model is close enough to the global optimum. Mostmethods proposed in the literature seem to agree on the l2-norm solution as a validcandidate (Ajo-Franklin et al., 2007; Portniaguine and Zhdanov, 2002; Sun and Li,2014). Current inversion strategies already rely on the assumption that a smoothmodel is a good approximation of the true solution. Applying a sparse lp-norm canthen segregate the most important model parameters, which in turn can reduce thecomplexity of the solution. This is mostly interesting if the true solution is knownto behave like a delta function.The basic IRLS algorithm becomes a two step process, as summarized in Ta-ble 4.1. During Phase-I, the algorithm finds a smooth solution with the l2-norm58regularization. The trade-off parameter β is monotonically reduced until the modelcan predict the data near the target misfit φ ∗d . The final l2-norm solution provides astarting model m(0) for the calculation of a weighting matrix R.In Phase-II of the IRLS method, model updates are computed iteratively byminimizing 4.17. This process is repeated until the algorithm converges to a stablesolution. I define a convergence criteria as:δφ (k)m =|φ (k)m −φ (k−1)m |φ (k)m×100% , (4.20)where the change in model norm falls below some pre-defined threshold, chosen tobe 2% in all my experiments. In the simplest form of the algorithm, the thresholdε and trade-off parameter β remain constant throughout Phase-II. For now, I willfollow the general consensus that ε has to be small, or near machine error (εp =εq= 1e−8) (Ajo-Franklin et al., 2007; Last and Kubik, 1983; Stocco et al., 2009). Iwill revisit this number in Section 4.3.4. In the method proposed by (Ajo-Franklinet al., 2007), a trade-off parameter β ∗ is chosen so that the initial IRLS solutionpredicts the data near the target misfit φ ∗d . This β∗ then remains constant throughoutthe iterative process.Table 4.1: IRLS Solver 1: Fix parametersPhase I - L2-norm iterations:{m(0) | φd ' φ ∗d}Adjust βmin φ(m)→m(0)Search:{β ∗ | φ (k)d ' φ ∗d}Phase II - Lp-norm iterations:{m(k) | δφ (k)m < 1%}Fix {β ∗,εp,εq}Update r(k)imin φ (k)→m(k)594.3.2 1-D synthetic exampleI demonstrate this simple IRLS implementation on a synthetic 1D problem similarto the problem used by Li and Oldenburg (1996). The model consists of a rectan-gular pulse and a Gaussian function on the interval [0 1] and discretized with 200uniform intervals. Data are generated from the equation:d j =∫ 10f j(x)m(x) dx, j = 0, ...,N , (4.21)where the kernel functions relating the model and the data are defined as:f j(z) = e− j x · cos(2pi jx) , (4.22)where x defines the distance along the x-axis.The model and kernel functions for j ∈ [1,30] are shown in Figure 4.4. Twopercent random noise is added to the data in order to simulate a true geophysicalexperiment. The data are weighted accordingly such that:φd = ‖Wd (F m−d)‖22 , (4.23)where the diagonal matrix Wd holds the estimated uncertainties associated witheach datum. Since the noise is assumed to be Gaussian and uncorrelated, the misfitfunction follows a chi-squared distribution with expected value of N, hence thetarget misfit.I here generalize the objective function presented in Chapter 3 to allow formultiple lp-norm regularization functions such that:φ(m) = φd+β[∫ 10ws∣∣m−mre f ∣∣p dl+∫ 10wx∣∣∣∣∂m∂x∣∣∣∣q dl] . (4.24)The only difference with the regularization function used in 3.6 is in the norms ap-plied on the model and model gradients | · |p and | · |q, where p and q can take anyvalue on the interval [0,2]. Since we are dealing with a 1-D problem, model gra-dients are only computed along the x-axis. I also apply a depth weighting function60Figure 4.4: (a) Synthetic 1D model made up of a rectangular pulse and a Gaus-sian function. (b) Kernel functions consisting of exponentially decaying cosinfunctions of the form f j(x) = e− j x ·cos(2pi jx). (c) Data generated from d= F m, with 5% random Gaussian noise added.to counter the decay of the kernel functions such that:ws = αs e jxwx = αx e jx .(4.25)Combining the IRLS approximation in (4.17) and the discrete objective function61presented in 3.8, the linear 1-D equation 4.24 becomes:φ(m) = ‖Wd (F m−d)‖22+β[‖Ws Rs (m−mref)‖22+‖Wx Rx Gx m‖22],(4.26)where the diagonal matrices Ws,Wx are the cell weights, and Gx is the spatialgradient operator presented in Chapter 3. The IRLS weights Rs and Rx are definedas:Rsii =[(mi(k−1))2+ ε2p](p/2−1)/2Rxii =[(∂m(k−1)i∂x)2+ ε2q](q/2−1)/2,(4.27)where εp,εq are the stabilizing parameters for the sparsity constraint applied on themodel and model gradient respectively.As explained in Section 3.1.1, the minimum norm solution is found where∂φ(m)∂m = 0. Taking the partial derivatives of 4.26 with respect to the model param-eters and setting the reference model to zero yields:(FTWTdWdF+β[RTs WTs Ws Rs+GTx RTx WTxWx Rx Gx])m= FTWTdWdd .(4.28)This linear system can be expressed as an overdetermined problem of the form: Wd F√βWs Rs√βWx Rx Gxm=Wd d00 . (4.29)Solving the least-squares problem 4.29 yields a model update at the kth iterationof the IRLS method. In this case, the left-hand side of 4.28 is linear with respectto m(k) and it forms a symmetric positive definite matrix, which can be solvedefficiently by the Conjugate Gradient descent method.For the initial phase of the IRLS method, a solution is found with the globallyconvex l2-norm regularization, in which case Rs and Rx reduce to the identity ma-trix. Figure 4.5(a) presents the inverted result obtained with the l2-norm, as well62as the convergence curve after achieving the target data misfit φ ∗d . As expectedfrom an l2-norm regularization, the solution is smooth and dispersed over the en-tire model domain.From there, the algorithm proceeds with a sequence of IRLS iterations withl0-norm penalty on the model gradient (αs = 0,q = 0). The goal is to use the reg-ularization to enforce a solution with sharp gradients and hopefully recover therectangular pulse. The IRLS weights are initiated with the smooth l2-norm solu-tion. The initial β ∗ is found by a search method where a solution to equation 4.29is computed for multiple trials. The iterative process is repeated until the inver-sion reaches the convergence criteria specified by equation 4.20, while keeping thetrade-off parameter β ∗ constant.As shown in Figure 4.5(b), the updated solution with the l0-norm penalty re-covers a blocky model with sharp edges. It is important to note that the final dataresidual φ (k)d is much larger than the target misfit (φ∗d = 30). Clearly the influenceof the regularization has overtaken the inversion process. The algorithm likely con-verged to some local minimum, far from the global optimizer previously found withthe smooth l2-norm. Since I know the true solution, I can measure the accuracy ofthe solution, or model error as:‖δm‖1 =nc∑i=1|m(k)i −m∗i | , (4.30)where m(k) is the solution found at the kth iteration. The final model error is largerthan the one found with the smooth l2-norm regularization, hence it is a poor esti-mate of the true solution.This example clearly illustrates some of the challenges related to non-convexobjective functions. We have so far only formulated the basic algorithm behindthe IRLS method, which has been used by many researchers in the past. Impor-tant details regarding the stability and predictability of the algorithm will now beaddressed.63Figure 4.5: (a) Recovered models from the smooth l2-norm regularization, and(bottom) measure of model misfit and model norm as a function of iterations.Both the rectangular pulse and Gaussian function are recovered at the right lo-cation along the x-axis, but the solution is smooth and dispersed over the entiremodel domain. (b) Solution obtained from the IRLS with spasity constraints onthe model gradients (q=0), using the l2-norm model (a) to initiate the IRLS steps.The algorithm uses a fixed threshold parameter (ε = 1e− 8) and fixed trade-offparameter β . The final solution is blocky, as expected from the norm chosen, butfails to achieve the target data misfit. Clearly the influence of the regularizationfunction has overtaken the minimization process.644.3.3 Regularization scalingFollowing Tikhonov’s approach, a solution to the inverse problem is found by pro-gressively reducing the influence of the regularization until reaching the target datamisfit φ ∗d . For strictly convex regularization functions, such as the l2-norm, the in-fluence of the regularization functions is scaled linearly by the trade-off parameterβ . For non-convex functions, in our case for q = 0, the iteration process involvesnon-linear transformations of the objective function, driven by the IRLS weightscomputed in (4.16). Those weights can change rapidly as ε→ 0, which directly im-pacts the influence of the model norm in a non-linear fashion. As demonstrated inFigure 4.5(b), special care must be taken in order to obtain a solution that satisfiesthe data while also being sparse.As a second strategy, I experiment with a brute-force approach where the op-timal trade-off parameter β is determined before each IRLS steps (Table 4.2). Asolution to 4.19 is computed several times for a range of β -values until a suitabletrade-off parameter is found. Figure 4.6(a) presents the model and convergencecurve following this strategy. The solution is blocky, as expected from the l0-normon model gradients, while also honoring the data within 1% of the target misfit.This iterative process is very expensive however, as it requires multiple β solvesper IRLS iteration. For large non-linear problems, this type of approach would becomputationally prohibitive.Table 4.2: IRLS Solver 2: β -searchPhase I - L2-norm iterations:{m(0),β (0) | φd ' φ ∗d}Adjust βmin φ(m)→m(0)Phase II - Lp-norm iterations:{m(k) | δφ (k)m < 1%}Fix {εp,εq}Update r(k)iSearch:{β (k) | φ (k)d ' φ ∗d}min φ (k)→m(k)65As a third strategy, I relax the requirement for the model to predict the dataexactly between each IRLS steps. The trade-off parameter is adjusted once aftereach IRLS step such that:β (k+1) = β (k) ∗ φ∗dφ (k)d. (4.31)I have found experimentally however that this type of posterior update is not suffi-cient to guarantee a smooth convergence. The initial IRLS update can change themodel substantially, moving away from the global l2-norm solution. In order topreserve the relative importance between misfit and regularization functions whilereducing the computational cost, I propose an iterative re-scaling of the model ob-jective function such that:φˆ (k)m = γ(k)(‖Ws Rs m‖22+‖Wx Rx Gx m‖22)γ(k) =φˆ (k−1)m‖Ws Rs m(k−1)‖22+‖Wx Rx Gx m(k−1)‖22,(4.32)where φ (k)m is the scaled model objective function at some kth iteration, and γ(k)is a scalar computed at the beginning of each IRLS iterations. Effectively we aresearching for model parameters that are close to the present solution. Since weare changing the rule by which the size of the model is measured, we attempt toaccount for it. I point out that for constant lp-norm regularizations, the scalingparameter γ(k) is equal to one. Table 4.3 summarizes the proposed method.As shown in Figure 4.6(b), the re-scaling procedure greatly reduces the totalnumber of CG solves. Even though nothing guarantees that non-convex norms willconverge to a global minimum, the re-scaling scheme proposed here forces localminima to be in the proximity of the l2-norm solution. In other words, the scalingparameter γ(k) helps preserve the character of all previous iterations, slowly con-verging to a new minimum, hopefully in the vicinity of the global l2-norm solution.66Figure 4.6: (Top) Recovered models from two different algorithms used to im-plement the IRLS method for q = 0 and a fix threshold parameter (ε = 1e− 8).(Bottom) Measure of model misfit and model norm as a function of iterations.(a) The first algorithm searches for an optimal trade-off parameter β (k) betweeneach IRLS step, requiring a solution to multiple sub-inversions. (b) The secondalgorithm only adjusts β (k) once after each IRLS step. A new scaling parameterγ(k) is added to smooth the transition between the IRLS updates. The algorithmrecovers a similar blocky model but is computationally cheaper, as indicated bythe total number of beta iterations and CG solves.67Table 4.3: IRLS Solver 3: Scaled regularizationPhase I: L2-norm iterations{m(0),β (0) | φd ' φ ∗d}Adjust βmin φ(m)→m(0)Phase II: Lp-norm iterations{m(k) | δφ (k)m < 1%, φd ' φ ∗d}Fix {εp,εq}Update r(k)iScale φˆ (k)mmin φ (k)→m(k)β (k+1) = β (k) ∗ φ∗dφ (k)d4.3.4 Threshold parameter εI have so far delayed providing details regarding the choice of threshold parameterε , which has been the subject of disagreement among researchers (Ajo-Franklinet al., 2007; Barbosa and Silva, 1994; Last and Kubik, 1983; Stocco et al., 2009;Sun and Li, 2014). Little has been said in the literature on how to determine aspecific value of ε .In the classic work of Last and Kubik (1983), and many other research papersafter, it has been suggested that the stabilizing parameter should be small (ε <10−8), or near machine error in order to approximate the lp-norm well. Otherresearchers, such as in Ajo-Franklin et al. (2007) have observed severe instabilitieswith such a small value, and found experimentally that ε should be between 10−4and 10−7. This may be due to the large spread in weights along the diagonalof R impacting the conditioning of the linear system described in equation (4.6).A large condition number can adversely affect the convergence rate of gradientdescent solvers.As an alternative to the regularization function used by Last and Kubik (1983),Gorodnitsky and Rao (1997) apply the IRLS weights directly to the sensitivity68matrix . The general objective function to be minimized takes the form:φ = ‖Fˆmˆ−d‖22+β‖xˆ(m)‖22 , (4.33)where the weighted forward model operator is written as:Fˆ= F diag[m(k−1)] . (4.34)The same technique was later revised by Portniaguine and Zhdanov (2002); Port-niaguine (1999) and coined Minimum Support (MS) functional. The method wasextended to penalties on the model gradients and named Mininum Gradient Sup-port (MGS) functional. This formulation is interesting as it eliminates the stabiliz-ing parameter ε . From a practical standpoint however, I have found issues whenused in concert with other sensitivity based weighting, which will be addressed inSection 4.4.1. Moreover, this type of penalty is less flexible than the general IRLSformulation as I will demonstrate in Section 4.4.It appears that the choice of ε is problem dependent and becomes a compro-mise between achieving the desired level of sparsity, while minimizing the numer-ical cost. Based on the method proposed by Chartrand (2007), I bring in a fourthalgorithm using an ε-cooling strategy. The goal is to progressively change thepenalty function, starting with a coarse approximation that resemble the l2-normfunction. Following the smooth inversion, the threshold ε is initialized as a largevalue (ε(0) ≈ max(x(m)) ∗ 10), then monotonically reduced between each IRLSstep such that:ε(k) =ε(0)2k. (4.35)The optimal threshold parameter ε is found after convergence of the algorithm asδφ (k)m → 0. In a third and final phase, the algorithm fixes ε(k) and continues re-adjusting the trade-off parameter β (k) until reaching the target misfit.Figure 4.7 presents the recovered model after convergence of the algorithm.The solution is blocky, as expected from an l0-norm penalty on the model gradi-ents. This example shows that ε can be much larger than machine error and stillaccomplish the same objective. By progressively reducing ε , changes in regular-ization between each IRLS step are reduced. The data residual remains close to69Table 4.4: IRLS Solver 4: ε-CoolingPhase I: L2-norm iterations{m(0),β (0) | φd ' φ ∗d}Adjust βmin φ(m)→m(0)Phase II: φm-iterations{m(k) | δφ (k)m < 1%}ε(k)p =ε(0)p2k , ε(k)p =ε(0)p2kUpdate r(k)iScale φˆ (k)mmin φ (k)→m(k)β (k+1) = β (k) ∗ φ∗dφ (k)dPhase III: φd-iterations{m(k) | φd ' φ ∗d}Fix {ε(k)p ,ε(k)q }Update r(k)iScale φˆ (k)mmin φ (k)→m(k)β (k+1) = β (k) ∗ φ∗dφ (k)dthe target misfit throughout the iteration process, indicative of a stable algorithm.Numerical experiments have shown that for large inverse problems, the coolingprocedure can make the algorithm substantially cheaper and more stable than witha fix and small ε approach.The method proposed above is advantageous as it does not require a specificchoice of threshold parameters εp and εq. It does depend however on the specifiedconvergence criteria δφ (k)m . In some cases, the algorithm may converge too quicklybefore reaching the desired level of sparsity. Alternatively, the algorithm may be-come overly expensive if the convergence criteria is too restrictive or the solutionoscillates around the minimum. Secondly, the stabilizing parameter ε can be in-terpreted as an e f f ective zero, penalizing specific ranges of model parameters. It70Figure 4.7: (Top) Recovered models from two different algorithms used to im-plement the IRLS method for q = 0 with cooling of the threshold parameter ε .(Bottom) Measure of model misfit and model norm as a function of iterations.Both algorithms adjust β and scaling parameter γ(k) between each IRLS itera-tion. (a) In the first case, the threshold parameter ε is monotonically reduceduntil reaching the convergence criteria. The solution is sparser than the previousalgorithm, even though ε is much larger than machine error. (b) In the secondcase, ε is decreased until reaching the target threshold ε∗. The solution is blocky,penalizing small model gradients.71may therefore be necessary to determine a minimum threshold value in order toguarantee convergence and penalize the right model values.Since the initialization of the IRLS requires a good approximation of the modelvia a least-squares solution, an estimate of the distribution of model parameters isavailable for analysis. The value of an optimal ε∗ can be chosen directly by theuser based on a priori information. Alternatively the choice can be based on thedistribution of model values. Figure 4.8(a) presents the distribution of model andmodel gradients obtained from the smooth l2-norm inversion. Both curves displaya sharp corner around which the model parameters rapidly change. Similarly, Zh-danov and Tolstaya (2004) suggest an L-curve based on the change in model normsuch that:sMS(ε) =mTRTs RsmsMGS(ε) =mTGTxRTxRxGxm ,(4.36)where the model norm sMS and model gradient norm sMGS are computed over arange of ε values as shown in 4.8(b). I found experimentally that both approacheswere valid but did not always yield a well defined corner. More research may beneeded to determine the most robust approach.As a final experiment, I invert the 1-D example using the point of maximumcurvature on the sMGS curve as a minimum threshold for the model gradients (ε∗q =2e−5) . Figure 4.7(b) presents the model after convergence. I note that the resultis more blocky than the one previously obtained with the convergence criteria. Theinversion is, computationally, slightly more expensive due to the larger number ofIRLS steps to reach the target threshold parameter ε∗q .72Figure 4.8: (a) Distribution of model parameters and (b) model gradients recov-ered from the smooth l2-norm regularization. Both curves show a sharp corneraround which the model functions vary rapidly. Similarly, a measure of (c) modelnorm sMS and (d) model gradient norm sMGS can be computed over a range of εvalues, yielding a similar L− curve. The point of maximum curvature can beused to determine the optimal e f f ective zero parameter εp and εq.734.4 Scaled-IRLS method (S-IRLS)The blocky solution found previously is expected from an l0-norm penalty on thegradient. But it is clearly not appropriate to recover the true synthetic model, whichis both smooth and sparse. Building upon the previous section, I explore differentcombinations of norms on the model and model gradients in order to shape thepenalty function. This function should allow any available a priori informationto be incorporated in the solution. My first attempt uses an l0-norm on the modelcombined with and l2-norm on the model gradients, or {p = 0,q = 2}. For thiscombination of norms, I would expect the solution to be both sparse and smooth,better approximating the width of the rectangular pulse and Gaussian anomaly.Figure 4.9(a) shows the solution after convergence of the IRLS algorithm. I hereidentify an important issue with the current IRLS method involving various normmeasures within the same objective function. The solution is clearly dominated bythe sparsity constraint. In this case, the l0-norm penalty on the model suppressesthe l2-norm penalty on the model gradient, yielding a strictly sparse solution with-out smoothness constraint.To understand the issue, it is useful to look at the linear system defined by 4.19.Recall that a solution is found along the gradient of the objective function, whichfrom (4.6) can be written explicitly as:∂φ(m)∂m=(FT d−FT F m)−β( gs(m)+gx(m) )gs(m) = RTs WTs Ws Rs mgx(m) =GTx RTx WTxWx Rx Gx m .(4.37)I divided (4.37) into three parts to highlight the different components of the totalgradient direction. The first term(FT d−FT F m) is related to the misfit functionand depends solely on the ability of the model to reproduce the data d. The secondterm β(gs(m)+ gx(m))is related to the regularization function, which dependson the IRLS weights Rs and Rx computed in 4.16.To illustrate the importance of scaling between the function gradient terms I74Figure 4.9: (a) (Top) Recovered model using a mixed-norm regularization forp= 0, q= 2, ε = 1e−3. (Bottom) Measure of model misfit and model norm as afunction of iterations. The inversion converges to a sparse solution without appar-ent penalty on the gradient, indicative of an imbalance in the regularization. (b)Recovered model and convergence curves after rescaling of the regularization,yielding a model that is both sparse and smooth as expected from the appliedmixed-norm .75form a slightly different 2-variable problem where:F=[0 1], d= 1. (4.38)Figure 4.10(a) displays contours along the surface formed by the l2-norm measureof data residual ‖F m−d‖22, where this time the null-space of F lies parallel to thex-axis. I simplify the problem by having cell size and weights all equal to one sothe objective function can be written as:φ(m) = φd+β [φs+φx]= ‖Fm−d‖22+β[‖m‖22+‖Gx m‖22].(4.39)By setting a small trade-off parameter (β = 1e− 3), I can focus on the solutionspace along the null-space of F marked by the red dash line, as shown in Figure4.10(e). The global minimum of the objective function lies at m = [0.5;1], wherethe partial gradients of ∂φs∂ m1 and∂φx∂ m1 have equal and opposite signs. Here the tworegularization functions use an l2-norm measure, hence the optimal solution is atmid-distance from their respective minimums. This simple example shows that theactual value of the individual norms do not matter, but rather the relative magnitudeof the gradients along the minimum of φd .The behavior of the l0-norm for small p and ε values can easily explain theresult obtain in Figure 4.9(a). Figure 4.11(a) compares the function gradients fordifferent norm penalties as a function of p values. In this case, the function gradi-ents gs measured with the l0-norm are always much larger than the function gradi-ents gx measured with the l2-norm, except at mi = {0,√1− ε2}. The solution cantherefore only be sparse with few non-zero model values required to satisfy φ ∗d .The penalty on model gradients had no influence on the solution.Just as I added the parameter γ(k) to balance the relative importance betweenthe misfit function and the regularization, I also want to control the relative scalingbetween the various regularization functions in terms of their gradients. In otherwords, I want to find a scaling parameter that makes the gradient of any lp-normfunction to intersect, elsewhere than at mi = {0,√1− ε2}.76Figure 4.10: Contour maps for (a) the misfit function φd , (b) the model norm φsand (c) the norm of model gradients φx. (d) The total objective function φ(m)has a global minimum located at m = (0.5,1.0) for a given small trade-off pa-rameter (β = 1e−3). The direction of update is shown for two starting modelsm(0) (black and white dot). (e) Section through the objective function along theminimum of φd . The global minimum occurs where the partial gradients of ∂φs∂ m1and ∂φx∂ m1 have equal and opposite signs.I here introduce a scaling parameter η :η = ε(1−p/2) , (4.40)where once again p denotes the lp-norm penalty and ε is a small value used to ap-proximate the norm. Adding the scaling parameter η to 4.13, the partial gradients77of the objective function become:∂φm∂m=ε(1−p/2)xi(xi2+ ε2)1−p/2 . (4.41)Figure 4.11: (a) Partial gradients of approximated lp-norm penalties for a fixstabilizing parameter ε = 1e−2. Gradients for p< 1 are consistently larger onthe interval [0< xi <√1− ε2], making it hard to combine multiple norm penal-ties within the same objective function. (b) Function gradients after applying ascale of ε(1−p/2), forcing each lp-norm to intersect at m=√ε .Looking at the gradient for xi =√ε , 4.41 simplifies to:∂φm∂m=ε1/2(1+ ε)1−p/2.Then assuming that ε 1:∂φm∂m≈ ε1/2 . (4.42)Hence after applying the scaling η , the gradient of the lp-norm is independent ofthe p-value exactly at m =√ε , and has also a gradient ∂φ∂m =√ε , as shown inFigure 4.3(c). Therefore,√ε is a critical point around which all norm penaltiescan influence the solution. In this manner, two penalty functions with different p-values can coexist within a regularization function and achieve different objectives.78Adding this final scaling to 4.32 yields a Scaled Iterative Re-weighted LeastSquares (S-IRLS) regularization function such that:φˆ (k)m = γ(k)(‖Ws Rˆs m‖22+‖Wx Rˆx Gx m‖22). (4.43)I define a Scaled-IRLS weighting matrix Rˆ such that:Rˆsii =√ηp[(mi(k−1))2+ ε2p](p/2−1)/2Rˆxii =√ηq[(∂m(k−1)i∂x)2+ ε2q](q/2−1)/2ηp = εp(1−p/2)ηq = εq(1−q/2) ,(4.44)where εp and εq are the stabilizing parameters for the lp-norm of the model andmodel gradients respectively. Note that for p = 2, the scaling parameters η isequal to one, hence we recover the the classic l2-norm regularization.Figure 4.9(b) presents the result after a re-scaling of the gradient descent. Boththe l0-norm penalty on the model and the l2-norm penalty on the model gradientsare represented, yielding a solution that is both sparse and smooth. I now have aflexible and robust regularization function. This allows us to explore a wide rangeof solutions, or lp-space of solutions, combining various norms on the model andmodel gradients.4.4.1 Cell-based weights (wr, wm)In Li and Oldenburg (1996), a depth weighting function is added to counter thenatural decay of potential fields. The rapid decay in sensitivity is an importantproblem in geophysics as most data sets are acquired from the surface. The sameidea can be used to incorporate any a priori information regarding the spatial distri-bution of model parameters. In this section, I investigate the effect of having suchcell-based weighting applied directly to the sensitivity matrix, compared to havingit applied to the regularization function as previously formulated in 3.8 .Following the weighted sensitivity formulation of Li and Oldenburg (1996),79the objective function with S-IRLS regularization takes the form:φ(mˆ) = ‖Wd (Fˆ mˆ−d)‖22+ γ(k)[‖Ws Rˆs mˆ‖22+‖Wx Rˆx Gx mˆ‖22], (4.45)where:Fˆ= FW−1rmˆ=Wr mWr = diag(e j z/(2pi)),where the matrixWr hold the inverse exponential describing the decay of the kernelfunctions presented in 4.22. Similar sensitivity based weighting is used in theMinimum Support functional of Portniaguine and Zhdanov (2002).Using both formulations, I invert the same 1-D example presented in Figure4.7and 4.9. In the first experiment, I apply a sparsity constraint on the model gradientsfor q = 0,αs = 0. Figure 4.12(a) compares the true and recovered models afterconvergence. I observe the solution obtained with the φ(mˆ) formulation is skewedin the direction of the sensitivity weighting. Similarly, Figure 4.12(b) shows therecovered model for the combination of p= 0, q= 2. While the recovered anomalyover the rectangular pulse is both smooth and sparse, I note that the solution tendsto favor a strictly sparse model further at depth over the Gaussian function. On theother end, the solution obtained with the φ(m) formulation remains smooth andsparse over both anomalies equally.Issues encountered with the φ(mˆ) formulation are due to the IRLS weightscomputed in 4.27. Written explicitly in terms of weighted model parameters, thelinearized IRLS norm from 4.14 can be written as:φ (k)mˆ =nc∑i=1(wixi)2[(wix(k−1)i )2+ ε2]1−p/2 . (4.46)80Factoring out the cell-based weight wi, equation 4.46 becomes:φ (k)mˆ =nc∑i=1wpi x2i[(x(k−1)i )2+ ε2/w2i]1−p/2 . (4.47)The important point to note is that the weights are now function of the p-value,while the threshold parameter ε depends on the cell-based weights. In most cases,the sensitivity weights increase with distance from the observation location, whichin turn reduce the threshold parameter ε and increase the influence of sparse norms.This explains the increase in sparsity constraint with distance as observed in Fig-ure 4.12. In order to preserve the character of cell-based weights, it is importantto isolate their effect from the lp-norm penalty. The same idea extends to any typeof cell-based weights (i.e. volumetric, a priori). In the following chapter, variousnorm measures will be used simultaneously on the model and model gradient. Thepredictability and stability of the algorithm depends greatly on the scaling betweenthe various components of the objective function.Figure 4.12: Recovered models for two different depth weighting formulations:(red) weighted sensitivity φ(mˆ), (black) weighted regularization φ(m). (a) Trueand recovered models using the φ(mˆ) and φ(m) formulations for penalty ap-plied on the model gradients for q = 0 and (b) for p = 0, q = 2. The weightedsensitivity formulation φ(mˆ) increases the influence the regularization functionwith distance along the x-axis, skewing the model towards the right.814.5 Mixed lp-norm regularizationI showcase the robustness and flexibility of the S-IRLS algorithm by inverting atotal of 441 models, corresponding to all the combinations of norms on the interval0 ≤ p ≤ 2 and 0 ≤ q ≤ 2, on a 0.1 increment. Figure 4.13(a) presents a contourmap of the final model errors ‖δm‖1 recovered after each inversion. The gradualvariation in model error is a good indicator of the predictability of the algorithm, assmall changes in the regularization yield equivalently small changes in the solution.For each inversion, the trade-off parameter β was adjusted to fit within±2% of thetarget data misfit for appropriate comparison of the results, as shown in 4.13(b). Inthis case, the optimal model is found with a combined norm of [p= 1.5, q= 0.4].In this particular case, the worst recovery was obtained with p= 0, q= 2, yieldinga sparse solution with few large oscillations. Even though this type of analysiswould be computationally prohibitive for large scale inverse problems, it illustratesthe flexibility of the S-IRLS method in designing an objective function reflectingspecific characteristics.Nine of the inverted models are shown in Figure 4.14 for comparison. Asexpected from the property of the IRLS functions, small norms on the gradientpromote blocky solutions. From right to left, the inversion result tends to favorright-angled anomalies. From bottom to top, smaller norms on the model parameterenforce sparse solutions with few non-zero values.4.5.1 Localized S-IRLSWhile I have so far implemented global regularization functions applied on theentire model, the S-IRLS method can easily be extended to problems made up ofseveral sub-regions with varying norm penalties. Also based on the IRLS method,Sun and Li (2014) introduce an automated process to extract structural informationfrom the data and impose either an l1 or l2-norm penalties to specific regions of themodel. While designing this objective function can be challenging and may stillrequire direct input from the user, their approach has great value as it increases theflexibility over traditional global penalty functions. The S-IRLS formulation hasthe potential to further generalize the work of Sun and Li (2014) in allowing theuse of any norm penalty in the range of 0≤ p≤ 2.82Figure 4.13: (a) Model error ‖m−m∗‖1 and (b) misfit function for the 441 in-verted models using a range of regularization with mixed-norm penalty on themodel for 0≤ p≤ 2 and on model gradients for 0≤ q≤ 2. (c) The largest modelerror (‖δm‖1) was obtained with the mixed-norm for p = 0, q = 2, comparedto (d) the optimal solution found with p= 1.5 and q= 0.4.83Figure 4.14: (a) Nine of the 441 inverted models for a range of mixed-normpenalties on the model and its gradient for 0≤ p≤ 2 and 0≤ q≤ 2.Suppose an N-dimensional model space m ∈Ω, divided in J sub-regions suchthat:S= {S1,S2, ...,SJ∣∣∣ S j ⊂Ω , j = 1, ...,J}Si∩S j = { /0∣∣∣ i, j = 1, ...,J , j 6= i} .84I can define a model of lp-norms by the union of each sub-domain:p=J⋃j=1S j p j , (4.48)where the scaler p j defines the lp-norm to be applied to the jth region of the modeldomain. Each sub-domain S j can have a distinct value for p and q, penalizingthe model values and model gradient differently than neighboring domains. Thetransition in lp-norm penalty can be smoothed after applying a linear transitionacross neighboring domains as defined by:p˜= A pA=1/2 1/2 0 . . . 0.... . . 00 . . . 0 1/2 1/2 , (4.49)where the transition between lp-norm is extended over several cells by applying theaveraging operator A. Allowing for multiple sub-domains requires only a slightchange to the S-IRLS method presented in Section 4.4. The general objectivefunction for a 3-D discretized problem becomes:φ(m) = φd+βγ(k)[‖Ws Rˆs (m−mref)‖22+ ∑i=x,y,z‖Wi Rˆi Gi m‖22], (4.50)such that the S-IRLS weights for the model and model gradients become:Rˆsii =√η˜pi[(mi(k−1))2+ ε2p](p˜i/2−1)/2Rˆxii =√η˜qxi[(∂m(k−1)i∂x)2+ ε2q](q˜xi/2−1)/2Rˆyii =√η˜qyi[(∂m(k−1)i∂y)2+ ε2q](q˜yi/2−1)/2Rˆzii =√η˜qzi[(∂m(k−1)i∂ z)2+ ε2q](q˜zi/2−1)/2.(4.51)85Figure 4.15: (Left) Improved solution for the 1-D problem after applying alocalized mixed-norm penalty, where the regularization is divided into tworegions with independent lp-norm regularization: (left) p = q = 0, (right)p = 1,q = 2. (Right) Convergence curves for the mixed-norm S-IRLS inver-sion.Figure 4.15 presents an optimal solution for the 1D problem, where I divide themodel space into two regions. The left half uses a p = 0,q = 0, whereas the righthalf imposes a p = 1 and q = 2. The mixed-norm penalty function enforces theright characteristic to each portion of the model. The inversion recovers a sparseand blocky model over the rectangular pulse while at the same time managing tomodel the smooth Gaussian anomaly. Consequently, the model error ‖δm‖ is muchsmaller than any of the previous attempts. The algorithm converges smoothly andrapidly to a stable solution.864.5.2 2-D exampleI test the S-IRLS algorithm on a 2-D synthetic model to further demonstrate theflexibility and robustness of the mixed-norm regularization. Analogous to the 1-Dproblem, the 2-D model consists of a square block and a smooth Gaussian anomaly,only this time the observations stations are located on both sides of the modeldomain. The kernel functions are exponentially decaying cosine functions of theform:e−ω r · sin(2piωr) , (4.52)as shown in Figure 4.16.Since we are now dealing with a 2-D problem, the regularization function in-volves a measure of model gradients in the x and z-direction. Most methods pro-posed in the literature, and presented in 3.9, involve computing gradients with finitedifference operators in orthogonal directions such that:∂mi∂x=m(i, j)−m(i−1, j)dx(i)∂mi∂ z=m(i, j)−m(i, j−1)dz( j),(4.53)where the indexes (i, j) represent the cells ordering in the x and y-direction respec-tively. A second option is to penalize the absolute model gradients such that:∇mi =[(∂m(k−1)i∂x)2+(∂m(k−1)i∂ z)2]1/2. (4.54)The IRLS weights then become:Rˆxii =√η˜qxi[(∇mi)2+ ε2q](q˜xi/2−1)/2Rˆzii =√η˜qzi[(∇mi)2+ ε2q](q˜zi/2−1)/2.(4.55)The same idea can easily be extended to three dimensions. I illustrate the dif-ference between both measures of model gradients on this 2-D problem. Figure4.17(a) presents the recovered model using an l1-norm penalty on the model gradi-87Figure 4.16: (a) Synthetic 2-D model made up of a square block and a smoothGaussian function. (b) Example of a kernel function for e−ω r · cos(2piωr) and(c) data generated from d= F m. Five percent random Gaussian noise is added.88ents from the finite difference formulation. This yields blocky right-angled anoma-lies. Alternatively, penalizing the absolute model gradients from (4.54) allows forsmooth corners as shown in 4.17(b). Both methods are valuable as they promotedifferent solutions, which in some cases may better reflect the geometry of the truemodel.To accommodate both the blocky and the smooth anomaly, I divide the modelspace into three zones. Figure 4.18 presents the different lp-norm zones definedover the two anomalies. I extend the procedure of sub-regions presented in Sec-tion 4.5.1 to this 2-D problem. I impose sparse model values with the l1-normover the background region in order to get a simple model. The norm on modelgradient is variable and reflects the expected characteristics of the true model. Ipurposefully chose regions larger than the extent of the known anomalies in orderto simulate a choice that could be made for a blind inversion. The goal is also tosee if the transition from an l2-norm to a l0-norm can be done smoothly withoutcreating artifacts. Figure 4.19(a) and (b) presents the recovered model after Phasel and Phase lll of the S-IRLS algorithm previously introduced in Table 4.4. Bothmodels fit the data within 2% of the target misfit φ ∗d The final mixed-norm modelnicely recovers the blocky and smooth Gaussian anomaly. The transitions betweenthe different lp-regions appear to be seamless.89Figure 4.17: (a) Recovered model for p = qx = qz = 1 penalizing finite differ-ence gradients in orthogonal directions, yielding right-angled anomalies. (b)Recovered model for the same norms but penalizing the absolute gradient ofthe model (|∇m|) recovering round edges.Figure 4.18: Distribution of lp-norm on the model and model gradients over the2-D model domain. The original boundary of each region was smoothed inorder to get a slow transition and reduce visible artifacts. Regions were chosento cover a larger area than the anomalies to simulate a blind-inversion.90Figure 4.19: (a) Smooth l2-norm solution used to initiate the IRLS iterations.(b) Recovered model using the mixed-norm regularization after seven S-IRLSiterations. The contour line (red) marks the value of εp, judged to be the effectivezero value of the model (mi ≤ 5e-2). Both models (a) and (b) fit the data within2% of the target misfit φ ∗d . (c) Dual plots showing the distribution of modelparameters and the gradient of the l0-norm penalty function gˆp(m) as a functionof S-IRLS iterations. High penalties are applied to progressively smaller modelvalues. The final model nicely recovers both the blocky and smooth Gaussiananomaly.914.5.3 3-D exampleAs a final test, I invert the same synthetic 3-D susceptibility example as presented inSection 3.2. The model consists of a folded anomaly with magnetic susceptibilityof 0.075 SI, arching around a discrete block with susceptibility of 0.05 SI. Thearc-shaped anomaly is dipping 20◦ towards the south.I impose sparsity constraints on the model gradients such that (p = 0, qx =2, qy = 2, qz = 2). Figure 4.20 shows the recovered model after five IRLS iter-ations. The solution is remarkably closer to the true solution, both in shape andmodel values compared to the smooth l2-norm solution. The inversion also hasan easier time reproducing the high frequency content from the data as shown inFigure 4.21. The inversion successfully recovers two distinct objects, no longerconnected at depth. Notice that both arms are well defined and extend along theentire length of the arc.As a second experiment, I impose a constraint on the model gradients in orderto recover a blocky model. Model gradients are measured with the standard finitedifference operators, expected to yield right-angled anomalies. Figures 4.22 and4.23 show the recovered model and predicted data for (p = 0, q = 1), yielding asparse and blocky solution. Once gain, compared to the smooth l2-norm solution,the correlated data residuals are substantially reduced. The recovered susceptibilityvalues within the block and arc shape are more uniform than with the smoothnessconstraint, ideal to get a bulk estimate in physical property.92Figure 4.20: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredsusceptibility model after five IRLS iterations for (p = 0, q = 2). The finalmodel is substantially closer to the true solution.Figure 4.21: Comparison between (a) observed and (b) predicted data from therecovered susceptibility model using compact norms for (p = 0, q = 2). (c)Normalized data residuals are within two standard deviations.93Figure 4.22: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredsusceptibility model after nine IRLS iterations for (p = 0, q = 1) . Sparsityconstraints on the model and model gradients yield a simple and blocky model.Figure 4.23: Comparison between (a) observed and (b) predicted data from therecovered susceptibility model using compact norms for (p = 0, q = 1). (c)Normalized data residuals are within two standard deviations.944.6 Case study - Tli Kwi Cho kimberlite complexAs a final test, I apply the mixed-norm S-IRLS algorithm on a magnetic data fromthe Tli Kwi Cho (TKC) diamondiferous kimberlite complex. The property is lo-cated in the Lac de Gras region, approximately 350 km northeast of Yellowknife,Northwest Territories, Canada. The TKC deposit was originally discovered by anairborne DIGHEM survey in 1992, including frequency-domain EM and magneticdata. Two kimberlite pipes, dubbed DO27 and DO18, were later confirmed bydrilling in 1993. Different volcaniclastic rock units suggests that the deposit wasformed over several events (Doyle et al., 1999). The pipes are intruding into oldergranitic rocks of the Archean Slave Province, known to have little to no magnetitecontent. A magnetic susceptibility contrast between the kimberlite and the countryrocks give rise to noticeable anomalies on the observed TMA data. Strong linearfeatures associated with the Mackenzie dyke swarms can be observed on the east-ern part of the survey. These dykes extend over the entire Lac de Gras region andgenerally strike 305◦N (Wilkinson et al., 2001). They are known to be near verticaland are typically between 20-50 m wide. This geological setting is therefore idealfor implementing a mixed-norm algorithm. My goal is to apply soft constraints onthe model values and model gradients in order to favor specific geometry.Prior to the inversion, I rotate the regional dataset by 30◦ counterclockwise toa local coordinate system in order to align the principal axis of the dykes parallelto the grid. The declination of the inducing field are also rotated by 30◦ to preservethe geometry of the problem. The region of interest is discretized into 25 m cubecells. In preparation for Phase II of the S-IRLS algorithm, a smooth solution isfound with the l2-norm regularization as shown in Figure 4.24(a). The inversionrecovers two discrete bodies corresponding to the known location of DO18 andDO27. Two linear anomalies striking north-south correspond with the Mackenziedyke swarms. A third anomaly intersects at right-angle running 90◦N. As expectedfrom the l2-norm regularization, the solution is smooth and edges are not clearlydefined. Note also that the highest recovered susceptibilities are strongly correlatedwith the observation locations. Dykes appear to break up between each survey line,which is clearly an inversion artifact due to changes in sensitivities.I want to modify the objective function in order to recover sharp and continuous95edges. In general terms, I want to enforce large gradients perpendicular to thestrike of the dykes, while imposing a sparsity constraint over the kimberlite pipes.I divide the region of interest into four sub-regions as shown in Figure 4.25:• S1: Sub-region covering the Mackenzie dyke swarms on the eastern edge ofthe survey. I impose an l0-norm on the model gradients in the xˆ-direction(qx = 0) to get sharp edges along the strike of the dykes, but smooth in theother directions (qy,qz, p= 2).• S2: Sub-region covering the dyke running perpendicular to the Mackenziedyke swarms. Similar to S1, I impose an l0-norm on the model gradients inthe yˆ-direction (qy = 0).• S3: Sub-region covering the kimberlite pipes DO27 and DO18. I imposean l1-norm on the model gradients (qx,qy,qz = 1) and an l0-norm on modelvalue (p= 0) in order to recover blocky and sparse anomalies.• S4: Background region covering the rest of the survey area. I impose an l2-norm on the model gradients (qx,qy,qz = 2) and an l1-norm on model value(p= 1) in order to recover a smooth and sparse model.Figure 4.24(b) presents the recovered model obtained after five iterations ofthe S-IRLS algorithm. In comparison with the smooth inversion, the edges of allthree dykes are better defined and extend vertically at depth. The inversion maysuggest the presence of a fourth dyke parallel to the Mackenzie dyke swarm, onthe eastern edge of the survey. Both DO18 and DO27 are recovered as compactbodies with higher susceptibility and well defined edges compared to the smoothl2-norm inversion. Because of the small l0-norm applied to model values, mostof the near surface anomalies are removed, yielding a more compact anomaly atdepth. This is also true away from the kimberlite pipe over the background region,where a simpler model is obtained. Predicted data for the smooth l2-norm modeland the final S-IRLS solution are presented in Figure 4.26, both within 5% of thetarget misfit.96Figure 4.24: (a) (Left) Horizontal section through the recovered susceptibilitymodel at 25 m depth below topography from the smooth l2-norm regulariza-tion. (Right) Iso-surface of susceptibility values around 0.002 SI. (b) Recov-ered model using the mixed-norm S-IRLS algorithm. Magnetic dykes are betterrecovered, imaged as continuous plates and extending vertically at depth. Sus-ceptibility values for DO-27 and DO-18 have increased, showing as compactvertical pipes.974.7 SummaryIn this chapter, I reviewed and experimented with several algorithms for the ap-proximation of lp-norm measures on the model and model gradients to promotesparse and blocky solutions. I proposed a Scaled-IRLS algorithm in order to stabi-lize current inversion strategies, while reducing the computational cost. My algo-rithm offers a strategy to determine an optimal value for the stabilizing parameterε , previously overlooked by other researchers. This robust method is further gen-eralized, allowing for mixed-norm penalties on the model and model gradients.Scalings are applied to the individual terms of the regularization function in orderto preserve their relative importance during the inverse process. A solution to thenon-linear inverse problem is found iteratively using Gauss-Newton steps.Varying the combination of lp-norm of the model and its gradients allows meto shape the objective function to exhibit specific characteristics. For differentgeological settings, there is likely a specific combination of norms that can bet-ter recover the true geometry of causative bodies. Preliminary results show greatpromise at reducing the non-uniqueness of current potential fields inversion codes.Because my formulation is general, the algorithm can be applied to a wide rangeof inverse problems.98Figure 4.25: Horizontal section through the mixed-norm models applied to foursub-regions with smooth transition across zones.99Figure 4.26: (a) Observed and predicted data over the TKC kimberlite complex.(b) Residuals between observed and predicted data normalized by the estimateduncertainties (10 nT). Both the smooth and mixed-norm inversions reproducethe data within four standard deviations.100Chapter 5Cooperative Magnetic Inversion(CMI)In Chapter 3, I reviewed three magnetic inversion algorithms introduced in the liter-ature. The magnetic susceptibility inversion is robust and computationally cheap,but it runs into limitations whenever the direction of magnetization is not paral-lel to the inducing field. The more general Magnetic Vector Inversion (MVI) canhandle any orientation of magnetization by solving directly for the magnetizationvector. This comes at the cost of increasing the number of unknowns in our largeunder-determined system of equations. The l2-norm regularization yields smoothsolutions that are often too complex for direct geological interpretation. Lastly, theinversion of magnetic amplitude data can provide a good estimate of the locationand strength of magnetized bodies. The method is less sensitive to the orientationof magnetization than an inversion carried out with an assumed direction of magne-tization. The main complication for amplitude inversion is in computing amplitudedata from the observed TMI, either via Fourier transform or by the equivalent-source. The inverse problem is also non-linear with respect to the model, whichincreases the computational cost and is less stable.While all of the above inversion methods have strength and weaknesses, theyeach bring complementary information. The susceptibility inversion can be usedto rapidly compute an equivalent-source layer and to generate amplitude data. Theamplitude data can then be fed into the amplitude inversion to locate magnetic101anomalies. In turn, imposing a geometric constraint on the MVI can greatly re-duce the non-uniqueness of the problem. All three algorithms share variations ofthe same system of equations. This is beneficial as calcuations for the sensitivitymatrix can be stored and repurposed at each step.In this chapter, I introduce a Cooperative Magnetic Inversion (CMI) method,which combines the equivalent-source, amplitude inversion and MVI method intoa single inversion algorithm. Using the S-IRLS method from Chapter 4, I impose asparsity constraint on the model and model gradients to further reduce the solutionspace. From an imaging stand point, compact and well defined magnetic anoma-lies can simplify the geological interpretation. From a practical aspect, having analgorithm that can manage all three codes reduces the number of manual steps thatwere previously required from the user.5.1 MethodologyIn Chapter 3, I inverted a synthetic model with complicated magnetization ori-entations. The model recovered from the MVI inversion was smooth with poorrecovery of the arc anomaly (Fig. 3.10). Using the CMI algorithm, I re-invert thesynthetic example using an l2-norm regularization. I will here breakdown each stepof the CMI algorithm, but I want to emphasize that the entire process is automatedand can all be done as a single workflow. I divide the algorithm in three parts asshown in Figure 5.1 in schematic form.Stage I of the algorithm is mainly a data preparation step, ahead of the am-plitude inversion and MVI. As in other inversion methods, the algorithm requiresTMA data and a mesh. The algorithm computes the large system of equationsneeded for sensitivity calculations as defined by 2.10. Following the sensitivitycalculations, the algorithm proceeds with an equivalent-source layer inversion asdescribed in Section3.5. Figure 5.2 presents the inverted effective susceptibilitylayer and predicted data at the station locations. Amplitude data are saved andpassed on to Stage II.In Stage II, the algorithm proceeds with the amplitude inversion from Sec-tion 3.4. Figure 5.3 shows the recovered smooth effective susceptibility model. Therecovered effective susceptibility model is then converted to sensitivity weighting102Figure 5.1: Schematic representation of the Cooperative Magnetic Inversion(CMI) algorithm. Input and output parameters are indicated by a dash arrow.matrix such that:Wmii =[ κei ∗0.9max(κe)+0.01]−1, (5.1)where Wm ∈Rnc×nc is a diagonal matrix of the normalized effective susceptibilitiesκe. A small number (1e-2) is added to assure that all values are between [1 100].This type of re-scaling was determined empirically to be robust. The sensitivityweighting matrix is saved and passed on to the MVI inversion.In Stage III, the regularization function is pre-multiplied by the effective sus-ceptibility weighting calculated in 5.1. The model objective function becomes:φ(m) = φd+βγ(k)[‖Wm Ws Rˆs (m−mref)‖22+ ∑i=x,y,z‖Wm Wi Rˆi Gi m‖22].(5.2)The model weighting matrix Wm imposes a high penalty on cells that received loweffective susceptibility from the amplitude inversion. Figure 5.5 presents the in-verted model after reaching the target data misfit. The added information from theamplitude inversion greatly improves the solution over the MVI alone. Althoughstill smoothly varying, the inversion manages to separate the arc and the blockanomaly, while also recovering the orientation of magnetization.103Figure 5.2: (a) Inverted equivalent-source layer and (b-f) predicted xˆ, yˆ, zˆ-component, TMI and magnetic amplitude data for the synthetic model.104Figure 5.3: (a) Iso-surface (0.002 SI) and (b) sections through the recoveredeffective susceptibility model. This effective susceptibility model is used to con-struct a weighting matrix to constrain the MVI.Figure 5.4: Comparison between (a) observed and (b) predicted data from therecovered effective susceptibility model. The inversion can predict most of thedata within one standard deviation.105Figure 5.5: (a) Iso-surface (0.005 SI) and (b) sections through the recoveredmagnetization model from the CMI algorithm (p = 2, q = 2) . The inversionrecovers both the arc and the block anomaly as distinct objects.Figure 5.6: Comparison between (a) observed and (b) predicted data from therecovered susceptibility model. (c) Normalized data residuals are within twostandard deviations.1065.2 CMI with S-IRLS regularizationAs a final experiment I impose sparsity constraints on the amplitude inversion for(p= 0, q= 1), in order to simplify the distribution of effective susceptibility. Thegoal is to get a simpler effective susceptibility model to further reduce the activespace of the MVI algorithm. It is a so f t constraint on the solution as I do notimpose a priori information to individual cells, but rather a general constraint onthe behavior of the solution. Figure 5.7 presents the recovered model after conver-gence.The lp-norm constraint considerably reduces the complexity of the model, al-though some artifacts remain at depth from the amplitude inversion. In Stage III,the compact effective susceptibility model is used to constrain the MVI result. Fig-ure 5.9 presents the final magnetization vector. The arc and the central block arerecovered at the right location and effective susceptibility values are near the truemodel. From a mineral exploration perspective, the model showed in 5.9 would beeasily interpreted. I note that most of the deep artifacts from the amplitude inver-sion have been removed. The final data residual is within one standard deviation,and no longer shows structure correlated with the magnetic anomaly.5.3 ConclusionIn this chapter, I implemented a Cooperative Magnetic Inversion for the inversionof magnetization vectors. The algorithm ties together three inversion algorithmspreviously introduced in the literature: the magnetic susceptibility, magnetic am-plitude and MVI. The recovered vector magnetization models are simpler and morerepresentative of the true solution. The inversion process is automated, hence re-ducing the overall work required, both numerically and from the end-user.107Figure 5.7: (a) Iso-surface (0.01 SI) and (b) sections through the recovered ef-fective susceptibility model from the amplitude inversion with sparsity constraintapplied (p = 0, q = 1). The lp-norm constraint considerably reduces the com-plexity of the model, although the model is still stretched vertically.Figure 5.8: Comparison between (a) observed and (b) predicted amplitude datafrom the recovered compact effective susceptibility model (p = 0, q = 1). (c)Normalized data residuals are within two standard deviations.108Figure 5.9: (a) Iso-surface (0.01 SI) and (b) sections through the recovered mag-netization model from the CMI algorithm. Compact norms (p= 0, q= 2) wereapplied during the amplitude inversion. The lp-norm constraint considerably re-duces the complexity of the model.Figure 5.10: Comparison between (a) observed and (b) predicted data from therecovered magnetization CMI model. (c) Normalized data residuals are withintwo standard deviations. Correlated residuals are no longer seen.109Chapter 6Case StudyIn this chapter, I apply the CMI algorithm on a 1992 aeromagnetic data set col-lected over the Ekati Property, Northwest Territories. Airborne magnetic surveyshave been an integral part of diamond exploration in the region since the early1990s (Pell, 1997). The Lac de Gras region has been particularly productive, host-ing two of the largest deposits found in Canada, the Ekati and Diavik mines. Over150 kimberlites have been identified on the Ekati Property (Carlson et al., 2015),but diamond grades are highly variable. It is estimated that less than 10% of the150 known kimberlites on the Ekati Property are of economic interest. While mostkimberlite pipes can easily be identified by their geophysical signature, estimatingthe economic potential early in the exploration stage remains challenging (Coop-ersmith et al., 2006).Kimberlite pipes are generally magnetic anomalies known to retain a strong re-manent component. There have been several studies specifically dedicated to mag-netic methods applied to kimberlite deposits. Previous geophysical studies reliedprimarily on data processing techniques to isolate exploration targets (Cowan et al.,2000). Direct inversion of magnetic data over remanently magnetized kimberlitepipes has long been challenging for reasons discussed in Chapter 3. In Shearer(2005), the magnetic amplitude inversion is implemented on a data set from theGalaxie trend, Nunavut. The location and geometry of negative magnetic anoma-lies were successfully inverted and attributed to hypabyssal kimberlite dykes, laterconfirmed by drilling. A few years later, Zhao (2012) inverted a subset of a large-110scale aeromagnetic survey over the Ekati Property. The inclination and amplitudeof magnetization were computed for the Grizzly and Leslie pipes and compared tothe 2D analysis done by Cheman (2006).In this chapter, I invert the large data set that was partially analyzed by Zhao(2012), over the Ekati Property. Using the CMI algorithm introduced in Chapter5, I extract regional information about the distribution of dyke swarms in relationto the known kimberlite deposits. I then proceed with deposit-scale inversionsover 17 of the known pipes and recover a bulk estimate of the magnetization. Therecovered orientations of magnetization are compared to physical property mea-surements from the literature. Finally, I compare the published radiometric datingon 11 of those pipes with the polarity predicted from the inversion and the globalgeomagnetic timescale.6.1 Ekati PropertyThe Ekati Property is located north of the Lac de Gras region, approximately 300km northeast of Yellowknife, Northwest Territories. The current mineral claimcovers an area of over 262,000 ha. Since discovery of the first kimberlite in 1991,over 150 additional pipes were identified from geophysical surveys, till samplingand drilling programs. Figure 6.1 shows the topography and locations of knownpipes over the region. The Panda and Beartooth pipes were actively mined betweenthe period of 1998 and 2010 and are now fully depleted. Three are still underproduction: Misery, Koala and Koala North. Since the beginning of operations in1999, over 58 million carats have been extracted from the various kimberlite pipes(Carlson et al., 2015). South of Ekati is the Diavik mine, which has also producedover 91 million carats between 2003 and 2014 (Yip and Thompson, 2015).6.1.1 GeologyThe Lac de Gras kimberlite field is located in the Slave craton, which consistsmainly of granitoid and meta-sediment rocks dated between 2.66-2.58 Ga-old (Creaseret al., 2004). The region is intruded by several phases of diabase dyke swarmsclearly visible on the airborne magnetic data. Table 6.1 summarizes the age andtrend of five important groups of dykes (Buchan et al., 2009). It has been argued111Figure 6.1: Topography and known kimberlites over the Lac de Gras region.The location of the main operations for the Ekati and Diavik mines are shown(red dot).that the dykes may serve as structural control for kimberlites (Wright, 1999).Kimberlite pipes in the Lac de Gras region are narrow, steeply sided invertedcones, with the majority of them covering less than 5 ha at the surface. The standardmodel for kimberlites in the Lac de Gras region is a mix of fragmented craterfacies and non-fragmented hypabyssal (HK) facies. The crater facies in the upperportion of the type can be further subdivided in various types of volcaniclastic(VK) and pyroclastic (PK) units (Pell, 1997). Found at greater depths, the HK unitis a coherent olivine-rich rock, and is also found as intrusive dykes and sills. Thecomposition and physical properties of kimberlites can vary greatly between pipes.Geophysically, it has been found that crater facies are generally associated withdensity and magnetic susceptibility lows. Weathering of the crater facies may also112Table 6.1: Dyke swarms of the Lac de Gras region listed in increasing age aspublished by Buchan et al. (2009)Name Orientation (◦) Age (Ma)Mackenzie 315 1267′305′ 305 -Lac de Gras 0 2027-2023MacKay 90 2210Malley 45 2230appear as resistivity lows due to the abundance of clay minerals, although they caneasily be confounded with fine grained glaciofluvial sediments in northern regions(Power and Hildes, 2007). The hypabyssal facies on the other hand are mostlyassociated with magnetic susceptibility highs and are prone to retaining remanentmagnetization.Fossil bearing xenoliths in some pipes are indicative of a sub-marine emplace-ment between the Late Cretaceous to Eocene time. Radiometric dating of 36 pipesindicates a wide range of emplacement times (Creaser et al., 2004), divided intofive temporally discrete episodes (Lockhart et al., 2004). Three of the most pro-ductive kimberlites at Ekati (Koala, Koala North and Panda pipes) were emplaced53 Ma ago during the Panda Age Array (PAA). Physical property measurementsof the rocks were also acquired over the Lac de Gras region and are summarizedin (Buchan et al., 2009; Enkin, 2003). It has been suggested that the orientationof magnetization could be used to determine the age of a kimberlite pipe, henceit could be used to estimate the economic potential of a deposit (Lockhart et al.,2004).6.1.2 Airborne magnetic dataIn 1993, a fixed-wing magnetic survey was flown over the Ekati Property as shownin Figure 6.2. The ′′Paul′s Lake′′ aeromagnetic data were downloaded from theNatural Resources Canada’s website. The data set consists of 91 survey lines ofTotal Magnetic Field (TMI) measurements, spaced 250 m apart, covering an areaof 23× 36 km. Observations were provided at a 50 m spacing along line, at an113average flight height of 120 m, as measured by radar altimeter. From historicalgeomagnetic data, the inducing field strength was 60,275 nT, [I: 27◦, D: 84◦].A Digital Elevation Model (DEM) was downloaded from the Canadian DigitalElevation Data collection (1:50,000), also provided by Natural Resources Canada.The DEM was used to convert the radar height to absolute elevation referenced tomean sea level.Figure 6.2: Topography, aeromagnetic survey and known kimberlites over theLac de Gras region. The location of the main operations for the Ekati and Diavikmines is marked with a red dot.6.2 Regional scale inversionThe inversion procedure follows the methodology presented in Chapter 5. The ob-jective of a regional scale inversion is to first gain knowledge about the general114distribution of magnetic anomalies. Moreover, this step is used to estimate uncer-tainties and to identify regions of interest for a deposit-scale inversion.Prior to the inversion, I remove a regional trend in the data in order to focusthe inversion on local anomalies. Small regions 1 km2 showing low variations inthe field were manually chosen as shown in Figure 6.3(a). Within each region,the median data value and horizontal location are recorded. Those values are thenused to create a first-order polynomial (Fig-6.3(b)). A regional trend of roughly250 nT trending up towards east is subtracted as shown in Figure 6.3(c). I assignfive percent uncertainties, plus 5 nT floor uncertainty to the data.Figure 6.3: Regional data (a) before and (c) after removal of a regional trend.The median value within selected regions (boxes) were used to compute (b) afirst-order polynomial trend.Table 6.2 presents some of the parameters used for the inversion. The modelspace is discretized by a 50 m cell size resolution. The depth extent of the meshcovers the top 1 km below topography, for a base mesh of 6 M cells. Inversion ofthe entire data set, approximately 50,000 observations, at this resolution remains115Table 6.2: Regional-scale inversion parameters.Core cell size 75×75×50 mBase mesh size 3.3 M cellsNumber of tiles 77M / tile ≈ 100k cellsN / tile ≈ 900 datap, qx, qy, qz 1, 1, 1, 1αs, αx, αy, αz 1.8e-4, 1, 1, 1Uncertainties 0.05*d + 5 nTcomputationally challenging for most computers. In order to reduce the compu-tational cost, I designed a tiled inversion procedure. Figure 6.4 presents the tilingconfiguration with respect to the data set. Neighboring regions overlap by 500 m inorder to share a minimum of two survey lines, guaranteeing good lateral continuity.The procedure is readily parallelized across multiple processors, each inverting asubset of the problem independently. The final magnetization model is interpolatedonto the base mesh using a weighted average scheme such that:mi =∑Tj=tm( j)i w( j)i∑Tj=t w( j)iw( j)i =[(xi− x( j)c )2+(yi− y( j)c )2]−1/2,(6.1)where the weights w( j)i are inverse distance between the ith cell and the centerlocation [xc;yc] of the jth tile. The number of tiles T overlapping at the ith cell isvariable depending on the horizontal location. The same weight is used for all cellsvertically at the same [xi;yi] location.Figure 6.5 presents a horizontal section at ≈50 m depth through the recoveredeffective susceptibility model. As expected from the observed data, several longand narrow anomalies are recovered. Based on the analysis by Buchan et al. (2009),I interpret the linear trends as regional dyke swarms belonging to the Mackenzie,Lac de Gras, MacKay and Malley groups. A large number of kimberlite pipes show116Figure 6.4: TMA data after removal of the regional signal and tiling configura-tion used for the regional inversion.as effective susceptibility highs, with few exceptions. The three Koala pipes, forexample, west of the Grizzly pipe, were not well recovered. Although importanteconomically, the anomalous magnetic response from the pipe is below the noiselevel. Proper imaging of the deposits would require a magnetic survey at lowerelevation to increase the noise to signal ratio.I note the clustering of kimberlite pipes near the intersection and along someof the dyke swarms. The apparent correlation between kimberlite pipes and dykesseems to reinforce the idea of structural controls as put forward by Wright (1999).While the topic is beyond the scope of this Masters thesis, future work in structuralgeology over the region may benefit from a regional scale inversion as presentedhere.Figure 6.6 compares the observed and predicted data from the tiled inversion.In the left panel, the predicted data from individual tiles were merged following6.1. The normalized residual data show horizontal striations due to variations inmagnetic data between adjacent lines, which could not be accounted for by theinversion. Overall, each invividual inversion could predict the data well within117the assigned uncertainties. On the other hand, the right panel shows the predicteddata from the final merged susceptibility model. I here notice that some of the lowfrequency content has been lost during the merging process, with higher residualsover the regions of low magnetic fields. This can be a problem when attempting toimage regional features and will require further research. In this case, we are onlyinterested in the dykes and kimberlite pipes, which for the most part were wellrecovered.Figure 6.5: Interpreted dykes (line) from the property-scale magnetization in-version and known kimberlite pipes location (circle). The 16 pipes chosen for adeposit-scale inversion are labeled.6.3 Deposit scale inversionFrom the property scale inversion, 16 sites were chosen based on the strength ofthe magnetic anomaly for a high resolution inversion. The goal of this section is togain knowledge about the bulk magnetization direction and compare those valuesto physical property measurements published in the literature.118Figure 6.6: Comparison between (a) observed and (b) merged predicted datafrom the recovered magnetization model over the Ekati Property. (c) Normal-ized residual data show horizontal striations due to variations in magnetic databetween adjacent lines, which could not be accounted for by the inversion. Mostof the predicted data from the individual tiles can fit the observed data within onestandard deviation. (d) Forward modeled data from the merged magnetizationmodel and (c) normalized data residual. It appears that a large portion of the thelow frequency content has been lost during the merging step. More research isrequired in order to preserve the long wavelength information.119For each inversion, a regional field was removed in order to focus the inversionon the local anomalies. Description of the regional field removal method can befound in Li and Oldenburg (1998). In summary, the following steps were taken:• Creation of a local mesh with large padding cells• Inversion of the data assuming induced magnetization• Removal of the local source within the core mesh. Only susceptible materialoutside the core region remain.• Forward modeling of the regional data at the observation location• Subtraction of the regional signal from the observed dataThe regional-removal step is important in regions where the regional field is largerthan the local anomaly. Figure 6.7 presents the magnetic data before and afterremoval of the regional field over the Misery pipe. Note that the original data aretrending up towards the east, with the largest magnetic field data occurring nearthe edge of the dataset. Consequently, the magnetic amplitude inversion wouldbe biased towards having a ring of effective susceptibility around the edge of thedataset, rather than the pipe itself. After removing the regional field, most of theanomalous fields come from a strong negative anomaly corresponding with thelocation of the Misery pipe.Following the regional field removal, each local tile is inverted with the CMIalgorithm on a 25 m cell size mesh. This resolution is required in order to properlymodel narrow anomalies in the range of 50 to 100 m in width. The process wasrepeated for the selected 16 kimberlite pipes. Table 6.3 summarizes the parametersused for the inversions. I imposed an l0-norm on the amplitude of magnetizationin order to focus the inversion within the core of the different kimberlites pipes.The direction of magnetization was allowed to vary smoothly in order to get anestimate of the variability, and hence to get a measure of uncertainty. A 27 pointaverage (3x3x3 cells) centered over the location of each pipe was used to computean effective susceptibility and direction as shown in Table 6.4. In the followingsection, the results are compared to those published in the literature.120Table 6.3: Deposit-scale inversion parameters.Core cell size 25×25×25 mBase mesh size 50 M cellsNumber of tiles 17M / tile ≈ 200k cellsN / tile ≈ 150 datap, qx, qy, qz 0, 2, 2, 2αs, αx, αy, αz 1.6e-3, 1, 1, 1Uncertainties 0.02*d + 5 nTTable 6.4: Average magnetization amplitude and direction from the local inver-sion over 16 known kimberlite pipes. Uncertainties were calculated from a three-cell cube standard deviation around the approximate location of each kimberlitepipe.Pipe ID X (m) Y (m) κe Incl (◦) Decl (◦)Arnie 532980 7176930 3.9e-2 ± 1.0e-2 -85 ± 4.2 355 ± 5.0Blackbear 519700 7176980 9.9e-4 ± 1.3e-4 -88 ± 1.9 211 ± 3.0Brent 530855 7174206 4.6e-3 ± 2.0e-3 -79 ± 3.1 299 ± 2.3Fox 515091 7169946 2.4e-2 ± 8.0e-3 66 ± 17.7 203 ±14.9Grizzly 521529 7177946 8.7e-2 ± 8.6e-3 -82 ± 1.3 100 ± 1.7Hawk 535550 7180934 2.0e-3 ± 1.6e-4 -84 ± 0.2 244 ± 0.1Jaeger 539575 7168970 2.3e-2 ± 6.6e-3 -73 ± 2.6 53 ± 2.8Leslie 515960 7173130 3.6e-1 ± 6.0e-2 65 ± 3.7 346 ± 5.2Mark 531202 7175715 4.5e-2 ± 7.0e-3 -83 ± 4.0 159 ± 5.4Nancy 515351 7180271 1.0e-2 ± 8.0e-4 -84 ± 2.1 11 ± 2.0Nora 514665 7167532 2.0e-3 ± 2.6e-4 -87 ± 2.1 287 ± 4.3Phoenix 540829 7162932 1.5e-4 ± 2.1e-5 -85 ± 1.8 238 ± 1.3Roger 521770 7174655 1.3e+0 ± 4.2e-1 -82 ± 5.3 160 ± 10.9Scorpion 531854 7167847 2.5e-2 ± 1.1e-2 68 ± 13.0 258 ± 23.6Zach 529812 7171900 3.4e-2 ± 8.8e-3 82 ± 7.9 185 ± 9.4Misery 539621 7159630 7.3e-3 ± 1.3e-3 -77 ± 5.0 172 ± 10.1121Figure 6.7: (a) Local data collected over the Misery pipe showing a local west-ern trend. (b) Regional field data are computed from a local inversion. Most ofthe signal comes from a dyke running north-south along the western edge of thelocal tile. (c) Residual data after regional field removal, showing a clear reverselymagnetized anomalie corresponding with the location of the Misery pipe.6.4 DiscussionMagnetization inversion at a regional scale reveals interesting patterns. Figure 6.8compares the induced and remanent components of magnetization over the prop-erty. The remanent component is noticeably stronger on most of the NE trendingMackenzie dykes, while the Lac de Gras and MacKay dykes seem to be mainlymagnetized along the inducing field. Several points of intersection between theLac de Gras, Mackenzie and Malley dykes in the south of survey [537,000 E ;7,176,000 N] show a strong remanent component. This result seems to agree withthe study of Buchan et al. (2009). The low resolution of this airborne survey doesnot support accurate modeling of narrow dykes on the order of 10 to 30 m width. Itcan however provide a relative time of emplacement when comparing the strength122and trend of effective susceptibility of intersecting dykes. Note the clear magneticoverprinting of a NW trending dyke onto a NE dyke near the Fox pipe [515,000 E; 7,170,000 N] and Leslie pipe [516,600 E ; 7,174,000 N]. This result fits with theidea of a younger Mackenzie dyke swarm (1270 Ma) crosscutting an older Lac deGras dyke (2020 Ma).Note that the NS trending dyke intercepting the Grizzly pipe has been recov-ered on its entire length with a relatively strong induced component. This is despitethe large footprint of a reversely magnetized body, which illustrates the strength ofthe CMI algorithm in its ability to distinguish between adjacent anomalies withvarying magnetic orientations.Several of the 16 pipes have been studied in the past, some of which haveborehole samples measured in laboratory. Table 6.5 summarizes the informationcollected from the literature. In general terms, my results seem to agree morewith Enkin (2003) and Cheman (2006) in terms of orientation compared to the re-sults of Zhao (2012). It is important to point out that the inversion results offera bulk estimate of magnetization within the pipe, which are known to be highlyheterogeneous. Laboratory measurements provided by Enkin (2003) showed largevariations in magnetic properties along individual boreholes. Moreover, larger un-certainties are expected from the recovered declination in cases where the magne-tization vector is nearly vertical (I > 80◦).Finally, I analyze the age of the kimberlite pipes in relation to the magnetizationdirection. Table 6.6 summarizes the age group as identified by Lockhart et al.(2004), as well as the predicted polarity from the geomagnetic timescale providedby Cande and Kent (1995). Figure 6.9 maps on a timescale the polarity of thevarious pipes with respect to Earth’s polarity reversal. It is interesting to note thatseveral of the predicted polarity directions from radiometric dating do not agreewith the recovered polarity. Those differences could be a result of uncertainties inthe Rb-Sr radiometric dating of ±1 Ma (2σ ).It is important to note the frequent changes in Earth’s polarity during the periodof 46 and 69 Ma, with at least 14 reversals identified. My results seem to suggestthat the age array MAA should be between 48 and 49 Ma in order to accountfor the recovered negative polarity. Similarly, the Leslie pipe should be eitheryounger than 51.7 or older than 52.7 to be normally magnetized. The Roger pipe123Figure 6.8: Horizontal sections through the recovered (a) induced and (b) re-manent magnetization model. Several pipes with strong remanence can easily beidentified as discrete circular anomalies.124Table 6.5: Summary of published results over selected kimberlite pipesSourceLeslie Misery GrizzlyI◦|D◦ I◦|D◦ I◦|D◦CMI 65 346 -77 172 -82 100Enkin (2003) 64 329 -74 168 ***Zhao (2012) 59 334 *** -58 316Cheman (2006) 82 *** *** -73 ***is particularly interesting as it is clearly reversely magnetized, and could have onlybeen formed during the short time window of 67.61 to 67.73 Ma.As suggested by Lockhart et al. (2004), the microdiamond abundance of kim-berlite pipes observed at Ekati may be correlated in time. It has been found how-ever that the magmatism can vary greatly between each eruption episode, evenif all episodes occur in close spatial proximity. Future research could look intocorrelations between magnetization strength and direction, as well as electricalproperties of the different pipes. Assuming that diamond grade and abundance arerelated to the mantle source, it can be expected that rock magnetic susceptibilityand electrical conductivity would be a better proxy for diamond content than ageof emplacement. High resolution magnetic and airborne EM inversions could leadto interesting insights into the Ekati property.125Figure 6.9: Age of 11 pipes from the Ekati region with respect to Earth’s polar-ity reversal.Table 6.6: Radiometric age and inverted magnetization direction for 11 pipes.Pipe IDAge Array AgePolarity ChronIncl DeclLockhart et al. (2004) (Ma) (+ , - ) (◦) (◦)BrentMAA47.1 C21n (+) -79.0 ± 3.1 299.0 ± 2.3Arnie 47.5 C21n (+) -85.9 ± 4.2 355.5 ± 5.0Mark 47.5 C21n (+) -83.9 ± 4.0 159.1 ± 5.4Hawk 48.0 (-) -84.3 ± 0.2 244.4 ± 0.1Grizzly 50.8 C23n.1n (+) -82.2 ± 1.3 100.1 ± 1.7LesliePAA52.1 (-) 65.4 ± 3.7 346.3 ± 5.2Zach 52.8 C24n.2n (+) 82.5 ± 7.9 185.8 ± 9.4RogerMSK67.6 C30n (+) -82.7 ± 5.3 160.1 ± 10.9Jaeger 69.1 (-) -73.2 ± 2.6 53.2 ± 2.8Nora 2106 - -87.4 ± 2.1 287.9 ± 4.3Misery 2480 - -77.6 ± 5.0 172.3 ± 10.1126Chapter 7ConclusionThe inversion of magnetic data affected by remanent magnetization is an activefield of research. Several inversion strategies have been proposed in the past. Themagnetic amplitude inversion of Shearer (2005) inverts for the distribution of ef-fective susceptibility from magnetic amplitude data. The algorithm is robust andcan define regions of high effective susceptibility, but no information is providedregarding the orientation of magnetization. On the other hand, the Magnetic Vec-tor Inversion (MVI) of Lelie`vre and Oldenburg (2009) has the ability to recover themagnetization direction over variable and complicated geology, but the algorithmrequires a priori information to constrain the solution. Moreover, most inversionalgorithms use l2-norm measures of model structure for regularization, yieldingsmooth and small models.In this thesis, I combined the magnetic amplitude inversion and the MVI algo-rithm into a Cooperative Magnetic Inversion (CMI) algorithm, in order to improvethe robustness and versatility of current magnetic inversion codes used separately. Ideveloped a flexible lp-norm regularization that allows for sparse and blocky mod-els. The regularization function can be applied on sub-regions of the model do-main, allowing for a smooth transition in norm penalties. Magnetic amplitude dataare computed by the Equivalent Source method adapted from Li and Oldenburg(2010). I tested the algorithm on a synthetic model, showing improvement over theMVI method alone, yielding a simpler and more compact magnetization model.This in turn can help differentiate between neighboring anomalies with variable127magnetization directions.Finally, I implemented the CMI algorithm on a large aeromagnetic survey ofthe Ekati Property. Information regarding the polarity and orientation of mag-netization for 16 kimberlite deposits are compared to values from the literature.Relative age of dyke swarms are inferred from the apparent overprinting of mag-netization. The inverted magnetization model yielded information about possibledistribution of reversed magnetization at a regional level.7.1 Future workThis thesis opens up several avenues for future research. While the CMI algorithmhas the potential to estimate the magnetization vector orientation, there is still aneed to differentiate between the induced and NRM components. In cases wherethe NRM is parallel to the present inducing field, the current method cannot dif-ferentiate between the induced and remanent part of magnetization. One optionwould be to use frequency-domain electromagnetic (FEM) data to estimate themagnetic susceptibility of rocks. FEM surveys are routinely used in mineral ex-ploration projects to identify conductive anomalies. Inversion of FEM data couldserve the dual purpose of providing susceptibility and conductivity values, as wellas additional geometrical constraints for the inversion. In the case of kimberlite ex-ploration, FEM data could also be used to characterize the electrical conductivity ofkimberlite pipes. Combined magnetization and conductivity information may helpdistinguish between different pipes and estimate the diamondiferous potential.The CMI algorithm could be implemented on other types of mineral deposits,as well as on global geophysical problems. The classic problem of ocean floor mag-netization would be an interesting subject to revisit with this inversion technique.Accurate modeling of the magnetization direction may provide important insightsinto the transient behavior of Earth’s magnetic field during polarity changes.From a general inverse problem standpoint, the mixed-norm regularization im-plemented by the S-IRLS method could be implemented in a broad variety of in-verse problems in geophysics. More specifically, the regularization may be appliedto constrain non-linear inverse problems such as encountered in electromagnetics.The issue regarding the vertical stretching experienced with the amplitude in-128version needs to be addressed. The closer the recovered effective susceptibility isto the true model, the more accurate the MVI solution will be. There might bea need to iterate between the amplitude and MVI steps in order to progressivelyrefine the magnetization model.The Cartesian formulation for the MVI is not well suited for compact normconstraints. Sparsity constraints on the individual components of magnetizationtend to yield polarized models along the {pˆ, sˆ, tˆ} directions. Future work shouldinvestigate the spherical formulation of the MVI, which would allow applicationof sparsity constraints directly on the amplitude of magnetization. In the sphericalformulation, the magnetization vector is parameterized as an amplitude and twoangles (Lelie`vre and Oldenburg, 2009). This will allow additional information tobe incorporated through a reference model with the constraint information com-ing from the amplitude of the magnetization or the Koenigsberger ratio of rocksamples.The current CMI code suffers from memory limitations due to the large size ofthe sensitivity matrix. Incorporating wavelet compression, as seen in Li and Old-enburg (2003), may greatly reduce the computational time and size of the problem.Some of these limitations may also be overcome by working in the differentialequation domain (Davis et al., 2013). While the tiled inversion procedure pro-posed in Chapter 6 was effective in reducing the overall size of the inverse prob-lem, there were several issues with the method. The robustness and accuracy ofthe final model depends on the distance of overlap between the tiles, as well as onthe merging strategy. A more efficient and rigorous approach would be to separatethe inversion mesh from the forward modeling mesh. This approach was recentlyproposed to solve large-scale EM problems (Haber and Schwartzbach, 2014; Yangand Oldenburg, 2014). The inversion code could handle large data sets discretizedon a fine mesh, would be highly parallelizable and could potentially eliminate theneed for compression.129BibliographyAjo-Franklin, J. B., B. J. Minsley, and T. M. Daley, 2007, Applying compactconstraints to differential traveltime tomography: Geophysics, 72, 67–75. →pages 56, 58, 59, 68Barbosa, V. C. F., and J. B. C. Silva, 1994, Generalized compact gravity inversion:Geophysics, 59, 57–68. → pages 56, 68Bhattacharyya, B. K., 1964, Magnetic anomalies due to prism-shaped bodies witharbitrary magnetization: Geophysics, 29, 517–531. → pages 43Blakely, R., 1996, Potential theory in gravity and magnetic applications:Cambridge University Press. → pages 3Buchan, K. L., A. N. LeCheminant, and O. van Breemen, 2009, Paleomagnetismand U-Pb geochronology of the Lac de Gras diabase dyke swarm, SlaveProvince Canada: Implications for relative drift of Slave and Superiorprovinces in the Paleoproterozoic: Canadian Journal of Earth Sciences, 46,361–379. (NRC Research Press). → pages viii, 31, 111, 113, 116, 122Campbell, W. H., 1997, Introduction to geomagnetic fields, second edition ed.:Cambridge University Press. → pages 2, 14Cande, S. C., and D. V. Kent, 1995, Revised calibration of the geomagneticpolarity timestime for the late cretaceous and cenozoic: Journal of GeophysicalResearch, 100, 6093–6095. → pages ix, 4, 123Carlson, J. A., P. J. Ravenscroft, C. Lavoie, and J. Cumming, 2015, Ekati diamondmine, Northwest Territories, Canada: Ni 43-101 technical report, DominionDiamond Coportation. → pages 110, 111Chartrand, R., 2007, Exact reconstruction of sparse signals via nonconvexminimization: IEEE Signal Processing Letters, 76, 167–188. → pages 56, 58,69Cheman, A., 2006, Inversion magnetique tridimensionnelle des anomalieskimberlitique isolees avec presence ou sans remanence: Msc a, EcolePolytechnique de Montreal. → pages 111, 123, 125Clark, D. A., 1991, Magnetic petrology of igneous intrusions: implications forexploration and magnetic interpretation: Exploration Geophysics, 30, no. 2,1305–26. → pages 3Coopersmith, H., J. Pell, and B. Scott Smith, 2006, The importance of kimberlitegeology in diamond deposit evaluation: A case study from the DO27/DO18kimberlite, NWT, Canada: Presented at the 2006 Kimberlite EmplacementWorkshop, Saskatoon. → pages 110Cowan, D. R., L. A. Tompkins, and S. Cowan, 2000, Screening kimberlitemagnetic anomalies in magnetic active areas: Exploration Geophysics, 31,66–72. → pages 5, 110Creaser, R. A., H. Grutter, J. Carlson, and B. Crawford, 2004, Macrocrystalphlogopite RbSr dates for the Ekati property kimberlites, Slave Province,Canada: Evidence for multiple intrusive episodes in the Paleocene and Eocene:Lithos, 76, 399–414. → pages 111, 113Dampney, C. N. G., 1969, The equivalent source technique: Geophysics, 34,39–53. → pages 44Dannemiller, N., and Y. Li, 2006, A new method for estimation of magnetizationdirection: Geophysics, 71, L69–L73. → pages 6Daubechies, I., R. DeVore, M. Fornasier, and C. S. Gunturk, 2010, Iterativelyre-weighted least squares minimization for sparse recovery: Communication onPure and Applied Mathematics, 63, 1–35. → pages 53Davis, K., E. Haber, and D. W. Oldenburg, 2013, Large-scale magnetic inversionusing differential equations and octrees: Presented at the 23rd GeophysicalConf. and Exhibition, Australian Soc. of Expl. Geophys. → pages 129Doyle, B. J., K. Kivi, and B. H. S. Smith, 1999, The Tli Kwi Cho (DO27 andDO18) Diamondiferous Kimberlite Complex, Northwest Territories, Canada:Proceedings of the VIIth International Kimberlite Conference, 194–204. →pages 95Dransfield, M., A. Christensen, and G. Liu, 2003, Airborne vector magneticmapping of remanently magnetized banded iron formation at Rocklea, WesternAustralia: Exploration Geophysics, 34, 93–96. → pages 5Ekblom, H., 1973, Calculation of linear best lp-approximation: BIT, 13, 292–300.→ pages 54Ellis, R. G., B. de Wet, and I. N. Macleod, 2012, Inversion of magnetic data fromremanent and induced sources: Presented at the 22th International GeophysicalConference and Exhibition, ASEG. → pages 7, 8Enkin, R., 2003, Paleomagnetic analysis of selected Ekati kimberlite.Unpublished report for BHP Billiton Diamonds. Geological Survey of Canada.→ pages 113, 123, 125Enkin, R. J., 2014, The rock physical property database of British Columbia, andthe distinct petrophysical signature of the Chilcotin basalts.: Canadian Journalof Earth Sciences, 51, 327–338. (NRC Research Press). → pages 5, 31131Farquharson, C. G., and D. W. Oldenburg, 1998, Non-linear inversion usinggeneral measures of data misfit and model structure: Geophysical JournalInternational, 134, 213–227. → pages 53Gorodnitsky, I. F., and B. D. Rao, 1997, Sparse singal reconstruction from limiteddata using FOCUSS: A re-weighted minimum norm algorithm: IEEETransactions on Signal Processing, 45, 600–616. → pages 68Haber, C., and C. Schwartzbach, 2014, Parallel inversion of large-scale airbornetime-domain electromagnetic data with multiple OcTree meshes: InverseProblems, 30, 1–28. → pages 129Kubota, R., and A. Uchiyama, 2005, Three-dimensional magnetization vectorinversion of a seamount: Earth Planets Space, 691–699. → pages 7Last, B. H., and K. Kubik, 1983, Compact gravity inversion: Geophysics, 48, no.8, 713–721. → pages 54, 56, 59, 68Lawson, C. L., 1961, Contribution to the theory of linear least maximumapproximation: PhD thesis, University of California, Los Angeles. → pages 54Lelie`vre, P., 2003, Forward modelling and inversion of geophysical magneticdata: Master’s thesis, University of British Columbia. → pages 2, 24Lelie`vre, P. G., 2009, Integrating geological and geophysical data throughadvanced constrained inversion: Phd thesis, The University of BritishColumbia. → pages 7, 8, 31, 34Lelie`vre, P. G., and D. W. Oldenburg, 2009, A 3D total magnetization inversionapplicable when significant, complicated remanence is present: Geophysics,74, 21–30. → pages 7, 8, 127, 129Li, Y., 2012, Recent advances in 3D generalized inversion of potential-field data:SEG Technical Program Expanded Abstracts, 1–7. → pages 6Li, Y., M. N. Nabighian, and D. Oldenburg, 2014, Using and equivalent sourcewith positivity for low-latitude reduction to the pole without stration:Geophysics, 79, J81–J90. → pages 44, 47Li, Y., and D. W. Oldenburg, 1996, 3-D inversion of magnetic data: Geophysics,61, 394–408. → pages 5, 8, 18, 25, 26, 53, 60, 79——–, 1998, Separation of regional and residual magnetic field data: Geophysics,63, 431–439. → pages 120——–, 2003, Fast inversion of large-scale magnetic data using wavelet transformsand a logarithmic barrier method: Geophysics Journal International, 152,251–265. → pages 24, 129——–, 2010, Rapid construction of equivalent sources using wavelets:Geophysics, 75, no. 3, L51–L59. → pages 7, 8, 44, 127Li, Y., S. E. Shearer, M. M. Haney, and N. Dannemiller, 2010, Comprehensiveapproaches to 3D inversion of magnetic data affected by remanentmagnetization: Geophysics, 75, L1–L11. → pages 5, 38132Liu, S., X. Hu, Y. Xi, T. Liu, and S. Xu, 2015, 2D sequential inversion of totalmagnitude and total magnetic anomaly data affected by remanentmagnetization: Geophysics, 80, K1–K12. → pages 8Lockhart, G., H. Grutter, and J. Carlson, 2004, Temporal, geomagnetic and relatedattributes of kimberlite magmatism at Ekati, Northwest Territories, Canada:Lithos, 77, 665–682. → pages 113, 123, 125, 126Nabighian, M. N., 1972, The analytic signal of two-dimensional magnetic bodieswith polygonal cross-section: Geophysics, 37, 507–517. → pages 6, 7, 38Nabighian, M. N., V. J. S. Grauch, R. O. Hansen, T. R. LaFehr, Y. Li, J. W. Peirce,J. D. Phillips, and M. E. Ruder, 2005, The historical development of themagnetic method in exploration: Geophysics, 70, no. 6, 33ND–61ND. →pages 1Pell, J. A., 1997, Kimberlites in the Slave Craton, Northwest Territories, Canada:Geosicence Canada, 24, 77–90. → pages 110, 112Phillips, J. D., 2003, Can we estimate total magnetization directions fromaeromagnetic data using Helbig’s formulas: Abstracts Week B, IUGG 2003,B261. → pages 6Pilkington, M., 1997, 3-D magnetic imaging using conjugate gradients:Geophysics, 62, 1132–1142. → pages 25Portniaguine, O., and M. S. Zhdanov, 2002, 3-D magnetic inversion with datacompression and image focusing: Geophysics, 67, 1532–1541. → pages 58,69, 80Portniaguine, O. N., 1999, Image focussing and data compression in the solutionof geophysical inverse problems: Phd thesis, Department of Geology andGeophysics, The University of Utah. → pages 54, 69Power, M., and D. Hildes, 2007, Geophysical strategies for kimberlite explorationin northern Canada: Proceedings of Exploration ’07: Fifth DecennialInternational Conference on Mineral Exploration, 1025–1031. → pages 113Roest, W. R., J. Verhoef, and M. Pilkington, 1992, Magnetic interpretation usingthe 3-D analytic signal: Geophysics, 57, 116–125. → pages 6Sharma, P. V., 1966, Rapid computation of magnetic anomalies anddemagnetization effects caused by bodies of arbitrary shape: Pure Appl.Geophys., 64, 89–109. → pages 11Shearer, S., 2005, Three-dimensional inversion of magnetic data in the presenceof remanent magnetization: Master’s thesis, Colorado School of Mines,Golden, Colorado. → pages 7, 8, 38, 110, 127Shewchuk, J. R., 1994, An introduction to the Conjugate Gradient method withoutthe agonizing pain: School of Computer Science. → pages 21Stocco, S., A. Godio, and L. Sambuelli, 2009, Modelling and compact inversionof magnetic data: A Matlab code: Computer and Geosciences, 35, 2111–2118.133→ pages 56, 59, 68Sun, J., and Y. Li, 2014, Adaptive lp inversion for simultaneous recovery of bothblocky and smooth feature in geophysical model: Geophysical JournalInternational, 197, 882–899. → pages 53, 58, 68, 82Tikhonov, A. N., and V. Y. Arsenin, 1977, Solution of Ill-posed Problems:Winston Press. → pages 17, 50van der Vorst, H. A., 2003, Iterative Krylov methods for large linear systems:Cambridge University Press. → pages 23Vine, F. J., and D. H. Matthews, 1963, Magnetic anomalies over the oceanicridges: Nature, 199, 947–949. → pages 1Vogel, C. R., 2002, Computational Methods for Inverse Problems (Frontiers inApplied Mathematics): Society for Industrial Mathematics. → pages 24Wilkinson, L., B. A. Kjarsgaard, A. N. LeCheminant, and J. Harris, 2001, Diabasedyke swarms in the Lac de Gras area, Northwest Territories, and theirsignificance to kimberlite exploration: initial results: Current Research2001-C8, Geological Survey of Canada. → pages 95Wright, K. J., 1999, Possible structural controls of kimberlites in the Lac de GrasRegion, Central Slave Province, Northwest Territories, Canada: Msc thesis,Queen’s University. → pages 112, 117Yang, D., and D. W. Oldenburg, 2014, 3-d inversion of airborne electromagneticdata parallelized and accelerated by local mesh and adaptive soundings:Geophysical Journal International, 196, 1492–1507. → pages 129Yip, C. G., and K. S. Thompson, 2015, Diavik diamond mine, NorthwestTerritories, Canada: Ni 43-101 tehnical report, Diavik Diamonf Mines Inc. →pages 111Zhao, P., 2012, Three dimensional inversion of magnetic survey data collectedover kimberlite pipes in presence of remanent magnetization: Msc. a thesis,Ecole Polytechnique de Montreal. → pages 110, 111, 123, 125Zhdanov, M., and E. Tolstaya, 2004, Minimum support nonlinear parametrizationin the solution of a 3d magnetotelluric inverse problem: Inverse Problems, 20,937–952. → pages 72134
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A cooperative magnetic inversion method with Lp-norm...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A cooperative magnetic inversion method with Lp-norm regularization Fournier, Dominique 2015
pdf
Page Metadata
Item Metadata
Title | A cooperative magnetic inversion method with Lp-norm regularization |
Creator |
Fournier, Dominique |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | In this master's thesis, I implement a Cooperative Magnetic Inversion (CMI) algorithm for the 3-D modeling of magnetic rocks at depth. While in most cases it is assumed that the magnetic response is purely induced, certain rocks have the ability to retain a permanent magnetic moment in any orientations, also known as remanent magnetization. The effect of remanence has long been recognized as an obstacle for the geological interpretation and modeling of magnetic data. My objective is to improve current magnetic inversion methods to recover simpler and better defined magnetization models. The CMI algorithm brings together three inversion techniques previously introduced in the literature. First, magnetic data are inverted for an equivalent-source layer, which is used to derive magnetic amplitude data. Next, amplitude data are inverted for an effective susceptibility model, providing information about the geometry and distribution of magnetized objects. Finally, the effective susceptibility model is used to constrain the Magnetic Vector Inversion (MVI), recovering the orientation and magnitude of magnetization. All three algorithms are formulated as regularized least-squares problems solved by the Gauss-Newton method. I further constrain the solution by imposing sparsity constraints on the model and model gradients via an approximated lp-norm penalty function. I elaborate a Scaled Iterative Re-weighted Least-Squares (S-IRLS) method, allowing for a stable and robust convergence of the algorithm while combining different lp-norms on the range 0 ≤ p ≤ 2. The goal is to reduce the complexity of magnetization models while also imposing geometrical constraints on the solution. As a final test, I implement the CMI algorithm on an airborne magnetic survey over the Ekaty Property, Northwest Territories. I formulate a tiled inversion scheme in order to reduce the computational cost and increase the level of parallelization. The final merged magnetization model provides insights into the distribution of dyke swarms and kimberlite pipes. Following the regional inversion, I focus my analysis on sixteen known kimberlite pipes. Magnetization vector inclinations are compared to the expected polarity of rocks inferred from radiometric dating. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-10-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0166794 |
URI | http://hdl.handle.net/2429/55089 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_november_fournier_dominique.pdf [ 22.18MB ]
- Metadata
- JSON: 24-1.0166794.json
- JSON-LD: 24-1.0166794-ld.json
- RDF/XML (Pretty): 24-1.0166794-rdf.xml
- RDF/JSON: 24-1.0166794-rdf.json
- Turtle: 24-1.0166794-turtle.txt
- N-Triples: 24-1.0166794-rdf-ntriples.txt
- Original Record: 24-1.0166794-source.json
- Full Text
- 24-1.0166794-fulltext.txt
- Citation
- 24-1.0166794.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-0166794/manifest