Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Advanced potential field data inversion with ℓp-norm regularization Fournier, Dominique 2019

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

Item Metadata


24-ubc_2019_november_fournier_dominique.pdf [ 53.81MB ]
JSON: 24-1.0380907.json
JSON-LD: 24-1.0380907-ld.json
RDF/XML (Pretty): 24-1.0380907-rdf.xml
RDF/JSON: 24-1.0380907-rdf.json
Turtle: 24-1.0380907-turtle.txt
N-Triples: 24-1.0380907-rdf-ntriples.txt
Original Record: 24-1.0380907-source.json
Full Text

Full Text

Advanced potential field data inversion with `p-normregularizationbyDominique FournierMSc, The University of British Columbia, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Geophysics)The University of British Columbia(Vancouver)September 2019c© Dominique Fournier, 2019The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitledAdvanced potential field data inversion with `p-norm regularizationsubmitted by Dominique Fournier in partial fulfillment of the requirementsfor the degree of Doctor of Philosophyin GeophysicsExamining Committee:Douglas W. Oldenburg, Earth, Ocean and Atmospheric SciencesSupervisorEldad Haber, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberPhil Austin, Earth, Ocean and Atmospheric SciencesUniversity ExaminerIan Mitchell, Computer ScienceUniversity ExaminerAdditional Supervisory Committee Member:Craig Hart, Earth, Ocean and Atmospheric SciencesSupervisory Committee MemberiiAbstractThe study of gravity and magnetism has a long history in Earth sciences and contin-ues to play an important role in exploration geophysics. Fundamental knowledgeabout our planet has been gained from the processing of potential field data at allscales: from large tectonic processes, down to the size of mineral deposits. At-tempting to model density and magnetization through the inverse process remainschallenging however as it is a highly non-unique problem. Conventional uncon-strained inversions generally yield either smooth models with poor edge definitionor sparse and compact models that are too simplistic to represent subtle features.In this doctoral dissertation, I provide three technical innovations to enhancethe capabilities of potential data inversion. First, I propose to explore a widerrange of solutions by independently varying the sparsity assumption imposed onthe model values and its gradients. At the core of this research is the use of mixed`p-norm regularization solved by a scaled Iterative Re-weighted Least Squares al-gorithm. I provide a path to extract dominant features from the solution space witha Principal Component Analysis and pattern recognition strategy.Secondly, I provide modification to the magnetic vector inversion in sphericalcoordinates to process magnetic field data affected by remanent magnetization. Itackle long-lasting issues related to the non-linear trigonometric coordinate trans-formation with an iterative sensitivity re-weighting strategy. The spherical formu-lation allows me to impose sparsity constraints on the amplitude and orientation ofthe magnetization vector independently. The algorithm can recover anomalies witha coherent magnetization direction.Lastly, I incorporate soft geological information in the inversion through the ro-tation of the objective function. Combined with sparsity assumptions, the rotatediiiregularization improves the imaging of orientated geological contacts. Rotationangles are either interpolated from structural measurements or inferred from thedata through a learning process. I demonstrate the benefits of structural constraintson gravity and magnetic datasets acquired over the Kevitsa Ni-Cu-PGE deposit,Finland. The recovered density model is compared to a seismic reflection profilefor validation and to complement the geological interpretation of the deposit. Ac-curate modelling of the magnetization vectors yields insights about past tectonicdeformations.ivLay SummaryThis research focuses on geophysical imaging of Earth’s crust based on the den-sity and magnetic properties of rocks. This is done through an inversion processwhereby gravity and magnetic field data collected above the surface are convertedto a 3D model representing the sub-surface. Due to the nature of the inverse prob-lem, there are infinite solutions that can satisfy the observed data. Conventionalinversion algorithms yield smooth models that are difficult to interpret within a ge-ological framework. I tackle these limitations with a flexible algorithm that allowsme to generate a suite of possible candidates with variable characteristics. The al-gorithm is implemented for both scalar density and magnetization vector models. Ithen provide a learning strategy to extract dominant features from this suite of so-lutions using Principal Component Analysis. Finally, the efficacy of my approachis demonstrated using gravity and magnetic data acquired over the Kevitsa nickeldeposit in Finland.vPrefaceThe research presented in this thesis was completed by myself, under the guidanceof my supervisor Professor Oldenburg. The algorithms and ideas brought forwardresulted in multiple articles, either published or in revision.The mixed norm algorithm presented in Chapter 3 and part of the learning algo-rithm presented in Chapter 6 resulted in the article “Inversion using spatially vari-able mixed `p-norm” published in Geophysical Journal International (Fournierand Oldenburg, 2019b). I am the primary author of this article with revisionsfrom Dr. Oldenburg. The content of the article had previously been presented at aconference (Fournier et al., 2016). The same algorithm has been used by other re-searchers and resulted in three research papers to which I am also co-author (Abediet al., 2018a,b; Miller et al., 2017).Improvements to the magnetic vector inversion in Chapter 4 will be featuredin the accepted paper “Sparse magnetic vector inversion in spherical coordinates:Application to the Kevitsa Ni-Cu-PGE magnetic anomaly, Finland” submitted toGeophysics (Fournier and Oldenburg, 2019b). The article was written by myselfunder the guidance of Dr. Oldenburg. I have also co-authored a research paper in-vestigating the geothermal resources at Mount Baker (Schermerhorn et al., 2017).Using the same methodology, I collaborated in the implementation of sparse vec-tor inversion applied to self-potential problems to map hydrothermal circulation atMount Tongariro, New Zealand (Miller et al., 2018).The material presented in Chapter 5 and part of Chapter 7 is currently in prepa-ration for publication. The article entitled “Sparse rotated objective function forstratigraphic constraints: Application to the Kevitsa Ni-Cu gravity anomaly, Fin-land” will be submitted within a few weeks.viAll programming work done in this thesis builds upon the open-source SimPEGlibrary as well as multiple packages from the Python ecosystem. Accreditation toopen-source algorithms have been made wherever necessary.viiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Geophysical processing for geological interpretation . . . 31.1.2 Constraining the inversion . . . . . . . . . . . . . . . . . 71.1.3 Exploring the solution space . . . . . . . . . . . . . . . . 81.1.4 Dealing with remanence . . . . . . . . . . . . . . . . . . 91.2 Research objectives . . . . . . . . . . . . . . . . . . . . . . . . . 101.3 Thesis arrangement . . . . . . . . . . . . . . . . . . . . . . . . . 112 Forward modeling of potential field data . . . . . . . . . . . . . . . 132.1 Discrete systems . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1.1 Choice of discretization . . . . . . . . . . . . . . . . . . 17viii2.2 Large scale problems . . . . . . . . . . . . . . . . . . . . . . . . 212.2.1 Numerical test . . . . . . . . . . . . . . . . . . . . . . . 233 Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.0.1 Synthetic gravity example . . . . . . . . . . . . . . . . . 313.1 General `p-norm regularization . . . . . . . . . . . . . . . . . . . 343.1.1 Synthetic 1D problem . . . . . . . . . . . . . . . . . . . 353.1.2 Iterative Re-weighted Least Squares algorithm . . . . . . 393.1.3 Case 1: `1-norm (ps = px = 1) . . . . . . . . . . . . . . . 423.1.4 Case 2: `0-norm (ps = px = 0) . . . . . . . . . . . . . . . 463.2 Mixed norm regularization . . . . . . . . . . . . . . . . . . . . . 473.2.1 Scaled-IRLS steps . . . . . . . . . . . . . . . . . . . . . 483.2.2 Threshold ε-parameter . . . . . . . . . . . . . . . . . . . 553.2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3 Exploring the model space . . . . . . . . . . . . . . . . . . . . . 613.3.1 Interpretation . . . . . . . . . . . . . . . . . . . . . . . . 624 Sparse Magnetic Vector Inversion . . . . . . . . . . . . . . . . . . . 664.1 Magnetic Vector Inversion - Cartesian parameters . . . . . . . . . 704.2 Magnetic Vector Inversion - Spherical parameters . . . . . . . . . 714.2.1 Scaled MVI-S algorithm . . . . . . . . . . . . . . . . . . 814.3 Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 825 Rotated `p-norm Regularization . . . . . . . . . . . . . . . . . . . . 845.1 Rotated objective function . . . . . . . . . . . . . . . . . . . . . 855.1.1 Diagonal derivatives . . . . . . . . . . . . . . . . . . . . 925.2 Synthetic fold model . . . . . . . . . . . . . . . . . . . . . . . . 945.2.1 `p-norm inversion . . . . . . . . . . . . . . . . . . . . . . 945.2.2 Directional `p-norms . . . . . . . . . . . . . . . . . . . . 965.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 986 Automated mixed-norm modeling . . . . . . . . . . . . . . . . . . . 1006.1 Selecting local p-parameters . . . . . . . . . . . . . . . . . . . . 1016.1.1 Model space inversions . . . . . . . . . . . . . . . . . . . 101ix6.1.2 Average PCA model . . . . . . . . . . . . . . . . . . . . 1046.1.3 Parameter extraction . . . . . . . . . . . . . . . . . . . . 1066.1.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . 1086.2 Dip and strike estimation . . . . . . . . . . . . . . . . . . . . . . 1096.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1167 Case Study - Kevitsa Ni-Cu-PGE deposit . . . . . . . . . . . . . . . 1177.0.1 Geological setting . . . . . . . . . . . . . . . . . . . . . 1177.1 Geophysical data . . . . . . . . . . . . . . . . . . . . . . . . . . 1207.1.1 Modeling objectives . . . . . . . . . . . . . . . . . . . . 1247.2 Gravity inversion . . . . . . . . . . . . . . . . . . . . . . . . . . 1247.3 Magnetic inversion . . . . . . . . . . . . . . . . . . . . . . . . . 1277.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1298 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1418.1 Limitations and future work . . . . . . . . . . . . . . . . . . . . 143Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154xList of TablesTable 2.1 Summary table for the forward modeling of a magnetized sphereusing a Tensor and Octree discretization . . . . . . . . . . . . 21Table 3.1 Norm values and proportionality ratio obtained for the 1D so-lution presented in Figure 3.5(c). A proportionality ratio ofλ∞≈ 1 indicates that the components of the regularization func-tion are both contributing significantly to the final solution. . . 39Table 3.2 IRLS algorithm in pseudo-code made of two stages: Stage 1Initialization with convex least-squares inversion, Stage 2 IRLSupdates with inner β -search steps. . . . . . . . . . . . . . . . . 43Table 3.3 Components of the objective function corresponding to the in-version results presented in Figure 3.12 for ps = px = 1. . . . 55Table 3.4 Inversion summary after convergence of the S-IRLS for variouscooling rates η as presented in Figure 3.16. Each inversion trialwas required to reach the target misfit φ ∗d and target ε∗ = 1e−6 . 60Table 7.1 Summary table grouping the various lithological units loggedfrom bore hole. Expected density and magnetic susceptibilitycontrasts are derived from Figure 7.2. . . . . . . . . . . . . . . 121Table 7.2 Intervals along boreholes KV200 and KV297 reporting signifi-cant remanent magnetization (Montonen, 2012) . . . . . . . . 124xiList of FiguresFigure 1.1 Cartoon illustrating the (left) gravity and (right) magnetic fieldsoriginating from dense and magnetized anomalies. Adaptedfrom (Blakely, 1996). . . . . . . . . . . . . . . . . . . . . . . 2Figure 1.2 Total magnetic map over northern Europe extracted from theglobal Earth Magnetic Anomaly Grid (2-arc-minute resolu-tion) compiled from satellite, ship and airborne measurements(NCEI, 2017). The case study used in this research, the KevitsaNi-Cu-PGE deposit in Finland, sits within a region of strongmagnetization. The image overlay in GoogleEarth wasgenerated with the open-sourced geosci.GeoToolkit pack-age. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 1.3 (a) Surface geology, (b) Topography (c) ground gravity and (d)airborne magnetic data acquired over the Kevitsa deposit. . . . 5Figure 1.4 Figure comparing three approaches used to model the Kevitsadeposit. (a) Geological interpretation derived from seismicand borehole data (borrowed from Koivisto et al. (2015)). (b)Vertical section through an unconstrained density inversion us-ing conventional smooth assumptions. (c) Geological surfacemodel guided by geophysical forward modelling using Gocad. 6Figure 2.1 Parameters describing the spatial relationship between an ob-servation point P and a rectangular prism C. . . . . . . . . . . 14xiiFigure 2.2 Simple representation of (a) tensor and (b) Octree meshes forthe organization of rectangular prisms within a core domain(red) and padding region over a square domain 4 x 4 m inwidth. The Tensor mesh can fill the space with 64 cells, com-pared to only 40 cells with the Octree mesh. . . . . . . . . . . 17Figure 2.3 Three refinement strategies used to discretize a Gaussian curvedefined by scattered points (red): (a) box, (b) radial and (c)surface refinement strategy from the SimPEG.discretizelibrary. Cells are coloured by their corresponding Octree level. 19Figure 2.4 (a) Discretization of a sphere defined by discrete points (red)using the conventional Tensor mesh with a core region andpadding cells. Octree meshes refined by (b) box, (c) radialand (d) surface methods from the SimPEG.discretizelibrary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20Figure 2.5 (a) Horizontal and (b) vertical sections through the syntheticdensity model used to test the mesh decoupling strategy. (c)The simulated data contain short and long wavelength infor-mation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 2.6 Example of mesh decoupling of the forward problem into nestedOctree tiles. (Middle) Local meshes are generated for a pre-defined number of tiles nested inside the global domain. (Bot-tom) The forward modelling operator for each tile is stored incompressed Zarr file format on solid-state memory. Memorychunks (grey) are accessed in parallel by the Dask library toperform the forward calculations. . . . . . . . . . . . . . . . 25Figure 2.7 Trade-off curves between the size of the forward problem (blue)and data residual (red) over a range of number of tiles. The to-tal size of the problem is calculated based on the sum of cellsin all the local meshes times the number of data. The optimalnumber of tiles would be at the point of intersection where boththe cost of forward calculations and the data residual changesignificantly. . . . . . . . . . . . . . . . . . . . . . . . . . . 26xiiiFigure 2.8 Simulated gravity data calculated from (a) the global modeland (b) the 12 forward tiled calculations. (c) Data residualsshow short wavelength discrepancies between adjacent tilesdue to interpolation effects. . . . . . . . . . . . . . . . . . . . 26Figure 3.1 (a) Vertical section through a 25 m cube with density ρ=0.2g/cc placed in a uniform zero density background. (b) Simu-lated gravity data responses on a 21× 21 survey grid placed5 m above the flat topography. (c) Gravity data with randomGaussian noise added, 10−3 mGal standard deviation. . . . . . 32Figure 3.2 (a) Vertical section through the inverted density model usingthe conventional `2-norm regularization, (b) predicted and (c)normalized data residual. Outline of the true model (red) isshown for reference. . . . . . . . . . . . . . . . . . . . . . . 33Figure 3.3 Approximated `p-norm using the Lawson measure (Lawson,1961) over a range of p-values and for a fixed threshold pa-rameter ε = 10−1. . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 3.4 Linear forward problem made up of: (a) an example kernelfunction ; (b) model; (c) observed data with assigned standarderrors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36Figure 3.5 Solution to the 1D inverse problem using (a) an `2-norm onthe model (αx = 0), (b) the `2-norm on model gradients (αs =0) and (c) combined regularization function (αs = 2500, αx =1). (d) Convergence curve comparing the misfit (φd) and theregularization (φm) as a function of iterations. (e) Comparativeplot for the relative contribution of the different components ofthe objective function measured in terms of maximum absolutegradient (‖gi‖∞) . . . . . . . . . . . . . . . . . . . . . . . . . 37xivFigure 3.6 (a) Two solutions using an `1-norm on the model: (blue) Sim-plex, and (black) IRLS method. (b) Solution obtained withthe approximated `1-norm (IRLS) penalty on model gradientsalone and (c) with the combined penalty functions (αs = 2500, αx =1). The calculated proportionality ratio λ∞ indicates that thecombined penalties is dominated by the φ 1s term. (d) Conver-gence curve and (e) maximum partial derivatives associatedwith the components of the objective function as a function ofiterations for the inversion in (c). The vertical dotted lines in-dicate the change in regularization from an `2-norm to `1-normmeasure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44Figure 3.7 (a) Solution obtained with the combined penalty functions αsφ 1s +αxφ 1x after re-adjustment of αs = 50, αx = 1. (b) Convergencecurve and (c) maximum partial derivatives associated with thecomponents of the objective function as a function of iteration. 46Figure 3.8 Solution to the 1D inverse problem using an approximate `0-norm (a) on the model, (b) on model gradients and (c) com-bined penalty functions using the IRLS algorithm (αs = 1, αx =1). All three solutions honor the data within the target misfit φ ∗d . 47Figure 3.9 (a) Recovered model and (b) convergence curves using theconventional IRLS method for ps = 0, px = 2 and a fixed thresh-old parameter ε = 10−3 (αs = αx = 1). (c) Trade-off param-eter and maximum gradients for the different components ofthe objective function. At the start of Stage 2 (iteration 6),the sudden increase in ‖gs‖∞ is matched with a decrease in β .Throughout the inversion, ‖gx‖∞ remains small in magnitude. 49Figure 3.10 Derivatives of the Lawson approximation over a range of modelvalues for (a) a fixed threshold parameter ε = 10−1 over arange of p values and for (b) a fixed p = 0 over a range ofε values. (c) Applying the γ-scaling to the gradients brings allmaximums to be equal irrespective of p and ε . . . . . . . . . 50xvFigure 3.11 (a) Recovered model and (b) convergence curves using theScaled-IRLS approach for ps = 0, px = 2 and a fixed thresh-old parameter ε = 1e−3 (αs = αx = 1). (c) Trade-off param-eter and maximum gradients for the different components ofthe objective function. The scaling procedures preserves theproportionality between∥∥gpss ∥∥∞ and ∥∥gpxx ∥∥∞ throughout the it-eration process. The trade-off β -parameter needed only to beadjusted slightly at the beginning of Stage 2. . . . . . . . . . 52Figure 3.12 Recovered 1D models for ps = px = 1 using different scalingstrategies and parameterization. (a) Solution previously shownin Figure 3.7 that uses the standard gradient measure (αs=50).(b) Solution obtained with the finite difference approach (αs =αx = 1). (c) Model recovered with a different parameterizationsuch that the right half of the domain has cells with half thesize. Small discrepancies between the three solutions can beattributed to slight differences in the iterative process. . . . . 55Figure 3.13 Recovered 1D models with variable threshold parameter on therange 10−5 < ε < 100 using a mixed-norm penalty functionφm = αsφ 0s +αxφ 2x . . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.14 (a)-(d) Recovered 1D models at different iteration steps (k).The value of ε (dash) is shown for reference, highlighting theidea of a progressive thresholding of model values. (e) Scaledpartial derivatives (γ2g0s ) as a function of the sorted model val-ues and at various iteration stages. . . . . . . . . . . . . . . . 58Figure 3.15 Partial derivatives of the objective function for iteration (top)k=15 and (bottom) k=55. The recovered models are shown inred for reference. . . . . . . . . . . . . . . . . . . . . . . . . 59Figure 3.16 (a-d) Recovered models and (e) convergence curves for theminimization of φd+βφpm with various cooling rates η but witha final ε∗ = 1e−6. . . . . . . . . . . . . . . . . . . . . . . . 59Figure 3.17 (a-i) Vertical section through a suite of density models recov-ered for varying `p-norm penalties applied on the model andmodel gradients for ps ∈ [0, 1, 2] and px = py = pz ∈ [0, 1, 2]. 63xviFigure 3.18 Iso-contour values for the 10th and 90th percentile of anoma-lous density calculated from the suite of models shown in Fig-ure 3.17. The outline of the target (red) is shown for reference.Contour lines tightly clustered indicate coherence between in-version trials. Negative anomalies (dash) appear to change sig-nificantly, to which I would assign lower confidence. . . . . . 64Figure 3.19 (a-i) Residual data map calculated from the suite of densitymodels for varying `p-norm penalties applied on the model andmodel gradients for ps ∈ [0, 1, 2] and px = py = pz ∈ [0, 1, 2]. 65Figure 4.1 (a) Vertical section through a 25 m cube with uniform mag-netization ~M [2.0 A/m, I: 45◦, D: 90◦]. (b) Simulated TMAdata response on a 21× 21 survey grid placed 15 m abovethe anomaly. (c) Magnetic data with random Gaussian noiseadded, 1 nT standard deviation. . . . . . . . . . . . . . . . . 68Figure 4.2 Vertical section through the recovered susceptibility model us-ing (a) smooth assumption (ps, px, py, pz = 2) and (b) sparse`p-norms to recover a compact model (ps, px, py, pz = 0).Both solutions show an anomaly with a false dip due to thewrong assumption of a vertical magnetization. . . . . . . . . . 69Figure 4.3 Vertical section through the recovered magnetization vectormodel using the Cartesian formulation with (a) smooth l2-normassumption and (b) sparsity constraints applied on all threeCartesian components(pis , pix , piy , piz = 0). Correspondingnormalized data residuals are shown in (b) and (d) . The trueposition and magnetization orientation of the block are shownin red for reference. . . . . . . . . . . . . . . . . . . . . . . . 72Figure 4.4 Measure of angle differences ∆θ converted to coterminal angle∆θˆ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74xviiFigure 4.5 Vertical section through the (a) starting model and (b) recov-ered magnetization vector model in Spherical coordinates. Thetrue position and magnetization orientation of the block areshown in red for reference. (c) Normalized data residuals showcorrelated signal. The inversion stopped after 15 iterations andwas enabled to further reduce the objective function. . . . . . 75Figure 4.6 Contour map for two objective functions and their gradients(arrows) for a two-parameter inverse problem solved in Carte-sian and polar coordinate systems. (a) The non-weighted (gray)problem yields a solution mC = [0.2, 0.4] that reflects the sizeof the forward coefficients. The sensitivity weighted (black)solution is more desirable as it predicts the data with equalmodel parameters m∗C = [0.33,0.33]. (b) Inversion steps per-formed in the non-linear polar coordinate system (solid red) arehighly oscillatory and do not reach the optimal solution m∗P =[0.47,0.79]. The same model updates are shown in Carte-sian coordinates (red dash) for reference. (green) Inversionsteps performed in polar coordinates with iterative sensitivityre-weighting and (blue) with the added Ω-scaling. . . . . . . . 77Figure 4.7 (a) Vertical section through the recovered magnetization vec-tor model and (b) normalized data residual using the sphericalformulation with sensitivity based weighting (ps = px = py =pz = 2). The true position and magnetization orientation of theblock are shown in red for reference. . . . . . . . . . . . . . 82Figure 4.8 (a) Vertical section through the recovered magnetization vectormodel in Spherical coordinates using sensitivity based weight-ing with smooth regularization (ps = px = py = pz = 2). MVI-C solution is used as a starting model. The true position andorientation of magnetization of the block are shown in red forreference. (c) Recovered vector model in Spherical coordi-nates using sparse and blocky assumptions on the amplitudeand angles of magnetization (pcs = pcx = pcy = pcz = 0). (b)and (d) Corresponding normalized data residuals. . . . . . . . 83xviiiFigure 5.1 Rotation parameters with respect to the Cartesian axes. . . . . 87Figure 5.2 Two-dimensional representation of finite difference operatorsusing a combination of (a) 3-cell (Li and Oldenburg, 2000),(b) 5-cell (Lelie`vre, 2009) and (c) the 7-cell scheme used forthe rotation of the objective function. Only the 7-cell strategyallows for the interaction between diagonal neighbors. . . . . 90Figure 5.3 Horizontal cross section through the recovered model with ro-tation of the objective function using (top) 5-point and (bot-tom) 7-point gradient operators. Blue and red arrows indicatethe rotation frame. (a, d) Recovered models with smooth as-sumptions (ps = px = py = pz = 2) stretched along the rotatex-axis at 45◦ (αx′ = 100). (b, e) Sparse solutions (ps = px =py = pz = 0) with rotation at 45◦. The target orientation of therotated block is marked by the black dashed line. (c, f) Sparsesolutions for a 30◦ rotation of the objective function. . . . . . 91Figure 5.4 Normalized data residuals calculated from the recovered mod-els presented in Figure 5.3(d), (e) and (f) respectively. . . . . . 93Figure 5.5 Synthetic density model made up of a folded layer with a N-S oriented axial plane dipping 75◦E. The intersection betweenthe fold axial plane and the cross-sections are marked by ared-dashed line: (left) topography draped section and (right)vertical cross-sections along (A-A’) and (B-B’). The hinge axisof the fold is plunging 30◦S. Two strike and dip measurementsare provided at the surface and used as structural constraints. . 95Figure 5.6 Synthetic gravity data generated 1 m above topography overthe synthetic density fold model. Gravity station locations(dot) and topographic contours (lines) are shown for reference.Structural data (dip, strike) are provided at two locations onopposite of the syncline. Dip direction of the fold axis (15◦ S)and limbs are shown with arrows . . . . . . . . . . . . . . . . 96xixFigure 5.7 Recovered smooth solution with `2-norm regularization of thesynthetic fold anomaly. The dip of the fold limbs are poorly re-covered and there is no indication of the limbs being connectedat depth. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Figure 5.8 Recovered density model from the mixed norm regularization(ps = 0, px,y,z = 2). Sparse assumptions helped in simplifyingthe density model over the smooth assumption. . . . . . . . . 97Figure 5.9 Interpolated dip and strike information used for the rotation ofthe objective function. The direction of the arrows define theestimated normal or the folded layer. . . . . . . . . . . . . . . 98Figure 5.10 Recovered density model from the mixed-norm regularization.Data locations (dots), topography and the outlines of the truemodel (black) are shown for reference. Rotated `p-norm ori-ented along the folded layer helped in resolving a continuousgeological unit at depth. . . . . . . . . . . . . . . . . . . . . 99Figure 6.1 (a) Synthetic 2D travel-time tomography problem made up of arectangular velocity low, and a smooth Gaussian velocity high.Contour lines (25th and 75th percentile) are shown in black forthe true position and shape of the velocity anomalies. Trans-mitters (red) and receivers (green) are positioned on oppositesides of the model domain. (b) Data map for the 143 line inte-gral data calculated for each transmitter-receiver pair. . . . . 102Figure 6.2 Suite of inverted models for various combinations of normsps, px = py ∈ [0,1,2] (αs = αx = αy = 1). Contour lines (25thand 75th percentile) are shown in black for the true positionand shape of the velocity anomalies. . . . . . . . . . . . . . . 103Figure 6.3 (a) to (i) Data residual contour map for the nine inversions pre-sented in Figure 6.2. Contour lines (±1 standard deviation) areshown in black. (j) Random noise added to the experiment. . . 104Figure 6.4 PCA vectors covering 75% of the model variances. . . . . . . 105xxFigure 6.5 (a) Averaged PCA model mA and (b) corresponding normalizedata residual map. Contour lines are shown in (a) for the trueposition and shape of the velocity anomalies. . . . . . . . . . 106Figure 6.6 Parametric representation of the nine inverted models derivedform the Canny edge detection algorithm. . . . . . . . . . . . 107Figure 6.7 Local p-values on the (a) the model norm φ pss and (b) modelgradients φ pxx , φpyy extracted from the solution space. . . . . . 107Figure 6.8 (a) The final SVMN inversion result and (b) normalized dataresidual map. Contour lines are shown in (a) for the true posi-tion and shape of the velocity anomalies. . . . . . . . . . . . 109Figure 6.9 Synthetic fold model introduced in Chapter 5 used here to testthe dip and strike estimation learning algorithm. . . . . . . . . 110Figure 6.10 Example of an automated azimuth estimation using the UBCCrest. (b) The image is first converted to a binary image us-ing the python routine ScikitImage.Skeleton. Imagemoment calculations are performed on moving windows todetermine the position and orientation of dominant features.(c) Vectors are then rotated by 90◦ to point along the normalof edges. The normal vectors are checked for consistency inpreparation for interpolation. . . . . . . . . . . . . . . . . . . 112Figure 6.11 (a) Image moment calculations over windowed regions of theinverted density model. Arrows are plotted for the principal di-rection (red) and normal (white) of the main density anomaly.(b) Normal vectors are interpolated in 3D and used to rotatethe objective function. . . . . . . . . . . . . . . . . . . . . . 113Figure 6.12 Suite of models with directional gradients for θ ∈ [−90◦, 75◦]. 114Figure 6.13 Sections through the average PCA model and rotation vectors(arrow) derived from the image moment analysis. . . . . . . . 115Figure 6.14 (a) Horizontal and (b, c) vertical sections through the recoveredmodel using rotation angles derived from the average PCAmodel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115xxiFigure 7.1 (a) Geological surface map of the Kevitsa-Satovaara intrusivecomplex and (b) 2D seismic reflection survey along with theE5 profile. Interpreted seismic reflectors and geological con-tact (black) from the interpretation of Koivisto et al. (2015) areshown for reference, as well as the outline of the gravity andmagnetic data sets analyzed in this chapter. . . . . . . . . . . 118Figure 7.2 Whisker plot of logged (a) density and (b) magnetic suscepti-bility measured along 279 boreholes. The coloured boxes havea width scaled by the calculated standard deviation and cen-tered on the mean value for all intercepts belonging to the samelithological classification. The black lines on either side definethe minimum and maximum values. The different lithologiesare colour coded and grouped based on relative age and simi-larities in physical properties. . . . . . . . . . . . . . . . . . 122Figure 7.3 (a) Ground gravity and (b) airborne magnetic data map ac-quired over the Kevitsa Ni-Cu-PGE deposit, Finland. Colorranges are histogram equalized and sun shaded to highlight de-tails. The location of the Kevitsa mine (red), geological bound-aries (black) and faults (dash) identified from surface mappingare shown for reference. . . . . . . . . . . . . . . . . . . . . 123Figure 7.4 (a) Horizontal section at ≈ 200 m depth below topographythrough the recovered density model after convergence of thealgorithm for ps = px = py = pz = 2. (b) Vertical section of thedensity model overlaid on the 2D seismic reflection line E5.Interpreted geological contacts (black) are shown for compar-ison. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Figure 7.5 Rotation parameters (normals) derived from (a) surface geol-ogy and (b) major reflectors detected on the seismic profile E5. 131Figure 7.6 Horizontal sections at≈ 200 m depth below topography throughthe density models recovered over a range of `p-norms appliedon the model and modeled derivatives for ps, px,y,z ∈ {0, 1, 2} 132xxiiFigure 7.7 Vertical sections overlaid on the 2D seismic reflection line E5through the density models recovered over a range of `p-normsapplied on the model and modeled derivatives for ps, px,y,z ∈{0, 1, 2} . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133Figure 7.8 Residual data maps for the nine gravity inversions. . . . . . . 134Figure 7.9 (a) Horizontal and (b) vertical section comparing iso-surfacesfor the 5th and 95th percentile of anomalous density values re-covered from the nine mixed norm inversions for ps, px,y,z ∈{0, 1, 2}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135Figure 7.10 (a) Horizontal and b) vertical sections through the recoveredsusceptibility model that ignores the effect of remanence. Litho-logical contacts (black) identified by Koivisto et al. (2015) areshown for reference. A large dome-shaped zero-susceptibilityanomaly is recovered at the center of the Kevitsa intrusion,likely due to remanence. (c) Residual data map shows strongcorrelation with the negative anomalies. . . . . . . . . . . . . 136Figure 7.11 Horizontal sections at ≈ 300 m below topography for a suitemodels using various sparsity assumptions put on the ampli-tude of magnetization for ps , px,y,z ∈ [0, : 2]. Norm measureson the magnetization angle are fixed to px,y,z = 0 in order topromote uniform magnetization . . . . . . . . . . . . . . . . 137Figure 7.12 Vertical sections through the recovered magnetization models.Effective susceptibilities are concentrated within the UDU andon the lower margin of the UPXO unit. . . . . . . . . . . . . 138Figure 7.13 Residual magnetic data maps for the nine MVI-S inversions. . 139Figure 7.14 (Top) Horizontal and (bottom) vertical sections through iso-contours of magnetization recovered from nine mixed `p-norminversions. Magnetization orientation (white) and standard de-viation on the angle (red) are shown. . . . . . . . . . . . . . . 140xxiiiAcknowledgmentsI would like to first thank my supervisor, Professor Doug Oldenburg, for givingme the opportunity to join the UBC-GIF group and for his tireless excitement forresearch. Your passion for science is contagious.I want to acknowledge the outstanding contribution of many colleagues fromthe open-source community. Thanks to Rowan Cockett, Lindsey Heagy and SeogiKang for establishing the SimPEG framework and pushing me to learn Python.Special mention goes to Joe Capriotti for his key update to the Octree mesh dis-cretization. The work presented in this thesis would not have been possible withoutyour contribution. Thank you to the rest of the UBC-GIF group for your support,kind words and friendship.I am grateful to the companies and staff that have supported research at UBC-GIF over all these years: Anglo American, Barrick, BHP, Cameco, CGI, Glencore,Newmont, Rio Tinto, Teck, Vale. Personal thanks to John McGaughey at MiraGeosciences for keeping on staff over all these years, as well as for lending mevaluable hardware and software.Thank you to the many public libraries across Canada where this thesis waswritten. Cities and towns from east to west (*with great pool nearby): Ottawa,Winnipeg*, Prince George, Fort Saint-James, Smithers*, Vancouver, Richmond*,Steveston, Nanaimo, Courtney, Alert Bay, Port Hardy.Last but not least, extra special thanks goes to my best friend, my partner incrimes and lover Chantelle Bellrichard. I am done being a student now. Thank youto sister-wife Candice Nichol for providing distractions to first-wife, allowing meto finish my thesis. So many more BBQ & sunsets at the beach to be enjoyed ...xxivChapter 1Introduction1.1 MotivationDensity and magnetic permeability are important physical properties used in min-eral exploration. Geologists routinely measure these properties on rock samples,along with textural and structural features, to target areas for further investigation.The vast majority of known mineral deposits have been discovered this way (Mar-joribanks, 2010). But the conventional approach to mineral exploration is limitedto shallow regions of the Earth with exposed bedrock. Large portions of the crustremain untouched, buried under thick sedimentary cover.In the absence of surface evidence, density and magnetic permeability can stillplay a role in mineral exploration as a source of geophysical signal. Variations inthe local gravity and magnetic fields are used to map geology and identify buriedstructures such as intrusions, faults, folds and alteration zones (Domzalski, 1966;Grant, 1984). The physical equations for gravity and magnetics are described withthe aid of Figure 1.1. For gravity we have the scalar gravitational potential:φg = G∫V1rρ(r) dV (1.1)The gravity field is obtained by taking the gradient of (1.1)~g = ∇φg (1.2)1which becomes~g =−G∫Vrˆr2ρ(r) dV . (1.3)Equation (1.3) describes the gravity field ~g( Nkg ) as observed at some arbitrary po-sition P due to an elementary volume dV with density ρ ( kgm3 ). The field scales asa function of radial distance r between the source and the observer, and by New-ton’s gravitational constant G (6.674× 1011 Nm2kg2 ). For the magnetic problem, weFigure 1.1: Cartoon illustrating the (left) gravity and (right) magnetic fieldsoriginating from dense and magnetized anomalies. Adapted from (Blakely,1996).follow the notation of Blakely (1996). We consider magnetic minerals as beingcomposed of small dipoles, each with their own magnetic dipole moment. Prefer-ential alignment of these dipoles gives rise to a magnetization per unit volume ~M( Am)~M =total magnetic moment∆V. (1.4)The scalar magnetic potential can than be defined as:φb(r) =µ04pi∫V∇1r· ~M(r) dV , (1.5)2Taking the gradient of (1.5) yields an expression for the magnetic flux density:~b = ∇φb (1.6)and thus~b(r) =µ04pi∫V∇∇1r· ~M(r) dV . (1.7)In this above equation,~b is the magnetic flux density (Wbm2 orNA m ) and µ0 is the mag-netic permeability of free space ( Hm orNA2 ). Both equation (1.1) and (1.5) describea source of geophysical signal that decays radially with distance and that scaleslinearly with changes in physical property contrasts. As a result, the gravity andmagnetic experiments are collectively referred to as potential field methods, andthere is a field of mathematics devoted to this subject.A large number of potential field datasets has been compiled over the years andmade available to the geoscientific community for analysis: from global satelliteand marine measurements, to airborne and ground surveys (Figure 1.2). Process-ing this wealth of data has encouraged collaboration between geologists and geo-physicists in their efforts to map the sub-surface. After nearly five decades of dataacquisition and interpretation, gravity and magnetic methods are still an area of ac-tive research in geophysics (Nabighian et al., 2005). This Ph.D. project contributesto that body of research.1.1.1 Geophysical processing for geological interpretationWhile the study of potential field data represents a common ground for geolo-gists and geophysicists to investigate the Earth, the two branches of geoscienceappear to diverge on the path taken for interpretation. The majority of geologi-cal studies have focused on data processing techniques for mapping structures andgeological domains in 2D. Fourier filters have been particularly popular, such asthe total derivatives (Cordell and Grauch, 2012), analytic signal (Macleod et al.,1993; Roest et al., 1992), tilt angle (Miller and Singh, 1994) or combination offilters (Sanchez et al., 2014; Verduzco et al., 2004). Euler deconvolution has beenused to estimate the depth of simple parametric shapes such as planes, dykes, andspherical bodies (Mushayandebvu et al., 2001; Reid et al., 1990; Thompson, 1982).3Figure 1.2: Total magnetic map over northern Europe extracted from the globalEarth Magnetic Anomaly Grid (2-arc-minute resolution) compiled from satellite,ship and airborne measurements (NCEI, 2017). The case study used in this re-search, the Kevitsa Ni-Cu-PGE deposit in Finland, sits within a region of strongmagnetization. The image overlay in GoogleEarth was generated with theopen-sourced geosci.GeoToolkit package.Similarly, the tilt-depth method (Salem et al., 2007)and Source Parameter Imaging(SPI) (Thurston and Smith, 1997) can provide an estimate for the depth and dipof elongated bodies in 2D (Phillips, 2010). A common limitation of these meth-ods is in dealing with data acquired over complex geology. It is challenging tointerpret overlapping signals from multiple sources with arbitrary shapes and posi-tions in 3D. It is even more challenging to make sense visually of dipolar magneticanomalies, especially at low latitudes.The Kevitsa Ni-Cu-PGE deposit shown in Figure 1.3 is an interesting casestudy to illustrate the difficulty of interpreting potential field data in a mineral ex-ploration context. The gravity and magnetic maps are both visually challenging tointerpret. Attempting to construct a 3D geological model of Kevitsa directly fromthe data maps would undoubtedly be a difficult task. Yet there is a lot of infor-mation that could be extracted to help unravel the geological history of the region.4Figure 1.3: (a) Surface geology, (b) Topography (c) ground gravity and (d) air-borne magnetic data acquired over the Kevitsa deposit.To this end, a large body of work from the geophysical community has focusedon inversion methodologies to model physical property contrasts in 3D. Amongthese methods, voxel-based inversions have received considerable attention. Oneof the main challenges posed by the inverse problem is that many models can fitthe observed data; this is referred to as non-uniqueness. Additional informationmust be provided to isolate geologically sensible solutions. The typical strategyis to assume smooth changes in physical properties. This strategy has become anindustry standard that is adopted by several commercial codes such as the UBC-GRAV3D and MAG3D, (Li and Oldenburg, 1996; Li. and Oldenburg, 1998).The conventional smooth physical property inversion framework is somewhatat odds with the geological mindset however which is geared towards defining dis-crete geological domains with sharp boundaries. This divide can once again beexemplified with sections through the Kevitsa deposit (Figure 1.4). After years ofexploration work in the area, geologists have identified and mapped over a dozen5lithologies along 2D sections. In an ideal case, we would like to use geophysicalinversion to complete the modelling work in 3D. Attempting to do so with currentmethodologies yields the model presented in Figure 1.4(b). The smooth densityanomaly is too difficult to interpret in this context as it shows poor correlation withthe conceptual model of Kevitsa. It is partially for this reason that the geologicalcommunity has been slow to adopt 3D inversion methods as part of their interpre-tation workflow. The apparent disconnect between the two fields of geoscienceremains a key issue to be addressed (Li et al., 2019).Figure 1.4: Figure comparing three approaches used to model the Kevitsa de-posit. (a) Geological interpretation derived from seismic and borehole data (bor-rowed from Koivisto et al. (2015)). (b) Vertical section through an unconstraineddensity inversion using conventional smooth assumptions. (c) Geological surfacemodel guided by geophysical forward modelling using Gocad.61.1.2 Constraining the inversionA number of inversion strategies have been proposed in the past to help bridgethe gap between geologists and geophysicists. For simple scenarios involving afew isolated targets, geophysical anomalies may be approximated by simple para-metric bodies such as ellipsoids (McMillan et al., 2016) and dipping dykes (Foss,2006; Krahenbuhl and Li, 2015). For most complex geological settings, full 3Dsurface modelling can be used to constrain the inversion (Bosch and McGaughey,2001; Fullagar et al., 2008). Figure 1.4(c) demonstrates such modelling efforts onthe Kevitsa intrusion. Building upon the seismic interpretation of Koivisto et al.(2015), borehole logs and surface mapping, I constructed a 3D surface model us-ing the Gocad-SKUA package. Generating and testing this model against thegeophysical data is a laborious process that requires skills, advanced software andthat can rapidly become intractable.On a more data-driven side, hard physical property constraints have been in-cluded in the inversion on a cell-by-cell basis (Phillips, 1996; Williams, 2008) orthrough clustering algorithms (Lin and Zhdanov, 2017; Sun and Li, 2013). Dip andstructural information have also been used to reinforce known trends; this has beenaccomplished by rotating the regularization functions with manual adjustment ofmodel gradients (Davis et al., 2012; Li and Oldenburg, 2000). A common limita-tion with these approaches is in building and testing different geological scenarios.Time and computational limitations often force practitioners to arbitrarily chooseone set of parameters for the regularization function and accept the outcome as a”best” model and have it used in subsequent interpretation.Few studies in geophysics have attempted to address non-uniqueness within alearning framework. The genetic algorithm of Wijns and Kowalczyk (2007) andsupervised learning of Haber and Tenorio (2003) can guide practitioners in thechoice of inversion parameters and they could serve as a starting point for advanced3D modelling work. But the efficacy of these approaches is once again impeded byassuming a smooth physical property distribution. The learning process is limitedby the ability to generate solutions that are significantly diverse.71.1.3 Exploring the solution spaceAs an alternative to smooth inversions, several studies have employed approxima-tions to `p-norms in order to recover compact anomalies. Methods such as theEkblom norm (p = 1) (Ekblom, 1973) and the Lawson approximation (Lawson,1961) (p = 0) favor impulsive models with fewer non-zero parameters. Severalpapers have demonstrated the efficacy of a sparsity assumption in reducing thecomplexity of physical property models (Ajo-Franklin et al., 2007; Barbosa andSilva, 1994; Blaschek et al., 2008; Chartrand, 2007; Last and Kubik, 1983; Port-niaguine and Zhdanov, 2002; Stocco et al., 2009). Likewise, a sparse solution formodel gradients can yield blocky (piece-wise constant) models (Daubechies et al.,2010; Farquharson and Oldenburg, 1998; Gorodnitsky and Rao, 1997; Li, 1993;Sun and Li, 2014).Sparsity assumptions have increased the flexibility of inversion algorithms forthe modelling of compact targets with sharp edges. The usual strategy is to defineglobal parameters applied over the entire model space. In most geological settings,however, we can expect to encounter a mix of features, either smooth and sparsein one region but elongated and blocky elsewhere. Determining exactly whichassumption to use in a particular geological setting remains largely user-driven.Sun and Li (2014) have made inroads in imposing variable sparsity assumptions,either the `2 or `1-norm measure of the model, to different regions of a 2D seismictomography problem. They demonstrate that this choice could be automated basedon information present in the data. Regions reacting favourably to the sparsityassumption are defined through a learning process.The next logical step, which I am taking in this research, is to further generalizethe use of sparse assumptions for p ∈ [0, 2] applied to both the model and modelgradients in 3D. I provide algorithmic details regarding the mixing of smooth,sparse and blocky assumptions. This general inversion framework can be used togenerate diverse solutions and thus help geologists in their interpretation. A suiteof models with broad characteristics can also form a basis for learning algorithms.81.1.4 Dealing with remanenceContrary to density, magnetization has the added complexity of being a vectorquantity defined by strength and orientation. In matter, the total magnetization,which is the magnetic moment per unit volume ~M in (1.7) can be separated in itsinduced and remanent component such that:~M = κ(~H0+ ~Hs)+ ~Mr , (1.8)where the magnetic susceptibility κ (SI) is the physical property describing theability of a rock to get magnetized under an applied field. In nature, this inducingfield has two components. The geomagnetic field ~H0 originates from Earth′s core.It is in most cases the dominant component. Secondary fields ~Hs are related tolocal magnetic anomalies. For highly susceptible material, the secondary fieldscan oppose the geomagnetic field direction and reduce the total magnetization. Itis also referred to as self-demagnetization effects. The remanent magnetization~Mrem is a permanent dipole moment that is preserved in the absence of an inducingfield.It has long been assumed that the induced component of magnetization wasdominant. The effect of remanence is often regarded as ‘noise‘ and simply ignoredby 2D filtering and 3D inversion methods. Recent studies have shown however thatremanent minerals, most often magnetite and members of the titanium-hematite se-ries, are commonly associated with mineral deposits such as diamondiferous kim-berlites, volcanic massive sulphides and porphyries (Enkin, 2014; Henkel, 1991).The orientation of remanent magnetization can make geologic interpretation morecomplicated and should not be ignored.Meanwhile, the same remanent component has been used extensively in pale-omagnetic studies and in that field, it has been regarded as geophysical ‘data‘. Anumber of researchers have used the permanent magnetization orientation to mapcontinental block rotation (Kissel and Laj, 1989; Norris and Black, 1961; Vine andMatthews, 1963), for fold and thrust belts reconstruction (Ramon et al., 2012; Vil-lalain et al., 2015) and in geochronology (Enkin, 2003; Henkel, 1991; Lockhartet al., 2004). While providing valuable information about Earth’s history, thesestudies have relied primarily on laboratory measurements performed on oriented9cores. The availability and cost to acquire the orientation of magnetization at apoint remains a limiting factor for paleomagnetic studies (Pueyo et al., 2016).Researchers have investigated ways to estimate the magnetization in order todeal with remanence. This includes search algorithms (Dannemiller and Li, 2006;Fedi et al., 1994), magnetic moment analysis (Helbig, 1963; Phillips, 2003), andinversion methods. The most common inversion approach uses simple paramet-ric shapes to approximate elongated and tabular bodies (Clark, 2014; Foss andMcKenzie, 2011; Fullagar and Pears, 2013; Pratt et al., 2014). For more complexsettings, magnetic data can be transformed into a quantity that is weakly sensi-tive to the orientation of magnetization, such as magnetic amplitude data, and in-verted for a scalar quantity, the effective susceptibility (Shearer, 2005). Lelievreand Oldenburg (2009) introduce the Magnetization Vector Inversion (MVI) for the3D modelling of magnetization parameters. Because a vector is sought, there arethree times the number of parameters compared to conventional inversion. This in-creases the non-uniqueness, and useful information from an inversion will requiremore complex regularization functions. Queitsch et al. (2019) have shown the ben-efits of inverting full magnetic gradient tensor data along with sparsity assumptionsimposed on the amplitude of magnetization. Other studies have explored the useof sparsity assumptions for the recovery of compact bodies, but imposing such as-sumptions on the orientation of magnetization remains difficult (Fournier, 2015;Liu et al., 2015; Zhu et al., 2015).1.2 Research objectivesThe challenges posed by the non-uniqueness of inverse problems and complica-tions posed by remanence form the motivation of my Ph.D. research. The accuracyof potential field inversions depends on our ability to recover the shape of anoma-lies and, in the case of magnetic inversion, to also accurately predict the magne-tization direction. I feel that more knowledge and confidence can be gained from3D potential field inversion if the algorithm can produce a suite of diverse, yetgeologically reasonable solutions with minimal input required from experts. Theinversion framework should allow for sharp and smooth geological contacts to fita broad range of geological scenarios supported by the geophysical data. I want10to extract subtle information present in the data to better define the geometry ofgeological bodies. More specifically, I want to answer the following questions:• Can I further generalize the use of sparse norms to generate diverse solutionsand determine dominant patterns?• Can sparse norms be used to constrain the magnetization vector inversion?• Can we improve the edge definition of geological domains?In this thesis, I will attempt to bridge the gap between geological and geophys-ical modelling with three technical innovations. First, I propose to vary sparsityassumptions imposed on the model and model gradients in a semi-automated fash-ion such that details about sparsity parameters can vary locally without the directinput from the user. At the core of this research is the implementation of a mixed`p-norms regularization. This algorithm will allow me to generate a suite of possi-ble candidates representing a broad range of characteristics. Secondly, I provide alearning strategy to extract dominant features from the solution space. This aspectis important to help practitioners in their interpretation. Thirdly, I improve the con-vergence of the non-linear MVI in spherical coordinates. When my improved MVIis combined with sparse norms, I can recover magnetized bodies with well-definededges and coherent magnetization orientation. The different components are testedon synthetic examples and finally implemented in the Kevitsa case study.1.3 Thesis arrangementThe chapters in this thesis are meant to build on each other. I progressively addthe tools needed for my analysis of the Kevitsa case study in Chapter 7. A numberof synthetic examples are provided along the way to showcase specific aspects ofinversion. In Chapter 2, I provide theoretical background information related tothe forward calculation of potential field data. I provide efficiency improvementsover conventional codes by using a nested OcTree mesh-decoupling strategy andout-of-core sensitivity storage.Chapter 3 introduces the theory needed for potential field inversion. I providenumerical details for the implementation of mixed `p-norm inversions. This chap-ter is an improvement of the work presented in my MSc thesis (Fournier et al.,112016). I use a 1D numerical example to justify changes in the algorithm. I show-case the flexibility of the algorithm on a synthetic 3D gravity problem.In Chapter 4, I review the work of Lelievre and Oldenburg (2009) on magneticvector inversion and provide improvements to the non-linear formulation in spheri-cal coordinates by using an iterative sensitivity re-weighting. Sparsity assumptionsintroduced in Chapter 3 are used to constrain the inversion for the recovery of com-pact magnetic anomalies with coherent magnetization orientation.Chapter 5 brings together advances in mixed norm inversion and the directionalobjective function of Li and Oldenburg (2000). I introduce a 7-cell gradient oper-ator to improve the symmetry of rotated bodies in 3D. I showcase the capability ofthe rotated norms for imaging a folded layer. Surface dip and strike data are inter-polated and used to enforce directionality for the recovery of a continuous layer in3D.In Chapter 6, I make inroads in exploring the model space with mixed `p-norminversions. Building upon Chapter 3, I propose a learning strategy to extract dom-inant features from a suite of models with different characteristics. I build an aver-age model with Principal Component Analysis (PCA) and extract local parameterswith edge detection and image moment algorithms. I extend the work presented inChapter 4 to estimate the dip and strike and geological units.In Chapter 7, I bring together all these components to process gravity and mag-netic data acquired over the Kevitsa deposit, Finland. I provide a suite of densityand magnetization models to characterize the 2.0 Ga old ultra-mafic intrusion andhost stratigraphy. I incorporate geological trends in the modelling using rotationangles derived from a seismic reflection profile and surface structural data. Anensemble of models is analyzed in relation to borehole physical property measure-ments and conceptual understanding of the deposit. I infer tectonic deformationfrom the recovered magnetization vector.I conclude the thesis with a summary chapter assessing the success of the pro-posed methodology, identifying weaknesses and providing avenues for future re-search.12Chapter 2Forward modeling of potentialfield dataMy goal is to characterize the sub-surface in terms of density and magnetization.In Chapter 1, I have introduced the basic physics describing the relation betweendensity and magnetization and their respective anomalous fields in integral form.In this chapter, I provide numerical details about the transformation from the con-tinuous space to a discrete linear system of equations. I also provide efficiencyimprovements over conventional implementations for the processing of large scaledata sets.2.1 Discrete systemsThe usual strategy is to represent the continuous Earth in terms of unit elementseach contributing to the total response observed at a given position in space. For thegravity problem, the integral in (1.3) can be evaluated analytically over a discreteprism with uniform density ρ . As derived by Nagy (1966), the integration givesrise to a linear system:g =TxTyTz ρ , (2.1)13where T relates the cell-centered density value to the observed gravity field g atsome location P(xP,yP,zP):Tx =−G(arctandy dzr dx+ log [ dz+ r ]+ log [ dy+ r ])∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTy =−G(arctandx dzr dy+ log [ dz+ r ]+ log [ dx+ r ])∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTz =−G(arctandy dyr dz+ log [ dy+ r ]+ log [ dy+ r ])∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLr = (dx2+dy2+dz2)1/2dx = (xP− x), dy = (yP− y), dz = (zP− z)(2.2)and G is Newton’s gravitational constant. Parameters needed to define the positionand shape of a unit cell are presented in Figure 2.1. Only the lower southwestL(xL,yL,zL) and upper northeast U(xU ,yU ,zU) corner coordinates are needed todefine the relative distance r(dx, dy, dz) between the observation point and thenodal limits. In this research I use right-handed Cartesian coordinate system suchthat the xˆ, yˆ and zˆ coordinate axes point along the east, north and vertical (up)direction respectively.Figure 2.1: Parameters describing the spatial relationship between an observa-tion point P and a rectangular prism C.14For most gravity field experiments only the vertical component of field gz ismeasured, such that (2.1) reduces togz = Tz ρ (2.3)Equation (2.3) defines the gravity response of a single rectangular prism as ob-served at a single position in space. I can augment equation (2.3) to describe agravity experiment conducted over a large volume of earth and at many observa-tion stationsgpre = Gρ (2.4)such that the linear forward operator G∈RN×M maps the contribution of M numberof prisms (ρ ∈RM), each contributing to the response measured over N observationlocations (gpre ∈ RN). There are many ways to organize the cells making up thisdiscrete model. In all the work presented in this thesis, I use an Octree-baseddiscretization. More details regarding this choice of parameterization are providedin the following section.Similarly for the magnetic response, the integral equation in (1.7) can be eval-uated analytically for a single prism (Sharma, 1966). This gives rise to a linearsystem of the formb = Tm , (2.5)where T is a dense 3-by-3 symmetric matrix describing the linear relation betweena prism with magnetization m = [Mx, My, Mz]T to the components of the fieldb = [bx, by, bz]T .T =µ04piTxx Txy TxzTxy Tyy TyzTxz Tyz Tzz (2.6)where µ0 is the magnetic permeability of free-space. It is important to note thatthe tensor T is a symmetric matrix with zero trace. Therefore only five of the nine15tensor components need to be calculated.Txx =−arctan dx dydx2+ r dz+dz2∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTyy =−arctan dx dydy2+ r dz+dz2∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTxy = log [ dz+ r ]∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTxz = log [ dy+ r ]∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzLTyz = log [ dx+ r ]∣∣∣∣xUxL∣∣∣∣yUyL∣∣∣∣zUzL(2.7)For most geophysical applications, we do not measure the vector field~b, butrather the Total Magnetic Intensity (TMI) of the field that includes both the ge-omagnetic and secondary fields. Since we are only interested in the anomalousresponse from rocks, the approximation is generally made thatbT MA =~b · Hˆ0−µ0‖~H0‖; (2.8)such that the Total Magnetic Anomaly bT MA is assumed to be small and parallel toEarth’s field ~H0. The simulated magnetic datum in (2.5) simplifies tobpre =[Hˆ>0 ·T]m (2.9)Just as for the gravity experiments, this system can be augmented for M number ofprisms (m ∈ R3M) and N observation locations (bpre ∈ RN)bpre = F m (2.10)such that F ∈RN×3M. This linear system has three times the number of parameterscompared to the gravity problem as magnetization is a vector property. As partof my contribution to the open-source community, both the gravity and magnetickernels have been added to the SimPEG.PF library.162.1.1 Choice of discretizationI have so far established the linear equations (2.4) and (2.10) that map the gravityand magnetic response for a collection of rectangular prisms making up a discreteEarth. I have yet to define how these cells are organized in a 3D space. Thisdecision will directly affect the size of the forward (and later inverse) calculationsas the size of the linear operators scale linearly with respect M. Approximating theEarth more efficiently will allow me to process large data sets.The usual strategy is to define a core region of interest with a fine discretiza-tion and surround it by coarser cells (padding) to absorb regional signals that maybe present in the data. Two options are available to organize rectangular prisms.The simplest implementation uses a regular grid, or tensor mesh shown in Fig-ure 2.2(a). Each unit element shares a face with 6 neighbours. Changes in cellsize propagate throughout the domain along the orthogonal direction. The use oftensor meshes has dominated the inversion literature due to the ease of storing andviewing uniformly gridded models.Figure 2.2: Simple representation of (a) tensor and (b) Octree meshes for theorganization of rectangular prisms within a core domain (red) and padding regionover a square domain 4 x 4 m in width. The Tensor mesh can fill the space with64 cells, compared to only 40 cells with the Octree mesh.17In an Octree mesh, unit cells are organized in a hierarchical structure as shownin Figure 2.2(b). The resolution of the grid is increased by dividing a parent cellinto 8 children (or 4 in 2D). This type of discretization offers the most flexibil-ity for increasing the resolution of the mesh in a specific region without affectingthe resolution of boundary cells. The main challenge is in generating a mesh thathonours the geometry of the problem in terms of data location, topography and ge-ological contacts. Over the course of my research, I contributed to the open-sourcecommunity with a suite of refinement functions to facilitate the creation process ofOctree meshes. As part of the SimPEG.discretize library, I implemented thefollowing three strategies:• Box refinement includes all cells intersected by a rectangular box containingthe input points.• Radial refinement is performed inside spheres centered on each input points.The radial distance is determined by the user.• Surface refinement is defined by a continuous Delaunay triangulation of theinput points (Barber et al., 1996). The Octree refinement is determined basedon the vertical distance between cell centers and the nearest triangle.Figure 2.3 compares the three refinement methods for the discretization of points(red) placed on a Gaussian surface. The box refinement is the simplest but alsothe least efficient strategy as it yields a uniform grid similar to the Tensor dis-cretization. For the radial refinement, I end up with a small number of cellsconcentrated around the input points. It is an optimal refinement for scattered ob-servation points. The surface is well suited to describe continuous features suchas topography and geological contacts. It gives, in this case, the most accuraterepresentation of the Gaussian surface.Forward simulation testI demonstrate the benefit of an Octree discretization by forward modelling the re-sponse of a 1 m sphere shown in Figure 2.4. I want to compare the numerical costto perform forward simulations using the standard Tensor discretization and an Oc-tree mesh with surface refinement. Octree refinement uses scatter points placed18Figure 2.3: Three refinement strategies used to discretize a Gaussian curve de-fined by scattered points (red): (a) box, (b) radial and (c) surface refine-ment strategy from the SimPEG.discretize library. Cells are coloured bytheir corresponding Octree level.on the outer surface of the sphere. Using equation (2.5), I calculate the magneticresponse over an 11-by-11 grid of observations placed 1 m above the anomaly. I re-peat the forward simulation over a range of cell sizes. Only cells inside the sphereare considered. The forward calculations are compared to the analytical verticalresponse (banaz ) for a vertical magnetization:banaz =µ04pi[−Mr3+3(M · r)rr5]· zˆ (2.11)where M = [0 xˆ ,0 yˆ, 1 zˆ] and r defines the vector between the center of the sphereand the observation location. Table 2.1 summarizes the forward calculations interms of total number of cells, run time and data residual between the simulationand the analytic solution.φd =121∑i=1(bi−banai )2 (2.12)All calculations were performed on a single thread 2.4 GHz Intel processor.For the Tensor mesh, a reduction in core cell size increases the total number ofcells by a factor 8. For the same reduction in cell size using the surface refinementincreases the number of cells by a factor 4. This is anticipated as only cells atthe surface of the sphere decrease in size compared to a full volume refinement in19Figure 2.4: (a) Discretization of a sphere defined by discrete points (red) us-ing the conventional Tensor mesh with a core region and padding cells. Octreemeshes refined by (b) box, (c) radial and (d) surface methods from theSimPEG.discretize library.the Tensor mesh. Both discretizations can reproduce the analytical response withroughly the same accuracy. Including padding cells to this problem would furtherincrease the efficiency gap between the two discretization methods as the Octreemesh can rapidly increase the cell size with little influence from the discretizationin the core region.20Tensor Octree (Surface)Cell size (m) # Cells Time (s) φd # Cells Time (s) φd1.00e-01 544 1.8e-01 4.2e+02 496 1.6e-01 7.7e+025.00e-02 4196 9.6e-01 4.7e+01 2516 5.1e-01 7.1e+012.50e-02 33478 6.5e+00 4.2e+00 10836 2.1e+00 5.2e+001.25e-02 268080 6.5e+01 2.0e+00 47276 1.0e+01 2.7e+00Table 2.1: Summary table for the forward modeling of a magnetized sphere us-ing a Tensor and Octree discretization2.2 Large scale problemsThe amount of memory needed to store the dense linear operators for gravity andmagnetics is a limiting factor for the simulation of potential field data with theintegral formulation. As prescribed in (2.10), the size of the problem is linearlyscaled by the number of model parameters M and the number of the data N. Formoderate size problems encountered in exploration geophysics (M≈ 106, N ≈ 104)the memory requirement to store a dense M×N matrix can exceed hundreds ofgigabytes. While modern supercomputers can handle problems of this size, it isstill out of reach for common desktop computers.A number of strategies have been proposed in the past to reduce the size ofthe problem. Compression methods, either in the Fourier (Pilkington, 1997) orwavelet domain (Li and Oldenburg, 2003), have proven successful in reducing thesize of the forward operators. On the downside, compression methods are proneto introducing modelling artifacts, which can be hard to differentiate from trueanomalies. This is especially an issue for the magnetic vector inversion exploredin Chapter 3.Another group of methods takes advantage of the rapid decay of geophysicalsignals to reduce the size of individual forward simulations. The assumption ismade that model parameters in the far-field of measurements contribute little to thetotal response and can, therefore, be approximated with fewer parameters. In thisregard, the concept is loosely related to the Fast Multipole Method (Engheta et al.,1992). This is especially valid for airborne EM experiments where cells beyondroughly 10 times the flight height of airborne EM systems have a negligible impacton the measured response (Reid and Pfaffling, 2006). In its simplest implemen-21tation, Cox et al. (2010) define a footprint approach such that model parametersoutside a pre-defined tolerance are simply ignored.More recently, Yang and Oldenburg (2014) introduced a mesh decoupling ap-proach. Local meshes are designed based on the source location and decay rateof the geophysical signal. Cell-centered conductivity values are homogenized us-ing volumetric-weighted averaging. Similarly, Haber and Schwartzbach (2014)designed a nested Octree mesh strategy. The core region of local meshes is co-located with cells in the global mesh and cover the same lateral extent. Physicalproperty values are homogenized from small cells in the global mesh to larger (orequal) cells in the local meshes. This strategy assures that far-field features can stillcontribute to the total response. The accuracy of individual forward simulationsdepends primarily on the interpolation scheme used to transfer physical propertiesfrom the global to the local meshes.In this research, I employ the methodology of Haber and Schwartzbach (2014).I want to divide the full data set into subsets (or tiles), each associated with a localmesh. Physical property values are transferred to local meshes using a volumet-ric weighted average. In order to further minimize the memory footprint of thesimulation, I make use of the parallel Dask package (Dask, 2016). The librarywas designed by the data science community for out-of-core processing of largedatasets (Rocklin, 2015). The zarr file format allows the storage of dense arrayson a solid-state drive (SSD) in compressed memory chunks. Out-of-core storageis appealing as it reduces the amount of Random-Access Memory (RAM) neededto store the large forward operators. Read and write operations are done in par-allel. The runtime depends largely on the hardware used in terms of processorspeed, number of processors and communication speed between the workers andthe solid-state memory. The size of individual memory chunks can be adjusted tooptimize the processing speed depending on the resources available.In the case of potential field data, the source of geophysical signals can span abroad range of wavelengths, from tens of meters (local) to hundreds of kilometres(regional). The memory footprint of a problem can be reduced considerably if wemanage to approximate far-field features with fewer cells. In the most extremecase, each observation point could have its own local mesh. Choosing the numberof tiles becomes a trade-off between accurately capture the heterogeneity of the22model for forward modeling while minimizing traffic between the parallel workersand the SSD memory.2.2.1 Numerical testI investigate this trade-off between accuracy and efficiency with a numerical exam-ple. I create a synthetic gravity model shown in Figure 2.5. The model containsboth short wavelength information from a small block anomaly and long wave-length signal from large domains in the north and southeast quadrants. From equa-tion (2.4), I compute gravity data on 21×21 grid placed 1 m above a flat topogra-phy. Forward simulation of the global model required 4.3 Mb of memory and wasperformed in parallel in 6.7 s with 4 processors (2.4 GHz). It is a relatively smallexample compared to industry standards, but similar trends are to be expected onlarger scale problems.I will assess the trade-off between efficiency and numerical accuracy by per-forming a series of five forward simulations using a range of tiling patterns. I breakup the dataset into 2, 4, 9, 12 and 16 square tiles. For each tiling experiment, I cal-culate the total memory footprint (Gb) and the computation time (s) to perform theforward simulations. I also compute the data residual between the global simula-tion (single mesh) and the tiled simulation using the `2-norm measure in equation(2.12).As seen in Figure 2.7, the memory needed to perform the forward calculationsrapidly decreases but eventually levels-off as the number of tiles increases. Thisis expected as local meshes necessitate a minimum number of cells to fill out theglobal domain and a minimum number of small padding cells near the edge of thetiles to reduce interpolation artifacts. The reduction in problem size is inverselycorrelated with an increase in computational time as communication between theworkers and the SSD memory becomes a bottleneck.I also note an increase in data residual as a function of tile size. Figure 2.8 com-pares the simulated data from the global mesh (single tile) to the combined forwardsimulation calculated with 12 tiles (last experiment). From the residual map, it ispossible to distinguish short wavelength discrepancies between adjacent forwardsimulations (tiles). These residuals are primarily due to the homogenization of23Figure 2.5: (a) Horizontal and (b) vertical sections through the synthetic densitymodel used to test the mesh decoupling strategy. (c) The simulated data containshort and long wavelength information.anomalies near the edges of tiles, which in turn is a function of the wavelengthinformation contained in the data. As a first pass, I establish experimentally theappropriate padding distance based on the estimated data uncertainties, such thatthe maximum residual falls below the experimental error. For this experiment, themeshing artifacts are at most 2% of the data amplitude (or 0.004 mGal), whichI achieved with a minimum padding distance of 4 cells per Octree level. Deter-mining an optimal padding distance as a function of the geophysical signal wouldwarrant further research.24Figure 2.6: Example of mesh decoupling of the forward problem into nestedOctree tiles. (Middle) Local meshes are generated for a pre-defined number oftiles nested inside the global domain. (Bottom) The forward modelling operatorfor each tile is stored in compressed Zarr file format on solid-state memory.Memory chunks (grey) are accessed in parallel by the Dask library to performthe forward calculations.25Figure 2.7: Trade-off curves between the size of the forward problem (blue) anddata residual (red) over a range of number of tiles. The total size of the problemis calculated based on the sum of cells in all the local meshes times the numberof data. The optimal number of tiles would be at the point of intersection whereboth the cost of forward calculations and the data residual change significantly.Figure 2.8: Simulated gravity data calculated from (a) the global model and (b)the 12 forward tiled calculations. (c) Data residuals show short wavelength dis-crepancies between adjacent tiles due to interpolation effects.26Chapter 3Inverse ProblemIn Chapter 2, I defined the linear equations relating density and magnetization togravity and magnetic data. I now review the theory needed to solve the inverseproblem, such that I can recover a 3D representation of the subsurface from theobserved data. A key issue related to the inverse problem is that there is an infinitenumber of possible models that can satisfy the data. Field measurements are gener-ally acquired from the surface resulting in the inverse problem to be ill-posed. Thepresence of experimental noise further complicates the problem. To circumventthese issues, the inversion is often formulated as an optimization problem of formminmφ(m) = φd +βφmsubject to φd ≤ φ ∗d .(3.1)where φd is the misfit functionφd =N∑i=1(dpredi −dobsiσi)2, (3.2)that measures the residuals between the observed and predicted data dpre, normal-ized by the estimated data uncertainties σ .The regularization function φm, or model objective function, serves as a vehicleto introduce a priori information in the inversion. Several regularization strategieshave been developed over the last decades such that the solution remains geologi-27cally plausible. I focus on the generic `p-norm regularization of the formφm = ∑r=s,x,y,zαr∫Vw(m)| fr(m)|p j dV . (3.3)The functions fr can take many forms but most often have beenfs = m−mre f , fx = dmdx , fy =dmdy, fz =dmdz. (3.4)Thus fs(m) measures the deviation from a reference model mre f and fx(m), fy(m)and fz(m) measure the roughness of the model m along orthogonal directions in3D. This optimization problem has multiple terms scaled by hyper-parameters. Thefirst parameter is β which controls the balance between misfit and regularization.It is assumed that a value can be found such that the target misfit is reached. Theα’s are constants that control the relative influence of the different regularizationfunctions. A larger α-value increases the focus of the optimization on the corre-sponding penalty function. User-defined weights w(m) are used to incorporate anytype of a priori information that may be available to guide the solution.Most often the `2-norms measure has been used giving rise to a discrete linearsystem of formφm = αsφs+αxφx+αyφy+αzφz= ∑r=s,x,y,zαr‖Wr Vr Gr (m−mre f )‖22 , (3.5)where φs measures the deviation of the discrete model m from a reference modelmre f and φx, φy and φz measure the roughness of the model along Cartesian di-rections. The reference model is sometimes omitted in the roughness terms. Thematrices Gx, Gy, and Gz are discrete gradient operators. For the smallness compo-nent, Gs reduces to the identity matrix. Volumes of integration resulting from theevaluation of (3.3) are applied through diagonal matrices such thatVs = diag[v1/2]. (3.6)where v holds the discrete volume elements corresponding to each cell. For the28model derivative terms, the volume of integration is computed over cell interfacessuch thatVx = diag[(AFxC v)1/2]. (3.7)where the matrix AFxC averages the cell-centered volumes to cell faces. Similaraveraging is performed along the orthogonal y and z-direction. Diagonal matricesWr hold user-defined weights. More details about these weights are provided in thefollowing section. Lastly, the α parameters control the relative importance given toindividual components of the regularization. What is sometimes sought, at least asa first pass, is that if φs, φx, φy, φz have about the same numerical value then they arecontributing equally. Dimensional analysis shows that for a uniform discretizationh:[φx][φs]= [h]−2 . (3.8)The common approach is to set αs accordingly in order to scale the components ofthe regularization function.The usual strategy to solve (3.1) is through a gradient descent algorithm, suchas a Gauss-Newton approach, where we attempt to find a solution that has zerogradientsg = ∇mφ(m) = ∇mφd +β[αs∇mφs+αx∇mφx+αy∇mφy+αz∇mφz]= 0 . (3.9)where ∇m stands for the partial derivatives of the function with respect to the dis-crete parameterization m. A solution to (3.9) can readily be calculated by gradientdescent methods (Hestenes and Stiefel, 1952; Nocedal and Wright, 1999).A large number of studies have made use of this formulation to incorporate avariety of a priori information: physical property data from rock and core samples(Lelie`vre et al., 2009), structural knowledge (Lelie`vre, 2009; Li and Oldenburg,2000) and advanced 3D geological modeling ((Bosch and McGaughey, 2001; Ful-lagar et al., 2008; Phillips, 1996; Williams, 2008). Although successful in identi-fying imaging anomalies at depth, penalty functions that rely on `2-norm measureshave a limited range of possible outcomes. The models tend to be smooth and diffi-cult to interpret in relation to known geological domains with discrete boundaries.29Moreover, substantial modelling work is generally required by experts to manuallyrefine these constraints in order to test different geological scenarios.Sensitivity weightingThe weighting matrices Wr introduced in (3.5) can take many forms dependingon the type of a priori information that may be available. For potential fieldsproblems, a sensitivity weighting function is generally used to counteract the rapiddecay of the geophysical signal as a function of distance. In the work of Li andOldenburg (1996), a distance weighting approximation is employed and fixed atthe onset. In this thesis, I resort to an iterative re-weighting strategy based on thesensitivity of a given inverse problemJ =∂ F [m]∂m, (3.10)where J, also referred to as the Jacobian matrix, holds the partial directives of theforward problem F [m] with respect to m. Adapted from Haber et al. (1997), Iformulate the sensitivity-based weighting function:Ws = diag[[wmax(w)]1/2]w j =[N∑i=1Ji j2+δ]1/2/v j ,(3.11)where w measures the sum square of the columns of the Jacobian, normalized bythe cell volumes. The constant δ is a small value (near machine precision) addedto avoid singularity. The weights are normalized by the maximum value such thatthe range of weights are bounded between [0, 1]. The same sensitivity weights canalso be applied to the model derivative terms using a cell averaging operation suchthatWx = diag[(AFxC[wmax(w)])1/2], (3.12)similar to the averaging of volumes used in (3.7). This sensitivity weighting strat-egy is general and adaptable to any inverse problems where the sensitivity matrix30can be calculated explicitly. While the initial purpose of the sensitivity weightingfunction of Li and Oldenburg (1996) is to simply counteract the decay of potentialfields, I will show numerically in Chapter 4 that the iterative re-scaling processcan also be beneficial in improving the convergence of gradient methods applied tonon-linear inverse problems.3.0.1 Synthetic gravity exampleAs an entry point to the inverse problem, I proceed with a simple synthetic gravityexample. I define a volume of interest 600 m wide by 300 m deep, over which Iplace a uniform survey grid of 21 x 21 stations placed 5 m above a flat topography.The core region directly below the survey grid is discretized at a 5 m resolution asshown in Figure 3.1(a) Within the core region, I build a geophysical target made upof a single dense cube, 25 m in width. I set the density contrast of the prism to 0.2g/cc in a uniform zero background. From (2.3), I simulate the vertical gravity fieldof the block and add random Gaussian noise with 10−3 mGal standard deviation.Figures 3.1(b) and (c) display the simulated data and ‘noisy’ observations gobsz usedin the inversion. I revisit this example in Chapter 4 to demonstrate the inversionprocess on magnetic data.From the noisy data I will attempt to recover the block anomaly by the inverseprocess. The objective function to be minimized takes in this case the form:minmφ(m) = ‖G ρ−dobs‖22+β ∑r=s,x,y,zαr‖WrVr Gr ρ‖22subject to φd ≤ φ ∗d(3.13)where I set ρre f = 0. Since (3.13) is linear with respect to the density contrastmodel ρ , I can solve it uniquely for a fixed trade-off parameter β . I repeat thisprocess for variable β values until a find a solution that satisfies φd ≈ N. Fig-ure 3.2(a) presents a vertical section through the recovered density model. Thedensity anomaly is imaged at roughly the right position, but the edges of the blockare poorly defined. As normally obtained with `2-norm penalties, the solution issmooth and density values remain near the zero reference model. Hence the needto explore other regularization functions that can better resolve compact objects.31Figure 3.1: (a) Vertical section through a 25 m cube with density ρ=0.2 g/ccplaced in a uniform zero density background. (b) Simulated gravity data re-sponses on a 21× 21 survey grid placed 5 m above the flat topography. (c)Gravity data with random Gaussian noise added, 10−3 mGal standard deviation.32Figure 3.2: (a) Vertical section through the inverted density model using theconventional `2-norm regularization, (b) predicted and (c) normalized data resid-ual. Outline of the true model (red) is shown for reference.333.1 General `p-norm regularizationAlternatively, researchers have explored the use of non-`2 measures to promote therecovery of compact anomalies. Approximations to `1-norm such as the Hubernorm (Huber, 1964)∑i|mi|p ≈∑i{m2i , |mi| ≤ ε,2ε|mi|− ε2, |mi|> ε,}and the Ekblom norm (Ekblom, 1973):∑i|mi|p ≈∑i(m2i + ε2)p/2(3.14)have received considerable attention in geophysical inversion and signal processing(Daubechies et al., 2010; Farquharson and Oldenburg, 1998; Gorodnitsky and Rao,1997; Li, 1993; Sun and Li, 2014). Likewise, the Lawson’s measure (Lawson,1961)∑i|mi|p ≈∑imi2(mi2+ ε2)1−p/2, (3.15)has been proposed to approximate `0-norm and it has proven useful in generatingminimum support models. This formulation has received considerable attentionin the literature. (Ajo-Franklin et al., 2007; Barbosa and Silva, 1994; Last andKubik, 1983; Portniaguine, 1999). Figure 3.3 compares the `p-norms with theLawson approximation over a range of model values. As ε→ 0, the approximationapproaches the `p-norm on the complete interval p ∈ [0,2]. While (3.15) would intheory permit us to explore a wide range of solutions for 0 ≤ p ≤ 2, its numericalimplementation remains challenging. Most algorithms have been limited to the `0,`1, and `2-norm measure applied evenly to all components of the model objectivefunction.Recent efforts by Sun and Li (2014) have shown promise in further exploringthe model space by varying `p-norm measures locally. They divided the inversiondomain into regions reacting favourably to either the `1 or `2-norm regularization.The automated process could adapt to complex geological scenarios where both34Figure 3.3: Approximated `p-norm using the Lawson measure (Lawson, 1961)over a range of p-values and for a fixed threshold parameter ε = 10−1.smooth and blocky anomalies are present. Building upon the work I introduced inmy M.Sc. Thesis (Fournier, 2015), I want to extend the work of Sun and Li (2014)and further generalized the mixed norm inversion for p ∈ [0 2].3.1.1 Synthetic 1D problemTo develop my methodology it suffices to work with a simple test example. InFigure 3.4(a) I present a synthetic 1D model made up of a boxcar anomaly. Theregion is divided into 50 uniform cells distributed along the interval [0≤ x ≤ 1]. I35Figure 3.4: Linear forward problem made up of: (a) an example kernel function; (b) model; (c) observed data with assigned standard errors.define a synthetic geophysical experiment such that the data (dobs) aredobs = F mtrue+ e , (3.16)The kernel coefficients Fi j are sampled from a standard normal distribution of pos-itive values multiplied by the discretization intervals. Choosing a stochastic kernelfunction for a linear inverse problem is unusual. Smooth functions are usually em-ployed (polynomials, decaying exponentials, sinusoidals), but my choice will serveto highlight the effects of various regularization functions. I generate 10 data, soF ∈ RN×M where M = 50 and N = 10. Random Gaussian error e (σ=0.025) isadded to simulate noise (Fig. 3.4(c)).To begin my analysis, I invert my synthetic dataset with two simple regulariza-tion functions. For this 1D problem, the objective function takes the form:minmφ(m) = ‖F m−dobs‖22+β ∑r=s,xαr‖WrVr Gr m‖22subject to φd ≤ φ ∗d(3.17)To simplify the analysis, I set all weighting terms to unity (Wr = I) and the ref-erence model to zero. Figure 3.5(a) presents the recovered model after reachingthe target misfit (φ ∗d = N) using the smallness term alone (αx = 0). The solutionexhibits high variability similar to the stochastic kernel function, but model param-36Figure 3.5: Solution to the 1D inverse problem using (a) an `2-norm on themodel (αx = 0), (b) the `2-norm on model gradients (αs = 0) and (c) combinedregularization function (αs = 2500, αx = 1). (d) Convergence curve comparingthe misfit (φd) and the regularization (φm) as a function of iterations. (e) Compar-ative plot for the relative contribution of the different components of the objectivefunction measured in terms of maximum absolute gradient (‖gi‖∞)eters remain near the implied zero reference value. Next, I invert the data using themodel gradient term (αs = 0); this yields the smoother model presented in 3.5(b).The solution shows less spatial variability and the horizontal position of the boxcaranomaly is better located.Next, I combine both regularization functions so that the solution remains closeto the implied zero reference value and it is smooth. I need to determine the lengthscale weighting proposed in (3.8). For my problem h = 0.02 and hence following37the relationship established in (3.8) I set αx = 1 and αs = 2500. Inverting withthese parameter values yields the model in 3.5(c). Convergence curves presentedin Figure 3.5(d) show the evolution of φd and φm as a function of iteration. Asβ decreases (shown in Figure 3.5(e)), the misfit, φd progressively decreases whilemodel complexity, indicated by φm, progressively increases.Visually, the solution 3.5(c) exhibits characteristics of remaining near the zeroreference value while also attempting to be smooth. Numerical evaluation of thetwo components of the regularization function, presented in Table 3.1, show thatαsφs = 3.31 and αxφx = 1.13. This might suggest that both φs and φx are roughlyequal in importance.Rather than working with global norms, in this study, I propose to quantifythe relative importance of the terms in the regularization function based on theirpartial derivatives, or gradients. From (3.9) I expect to find an optimal solutionwhere the sum of the gradients vanishes, either because all components are equalto zero, or because multiple gradients have opposite signs. To quantify the size ofthe gradients I use the infinity-norm‖gr‖∞ = ‖∇mφr‖∞ , (3.18)corresponding to the maximum absolute value of gr. The ‖gr‖∞ metric is appealingfor a few reasons: (a) it is directly linked to the minimization process because Iuse gradient descent methods, (b) it does not depend on the dimension M of theparameter space as do other measures that involve a sum of components of thevector, (c) the theoretical maximum can be calculated analytically for any given`p-norm function. These properties will become useful in the following sectionwhen I attempt to balance different norm penalties applied on a cell-by-cell basis.Figure 3.5(e) compares ‖gd‖∞, αs‖gs‖∞ and αx‖gx‖∞ over the iterative process.I note that, under the current α-scaling strategy proposed in (3.8), the individualpartial derivatives for φs and φx also appear to be proportional in magnitude. Toquantify this I define a proportionality ratio:λ∞ =αs ‖gs‖∞αx ‖gx‖∞(3.19)38αs αx αs φs αx φx λ∞2500 1 3.31 1.13 1.23Table 3.1: Norm values and proportionality ratio obtained for the 1D solutionpresented in Figure 3.5(c). A proportionality ratio of λ∞ ≈ 1 indicates that thecomponents of the regularization function are both contributing significantly tothe final solution.I shall use λ∞ as an indicator to evaluate the relative influence of (any) two terms inthe regularization function. For my example λ∞ = 1.23, from which I infer that φsand φx are contributing nearly equally to the solution (Table 3.1). As I further gen-eralize the regularization function for arbitrary `p-norm measures, I will attempt tomaintain this proportionality ratio (λ∞ ≈ 1) between competing functions so thatmy modeling objectives are preserved throughout the inversion process.3.1.2 Iterative Re-weighted Least Squares algorithmSolutions obtained with `2-norm regularization functions provided some insightabout the sought model but better representations can be obtained by employinggeneral `p-norms:φ ps =∑i|mi|p (3.20)My main focus is in the regularization function in (3.3) which I approximate withthe Lawson norm such thatφm = ∑r=s,xαr∫Vw(m)fr(m)2(fr(m)2+ ε2)1−pr/2 dV , (3.21)This measure makes the inverse problem non-linear with respect to the model. Thecommon strategy is to solve the inverse problem through an Iterative ReweightedLeast-Squares (IRLS) approach such that (3.21) is expressed as a weighted least-squares problem. The denominator is evaluated for model parameters obtainedfrom the most recent iteration such thatφ prr =M∑i=1wivifri(m)2[( fri(m(k−1)))2+ ε2]1−pr/2 (3.22)39where m(k−1)i are model parameters obtained at a previous iteration. The integralcorresponding to the smallest model component can be written as:φ pss =M∑i=1wsivsim2i((m(k−1)i )2+ ε2)1−ps/2 (3.23)In (3.23) I have explicitly written the objective function as φ pss to indicate that Iam evaluating a smallest model component with an `p-norm with p = ps. Thisapproximation of the `p-norm can be implemented within the same least-squaresframework used in (3.5) such that:φ pss = ‖WsVs Rs m‖22 . (3.24)where the IRLS weights Rs are defined asRs = diag [rs]1/2rsi =((mi(k−1))2+ ε2)ps/2−1.(3.25)Carrying out the same procedure on the measure of model derivatives yieldsφ pxx =M−1∑i=1wxivxi(mi+1−mihi)2[(m(k−1)i+1 −m(k−1)ihi)2+ ε2]1−px/2 (3.26)where hi defines the cell-center distance between neighboring model parameters.Equation (3.26) can also be expressed in linear form asφ pxx = ‖WxVx Rx Gx m‖22 , (3.27)where the gradient operator and the corresponding IRLS weights are calculated byGx =−h−11 h−11 0 . . . 00. . . . . . . . ........ . . 0 −h−1M−1 h−1M−1 . (3.28)40andRx = diag [rx]1/2rxi =(m(k−1)i+1 −m(k−1)ihi)2+ ε2px/2−1 . (3.29)respectively. The final regularization function is thusφ pm = αs‖WsVs Rs m‖22+αx‖WxVx Rx Gx m‖22 . (3.30)The core IRLS procedure described in Table 3.2 involves two main stages:1. Stage 1 solves the inverse problem using `2-norms presented in (3.5). Theassumption is made that the globally convex `2-norm regularized inversionis a good approximation of the true solution and it is used to form the initialIRLS weights defined in (3.25). The β parameter is controlled by a coolingschedule that starts with a high value and is successively decreased untilφd ≈ φ ∗d .2. Stage 2 starts from the solution obtained in Stage 1 and solves the inverseproblem iteratively using the regularization in (3.30) and a standard Gauss-Newton procedure. A gradient descent direction δm is found by solvingH δm = g (3.31)where H is the approximate Hessian and g is the gradient of the objectivefunction. I use the Conjugate Gradient method (Hestenes and Stiefel, 1952)to solve this system.The model update at the kth iteration ism = m(k−1)+αδm (3.32)where the step length α is found by a line-search back-stepping method (Nocedaland Wright, 1999). Gradient steps are only performed if the data misfit remains41within the user-defined tolerance ηφd .|φd−φ ∗d |φ ∗d≤ ηφd (3.33)If outside the tolerance, the algorithm repeats the Gauss-Newton calculation withthe previous m(k−1) and a different β -value, either lower or higher depending onthe achieved φd . This β -search step is an important component in the workflowwhen the minimization switches between an l2 to an lp objective function becauseφ pm can vary markedly. This can force a change of β by a few orders of magnitudein some cases. Once an appropriate β has been found such that (3.33) is respected,the model update m(k) is accepted and used for the next iteration cycle. The IRLSprocess continues until the change in regularization falls below some pre-definedtolerance ηφm|φ (k−1)m −φ (k)m |φ (k)m< ηφm (3.34)I set to ηφm = 10−5 (0.01% change) in all my experiments. Using the above algo-rithm I now explore specific inversions for a fixed ε = 10−3 and uniform norms,with p = 1 and p = 0, applied on the model and model gradients.3.1.3 Case 1: `1-norm (ps = px = 1)I first explore the end member of convex functions for ps = px = 1 for whichoptimality can be guaranteed (Daubechies et al., 2010; Osborne, 1985). Using theprocedure prescribed in Table (3.2), I invert the 1D problem with three differentregularization functions: (a) l1-norm measure of the model (αx = 0), (b) l1-normmeasure of the model gradients (αs = 0) and (c) for the combined penalties usingαs = 2500, αx = 1, which I previously used for the l2-norm inversion.As shown in Figure 3.6(a), the first inversion is successful in recovering asparse solution. From Linear Programming (LP) theory, the expected optimal so-lution would have as many non-zero parameters as there are linearly independentconstraints or 10 values in this case. For comparison, I solve the LP problem by theSimplex routine from the open-source library Scipy.Optimization.linprog(Jones et al., 2001). Figure 3.6(a) compares both solutions and shows that my im-42Stage 1: Initialization (φ 2m)minmφd +βφ 2ms.t. φd = φ ∗dβ (0), m(0), R(0), φ (0)mStage 2: IRLS (φ pm)while |φ(k−1)m −φ (k)m |φ (k)m> ηφmdo β -Searchk := k+1β (k), m(k), R(k)β -SearchSolve H δm = gm = m(k−1)+αδmif |φd−φ∗d |φ∗d> ηφdadjust β , re-doelsecontinueTable 3.2: IRLS algorithm in pseudo-code made of two stages: Stage 1 Ini-tialization with convex least-squares inversion, Stage 2 IRLS updates with innerβ -search steps.plementation of IRLS for l1-norm yields a solution in close agreement with theSimplex routine. A better approximation could be obtained (not shown here) bylowering the threshold parameter ε . I will examine this aspect of the algorithm inthe following section. Figure 3.6(b) presents the solution for the l1-norm appliedto the model gradients. The final solution is blockier and the general shape of theboxcar model has been improved.Lastly, the solution obtained with the combined l1-norm regularization on themodel and model derivative is shown in Figure 3.6(c); it is similar to that in Fig-ure 3.6(a). This shows that the smallest model component has dominated the so-lution. This is quantified by the evaluated proportionality ratio λ∞ = 50; settingαs = 2500 is too large. To understand this result, I can factor out a base cell length43Figure 3.6: (a) Two solutions using an `1-norm on the model: (blue) Simplex,and (black) IRLS method. (b) Solution obtained with the approximated `1-norm(IRLS) penalty on model gradients alone and (c) with the combined penalty func-tions (αs = 2500, αx = 1). The calculated proportionality ratio λ∞ indicates thatthe combined penalties is dominated by the φ 1s term. (d) Convergence curve and(e) maximum partial derivatives associated with the components of the objectivefunction as a function of iterations for the inversion in (c). The vertical dottedlines indicate the change in regularization from an `2-norm to `1-norm measure.44h from (3.26) such thatφ pxx =M−1∑i=1wxivxih−2(mi+1−mihˆi)2[h−2((m(k−1)i+1 −m(k−1)ihˆi)2+h2ε2)]1−px/2=M−1∑i=1wxivxih−2(mi+1−mihˆi)2hpx−2[(m(k−1)i+1 −m(k−1)ihˆi)2+h2ε2]1−px/2= h−pxM−1∑i=1wxivxi(mi+1−mihˆi)2[(m(k−1)i+1 −m(k−1)ihˆi)2+h2ε2]1−px/2(3.35)where h = min(h) represents the core discretization length and hˆi = hi/h. For auniform grid, hˆi simply reduces to unity everywhere. Comparing this expression to(3.22) clearly shows a difference in scales between φ ps and φ px , previously fixed in(3.8), that now depends on the chosen p-value such that[φ px ][φ ps ]= [h]−px . (3.36)It also highlights a dependency between the chosen threshold parameter ε and thediscretization length h. I will revisit this parameter in later sections.In accordance to this new relationship, I can re-adjust the importance of φ ps bysetting αs = 50, where p= 1 and h= 0.02. After applying this change I recover themodel presented in Figure 3.7(a). The combined assumption of a piece-wise con-tinuous and sparse model yields a solution that closely resembles the true boxcarmodel. The recovery of mtrue has remarkably improved compared to the `2-normsolutions (Fig. 3.5), and this demonstrates the power of customizable objectivefunctions. It is important to notice that the re-adjustment of αs has brought thepartial derivatives of φ pss and φ pxx to a comparable level, with a final proportionalityratio λ∞ = 1.01. Even though I have changed the norms during the inversion, thecontribution of both penalty functions has remained at a comparable level during45Figure 3.7: (a) Solution obtained with the combined penalty functions αsφ 1s +αxφ 1x after re-adjustment of αs = 50, αx = 1. (b) Convergence curve and (c)maximum partial derivatives associated with the components of the objectivefunction as a function of iteration.the transition between Stage 1 and 2 of the algorithm (Fig 3.7(c)).3.1.4 Case 2: `0-norm (ps = px = 0)The main advantage of the IRLS formulation is that it permits, in theory, approx-imating any norm including the non-linear approximation for p < 1. The goal isto potentially recover a model with even fewer non-zero parameters than that ob-tained by solving the problem with p= 1. The IRLS formulation for p= 0 has beenimplemented for various geophysical problems under different names: such as thecompact inversion (Last and Kubik, 1983), minimum support functional (Portni-aguine and Zhdanov, 2002), and others (Ajo-Franklin et al., 2007; Barbosa andSilva, 1994; Blaschek et al., 2008; Chartrand, 2007; Stocco et al., 2009).Following the same IRLS methodology as described in Table 3.2, I invert thesynthetic 1D problem with three assumptions: (a) `0-norm applied on the model(αx = 0), (b) `0 on model gradients (αs = 0) and combined penalties (αs = 1, αx =1). Figure 3.8 presents the solutions for all three cases. I note that in the first case,(a), the approximate `0-norm inversion recovers a sparser solution than obtainedwith the `1-norm; there are only eight non-zeros parameters. Similarly for case(b), I recover a model with fewer changes in model values. Finally in case (c) the46Figure 3.8: Solution to the 1D inverse problem using an approximate `0-norm(a) on the model, (b) on model gradients and (c) combined penalty functionsusing the IRLS algorithm (αs = 1, αx = 1). All three solutions honor the datawithin the target misfit φ ∗d .solution obtained with the combined `0-norm penalties matches almost perfectlythe true boxcar anomaly. The final proportionality ratio as calculated from (3.19)indicates once again a good balance between the penalty functions (λ∞ = 1.13).To summarize this section, I have now recovered nine models using different`p-norm penalties applied on the model and model gradient. All solutions pre-sented in Figure 3.5, 3.6 and 3.8 can reproduce the data within the predefined datatolerance (φ (k)d ≈ N). Without prior knowledge about the true signal, all these so-lutions would be valid candidates to explain the observed geophysical data.3.2 Mixed norm regularizationWhile I was successful in recovering a solution that closely resembles the boxcarmodel, the same penalty functions might not be appropriate for other models, suchas compact targets with smooth edges. I thus explore a broader range of solutionsby using the Lawson approximation in (3.23) for any combination of norms on therange 0≤ p≤ 2.The idea of combining different norm measures for the simultaneous recoveryof smooth and compact features has partially been explored by Sun and Li (2014)on a 2D seismic tomography problem. They demonstrated the benefits of dividingmodel space into regions with different `p-norm penalties. The choice of norms47was limited to be either l1 or l2. Little has been published however on the inde-pendent mixing of model and gradient norms on the range p ∈ [0,2], although thisproblem was initially addressed in (Fournier, 2015). I now apply my algorithm tominimizeφ pm = αsφ0s +αxφ2x . (3.37)where once again the superscript indicates the `p-norm measure used in each func-tion. Based upon my previous work, I expect the solution to be sparse, in termsof non-zero model parameters, while smooth with respect to the model gradients.Unfortunately, following the current IRLS strategy, I recover the model presentedin Figure 3.9(a). The anomaly is concentrated near the boxcar but appears to bedominated by φ 0s . There seems to be only marginal influence from φ 2x . Comparingthe partial derivatives of the objective function confirms this. After convergencethe calculated proportionality ratio is λ∞ = 159. This is a significant change fromthe end of Stage 1 where λ∞ ≈ 1. Clearly, iteration 6, at Stage 2 of the IRLS, tookthe solution away from the proportionality condition (Fig 3.9(c)). I hypothesizethat a more desirable solution could be obtained if proportionality was preservedamong the components of the objective function throughout the IRLS process. Inthe following sections, I provide an important modification to the standard IRLSalgorithm to achieve this goal.3.2.1 Scaled-IRLS stepsSince the inverse problem is solved using gradients of the composite objectivefunction φ(m), the relative magnitude of the individual gradients is a driving forcein controlling the iteration step in (3.31). I want to ensure that each penalty term inthe objective function is playing a significant role. Taking the partial derivatives of(3.24) with respect to m yields:gps =∂φ ps∂m= R>s V>s W>s WsVsRsm (3.38)where I purposely omitted a factor 2 from the differentiations of the `2-norm as itgets absorbed by the zero right-end side of (3.9). From Figure 3.10(a), I note thatthe magnitude of the derivatives increases rapidly for small p values as mi → 0.48Figure 3.9: (a) Recovered model and (b) convergence curves using the conven-tional IRLS method for ps = 0, px = 2 and a fixed threshold parameter ε = 10−3(αs = αx = 1). (c) Trade-off parameter and maximum gradients for the differ-ent components of the objective function. At the start of Stage 2 (iteration 6),the sudden increase in ‖gs‖∞ is matched with a decrease in β . Throughout theinversion, ‖gx‖∞ remains small in magnitude.This trend is accentuated for small ε values as demonstrated in Figure 3.10(b) forp= 0. The magnitude of derivatives for p< 2 increase rapidly as m→ 0 and ε→ 0.This results in gradient steps in equation (3.9) that are dominated by sparse norms.This property of the IRLS approximation is important because, when attempting tocombine different norm penalties within the same objective function, there will be asystematic bias towards small `p-norm penalties. To circumvent this bias I proposethe re-scale the contribution of each regularization function during the iterativeprocess in order to preserve proportionality. I define the following gradient-basedscalingγ =[‖g2‖∞‖gp‖∞]1/2. (3.39)By using this scaling strategy I can equalize the influence of each regularizationfunction. I can evaluate the theoretical maximum of gpr by taking the deriva-tive of (3.21) in terms of f (m), followed by a second derivative after substitutingf (m(k−1)) for f (m) and setting the expression to zero such that( f (m)2+ ε2)p/2−1+(p−2) f (m)2( f (m)2+ ε2)p/2−2 = 0 (3.40)49Figure 3.10: Derivatives of the Lawson approximation over a range of modelvalues for (a) a fixed threshold parameter ε = 10−1 over a range of p values andfor (b) a fixed p = 0 over a range of ε values. (c) Applying the γ-scaling to thegradients brings all maximums to be equal irrespective of p and ε .The maximum gradient of the Lawson approximation occurs at f (m)∗f (m)∗ =∞ or f (m)max, p≥ 1ε√1−p , p < 1 ,(3.41)from which I can calculate ‖gp‖∞ by substituting f (m)∗ into (3.22). I note that forp < 1, the maximum gradient does not depend on f (m) but only on the chosen pand ε value. Figure 3.10(c) presents the derivatives of different approximated `p-norms after applying the corresponding γ-scale. The role of γs is to reference thepartial derivatives of the approximated `p-norms to the derivatives of its `2-normmeasure. This re-scaling is done for two reasons. First, at the transition betweenStage 1 and 2, it preserves the balance between the misfit and regularization termsand thus no large adjustment in the trade-off parameter β is needed. Secondly, thescaling based on the gradients guarantees that two penalties can co-exist and impactthe solution at every step of the IRLS, regardless of the chosen {p,ε}-values or theamplitude of f (m).50I therefore define Scaled-IRLS weights such that (3.25) become:Rˆs = γs diag[((m(k−1))2+ ε2)ps/2−1]1/2(3.42)where the scaling parameter γs is:γs =[ ∥∥g2s∥∥∞∥∥gpss ∥∥∞]1/2(3.43)Two options are possible to compute the γ-scalings: (a) take the maximum absolutegradient directly from the gradient values in (3.38), or (b) calculate∥∥∥gpj∥∥∥∞ analyt-ically from (3.41). I have found that Option 2 is more stable since it is based upona theoretical maximum of the gradient and not on a particular realization of thatmaximum that arises from the distribution of values in the current model mk.The outcome of the re-scaling strategy is shown in Figure 3.11(a). The solu-tion seems to have my desired properties of being sparse in terms of the numberof non-zero model values and the model has smooth edges. The maximum partialderivatives, shown in Figure 3.11(c), confirm that the scaling strategy was success-ful in balancing the impact of the two components of the regularization. This isquantified by the calculated proportionality ratio λ∞ = 0.7. It is an improvementover the previous solution with a ratio of 150 (Figure 3.9), but it appears that thealgorithm has reached a steady state solution with slightly more influenced fromφs. In the following section I provide a strategy to better preserve proportionalitybetween each model update through a cooling strategy.Scaled model derivativesApplying the same scaling strategy to the model derivative requires additional careas the measure is also dependent on length scales. Following the strategy estab-lished in (3.35), I can factor out the length scales from the model derivative term51Figure 3.11: (a) Recovered model and (b) convergence curves using the Scaled-IRLS approach for ps = 0, px = 2 and a fixed threshold parameter ε = 1e− 3(αs = αx = 1). (c) Trade-off parameter and maximum gradients for the differentcomponents of the objective function. The scaling procedures preserves theproportionality between∥∥gpss ∥∥∞ and ∥∥gpxx ∥∥∞ throughout the iteration process.The trade-off β -parameter needed only to be adjusted slightly at the beginningof Stage 2.such that the partial derivative of φ px with respect to m can be written asgpx =∂φ px∂m= h−px gˆpx= h−pxD>x Rˆ>x V>x W>x WxVxRˆxDxm(3.44)whereDx =−hˆ−11 hˆ−11 0 . . . 00. . . . . . . . ........ . . 0 −hˆ−1M−1 hˆ−1M−1 . (3.45)measures the model derivatives over normalized length scales hˆi. The IRLS weightsin (3.29) becomeRˆx = diag [rˆx]1/2rˆxi =(m(k−1)i+1 −m(k−1)ihˆi)2+h2ε2px/2−1 . (3.46)52Since the maximum of gpx scales linearly with the base length scale h−px :‖gpx‖∞ = h−p‖gˆpx‖∞ , (3.47)I can express the γx-scaling asγx =[hp−2∥∥gˆ2x∥∥∞∥∥gˆpx∥∥∞]1/2. (3.48)By simply factoring out the base length h, I recover a multiplication factor hp−2that is closely related to the αs scaling strategy specified in (3.36). If I use theγx-scaling in (3.48) and the constant αs scaling prescribed in (3.8) I get thatφm = h−2φs+hp−2φx . (3.49)Conversely, if I set the α parameters based on the strategy put forward in (3.36) Iget thatφm = h−pφs+φx . (3.50)Expression (3.50) and (3.49) are related by a multiplication factor hp−2. Since thepartial derivatives of φm are invariant with respect to a global constant, I can expectto get the same solution with either approach, albeit a re-adjustment of the tradeoffβ parameter.Rather than choosing one approach over the other, and in order to simplify thedefinition of the hyper-parameters α , I propose to remove the constant factor hfrom the regularization function altogether. I can do this by directly evaluatingthe model derivatives with the finite difference operator described in (3.45). Thissimple change brings both φ pss and φ pxx to be dimensionally equivalent such that[φˆ px ][φ ps ]= 1 . (3.51)where φˆ px denotes the measure of model derivative using the finite difference ap-53proach. The scaled IRLS weights becomeRˆx = γx diag[((Dxm(k−1))2+ ε2)px/2−1]1/2, (3.52)whereγx =[ ∥∥gˆ2x∥∥∞∥∥gˆpxx ∥∥∞]1/2. (3.53)The Scaled-IRLS regulatization becomes:φ pm = αs‖WsVs Rˆs m‖22+αx‖WxVx Rˆx Dx m‖22 . (3.54)with default parameters αs = αx = 1.In order to demonstrate that this change in length scale does not change theoverall objective function, I proceed with two inversions. I repeat the experimentpresented in Section 3.1.3 for ps = px = 1, also shown in Figure 3.12(a) for com-parison. First, I use the regularization function in (3.54) with uniform scaling withthe same uniform discretization (αs = αx = 1). The recovered model shown inFigure 3.12(b) is almost identical to the previous solution. Small discrepancies be-tween the two models can be attributed to slight differences in the iterative process.As presented in Table 3.3, the global objective function φ(m) remains unchangedwith both approaches. Changes in scale between individual components of the ob-jective function (φd , φps , φ px ) are absorbed by their respective hyper-parameter (β ,αs, αx).The second experiment tests the case of a non-uniform discretization. Using thesame noisy data, I refine the mesh on the right-half of the domain by dividing thecell size by a factor 2 (hi = 0.01). I also re-adjust the kernel functions Fi j withinthe refined region such that each random coefficient is sampled twice but over asmaller cell length. The final inversion mesh contains 75 parameters, compared to50 parameters in the previous experiment. After convergence of the algorithm Irecover the model presented in Figure 3.12(c). Once again, the final model and thecalculated components of the objective function (Table 3.3) are almost identical tothe previous experiments. This demonstrates that the normalization of the modelderivatives does not change the global objective function, and that the solution54Figure 3.12: Recovered 1D models for ps = px = 1 using different scalingstrategies and parameterization. (a) Solution previously shown in Figure 3.7that uses the standard gradient measure (αs=50). (b) Solution obtained with thefinite difference approach (αs = αx = 1). (c) Model recovered with a differentparameterization such that the right half of the domain has cells with half thesize. Small discrepancies between the three solutions can be attributed to slightdifferences in the iterative process.φd β αs αx βαsφs βαxφxFigure 3.12(a) 4.87 65.13 50.0 1.0 56.18 8.60Figure 3.12(b) 5.13 3381.87 1.0 1.0 58.33 8.81Figure 3.12(c) 4.93 3236.60 1.0 1.0 55.81 8.59Table 3.3: Components of the objective function corresponding to the inversionresults presented in Figure 3.12 for ps = px = 1.remains independent on the choice of discretization. For clarity, I will use φ px todenote the measure of model derivatives using the finite difference operator.3.2.2 Threshold ε-parameterWhile I have improved the flexibility of the IRLS algorithm, I have yet to addressthe threshold ε-parameter which has been held fixed. The choice of thresholdparameters remains a subject of disagreement among researchers. In the earlywork of Last and Kubik (1983), it was suggested that the threshold value should besmall or near machine error (ε < 10−8) in order to best approximate the `p-norm.The same strategy was later adopted by others (Barbosa and Silva, 1994; Stoccoet al., 2009). Other researchers, such as in Ajo-Franklin et al. (2007) observed55instabilities with small values, and opted for a wider range (10−4 < ε < 10−7).More recently, Sun and Li (2014) proposed an ε-search phase to highlight re-gions reacting favourably to the sparsity constraints. A final inversion step was thencarried out with a fixed threshold value (ε  1e− 4). A similar strategy has alsobeen proposed by Zhdanov and Tolstaya (2004) after selecting an optimal point ona trade-off curve.Selecting an appropriate threshold value becomes more complicated when com-bining different penalty functions. I not only need to contend with the range ofmodel values, but also with the relative influence of the components of a non-linear regularization function. To illustrate the challenge, I invert the 1D exampleagain with the mixed norm penalty function (φm = αsφ 0s +αxφ 2x ) but this time overa range of threshold values (100 < ε < 10−5). The resulting models are shown inFigure 3.13. I identify the following trends:• For large values (ε > 10−1), no sparsity is achieved and the model resemblesthe solution previously obtained with φm = αsφ 2s +αxφ 2x .• With small values (ε < 10−4), φ pss appears to have little influence on the so-lution and the model resembles the solution obtained with smooth penaltiesαs = 0. The proportionality ratio λ∞ 1 confirms this bias towards φ 2x .• The mid-range values (ε−1 < ε < 10−3) show the most significant variabilityin the solutions with an achieved proportionality ratio λ∞ ≈ 1.From this numerical experiment, there appears to be an optimal ε-parameterin the mid-range (ε−1 < ε < 10−3) that can promote both a sparse and smoothsolution. Ideally, I want to automate the selection process.In this study, I opt for a cooling strategy. Threshold value ε is initialized at alarge value then monotonically reduced such that:ε(k) =‖m(0)‖∞ηk, (3.55)In (3.55), η is a user-defined cooling rate constant and ‖ f j(m)(0)‖∞ denotes thelargest function value obtained at the end of Stage 1 of the algorithm. At thestart of Stage 2, the Lawson approximation with large ε is effectively an `2-norm.56Figure 3.13: Recovered 1D models with variable threshold parameter on therange 10−5 < ε < 100 using a mixed-norm penalty function φm = αsφ 0s +αxφ 2x .Thus there is only a small change in regularization between Stages 1 and 2 of thealgorithm. This is desired since the `p-norm regularization is highly non-linear andI want to reduce the risk of moving away from my initial proportionality conditions.I proceed with an inversion with a cooling rate η = 1.25. Recovered modelsas a function of iterations are shown in Figure 3.14(a) to (d). As the number ofiterations increases and ε → 0, the emphasis of the sparse penalties (γ2s g0s ) sweepsthrough the range of model values, progressively focusing on smaller model pa-rameters. Figure 3.14(e) plots the scaled gradients as a function of absolute modelvalues obtained at iteration k=10, 15, 24 and 55. This plot can be compared to Fig-ure 3.10(c). The gradients associated with sparse penalties (g0s ) force small modelvalues (m ≈ ε) towards the reference (mre f = 0). Large model values (m >> ε)are free to increase unless penalized by the other competing functions (φ pxx , φd).Figure 3.15 illustrates the evolution of penalties by plotting the partial derivativesof the objective function for iteration k=15 (early stage) and iteration k=55 (latestage). As ε → 0, small model values are primarily influenced by ∂φd∂m and βαs ∂φs∂m ,57Figure 3.14: (a)-(d) Recovered 1D models at different iteration steps (k). Thevalue of ε (dash) is shown for reference, highlighting the idea of a progressivethresholding of model values. (e) Scaled partial derivatives (γ2g0s ) as a functionof the sorted model values and at various iteration stages.while large model values are influenced by ∂φd∂m and βαx∂φx∂m .From a user standpoint, the cooling strategy is attractive as it eliminates therequirement to predetermine an optimal ε threshold values and instead relies on acooling rate η . To investigate the impact of the cooling rate on the solution, I solvethe inverse problem (ps = 0, px = 2) for various cooling rates η = [6, 3, 1.5, 1.125].For each independent trial, the iteration process continues until reaching the addi-tional criteria that ε(k) = 10−6 (near machine single precision). Recovered solu-tions for the 4 cooling rates are presented in Figure 3.16(a-d). I note importantdifferences between the solutions. A summary of the inversions is provided inTable 3.4.Figure 3.16(e) displays convergence curves for each inversion trial. The finalnorm φ pm is evaluated by using the expression (3.30), that is, all γ scaling parame-58Figure 3.15: Partial derivatives of the objective function for iteration (top) k=15and (bottom) k=55. The recovered models are shown in red for reference.Figure 3.16: (a-d) Recovered models and (e) convergence curves for the mini-mization of φd +βφpm with various cooling rates η but with a final ε∗ = 1e−6.59η # Iterations φ pm ε(k) φ(k)d λ1.125 77 12.6 1e-06 9.2 0.971.5 39 17.8 1e-06 9.1 1.273 24 24.1 1e-06 8.5 0.686 23 31.1 1e-06 8.9 0.54Table 3.4: Inversion summary after convergence of the S-IRLS for various cool-ing rates η as presented in Figure 3.16. Each inversion trial was required to reachthe target misfit φ ∗d and target ε∗ = 1e−6 .ters have been removed. Thus I am comparing my ability to minimize the globalobjective function even though the algorithm has used variable scalings to reachthat result. This promotes insight into the robustness of the final model as a func-tion of the cooling rate. I note that a lower model norm can be obtained by coolingslowly. Cooling at an even slower rate (not shown here) has reaffirmed this. Arate of η = 1.025 yielded φm = 12.5 in 155 iterations. It appears that φ pm → 12.5as η → 1. I found experimentally that for η ≈ 1.25 generally yielded an optimaltrade-off between computational time (number of iterations) and convergence to asuitable solution.3.2.3 SummaryMy goal is to solve an inverse problem where the regularization function is com-posed of multiple terms, each defined as an `p-norm premultiplied by a scalingparameter. The scaling parameters, α’s, are used to control how much each com-ponent contributes to the final solution. The relative influence of these componentsis quantified by evaluating the proportionality ratio λ∞ (3.19). If two componentscontribute equally, then λ∞ should be close to unity. Unfortunately, when the com-ponents of the regularization include model and gradient terms, the scaling is af-fected by the cell size chosen for discretization. To simplify the implementation Inormalize the length scales used in the measure of model derivatives. This makesφ pss and φ pxx dimensionally equivalent. The default values for obtaining equal con-tributions are thus αs = αx = 1 for all combinations of `p-norms on the gradients.I solve the inverse problem by replacing the `p-norms with their Lawson norm60approximations. Thus I search for the model that minimizesφ(m) =‖Wd ∆d‖22+βαsM∑i=1wsivsim2i(m2i + ε2)1−ps/2+βαxM−1∑i=1wxivxi(mi+1−mihˆi)2[(mi+1−mihˆi)2+ ε2]1−px/2Vi ,(3.56)for arbitrarily small ε value. I solve the inverse problem using a two-stage ap-proach. I first find a solution for the `2-norm problem and then I change the objec-tive functions to their final desired `p-norms and solve the optimization problemusing IRLS. To keep stability in the iterative process I successively rescale theIRLS weights R. Thus at each iteration I solve a locally convex problemφ(m(k)) = ‖Wd ∆d‖22+β(αs‖Ws Vs Rˆs m‖22+αx‖WxVx Rˆx Dx m‖22),(3.57)Although the local minimization problems involve scaled gradients, the finaldesired solution is that which minimizes (3.56) such that all components contributeequally. My ability to achieve this goal depends upon the value of ε and the chosencooling rate. I find that the best (i.e. minimum norm) solution is obtained whenε is cooled slowly to a final small value. If the cooling is too fast then I obtain asubstandard solution in which λ∞ is not close to unity and my modelling objectivesare not satisfied. Slower cooling and the frequent re-scaling of the gradients keepsthe proportionality ratio near unity.3.3 Exploring the model spaceThe smooth density model presented in Figure 3.2 was a poor approximation of acompact block, but it is one of many possible solutions. Now that I have developedan algorithm that can combine multiple regularization functions with different `p-norm measure, I can explore the model space by generating a suite of solutionsthat have variable characteristics. I will demonstrate this on the synthetic gravity61example shown in Figure 3.1. The function to be minimized for the 3D gravityproblem becomesminρφ(ρ) = ‖G ρ−dobs‖22+β ∑r=s,x,y,zαr‖WrVr Rˆr Dr ρ‖22s.t. φd ≤ φ ∗d(3.58)I carry out eight additional inversions. I use a combination of norms on a rangeof ps, p[x,y,z],∈ [0,1,2] values. I set px = py = pz in all cases. The solutions,nine models in total, are presented in Figure 3.17. All models have a final misfitφ ∗d ≈ 441 and use the same `2-norm solution to initiate the IRLS steps. I make thefollowing general observations. There is a progressive transition from a smoothmodel (upper left) to a blocky solution (lower right) as ps, px, py and pz decrease.The top of the density high is most often recovered at 10 m depth. Away fromthe anomalous region the density is relatively smooth and close to the backgroundreference model of 0 g/cc. There is also a clear trend in the data misfit map suchthat the correlated residual decreases as ps, px,py, pz→ InterpretationAccessing a range of solutions is important to assess the stability of different fea-tures and to avoid over-interpreting one specific realization. The next step requiresto compare this ensemble of models and make a geological interpretation. In Chap-ter 6 I provide a more evolved methodology to extract local parameters, but fornow, I will compare the solutions visually. Figure3.18 presents an overlay for the10th and 90th percentiles anomalous densities calculated from the suite of modelsshown in Figure 3.17. I can assess the robustness of features by comparing theiso-contour lines of each model: tight clustering of the contours indicates that sev-eral models agree on the position of an edge, while a large spread indicates highvariability. At the center of the model, I note that the top of the anomaly (solid)is highly correlated among models, but less so the bottom limit. This is expectedas the resolution of the survey decreases with depth. Meanwhile, on the edgesof the domain, the shape and sign of density contrasts vary substantially. Fromthis simple analysis, I would assign high confidence on the horizontal and top of62Figure 3.17: (a-i) Vertical section through a suite of density models recoveredfor varying `p-norm penalties applied on the model and model gradients forps ∈ [0, 1, 2] and px = py = pz ∈ [0, 1, 2].the positive density anomaly, and low confidence on features on either side of theinversion domain.The normalized data residual maps for each inversion are shown in Figure 3.19.The decrease in correlated residual observed on the misfit maps (Figure 3.18) isalso an important aspect to consider. While all the inversions have achieved theglobal target misfit, only after applying the proper constraints (blocky and compact)that the inversion was able to predict the short wavelength information presentin the data. This is an important aspect of this research as it further stresses theimportance of exploring a range of solutions with broadly different characteristicssuch that more subtle features can be extracted. It is also a motivation to automatethe search for suitable inversion parameters that can better resolve the data.63Figure 3.18: Iso-contour values for the 10th and 90th percentile of anomalousdensity calculated from the suite of models shown in Figure 3.17. The outlineof the target (red) is shown for reference. Contour lines tightly clustered indicatecoherence between inversion trials. Negative anomalies (dash) appear to changesignificantly, to which I would assign lower confidence.64Figure 3.19: (a-i) Residual data map calculated from the suite of density modelsfor varying `p-norm penalties applied on the model and model gradients forps ∈ [0, 1, 2] and px = py = pz ∈ [0, 1, 2].65Chapter 4Sparse Magnetic Vector InversionI have so far showcased the benefits of employing a mixed-norm inversion to solvethe gravity problem. In this chapter, I review and improve inversion methods forthe modeling of magnetization parameters. Contrary to the scalar density, magneti-zation is a vector quantity defined by a magnitude and an orientation as previouslydefined in (1.8)~M = κ(~H0+ ~Hs)+ ~Mr . (4.1)Researchers have investigated ways to extract information about magnetizationwith inversion methods, but this process remains challenging as there are threetimes the number of unknowns over conventional problems. To simplify the in-verse problem, the assumption is often made that Earths inducing field is muchlarger than the other component such that secondary fields and the presence of re-manence are ignored (~Mr = ~Hs = 0). This assumption has dominated the inversionliterature over the past 30 years (Li and Oldenburg, 1996; Pilkington, 1997). Underthis assumption, the definition of magnetization (1.8) simplifies to~M = κ~H0 ,66which gives rise to a linear system relating N data dpre to the M discrete modelcells of magnetic susceptibility κdpre = F κdpre ∈ RN ,F ∈ RN×M,κ ∈ RM .(4.2)To illustrate potential issues with this induced assumption I revisit the syntheticblock example in terms of magnetic properties. I set the magnetic susceptibility(κ) of the block to 0.035 SI and the vertical inducing flux to ~B0 [50,000 nT, I: 90◦,D: 0◦]. I also add a remanent component equal in magnitude and pointing alongthe x-axis ~Mr [1.4 A/m, I: 0◦, D: 90◦]. This results in a total magnetization ~M[2.0 A/m, I: 45◦, D: 90◦]. Using the linear relationship presented in (2.10), I cansimulate magnetic data presented in 4.1(b) on which I add random Gaussian noiseof 1 nT standard deviation to simulate field conditions. I will attempt to recoverthe magnetized block from the noisy data shown in Figure 4.1(c).I begin with the conventional smooth assumptions (ps, px, py, pz = 2) andattempt to recover the position and shape of the magnetic block. Since (4.2) islinear with respect to κ , I can use the same inversion methodology established inChapter 3. The objective function to be minimized becomesminκφ(κ) = ‖F κ−dobs‖22+β ∑r=s,x,y,zαr‖WrVr Rr Dr κ‖22s.t. φd ≤ φ ∗d(4.3)As I am dealing with strictly positive magnetic susceptibility κ , I impose boundconstraints by the projected gradient method (Vogel, 2002). Model parameters κithat become negative are set to zero and ignored for the following Gauss-Newtonstep in (3.32). After reaching the target misfit criterion in equation (3.33), I re-cover the susceptibility model shown in Figure 4.2(a). I note that the position ofthe susceptibility anomaly is shifted to the side of the true block and appears todip at 45◦ angle. This is due to the large negative data lobe introduced by the re-manent component that I have purposefully ignored. Attempting to improve thesolution by solving for a sparse model (ps, px, py, pz = 0) yields the solution pre-sented in Figure 4.2(c). The magnetic anomaly is imaged at the right depth and67Figure 4.1: (a) Vertical section through a 25 m cube with uniform magnetization~M [2.0 A/m, I: 45◦, D: 90◦]. (b) Simulated TMA data response on a 21× 21survey grid placed 15 m above the anomaly. (c) Magnetic data with randomGaussian noise added, 1 nT standard deviation.the vertical extent is better recovered, but position and shape of the anomaly havenot improved. It is also important to note the correlated negative data residuals inFigure 4.2(b) and (d). The inversion struggled to reproduce the negative portion of68the magnetic data using strictly positive susceptibility values.Figure 4.2: Vertical section through the recovered susceptibility model using(a) smooth assumption (ps, px, py, pz = 2) and (b) sparse `p-norms to recovera compact model (ps, px, py, pz = 0). Both solutions show an anomaly with afalse dip due to the wrong assumption of a vertical magnetization.The presence of remanence has long been recognized as an obstacle for thegeological interpretation of magnetic data. Several types of mineral deposits areassociated with remanent magnetization such as diamondiferous kimberlites, vol-canic massive sulphides and porphyries (Enkin, 2014; Henkel, 1991). In a miningexploration context, having the wrong image could result in false drilling targets:costly both in time, resources and confidence in geophysical methods. Hence theneed for a more robust algorithm that does not require knowledge about the ori-entation of magnetization. In the following sections, I revisit the Magnetic VectorInversion introduced in the Ph.D. thesis of Lelie`vre (2009). I address numerical69limitations encountered with the spherical formulation. I also leverage advance-ments made in Chapter 3 and impose sparsity penalties on the magnetization di-rection and amplitude independently. This allows to recover compact bodies withcoherent magnetization direction.4.1 Magnetic Vector Inversion - Cartesian parametersRather than assuming a purely susceptible response, Lelievre and Oldenburg (2009)proposed a strategy to directly invert for the magnetization vector ~M. Re-writingthe discrete system (4.2) in terms of vector magnetization:dpre = FCmC= [Fu Fv Fw][ κuκvκw]Fu, Fv, Fw ∈ RN×M(4.4)where the model mC = [κu,κv,κw]> describes the strength of magnetization alongthe Cartesian directions in terms of an effective susceptibility parameter,κe =~M‖~H0‖, (4.5)The effective susceptibility scales the amplitude of magnetization with respect tothe inducing field strength ‖~H0‖. The forward relationship (4.4) is still a linearsystem of equations but it has three times the number of unknown parameters com-pared to the susceptibility assumption (mC ∈ R3M). The objective function (3.58)becomesminmCφ(mC) = ‖FCmC−dobs‖22+β ∑c=u,v,w∑r=s,x,y,zαcr‖Wcr Vcr Rcr Dcr Pc mC‖22 ,(4.6)where the projection matrices Pc select individual Cartesian component of the vec-tor model mC. The regularization function is made up of twelve terms. Differentnorm measures can be applied to each Cartesian component independently.Keeping the same inversion methodology and smooth assumptions (pcs , pcx ,70pcy , pcz = 2), I recover the magnetization model presented in Figure 4.3(a). Thissolution is an improvement over the susceptibility model; the bulk magnetization isrecovered at the right position. I note however that the solution is distributed overa large volume and there is a broad distribution in magnetization direction insideand around the block.In order to reduce the complexity of the solution, I once again resort to `p-normmeasures (pcs , pcx , pcy , pcz = 0). As shown in Figure 4.3(b), the recovery of theblock has clearly improved. It is important to point out however that the magneti-zation vectors are pointing along the Cartesian directions and the final anomaly isslightly wider than the true model. In the Cartesian formulation, both the directionand strength of magnetizations are coupled in the vector components and thereforethe method lacks flexibility in recovering sparse vector along arbitrary orientations.This was the main motivation behind recent research investigating advanced reg-ularization methodologies and cooperative approaches (Fournier, 2015; Liu et al.,2015; Zhu et al., 2015). I will attempt to improve on this solution by decoupling thestrength and direction of the magnetization vector with the spherical formulation.4.2 Magnetic Vector Inversion - Spherical parametersAs an alternative to the Cartesian formulation, Lelievre and Oldenburg (2009) alsoproposed the vector inversion in a spherical coordinate system. The conversionbetween Cartesian to spherical system follows the relation:κu =ρ cos(θ) cos(φ)κv =ρ cos(θ) sin(φ)κw =ρ sin(θ)(4.7)where the magnetization vector is defined by parameters of amplitude (ρ) and twoangles (θ , φ ). Performing the forward calculations using spherical parameters mScan expressed asF[mS] = FCρ cos(θ) cos(φ)ρ cos(θ) sin(φ)ρ sin(θ) (4.8)71Figure 4.3: Vertical section through the recovered magnetization vector modelusing the Cartesian formulation with (a) smooth l2-norm assumption and (b) spar-sity constraints applied on all three Cartesian components(pis , pix , piy , piz = 0).Corresponding normalized data residuals are shown in (b) and (d) . The true po-sition and magnetization orientation of the block are shown in red for reference.The spherical formulation separates the magnitude and orientation of magnetiza-tion vector; this has two advantages. First, physical property constraints can easilybe incorporated in the inversion. Magnetization measurements performed on rocksamples are often provided in term of susceptibility κ and Koenigsberger ratio QQ =~Mr‖κ~H‖ (4.9)measuring the ratio between the induced and remanent component of magnetiza-tion. This quantity does not contain directionality but it provides information about72the amplitude of magnetization. Magnetization direction may also be used in theinversion when laboratory measurements are performed on oriented cores. Ref-erence angles can than be introduced as constraints. The second advantage of thespherical coordinate system is that it allows for sparsity assumptions to be imposedon the magnitude and orientation independently. This has the potential of resolv-ing compact bodies with uniform magnetization direction in any orientation notrestricted to the Cartesian axes.Despite its obvious advantages, the MVI-S method has received little atten-tion in the literature because the non-linear transformation between Cartesian andspherical coordinates greatly complicates the inverse problem. I demonstrate chal-lenges encountered with the spherical coordinate system on my synthetic problem.Taking the partial derivatives of (4.4) as a function of mS(ρ,θ ,φ) yields:J =∂F[mS]∂mS=∂F[mS]∂mC∂mC∂mS(4.10)where ∂mC∂mS involves partial derivatives of trigonometric functions prescribed in(4.7). The sensitivity matrix can be linearized before each Gauss-Newton stepas:J = FC S (4.11)where the matrix S holds the partial derivativesS =cosθ cosφ −ρ sinθ cosφ −ρ cosθ sinφcosθ sinφ −ρ sinθ sinφ ρ cosθ cosφsinθ ρ cosθ 0 (4.12)The regularization function becomesφm = ∑c=ρ,θ ,φ∑r=s,x,y,zαcr‖Wcr Vcr Rcr Dcr Pc mS‖22 , (4.13)I also take three additional precautions to deal with the spherical parameterization.First, a zero reference value θ re f and φ re f would imply a magnetization directionpointing along the x-axis. In the absence of physical constraints, and since I do notwant to assume a specific orientation, I set αθs = αφs = 0 in all my experiments.73Thus the regularization only penalizes the change in angle between neighboringcells. Secondly, I perform an adjustment to the measure of angle differences to dealwith the discontinuity at θ =−pi, pi in order to prevent over penalizing angles thatdescribe a similar orientation along the x-axis. I convert the difference in anglesbetween adjacent cells to the coterminal angles as shown in Figure 4.4. Thirdly,rather than resorting to a projected gradient approach, I proceed with a doubletransformation (mS → mC → mS) such that spherical parameter remain boundedas ρ ∈ [0, ∞], θ ∈ [−pi, pi] and φ ∈ [−pi/2, pi/2].Figure 4.4: Measure of angle differences ∆θ converted to coterminal angle ∆θˆ .To demonstrate the difficulties encountered with the MVI-S formulation, I in-vert the synthetic TMI data with a starting magnetization model pointing upwardm(0)S (ρ = 10−2,θ = −45◦,φ = 0◦) as shown in Figure 4.5(a). This represents aworst-case scenario such that the assumed magnetization orientation is at 90◦ fromthe true orientation. After convergence of the algorithm, I recover the model shownin Figure 4.5(b). The solution is a poor representation of the true magnetization.Model updates were forced to stop before reaching the target data misfit as theinversion was likely trapped in a local minimum. I note that most of the modelupdates were performed on the amplitude ρ , with only marginal changes on theangle of magnetization. Similar behaviors have been documented by Lelievre andOldenburg (2009) and later by Liu et al. (2017). Poor convergence was attributedto an imbalance between the model parameters. Before attempting to implementmore advanced constraints, I make inroads in improving the convergence of the74non-linear MVI-S formulation.Figure 4.5: Vertical section through the (a) starting model and (b) recoveredmagnetization vector model in Spherical coordinates. The true position and mag-netization orientation of the block are shown in red for reference. (c) Normalizeddata residuals show correlated signal. The inversion stopped after 15 iterationsand was enabled to further reduce the objective function.To gain some insight about the issues encountered with the MVI-S formulation,I consider a simpler two-parameter linear problem of the formmx+ 2∗my = 1 , (4.14)which I can express in matrix form asFC mC = dobs ,FC = [1 2], mC =[mxmy], dobs = 1(4.15)75This defines an under-determined linear system of equations. Just as I did forthe magnetic inverse problem, I can isolate a solution by minimizing an objectivefunction of the formφ(m) = ‖FC mC−dobs‖22+βC‖mC‖22 , (4.16)Figure 4.6(a) displays a contour map of the objective function along with its gra-dients. Following the same methodology as in (3.9), I find a solution such that thegradient of the objective function φ(m) vanish∂φ∂m= g =(F>C FC +βCI)mC−F>C dobs = 0 (4.17)where I is the identity matrix. The factor 2 from the derivative of the `2-normis absorbed by the zero right-hand side. After determining a trade-off parameterβC such that φd ≤ 1e− 3, I recover the Cartesian model mC = [0.2,0.4]. It isthe solution that minimizes the distance (evaluated with the `2-norm) between theorigin and the solution space of F. I note that the relative magnitudes of modelparameters in mC reflect the size of the forward coefficients in F.As previously discussed for the gravity and magnetic problems, the smallestsolution is often not satisfactory as it is strongly influenced by the physics of theexperiment. From equation (3.11) in Chapter 3, I can introduce sensitivity basedweights to counteract this bias towards a large my valueWC = diag[[wCmax(wC)]1/2]wC j =[N∑i=1Fi j2]1/2,(4.18)where WC holds sensitivity weights added to the regularization (wC = [1, 2]>). Thenew weighted objective function becomesφ(m) = ‖FC mC−dobs‖22+βC‖WCmC‖22 , (4.19)76Figure 4.6: Contour map for two objective functions and their gradients (ar-rows) for a two-parameter inverse problem solved in Cartesian and polar coor-dinate systems. (a) The non-weighted (gray) problem yields a solution mC =[0.2, 0.4] that reflects the size of the forward coefficients. The sensitivityweighted (black) solution is more desirable as it predicts the data with equalmodel parameters m∗C = [0.33,0.33]. (b) Inversion steps performed in the non-linear polar coordinate system (solid red) are highly oscillatory and do not reachthe optimal solution m∗P = [0.47,0.79]. The same model updates are shown inCartesian coordinates (red dash) for reference. (green) Inversion steps performedin polar coordinates with iterative sensitivity re-weighting and (blue) with theadded Ω-scaling.and a weighted gradientgC = F>C FCmC +βCW>C WCmC−F>C dobs (4.20)After determining the appropriate β ∗C , I get the solution m∗C = [0.33,0.33] markedwith a black circle in Figure 4.6(a). I have reached the only solution with equalcontribution from both model parameters that also predicts the data within the tol-erance.Alternatively, I can attempt to solve the same problem in a polar coordinate77system under the transformationmP = [ρ, θ ]>mx = ρ cosθmy = ρ sinθ ,(4.21)where the polar model mP is defined by a radius ρ and an angle θ . This is analo-gous to the spherical transformation performed for the MVI-S formulation in (4.7).The objective function to be minimized becomesφ(mP) = ‖F[mP]−dobs‖22+βP‖WCmP‖22 (4.22)The inverse problem is now non-linear with respect to the polar model so I solve ititeratively with the standard Gauss-Newton procedure described in equation (3.31).The partial derivatives of the forward mapping with respect to the polar coordinatesare calculated byJ =∂F[mP]∂mP=∂F[mP]∂mC∂mC∂mP= FCS (4.23)where the matrix S holds the partial derivatives of the model with respect to thepolar parametersS =[cosθ −ρ sinθsinθ ρ cosθ]. (4.24)The gradient of the objective function becomesg = S>F>C F[mP]+βPW>C WCmP−S>F>C dobs (4.25)A trade-off parameter βP is determined through the cooling schedule establishedin Chapter 3. The inversion is terminated once the data misfit and change in modelnorm fall below the tolerances ηφd and ηφm defined in equation (3.33) and (3.34)respectively.Since m∗C is a desirable model, I would like to be able to recover a similar78solution in polar parameters (m∗P = [0.47,0.76]). Unfortunately, as shown in Fig-ure 4.6(b), the minimization process performed in polar coordinates converges to adifferent solution (mP [ρ = 0.67, θ = 0.26]) and the iterations steps are oscillatory.I display the equivalent iterations (red dash) in the Cartesian space for comparison(mCP [mx = 0.65, my = 0.17]).My main goal is to recover the solution m∗P, and I want to reach this solutionwith only a few model updates. To understand the discrepancy between the twoformulations, I compare their respective gradients for a given starting model m(0)Cand its equivalent polar model m(0)P . In Cartesian coordinates, the gradient directionisg(0)C = F>C FCm(0)C +β∗CW>C WCm(0)C −F>C dobs (4.26)assuming that I know the optimal trade-off parameter β ∗C . I can convert these gradi-ents to polar coordinate by multiplying (4.26) with the matrix of partial derivativesS such thatgPC = S>[F>C FCm(0)C +β∗CW>C WCm(0)C −F>C dobs](4.27)I want to compare this gradient to the gradient calculated in polar coordinates tofind some equivalence between the two systemsgPC ' S>F>C F[m(0)P ]+βPW>C WCm(0)P −S>F>C dobs . (4.28)I can simplify both sides of equation (4.28) by noting that FCm(0)C = F[m(0)P ] .Equation 4.28 and 4.29 are therefore the same ifβ ∗CS>W>C WCm(0)C ' βPW>C WCm(0)P . (4.29)I would like both sides to be roughly equal such that the gradient direction in polarspace resemble the gradient direction calculated in the Cartesian space. First, I notefrom equation (4.29) that the transformation matrix S is missing from the right-hand side. The current sensitivity weights were designed to compensate for fixed79experimental bias but did not account for a constantly varying sensitivity matrix J.I address this shortcoming by resorting to an iterative update of sensitivity weightsWP = diag[[wPmax(wP)]1/2]wPj =[N∑i=1Ji j2]1/2,(4.30)where the weights wP are updated between each Gauss-Newton step. Invertingonce again the non-linear problem with the iterative scaling strategies (green) Irecover the model mP[ρ = 0.53, θ = 0.53] (Fig. 4.6(b)). The solution has equalparameters of ρ and θ , and I reached this solution in a few iterations. In mostapplications, however, obtaining proportionality between the magnitude and angleof the vector is not meaningful. Converted to Cartesian space mCP[mx = 0.46, my =0.27], I note that the solution is still different from m∗C.To understand this result, I now examine equation (4.29) in terms of the sizeof the model parameters in mP. I have used a regularization function to penalizetwo parameters with different units: the radius ρ ∈ [0,∞] in lengths and angleθ ∈ [−pi, pi] in radians. The range of values spanned by these parameters differ inscale as depicted by the aspect ratio of Figure 4.6(b). Comparing the largest changein model values (gradient) in relation to the Cartesian space mC, it is easy to showthat‖gρ‖∞ ∝ ‖gx‖∞+‖gy‖∞ (4.31)such that a change in the radius ρ is proportional in magnitude to a change incomponents in Cartesian space. The same relation does not hold for the angleparameter as an equivalent change in Cartesian parameters can be achieved with arotation ∆θ = pi/2 independent of length. In order to scale the gradient steps takenalong different dimensions, I define a proportionality relationω =‖gρ‖∞‖gθ‖∞ = ‖ρ‖∞2pi, (4.32)80and scale the regularization asWˆP = diag([1 ω](1/2))WP (4.33)My new scaled objective function becomesφ(mP) = ‖F[mP]−dobs‖22+βP‖WˆPmP‖22 (4.34)Minimizing this function I get the model (blue) presented in Figure 4.6 (mCP[m1 =0.47, m2 = 0.74]). Converted to Cartesian space, this solution mCP[m1 = 0.35, m2 =0.32] closely matches m∗C.4.2.1 Scaled MVI-S algorithmNow that I have demonstrated the benefit of an iterative sensitivity re-weightingof the regularization, I re-visit my synthetic magnetic problem. I invert the TMIdata once more with the starting model oriented at 90◦ from the true magnetizationdirection (m(0)P [ρ = 10−2, θ = −45◦, φ = 0◦]) with smooth assumptions (ps, px,py, pz=2). The recovered magnetization model obtained after convergence of thescaled MVIS-S algorithm is presented in Figure 4.7(a). The inversion took 15iterations to converge. I note close similarities with the MVI-C solution presentedin 4.3(a), with the bulk magnetization centered around the position of the block.This is a clear improvement over the solution previously shown in Figure 4.5(c).From a practical standpoint, I have found that it is more efficient to initializethe MVI-S algorithm with the Cartesian solution. The linear MVI-C approachallows me to rapidly find a model that fits the observed data and it provides a goodapproximation to the MVI-S solution. I invert the data once more using the smoothCartesian solution in 4.3(a) as a starting model. Figure 4.8(a) shows the recoveredsolution obtained after a single iteration of MVI-S.Having achieved a stable and reasonable solution with the `2-norm, I can nowconsider applying sparsity constraints to recover a block with a coherent magneti-zation direction. I vary the regularization measure on the amplitude, derivative ofamplitude and derivatives of angles uniformly such that (pcs , pcx , pcy pcz = 0) wherec ∈ [ρ,θ ,φ ]. Figure 4.8(c) presents a section through the magnetic vector model.81Figure 4.7: (a) Vertical section through the recovered magnetization vectormodel and (b) normalized data residual using the spherical formulation with sen-sitivity based weighting (ps = px = py = pz = 2). The true position and magne-tization orientation of the block are shown in red for reference.The shape of the anomaly matches that of the magnetic block, with a magnetizationdirection uniformly orientated at 45◦ inclination. The uniform effective suscepti-bility recovered inside the block matches exactly the true model with κe = 0.05SI. The normalized data residual shows no apparent correlated signal (Fig. 4.8 d).This simple example increases my confidence in my ability to accurately recoverthe magnetization of geological bodies in 3D.4.3 SynthesisIn this section, I introduced an iterative sensitivity re-weighting strategy to improvethe convergence of non-linear inverse problems. The iterative re-scaling of theregularization function associated with the amplitude and angles of magnetizationwas crucial in order to get a stable convergence of the MVI-S algorithm. Smootherand more robust solutions allowed me to apply compact norms on the three modelparameters independently. This can greatly simplify the solution compared to thatfrom the conventional MVI-C formulation.Despite this improvement, the MVI problem remains fundamentally under-determined. Incorporating a priori information, either through model constraints82Figure 4.8: (a) Vertical section through the recovered magnetization vectormodel in Spherical coordinates using sensitivity based weighting with smoothregularization (ps = px = py = pz = 2). MVI-C solution is used as a startingmodel. The true position and orientation of magnetization of the block are shownin red for reference. (c) Recovered vector model in Spherical coordinates usingsparse and blocky assumptions on the amplitude and angles of magnetization(pcs = pcx = pcy = pcz = 0). (b) and (d) Corresponding normalized data residu-als.or joint physical properties, remains important to accurately represent the geology.In the following section, I introduce structural information to constrain the shapeof the recovered magnetic bodies.83Chapter 5Rotated `p-norm RegularizationIn Chapter 3, I introduced a mixed-norm regularization function that could be usedto recover a suite of models that had a broad range of characteristics. This wasaccomplished by altering the norms on the model and model derivatives indepen-dently. I have shown the value in interpreting an ensemble of models to assess therobustness of certain features. When inverting gravity data arising from a buriedprism, I managed to recover an anomaly at the right depth and with sharp edges.In this case, the prism was aligned with the Cartesian frame. In most cases, how-ever, the shape of geophysical anomalies can take many forms. Dip and strikesof geological contacts rarely occur at right angles and aligned with the geographicgrid. I, therefore, need a more general strategy that is adaptable to any geologicalscenarios.Some research has addressed this issue by imposing directionality in the inver-sion. Barbosa and Silva (2006) introduce a regularization function to force anoma-lies to be closer to geometric elements in 2D. Their approach can recover compactand elongated bodies, but it requires direct input from the user to define the loca-tion of different anomalies. In Li and Oldenburg (2000), directionality is enforcedby rotating the objective function along preferential axes. The method was laterrevised by Lelie`vre (2009) and Davis et al. (2012) who advocated the use of for-ward and backward difference operators to improve symmetry. While successfulin recovering elongated bodies, their approach still requires manual adjustment ofinversion parameters to control stretching factors along different orientations. As84a natural extension to their work, I am proposing in this chapter to combine themixed norm regularization function with the rotation strategy for the recovery ofsharp and oriented edges.5.1 Rotated objective functionI begin by reviewing the strategy presented in Li and Oldenburg (2000) for therotation of the objective function using the conventional `2-norm regularizationφm = αs∫Vm2 dV + ∑r=x,y,zαr∫V(dmdr)2dV . (5.1)Model gradients are traditionally calculated along the Cartesian axes x, y and z.Evaluating the integral over the discrete domain gives rise to linear penalty func-tions introduced in equation (3.5)φs = αs‖WsVs (m−mre f )‖22+ ∑r=x,y,zαr‖WrVrDr m‖22 , (5.2)where W matrices contain sensitivity or volumetric based weighting and weight-ing parameters set by the user. As explained in Section 3.1.3, I use finite differenceoperators D rather than the conventional gradient operators G to simplify the imple-mentation. The scaling parameters αr control the stretching along the three axes.This formulation is restrictive however as the dip and strike of geological bound-aries can occur at any orientation. To address this shortcoming, Li and Oldenburg(2000) introduced rotated gradient measures of the form∂m∂x′= ωxx∂m∂x+ωxy∂m∂y+ωxz∂m∂ z∂m∂y′= ωyx∂m∂x+ωyy∂m∂y+ωyz∂m∂ z∂m∂ z′= ωzx∂m∂x+ωzy∂m∂y+ωzz∂m∂ z(5.3)85Rotation coefficients ω are defined by the triple rotation around the Cartesian axessuch thatΩ=ΩzΩxΩy =ωxx ωxy ωxzωyx ωyy ωyzωzx ωzy ωzz= cosφ cosψ− sinφ cosθ sinψ sinφ cosψ+ cosφ cosθ sinψ sinθ sinψ−sinφ sinθ cosφ sinθ −cosθ−cosφ sinψ− sinφ cosθ cosψ −sinφ sinψ+ cosφ cosθ cosψ sinθ cosψ(5.4)where the matrix Ωy defines the rotation angle ψ counterclockwise around the y-axis (dip)Ωx =1 0 00 cos(ψ) sin(ψ)0 −sin(ψ) cos(ψ) (5.5)followed by Ωx defining the rotation angle θ around the x-axis (plunge)Ωy = cos(θ) 0 sin(θ)0 1 0−sin(θ) 0 cos(θ) (5.6)and a final Ωz defining the rotation φ around the z-axis (strike)Ωz = cos(φ) sin(φ) 0−sin(φ) cos(φ) 00 0 1 . (5.7)Figure 5.1 depicts the rotation angles with respect to the Cartesian system. Theregularization functions along the rotated axis x′ is defined by substituting (5.3)into (5.1)φx′ = αx∫V(ωxxdmdx+ωxydmdy+ωxzdmdz)2dV . (5.8)86Figure 5.1: Rotation parameters with respect to the Cartesian axes.which after expanding can be written in discrete form as:φx′ = αxm>(ω2xxD>x W>x WxDx+ω2xyD>y W>y WyDy+ω2xzD>z W>z WzDz+2∗ωxxωxyD>x W>x WyDy+2∗ωxxωxzD>x W>x WzDz+2∗ωxyωxzD>y W>y WzDz)m(5.9)where I absorbed the volumes of integration Vr in the weighting terms Wr forclarity. Derivation from the integral to discrete form of (5.9) can be found in Liand Oldenburg (2000). It is important to point out that the difference operators in(5.9) require square matrices such that Dx, Dy, Dz ∈RM×M. This is done by addinga backward difference to the edge cells such that the divergence operator in (3.28)becomesDx =−1 1 0 . . . 00. . . . . . . . ........ . . 0 −1 1... . . . 0 1 −1 . (5.10)87Carrying out the same procedure for φy′ and φz′ , I obtain the rotated regularizationfunction φm′ made up of seven terms:φm′ = m>W>s Wsm+m>(D>x W>x BxxWxDx+D>y W>y ByyWyDy+D>z W>z BzzWzDz+2∗D>x W>x BxyWyDy+2∗D>x W>x BxzWzDz+2∗D>y W>y ByzWzDz)m(5.11)Rotation coefficients ω and α parameters corresponding to each term are col-lected and added to form the different diagonal matrices B. Rotation parameterscan be defined on a cell by cell basis.Lelie`vre (2009) found experimentally that this formulation resulted in asym-metric solutions due to the cancellation of rotation parameters. Checkerboard pat-terns were observed for rotation angles at 45◦. This was addressed with an aver-aged combination of forward and backward difference operators. Equation (5.11)becomesφm′ = m>W>s Wsm+m>188∑i=1(D>x,iW>x BxxWxDx,i+D>y,iW>y ByyWyDy,i+D>z,iW>z BzzWzDz,i+2∗D>x,iW>x BxyWyDy,i+2∗D>x,iW>x BxzWzDz,i+2∗D>y,iW>y ByzWzDz,i)m(5.12)such that the ith components correspond to all the combination of forward andbackward difference operators. Figure 5.2(a) and (b) compares the different gra-dient scheme along the xy-plane. For 3D problems, this formulation results in aregularization function containing 49 terms in total.I showcase the effect of rotation with three inversions applied to the density88ID Dx,i Dy,i Dz,i1 Forward Forward Forward2 Backward Forward Forward3 Forward Backward Forward4 Forward Forward Backward5 Forward Backward Backward6 Backward Forward Backward7 Backward Backward Forward8 Backward Backward Backwardblock model introduced in Chapter 2. I will take advantage of the non-uniquenessof potential fields to accentuate trends in the model, even though the true blockanomaly is a single prism orientated parallel to the Cartesian axes. In the firstinversion, I use the conventional `2-norm regularization for a rotation angle ofφ = 45◦. Since I am dealing with smooth assumption, I must impose a stretchingfactor along one of the rotated axes (αx = 100, αs = αy = αz = 1). As expectedthe recovered density anomaly is stretched and symmetric about the rotated axis x′(Figure5.3(a)). I have managed to replicate the implementation of Lelie`vre (2009).For the second inversion, I attempt to impose sparse assumptions (ps = px =py = pz = 0) to recover a compact block anomaly oriented at 45◦. As a first pass,I attempt to directly incorporate the IRLS weights into the rotated regularizationfunction such thatφm = m>W>s R>s RsWsρ+m>188∑i=1(D>x,iR>xxW>x BxxWxRxxDx,i+D>y,iR>yyW>y ByyWyRyyDy,i+D>z,iR>zzW>z BzzWzRzzDz,i+2∗D>x,iR>xyW>x BxyWyRxyDy,i+2∗D>x,iR>xzW>x BxzWzRxzDz,i+2∗D>y,iR>yzW>y ByzWzRyzDz,i)m(5.13)89Figure 5.2: Two-dimensional representation of finite difference operators usinga combination of (a) 3-cell (Li and Oldenburg, 2000), (b) 5-cell (Lelie`vre, 2009)and (c) the 7-cell scheme used for the rotation of the objective function. Only the7-cell strategy allows for the interaction between diagonal neighbors.The IRLS weights for the model derivatives are defined asRi j = diag[(|(Gi m(k−1)) (G j m(k−1))|+ ε2i j)pi j/2−1]1/2. (5.14)where denotes the Hadamard (element-wise) product between the model deriva-tives. Figure 5.3(b) shows the recovered density models after convergence of the90IRLS. While the general orientation of the anomaly appears to be aligned at 45◦,edges are poorly defined and the solution is not exactly symmetric about the x′ axis.Similarly for the third inversion, I rotate the objective function by a smaller angleof 30◦. The solution in Figure 5.3(c) shows little difference from the 45◦ rotation.This is primarily due to the averaging over the 8 combination of backward and for-ward measures. The regularization function in (5.13) also comes with a significantincrease in computational cost as it involves 16 times more operations comparedto the conventional Cartesian approach.Figure 5.3: Horizontal cross section through the recovered model with rotationof the objective function using (top) 5-point and (bottom) 7-point gradient oper-ators. Blue and red arrows indicate the rotation frame. (a, d) Recovered modelswith smooth assumptions (ps = px = py = pz = 2) stretched along the rotate x-axis at 45◦ (αx′ = 100). (b, e) Sparse solutions (ps = px = py = pz = 0) withrotation at 45◦. The target orientation of the rotated block is marked by the blackdashed line. (c, f) Sparse solutions for a 30◦ rotation of the objective function.915.1.1 Diagonal derivativesBuilding upon the previous results, I will attempt to improve the recovery of ori-ented bodies by employing a new derivative scheme. Rather than adding and aver-aging the contribution of gradient operators along the Cartesian axes, I propose toalso measure the interaction with diagonal neighbours. For clarity, I first examinethe rotated finite difference operator Dx′ along the xy-plane (rotation about zˆ). Fora given cell mi, there are 8 cells that either share a node or a face. Given a rotationangle θ , I need to determine which of the 8 neighbours participate in the finite dif-ference. I determine the contribution of each cell based on the intersecting volumewith a test cell mx′ as shown in Figure 5.2(c). To determine the position of mx′ , Icalculate the displacement of its nodes about the center of miMx′ = R(θ)1 10 00 0+xU xLyU yLzU zL (5.15)where U and L define the lower southwest and upper northeast node locations ofthe target cell mi. The rotated finite difference is calculated by a weighted averagedx′ = mi− ∑8i=1 vimi∑8i=1 vi(5.16)where the weights vi are partial volume intersected with mx′ as presented in theAppendix. The same calculation can be carried out for rotated test cells my′ andmz′ and rotation angle ψ and φ . To improve symmetry, this process is repeatedfor both a forward and backward difference. These calculations are performed foreach cell and their neighbors in our 3D domain, giving rise to the regularizationfunction:φm =‖αsWs Vs Rs (ρ−ρre f )‖22+∑r=x′,y′,z′‖αrWr Vr Rr DFr m‖22+‖αrWr Vr Rr DBr m‖22 ,(5.17)where DFr and DBr denotes the rotated forward and backward difference operators.The upfront cost to build the finite difference operators is higher than for simple92Cartesian operators but the regularization function has only seven terms in total.This translates to about 6 times fewer operations needed during the inversion pro-cess compared to the methodology proposed by Lelie`vre (2009).For comparison, I re-invert the synthetic data using this new regularizationfunction. Models are presented in Figure 5.3 for (d) 45◦ rotation with smooth `2-norm (αx′ = 100), (e) 45◦ rotation with `p-norm (ps = px = py = pz = 0) and (f)30◦ rotation. In the case of the smooth inversion, both formulations give equivalentresults with anomalous density stretched along the rotated x′-axis. I note a signifi-cant improvement however in the geometry of the rotated sparse models, with goodsymmetry and better edge definition along the rotation angle.I have achieved my initial goal of recovering rotated anomalies, even thoughin this specific case the true model is a simple prism oriented along the Cartesianaxes. Figure 5.4 presents the normalized data residual calculated from the last threeexperiments. It is important to note that all solutions fit the observed data withinthe tolerance φd−φ ∗d ≤ δd, yet large correlated residuals are clearly visible. I haveforced the solution to exhibit characteristics that are not consistent with the truesolution. This reinforces the importance of designing an objective function thatcan satisfy subtle information that may be present in the data in order to betterrepresent the local geology.Figure 5.4: Normalized data residuals calculated from the recovered modelspresented in Figure 5.3(d), (e) and (f) respectively.935.2 Synthetic fold modelI have so far demonstrated the use of rotated sparse norms for the recovery oforiented blocks by taking advantage of the non-uniqueness of the inverse problem.The true value of this strategy is to promote the recovery of oriented edges to betterrepresent geological units of arbitrary shapes. I demonstrate this by simulating amore realistic scenario. I generate a synthetic fold model made up of a layer withdensity contrast of 0.1 g/cc placed in uniform background as shown in Figure 5.5.The axial plane of the syncline strikes NS and dips 75◦E. The fold axis plunges 15◦S. The intersection of the axial plane with the various cross-sections are shown asred dashed lines in Figure 5.5. I divide the ground into an OcTree mesh with corecells that are 20 m meters in width.I simulate gravity data at 560 stations placed one meter above the topography.The locations are randomly sampled from a grid to emulate a field survey withunevenly spaced data. From the linear relation introduced in (2.3), the densityanomaly of the syncline gives rise to the gravity response presented in Figure 5.6.I will attempt to use the rotated mixed-norm regularization to recover the foldedlayer.5.2.1 `p-norm inversionI begin exploring model space by first generating a density model that is smoothand has minimum structure. Using the regularization function in (5.17) for ps =px = py = pz = 2, I recover the density model presented in Figure 5.7. As expectedfrom the smooth regularization function, the thickness of the layer near the surfaceis not clearly defined; there is a wide transition from a positive to a negative densitycontrast. More importantly, the dip and extent of the fold limbs at depth are poorlyresolved. Without prior knowledge, geological interpretation of this result wouldbe difficult with low confidence about the geometry of the fold.I follow up with a mixed norm inversion and attempt to reduce the complex-ity of the solution, with sparsity assumptions applied on the model (ps = 0) andsmooth model gradients (px,y,z = 2). Figure 5.8 presents a vertical section at thecenter of the recovered density model. I note that this solution is an improvementover the smooth model for at least three reasons.94Figure 5.5: Synthetic density model made up of a folded layer with a N-S ori-ented axial plane dipping 75◦E. The intersection between the fold axial plane andthe cross-sections are marked by a red-dashed line: (left) topography draped sec-tion and (right) vertical cross-sections along (A-A’) and (B-B’). The hinge axisof the fold is plunging 30◦S. Two strike and dip measurements are provided atthe surface and used as structural constraints.1. The density contrasts values are mostly positive with the maximum densityapproaching the true density of 0.1 g/cc. The near-surface trace of the foldappears thinner with fewer negative density contrasts seen along the hori-zontally draped section.2. The fold is imaged as a continuous body connected at depth along the sectionA-A’3. The dip angle along the fold axis (B-B’) is improved.While I have made some progress in imaging a continuous layer connected atdepth, there is still space for improvement about the shape and thickness of thefold. This is a motivation for incorporating geological information in the inver-sion. Rather than building a 3D model to constrain the solution (the usual path),I will attempt to drive the solution with so f t directional constraints in the form ofdip and strike measurements and let the data determine the position and shape ofthe fold.95Figure 5.6: Synthetic gravity data generated 1 m above topography over thesynthetic density fold model. Gravity station locations (dot) and topographiccontours (lines) are shown for reference. Structural data (dip, strike) are pro-vided at two locations on opposite of the syncline. Dip direction of the fold axis(15◦ S) and limbs are shown with arrows .5.2.2 Directional `p-normsMy goal is to introduce stratigraphic constraints in order to recover a continuouslayer connected at depth. Apart from the gravity data, structural data are oftenavailable in the form of dip and strike information measured on outcrops. I willinclude three points of structural data: one on each limb and one on the fold axis(Fig. 5.6). I will use this information to build a rotated objective function as pre-scribed in (5.17).First I need to extrapolate the structural data to the full 3D domain. I use aminimum curvature approach (Briggs, 1974). The interpolated normal vectors are96Figure 5.7: Recovered smooth solution with `2-norm regularization of the syn-thetic fold anomaly. The dip of the fold limbs are poorly recovered and there isno indication of the limbs being connected at depth.Figure 5.8: Recovered density model from the mixed norm regularization (ps =0, px,y,z = 2). Sparse assumptions helped in simplifying the density model overthe smooth assumption.shown in Figure 5.9. The orientation is then used to rotate the derivative terms ona cell-by-cell basis such that Dx′ points along the normal of the fold, while Dy′ andDz′ are parallel to the stratigraphy.97Figure 5.9: Interpolated dip and strike information used for the rotation of theobjective function. The direction of the arrows define the estimated normal orthe folded layer.I then invert the gravity data with sparse rotated norms. To re-enforce lateralcontinuity, I lower the norm applied to model derivatives parallel to the stratigraphy(px′ = 2, py′ , pz′ = 0). Keeping smooth penalties on the normal component (Dx′)helps for gradual readjustment of the fold position. Sparsity on the model valuesare also used (ps = 0). Figure 5.10 presents sections through the final model. Thedip and lateral continuity of the layer has improved remarkably, which has beenachieved with only three-point constraints at the surface5.3 SummaryIn this chapter, I have successfully combined sparse norms and rotation of the ob-jective function for the modelling of geological bodies with oriented edges. I haveshown that it is possible to use point structural measurements, in the form of strikeand dip angles, to constrain the model. The recovered fold model closely matchedthe true solution even though I did not specify the position of the anomaly or itsphysical property contrast. In this regard, sparse rotated norms can be seen as a softgeological constraint that requires only a general understanding of the geometry ofthe problem. The resulting model, or suite of models, could subsequently be used98Figure 5.10: Recovered density model from the mixed-norm regularization.Data locations (dots), topography and the outlines of the true model (black) areshown for reference. Rotated `p-norm oriented along the folded layer helped inresolving a continuous geological unit at a reference for more advanced modelling efforts.The main challenge with the methodology presented so far is in generatingrotation parameters in 3D. Simple extrapolation of scarce structural measurementmay not accurately capture complex geological settings involving multiple foldingand faulting events. For such complicated cases, such as the Kevitsa deposit intro-duced in Chapter 1, some level of 3D modelling and interpretation may be requiredby experts. Surface mapping and other types of geophysical data may be availableto infer strike and dip information. In the following chapter, I will investigate waysto learn about preferential orientations using geophysical inversion.99Chapter 6Automated mixed-norm modelingI have so far developed a methodology that can generate a suite of models thatadequately fit the data and, at the same time, sample model space in a consistentfashion. This flexibility promotes two objectives. First, inverting for diverse mod-els provides immediate insight about non-uniqueness and that can prevent prac-titioners from over-interpreting a single inversion result. This is one of the mostimportant aspects of my work.The next objective is more challenging. I want to use different combinationsof norms on different parts of the model space so that I obtain solutions that aregeologically informative. As exemplified with the residual shown in Figure 3.19,certain assumptions may be better suited to fit portions of the data. However, theselection of different norms to be used in different portions of the model can rapidlybecome overwhelming for large problems in complex geological settings. A semi-automated approach that can capture the heterogeneity of the Earth, with minimalintervention from an expert, is needed. I am proposing the following strategy:1. Run a suite of inversions over a range of p parameters applied to the modelnorm and its spatial gradients. I will assume that the suit of models is arepresentative sampling of diverse solutions from the `p-space.2. Form an average model (mA) that captures most of the variability in thesolution space3. Extract p parameters or rotation angles for windowed portions of model100space based on the correlation between the average mA and individual mod-els.4. Carry out a final SVMN inversion with local parameters.There are numerous ways to implement the above strategy. Here, I use Princi-pal Component Analysis. I illustrate my approach with a 2D seismic tomographyexample and on the 3D fold example introduced in Chapter 5.6.1 Selecting local p-parametersI will first demonstrate how it is possible to extract local p parameters for the simul-taneous recovery of smooth and compact features as introduced in Chapter 3. Todemonstrate my methodology I use the 2D travel-time tomography example pre-sented in Sun and Li (2014). The synthetic velocity model shown in Figure 6.1(a)is made up of a smooth velocity high and a blocky velocity low centered at a 1,000m along the x-axis; the background velocity is 2000 m/s. Contour lines, corre-sponding to the 25th and 75th percentile values for both anomalies, are also plottedfor reference. An array of 11 transmitters and 13 receivers is positioned on eitherside of the model which is discretized into 32×64 square 25 m cells. First arrivaldata are calculated for each transmitter-receiver pair by taking the line integral ofslowness (reciprocal of velocity) along each ray path; this yields a total of 143observations. Gaussian noise of 5% is added to the simulated data shown in Fig-ure 6.1(b). I will first attempt to invert this data with mixed-norm penalties that areapplied to the entire model.6.1.1 Model space inversionsHaving dealt with scaling issues between `p-norm measures, I can now exploit thefull flexibility of the regularization function in (3.3). My goal is to generate a suiteof models under a broad range of assumptions. For this 2D problem, the objectivefunction to be minimized takes the form:minmφ(m) = ‖F m−dobs‖22+β ∑r=s,x,yαr‖Wr Vr Rr Dr m‖22s.t. φd ≤ φ ∗d(6.1)101Figure 6.1: (a) Synthetic 2D travel-time tomography problem made up of a rect-angular velocity low, and a smooth Gaussian velocity high. Contour lines (25thand 75th percentile) are shown in black for the true position and shape of thevelocity anomalies. Transmitters (red) and receivers (green) are positioned onopposite sides of the model domain. (b) Data map for the 143 line integral datacalculated for each transmitter-receiver pair.To demonstrate the flexibility of S-IRLS algorithm, I carry out nine inversions,using a combination of norms on a range of ps, px, py ∈ [0,1,2] values (px = py inall cases). The solutions, nine in total, are presented in Figure 6.2. All models havea final misfit φ ∗d ≈ 143 and use the same `2-norm solution to initiate the S-IRLSsteps. I make the following observations:• The two velocity anomalies centered at 1000 m on the x-axis are dominantfeatures.• Anomalies are generally stretched along the ray path due to the experimentalbias, even though I have attempted to compensate for it with the sensitivityweighting Wr. It is less pronounced for p≤ 1.• There is a progressive transition from a smooth model (upper left) to a blockysolution (lower right) as ps, px and py decrease.• The upper body (velocity low) appears to be most often represented as ablocky body with sharp edges.• The lower anomaly (velocity high) tends to be more smooth.102Figure 6.2: Suite of inverted models for various combinations of norms ps, px =py ∈ [0,1,2] (αs = αx = αy = 1). Contour lines (25th and 75th percentile) areshown in black for the true position and shape of the velocity anomalies.• Away from the anomalous regions the velocity is relatively smooth and closeto the background reference model of 2000 m/s.Figure 6.3 presents the data residual maps predicted by these nine inversions.The random noise originally added to the data is shown in Figure 6.3(j) for com-parison. I note that the strongest correlated residuals appear as pant-legs, whichcorresponds to ray paths travelling through the velocity anomalies. The residualtrend associated with the blocky velocity low (low Tx# to high Rx#) is generallymore obvious, which indicates that the experiment may be more sensitive to uni-form velocity anomalies with sharp edges as it is the source of coherent signal.This pant-leg vanishes as p, q→ 0. Meanwhile, the residual trend associated withthe smooth velocity high (low Rx# to high Tx#) diminishes for ps = px = py = 1.This change in data residual is a strong indicator that a specific combination ofnorms may be better suited to represent individual velocity anomalies.103Figure 6.3: (a) to (i) Data residual contour map for the nine inversions presentedin Figure 6.2. Contour lines (±1 standard deviation) are shown in black. (j)Random noise added to the experiment.6.1.2 Average PCA modelI assume that the suite of models presented in Figure 6.2 contains sufficient vari-ability to be representative of my solution space and I wish to use these to extractthe main features. There are numerous approaches to accomplish this and they varyin complexity from simple averaging to advanced machine learning algorithms. Iuse a Principal Component Analysis (Hotelling, 1933; Pearson, 1901). Consid-ering each model in Figure 6.2 as a data vector, the principal components of mysolution space can be written as:[m1 , ... ,mnM]≈ A W . (6.2)such that PCA vectors along the columns of A ∈ RnM×nV contains a subset of nVeigenvectors spanning the model space. These vectors encode the principal sourceof variation across the nine models recovered. The corresponding weights W ∈104Figure 6.4: PCA vectors covering 75% of the model variances.RnV×nM, also known as loadings, are scalars relating to how each principal com-ponent contributes to each model. I use the PCA algorithm from the open-sourcePython library Scikit-Learn.decomposition.PCA (Pedregosa et al., 2011).The number of principal components to be used in the analysis is determined by theexperimenter. Figure 6.4 presents the four largest principal components, coveringin this case over 75% of the variance present in the model space.Next, I generate a representative model by computing a weighted averagedmodel based on the positive PCA loadings such that:mA =∑nVi=1∑nMj=1Wi jmi∑nVi=1∑nMj=1Wi j. (6.3)The average model mA is presented in Figure 6.5. I note the close resemblance withthe true model. The background is fairly uniform near 2000 m/s; the low velocityanomaly appears to be a block with a velocity near 1900 m/s, and the high velocityanomaly appears is a smooth feature with a maximum velocity near 2100 m/s.Since mA is not the result of inversion, I have no guarantee that this model cansatisfy the observed data. The corresponding normalized data residual map shows105Figure 6.5: (a) Averaged PCA model mA and (b) corresponding normalize dataresidual map. Contour lines are shown in (a) for the true position and shape ofthe velocity anomalies.some level of correlation with the anomalies (Figure 6.5(b)). My next objective isto identify inversion parameters that will allow me to fit both the data and featureshighlighted by the average model.6.1.3 Parameter extractionNext, I want to extract optimal inversion parameters on a cell-by-cell basis in orderto best describe local features. To do so, I resort to a pattern recognition approach.In order to remove biases towards extreme model values, I transform my modelspace into a simpler parametric representation. I use a Canny edge detection al-gorithm from the open-source library Scikit-Image.feature.Canny (Pe-dregosa et al., 2011). Figure 6.6 shows the parametric edges extracted from all nineinversions.From this simplified representation of each model, I perform a moving windowcorrelation rmimP between the average PCA model mA and each of the mi solutions:rmimP =∑nj=1(mi j− m¯i)(mP j− m¯P)√∑nj=1(mi j− m¯i)2∑nj=1(mP− m¯P j)2(6.4)where m¯i and m¯P are the average model and PCA model values inside the window106Figure 6.6: Parametric representation of the nine inverted models derived formthe Canny edge detection algorithm.Figure 6.7: Local p-values on the (a) the model norm φ pss and (b) model gradi-ents φ pxx , φpyy extracted from the solution space.denoted by the subscript j. The parameters ps, px, py associated with the highestcorrelation are used in a weighted average as defined in (6.3). The process is re-peated over the entire model space. For my example I use 20×20 pixels window.The recovered ps and px, py values are presented in Figure 6.7(a) and (b) respec-tively. I note that the norm on the model gradients is generally larger in the bottomregion of the model corresponding to the location of the smooth positive anomaly.107I now define Scaled-IRLS weights on a cell-by-cell basis as:Rˆs = diag[γs((m(k−1)−mmre f )2+ ε2s)◦ps/2−1]1/2,Rˆx = diag[γx((Dx (m(k−1)−mmre f ))2+ ε2x)px/2−1]1/2,Rˆy = diag[γy((Dy (m(k−1)−mmre f ))2+ ε2y)py/2−1]1/2,(6.5)where ps, px and py define the element-wise Hadamard power and γs, γxand γy are the element-wise scaling multiplications defined in (3.52). In thisfashion, the approximated mixed norm regularization can vary on a cell-by-cellbasis.Finally, I proceed with my SVMN inversion. The data are now inverted us-ing the extracted local parameters, applied on a cell by cell basis, and the result isshown in Figure 6.8(a). There is good correspondence with the true model: boththe low-velocity blocky anomaly and the high-velocity smooth anomaly are im-aged at the right location, with the right shape and near their respective seismicvelocity. The normalized data residual in Figure 6.8(b) do not clearly show pant-legs associated with the targets. This is a good indication that most of the importantstructures have been recovered.6.1.4 SummaryThis example shows that, as long as the data are sufficiently sensitive to the geom-etry of causative bodies, it is possible to extract preferential inversion parameters.In the original study of Sun and Li (2014), inversion parameters (either `2 or `1-norm on the model gradients) were extracted based on the data residual. In thisstudy, I proposed a different path using modelling trends identified in the solutionspace. I would argue that the modelling route is a more robust approach and easilyapplicable to other geophysical methods (EM, seismic) as complex signals frommultiple sources can be unravelled by the inverse process.It is also important to note that I have achieved my final result without direct in-put from the user, other than setting tuning parameters used in the pattern recogni-108Figure 6.8: (a) The final SVMN inversion result and (b) normalized data resid-ual map. Contour lines are shown in (a) for the true position and shape of thevelocity anomalies.tion phase. The modelling process is therefore completely replicable. The learningprocess remains sensitive however to choices made by the user. The size and shapeof the averaging window is an important aspect that can significantly impact theoutcome. Choosing a window that is too large can lump smaller features together,while a window that is too small might introduce unwanted variability. The workpresented here does provide however an interesting basis for more advanced workin machine learning.6.2 Dip and strike estimationIn Chapter 5, I have demonstrated the benefits of using rotated sparse gradient toaccentuate geological trends and better image continuous bodies at depth. I haveused the fold model shown in Figure 6.9 to showcase the use of surface structuraldata. In most greenfield exploration settings however such data may either not beavailable or too scarce for accurate interpolation. Practitioners would benefit frombeing able to estimate the strike and dip of geological features in a semi-automatedway. In this section, I propose to repeat the rotated sparse norm experiment througha learning process. I will make use of a similar strategy elaborated in the previoussection such that preferential orientations are extracted from the solution spacewith a pattern recognition approach.109Figure 6.9: Synthetic fold model introduced in Chapter 5 used here to test thedip and strike estimation learning algorithm.I propose to estimate the angle of rotation in two stages: first along a horizontalplane to determine the strike, then on vertical sections to estimate the dips. I createa workflow divided into five steps:1. Run an unconstrained inversion with sparse and smooth assumptions (ps = 0,p[x,y,z] = 2)2. Estimate horizontal trends and rotate the objective function on the xy-plane3. Run a suite of inversions over a range of dip angles (ps = 1, px = 2, p[y,z]= 0.4. Form an average model (mA) and extract dip angles locally.5. Carry out a final SVMN inversion with local p parameters and rotation an-gles.I apply this procedure to recover the folded density layer used in Chapter 5 (Fig-ure 5.5). I will attempt to do this without any external input from the user otherthan general tuning parameters in Step 2 and 4.110Step 1: Unconstrained inversionPotential field data are most sensitive to lateral changes in density and magneticproperties, so it is generally easier to determine horizontal trends related to geol-ogy. I am using a mixed norm inversion (ps = 0, p[x,y,z] = 2) to get a first estimateof the geometry of the problem such that physical property contrasts are increasedwithout enforcing assumptions about directionality. While similar analysis couldbe performed directly on a gridded data map at a much smaller computational cost,the inversion route has proven experimentally to be more robust. It can deal withtopographic effects and it can isolate long wavelength trends from near surfaceanomalies. Since I already have inverted this synthetic example with sparsity as-sumptions in Chapter 5 (Figure 5.8), I can advance directly to Step 2.Step 2: Azimuth estimationFrom the recovered model, I want to extract directional information. I resort toan image moment algorithm to extract dominant patterns on the windowed por-tion of the model. I first demonstrate the procedure on the UBC Crest shown inFigure Density values are converted to a binary image using the python routineScikitImage.Skeleton to extract dominant feature (Lee et al., 1994). The2D image is reduced to one pixel wide representations by recursively identifyingand removing border pixels. The process is carried out until no pixels can beremoved without breaking the connectivity of a given feature.2.2: Position and orientation of dominant features are extracted with an imagemoment approach (Hu, 1962). For a set window of the image with pixel intensitydefined by I(x,y), the image moment is given bymi j =∑x∑yxiy jI(x,y) . (6.6)The angle defining the principal axis is calculated byθ =12arctan(2m11m20−m02)(6.7)111Figure 6.10: Example of an automated azimuth estimation using the UBC Crest.(b) The image is first converted to a binary image using the python routineScikitImage.Skeleton. Image moment calculations are performed onmoving windows to determine the position and orientation of dominant features.(c) Vectors are then rotated by 90◦ to point along the normal of edges. Thenormal vectors are checked for consistency in preparation for interpolation.and the center of mass of the image is given byxc =m10∑x∑y I(x,y)yc =m01∑x∑y I(x,y),(6.8)Calculations are only performed for windows that contain at least 5% of non-zeropixels. This process is repeated across the entire section by incrementally movingthe window location. Efficiency can be gained by displacing the window basedon the calculated orientation. This results in a collection of vectors pointing alongwith the main features, as shown in Figure 6.10(b).2.3: Orientation vectors are rotated by 90◦ to point normal to the edge of fea-tures (Figure 6.10(c)). Neighbouring normals are compared to each other for con-sistency, such that vectors pointing against the general trend are rotated by 180◦.This final step is important for a smooth interpolation of angles in 3D space. Incomplex geological settings, this step may require quality control checks from theuser.112Figure 6.11: (a) Image moment calculations over windowed regions of the in-verted density model. Arrows are plotted for the principal direction (red) andnormal (white) of the main density anomaly. (b) Normal vectors are interpo-lated in 3D and used to rotate the objective function.Carrying out the process described above on a horizontal section of the densitymodel, I recover the orientation vectors presented in Figure 6.11. In this case, Ichoose a square window 200 m in width that sweeps through a draped section ofthe density model at 40 m depth (two vertical cells). This choice depends on thediscretization and wavelength information present in the data. Going down a fewcells below topography generally offers a good compromise between damping thenear-surface variations while preserving major trends. Just as for the structuraldata in Chapter 5, I interpolate the normal vectors in 3D using minimum curvatureapproach (Briggs, 1974).Step 3: Model space inversionNow that I have determined the preferential orientation (normal) of geologicalbodies, I will attempt to extract preferential dip directions. I proceed with a se-ries of 12 inversions with rotated sparse norms over a range of rotation angles forθ ∈ [−90◦, 75◦] at 15◦. Figure 6.12 compares the recovered models along the EWsection B-B’.113Figure 6.12: Suite of models with directional gradients for θ ∈ [−90◦, 75◦].Step 4: Dip estimationFrom the twelve recovered models, I proceed with PCA and form an average modelas prescribed in Section 6.1.2. The calculated average model is shown in Fig-ure 6.13. From this average model, I extract dip information along vertical sec-tions using the image moment strategy presented in Step 2. The resulting anglesare extrapolated in 3D and used to rotate the objective function.Step 5: SVMN inversionIn the final step, I re-invert the dataset with rotated sparse norms. The final so-lution is shown in Figure 6.14. I note the good match between the solution andthe true position of the dense layer. I have obtained this solution without externalinput, other than the window size and threshold value used in the image momentestimation.114Figure 6.13: Sections through the average PCA model and rotation vectors (ar-row) derived from the image moment analysis.Figure 6.14: (a) Horizontal and (b, c) vertical sections through the recoveredmodel using rotation angles derived from the average PCA model.1156.3 SummaryLeveraging the mixed norm regularization function of Chapter 3 and pattern recog-nition algorithms, I have explored ways to extract inversion parameters from thesolution space. The assumption is made that the shape and position of robust fea-tures can be highlighted, as long as the data is sufficiently sensitive to their ge-ometry. I have shown that gravity surveys could be used to determine the strikeand dip of geological features with detectable edges. I have demonstrated this bysuccessfully imaging a continuous and folded density layer. This process was onlypossible because the edges of the anomaly were detectable from the surface. Simi-lar procedures would not have been successful in a strictly layered environment towhich potential fields are not sensitive.The rotation parameters are soft geological constraints that reinforce modellingtrends without having to specify the location of physical property contrasts. This isan appealing property of the algorithm for greenfield settings as the implementationrequires little input from practitioners and easily be automated.The methodology developed in this chapter is general and can easily be ex-tended to other geophysical methods, such as electromagnetic and seismic data.The rotation parameters could potentially be derived from multiple physical prop-erties. Preferential orientations could form a basis for joint or cooperative physicalproperty inversion.116Chapter 7Case Study - Kevitsa Ni-Cu-PGEdepositIn the previous chapters, I have provided technical developments to explore themodel space, interpret a suite of models and extract dominant inversion parame-ters. The different components were tested on synthetic examples that showcasedspecific aspects. In this final chapter, I demonstrate the practical implementationof these advances on gravity and magnetic data sets acquired over the Kevitsa Ni-Cu-PGE deposit. The complex geology of Kevitsa makes it an ideal candidate toshowcase my research.7.0.1 Geological settingThe Kevitsa deposit was discovered in the mid-1980’s through extensive explo-ration programs sponsored by the Geological Survey of Finland. The proven 160million tons of nickel is economically significant in the region as it is expected tobecome the largest mine in the country. Figure 7.1 presents a simplified geolog-ical map of the Kevitsa-Satovaara intrusive complex adapted from Koivisto et al.(2015). The ore deposit is hosted in a funnel-shaped ultra-mafic olivine pyroxen-ite (UPXO) unit dipping towards the southwest, dated to approximately 2054 Ma(Gregory et al., 2011; Mutanen, 1997). Directly adjacent to the southwest marginis a large gabbro (IGB) unit formed during a late phase of the intrusion. Similargabbro, found west of the Satovaara Fault, is likely related temporally. The body117Figure 7.1: (a) Geological surface map of the Kevitsa-Satovaara intrusive com-plex and (b) 2D seismic reflection survey along with the E5 profile. Inter-preted seismic reflectors and geological contact (black) from the interpretationof Koivisto et al. (2015) are shown for reference, as well as the outline of thegravity and magnetic data sets analyzed in this chapter.118intruded a mostly layered sequence of mafic volcanic (MVS) and carbonaceousphyllites (MPH) units. Discontinuous layers of arkose (ARK), arenite (ARN) andfelsic volcanic (FVS) units are interbedded within the volcanic sequence. Thecurrent folded configuration of the geology is due to large tectonic events that de-formed the host stratigraphy around the more competent UPXO block.It is believed that the disseminated sulphide mineralization within the UPXOwas precipitated from the dissolution of Ni-Cu rich proterozoic metasediments,referred to as black schist (MPHB) (Mutanen, 1997). Mineralogical and texturalchanges observed on core samples within the intrusion reveal the highly heteroge-neous nature of the deposit. A highly altered dunite raft (UDU), referred to as thecentral dunite, is located between the IGB and UPXO. It has likely been partiallyassimilated by the UPXO intrusion.Seismic surveys acquired over the deposit have been studied extensively andyielded the most complete picture of the deposit to date. Koivisto et al. (2015)interpreted the profiles along side drill hole data to produce a partial 3D surfacemodel of the deposit. The bottom limit of the intrusion is marked by a series ofstrong reflectors, interpreted as footwall dunite units and host stratigraphy (Fig-ure 7.1). Additional UDU units have been intercepted by two deep holes, but itis unclear if this is related to the central dunite or to the komatiite (UKO) layersfound in the host stratigraphy. The southern continuation of the Kevitsa intrusioncan clearly be seen along the E5 seismic profile, as a thick 1km non-reflective zone,but no bore hole is available to confirm the lateral extent.Petrophysical dataAlong with geophysical data, physical property measurements collected along 279bore holes were made available for analysis; they include density, magnetic suscep-tibility, conductivity, seismic velocity measurements. Figure 7.2 presents whiskerplots summarizing the measured density and susceptibility grouped by lithologicalunits. The interpretation of core logs is challenging due to the complex geologyand broad terminology used to describe the rocks. In order to focus my analysis, Igeneralize the rock classification into 10 main geological groups as summarized inTable 7.1. These groups were defined in terms of relative age and expected physical119property contrasts. I note a few obvious trends:• Intermediate and felsic volcanic rocks (VIO) show low density and moderateto low magnetic susceptibility• The komatiite (UKO) and dunite (UDU) show significantly higher suscepti-bility but low density• The intrusive rocks (UPX) and (IGB) are both dense and susceptibleOverall I should expect high physical property contrasts between the differentlithologies. While I will not directly incorporate this information in the inversion,it will be used to define research questions and for the final interpretation.7.1 Geophysical dataKevitsa is an interesting case study considering the large collection of data sets thatwere acquired and made available to researchers: bore hole petrophysical measure-ments (Montonen, 2012), seismic refraction (Koivisto et al., 2015), direct-currentresistivity and magnetotelluric (Le et al., 2016), airborne time-domain EM surveys(VTEM 2009, SkyTEM 2010). In this section, I focus on the ground gravity datasetand airborne magnetic data collected by the VTEM survey shown in Figure 7.3(a)and (b) respectively.From simple visual inspection, I note some obvious connections between thepotential field data and the surface geology:• High gravity and magnetic data correlated with the UPXO and hydrothermalBXH• Moderate response from the VMO• Weak fields over most of the MPHB units.• Magnetic high and gravity low over the komatiite (UKO),• Large negative magnetic anomaly within the UDU and near the southernedge of UPXO120Code Description Density SusceptibilityOVB Overburden Low MediumIGBIGB: GabbroIPG: PegmatiteIGBO: Olivine gabbroIDI: DioriteIDB: DiabaseIMO: Intrusive (mafic)IGBM: Magnetite gabbroMediumMediumMediumHighUPXUPXO: Olivine peroxiniteUOO: Ultramafic (Undiff.)MPE: MetaperidotiteUWB: WebsteriteMedium MediumUDUUDU: DuniteUPE: PeridotiteUSP: SerpentiniteLow HighUKO Komatiite Low HighBXHBXO (undiff.)BXHC: Hydrothermal (crackle)LowHighMediumSOOMAB: AlbititeMAM: AmphiboliteMQZ: QuartziteMSCSD: SchistLow LowMPHMPH: PhylliteMPHB: Black PhylliteMSCBK: Black SchistMHF: HornfelsHigh LowVMOVMO: Volcanic MaficVBA: BasaltVTUM: Volcanic tuffHigh LowVIOVIO: Volcanic IntermediateVTUI: Volcanic tuffVAN: AndesiteVOO: Volcanic (undiff.)Low MediumTable 7.1: Summary table grouping the various lithological units logged frombore hole. Expected density and magnetic susceptibility contrasts are derived fromFigure 7.2.121Figure 7.2: Whisker plot of logged (a) density and (b) magnetic susceptibilitymeasured along 279 boreholes. The coloured boxes have a width scaled by thecalculated standard deviation and centered on the mean value for all interceptsbelonging to the same lithological classification. The black lines on either sidedefine the minimum and maximum values. The different lithologies are colourcoded and grouped based on relative age and similarities in physical properties.The strong negative magnetic fields observed over the dunite unit are of particularinterest for this study, as it is likely related to remanent magnetization. Analy-sis from core samples by Montonen (2012) reported large Keonigsberger ratiosand reversed magnetization direction in the UDU unit as summarized in Table 7.2.It is important to note that large Koenigsberger ratios were also measured in the122Figure 7.3: (a) Ground gravity and (b) airborne magnetic data map acquiredover the Kevitsa Ni-Cu-PGE deposit, Finland. Color ranges are histogram equal-ized and sun shaded to highlight details. The location of the Kevitsa mine (red),geological boundaries (black) and faults (dash) identified from surface mappingare shown for reference.123Hole ID Interval κ (SI) Inc. (◦) QKV297 0-52.9 0.034 −42.4±18.9◦ [2,10]KV200 29.9 0.038 −50.9 5.4Table 7.2: Intervals along boreholes KV200 and KV297 reporting significantremanent magnetization (Montonen, 2012)lower UPXO unit, although susceptibility values remained small. In the absenceof oriented core, no magnetic declinations were provided. From forward mod-elling of magnetized sheets, Montonen (2012) estimated that a magnetized unitwith effective susceptibility 0.82 SI and orientated [I = −42.5◦,D = 240◦] couldbe responsible for the observed negative magnetic anomaly.7.1.1 Modeling objectivesDue to the complex and disseminated nature of the mineralization, substantial mod-elling efforts continue to be invested to better understand the intrusion. Most of thedrill holes made public are concentrated within the UPXO, which leaves a fewquestions unanswered1. Can potential fields be used to define the 3D geometry of the intrusion?2. Can we distinguish variations in physical properties within the UPXO relatedto mineralization?3. Can we characterize the central dunite and possible extension at depth?4. Can the magnetic response be used to infer tectonic deformation?7.2 Gravity inversionI first invert the ground gravity data set acquired over the Kevitsa deposit shownin Figure 7.3(a). I will attempt to model the 3D geometry of the intrusion andhost stratigraphy based on density contrasts. The survey consists of 28,860 gravitystations collected at 20 m intervals along east-west and north-south lines spaced at100 m. Data were terrain corrected with a reference density of 2.67 g/cc. I design124an OcTree mesh for the inversion with 25 m cells in the core region extending 2km at depth. After several inversion trials, data uncertainties are set to 0.1 mGal.As a first pass, I invert the data with the linear `2-norm regularization. Fig-ure 7.4(a) presents a section through the recovered density model at ≈ 200 mdepth. I also extract a vertical section along the seismic profile E5 in Figure 7.4(b).I note that the bulk of the density anomaly is located inside the upper portion ofthe UPXO. As expected from smooth assumptions, the edges of the anomaly arepoorly defined.Clearly the assumption of a smooth density distribution is not sufficient to es-tablish a clear correlation between the known lithology and density. Taking ad-vantage of the methodology developed in this thesis, I will attempt to improve thedefinition of geological domains. First I need to define rotation angles in 3D tobetter represent the folded geology of Kevitsa. In this case, I do not have to learnthe orientation as a I have access to geological interpretation at the surface and aseismic profile at depth. The folded geology of Kevitsa appears to be more or lesssymmetric radially around the center of the intrusion. I will, therefore, assume thatvector components extracted from the vertical section can be interpolated in 3D.Using the method of image moment introduced in Chapter 6, I extract orientationof major contacts over the study area. Figure 7.5 presents the orientation vectorsused for the rotation of the objective function.I then proceed with nine inversions over a range of mix-norm penalty functionsfor ps, py′,z′ ∈ {0, 1, 2}, px′ = 2. Only the normal component of the geologicalcontacts is fixed to the `2-norm. Horizontal and vertical sections through the ninerecovered density models are shown in Figure 7.6 and Figure 7.7 respectively.Residual data maps are shown in Figure 7.8. Correlated residuals indicate thatsome of the short short wavelength information is under fitted.In order to simplify the interpretation, I overlay iso-contour values of low andhigh density using the 5th and 95th percentile of density values for each of the nineindependent inversions. From Figure 7.9 I make a preliminary interpretation.• The recovered density within the UPXO unit appears to be highly variable.The known ore deposit appears to be sitting within a region of moderatedensity.125Figure 7.4: (a) Horizontal section at ≈ 200 m depth below topography throughthe recovered density model after convergence of the algorithm for ps = px =py = pz = 2. (b) Vertical section of the density model overlaid on the 2D seismicreflection line E5. Interpreted geological contacts (black) are shown for compar-ison.126• Density lows are strongly correlated with mapped schist (SOO) and phyllite(MPH) units on the outer perimeter of the intrusion. The same does not applyto the black schist unit (MPHB) found immediately south of Kevitsa. WhileMPH and MPHB appear to have been used interchangeably to build the ge-ological map, the recovered density appears to be a distinguishing property.The most important aspect of this result is the general trend observed for the mainintrusion that extends south-east below the gabbro. This result is interesting as Imanaged to recover models that are in good agreement with the seismic interpre-tation of Koivisto et al. (2015), even though no information about the depth of theintrusion was included in the inversion. So f t geological constraints, in terms ofgeneral trend, were sufficient to push the lateral extent of the density anomaly inaccordance with the gravity data.7.3 Magnetic inversionI follow my analysis with the processing of magnetic data collected during the2009 VTEM survey (Figure 7.3(b)). The inducing field parameters at the time ofacquisition were B0 = [A : 52,800 nT, I : 77.5◦,D : 12.2◦]. Similar to the gravityinversion, I design an OcTree mesh with 50 m core cells extending 2 km at depth.I determined experimentally a 10 nT uncertainty floor value to the data.I first invert the TMI data with the conventional smooth susceptibility assump-tion (ps = px = py = pz = 2) and ignore the effect of remanence. From sectionsthrough the recovered susceptibility model, presented in Figure 7.10(a) and (b), Inote discrepancies with the known geology:• In plan view, the arc shaped anomaly, SW of the deposit, is recovered outsidethe mapped UKO unit.• Along the E5 seismic section, no susceptibility anomaly is recovered overthe central dunite. This directly contradicts the core sample measurementsmade by Montonen (2012).• The shape and extent of the large anomaly correlates poorly with the UPXOunit interpreted by Koivisto et al. (2015) from seismic reflectors.127• The residual data map shows strong correlation with the location of strongnegative magnetic fields.I conclude that the magnetic response observed at Kevitsa cannot solely be at-tributed to an induced magnetization.To address issues posed by remanence, I proceed with the MVI-S algorithm.I perform a series of nine inversions with varying sparsity measures to assess thevariability in the magnetization model. Starting from a common `2-norm MVI-Cmodel, I sequentially vary the combination of norms applied to the amplitude andits derivatives for (pρs , pρxyz ∈ [0,2]) and I apply the same rotation parameters aspreviously used for the gravity inversions. In this regard, the sparsity and rotationparameters tie the density and magnetization inversions together. While not jointlyinverted, I expect some correlation between the two physical properties. Horizontaland vertical sections through the recovered nine magnetization models are shownin Figure 7.11 and 7.12 respectively. Residual data maps are shown in Figure 7.13for all nine inversions.To simplify the analysis, I superimpose the 90th percentile iso-value of ampli-tude for each of the nine models (Figure 7.14). I calculate an average magnetizationdirection (white) and standard deviation on the angle (red). I observe the following:• The known ore deposit appears to sit within a volume of low magnetization.• Parts of the central dunite unit appear to be reversely magnetized [κe =0.09 SI, I =−52◦±15◦, D= 246◦±5◦] This is in excellent agreement withthe laboratory results published by Montonen (2012)• The vertical magnetic anomaly between the IGB and UPXO unit, likely re-lated to the central dunite unit, appears to be plunging towards SE, poten-tially extending below the UPXO unit as hypothesized by Koivisto et al.(2015).• Strong magnetization recovered along the outer-shell of the ultra-mafic intru-sion appears to be pointing radially outward. Largest magnetization appearsto originate below the UPXO unit.128• Similar magnetization direction pointing normal to the arc-shaped peridotiteunit.The last two remarks are interesting for a few reasons. First, strong magneti-zation below the base of the ultra-mafic supports the presence of magnetic UDUlayers. Secondly, the orientation of magnetization pointing normal to the base ofthe unit may be indicative of past tectonic deformation. Under the assumption thatthe remanent magnetization component had been fairly uniform within the layeredUDU, UKO and UPXO unit at the time of formation, then the current radiatingmagnetization pattern would be caused by subsequent folding of the units. If thisis the case, then this is one of the first times that airborne magnetic data would havebeen used as a marker for tectonic deformation of elongated and folder geologicalunits.While my modelling of the central dunite unit agrees with published laboratorymeasurements, the cause for this reverse magnetization direction remains unclear.No other lithological units appear to share the same orientation. Re-magnetizationafter the emplacement of the ultra-mafic intrusion is unlikely as a similar reversedpolarity pattern would also be expected elsewhere at Kevitsa. I speculate that thedunite block could be related to the lower UDU unit, which would have been ro-tated to its current sub-vertical location.7.4 SummaryThe inversion of gravity and magnetic data over the Kevitsa deposit yielded valu-able information. First from the gravity inversion, I have imaged the UPXO atdepth and extending below the gabbro unit. Second, variations in density and mag-netization inside the UPXO unit is of interest as it appears to correlate well with theknown mineralization. This could potentially be used for future exploration work.Third, the average magnetization model confirms that the central dunite unitis associated with strong reversed magnetization oriented roughly [κe = 0.09 SI,I = −52◦±15◦, D = 246◦±5◦]. There is also a strong indication that it could beconnected to the lower dunite units, although the magnetization direction wouldrequire the central dunite to have undergone a rotation of almost 180◦.Potentially the most significant outcome of this case study is the recovered129magnetization pointing normal to the base of the UPXO unit. A similar radiatingpattern is observed in the komatiite unit east of the deposit. If confirmed by labo-ratory measurements, this result would be one of the first cases of paleomagnetismbased on the inversion of airborne magnetic data over folded geology.130Figure 7.5: Rotation parameters (normals) derived from (a) surface geology and(b) major reflectors detected on the seismic profile E5.131Figure 7.6: Horizontal sections at ≈ 200 m depth below topography throughthe density models recovered over a range of `p-norms applied on the model andmodeled derivatives for ps, px,y,z ∈ {0, 1, 2}132Figure 7.7: Vertical sections overlaid on the 2D seismic reflection line E5through the density models recovered over a range of `p-norms applied on themodel and modeled derivatives for ps, px,y,z ∈ {0, 1, 2}133Figure 7.8: Residual data maps for the nine gravity inversions.134Figure 7.9: (a) Horizontal and (b) vertical section comparing iso-surfaces forthe 5th and 95th percentile of anomalous density values recovered from the ninemixed norm inversions for ps, px,y,z ∈ {0, 1, 2}.135Figure 7.10: (a) Horizontal and b) vertical sections through the recovered sus-ceptibility model that ignores the effect of remanence. Lithological contacts(black) identified by Koivisto et al. (2015) are shown for reference. A largedome-shaped zero-susceptibility anomaly is recovered at the center of the Ke-vitsa intrusion, likely due to remanence. (c) Residual data map shows strongcorrelation with the negative anomalies.136Figure 7.11: Horizontal sections at≈ 300 m below topography for a suite mod-els using various sparsity assumptions put on the amplitude of magnetizationfor ps , px,y,z ∈ [0, : 2]. Norm measures on the magnetization angle are fixed topx,y,z = 0 in order to promote uniform magnetization137Figure 7.12: Vertical sections through the recovered magnetization models. Ef-fective susceptibilities are concentrated within the UDU and on the lower mar-gin of the UPXO unit.138Figure 7.13: Residual magnetic data maps for the nine MVI-S inversions.139Figure 7.14: (Top) Horizontal and (bottom) vertical sections through iso-contours of magnetization recovered from nine mixed `p-norm inversions. Mag-netization orientation (white) and standard deviation on the angle (red) areshown.140Chapter 8ConclusionThe overarching goal of this research thesis was to facilitate the interpretation ofpotential field data acquired over complex geology and to extract as much infor-mation out of the data through a semi-automated learning process. These goalwas motivated by technical limitations encountered with conventional inversionmethodologies. Geophysical inverse problems are inherently non-unique and thecharacter of the solution depends on assumptions set by the user. Smooth physicalproperty inversions yield models that poorly represent sharp geological contact,while sparsity assumptions generally yield simplistic blocky anomalies. The rangeof possible solutions available for interpretation has generally been fairly limited.More experienced users may be able to improve the interpretation with geologicalconstraints, but building and testing different scenarios remains a laborious processthat is difficult to track. In this regard, the methodology presented in this thesis sitsomewhere between the blind unconstrained inversion and the expert-driven geo-logical inversion. I have developed a methodology and software that can generatea suite of models that honor so f t geological assumptions.In Chapter 2, I reviewed the numerical implementation of gravity and magneticforward modeling in integral form. I tackle numerical limitations associated withthe storage and manipulation of large dense matrices. The memory footprint of po-tential fields problem is reduced at two levels. First, I borrow the mesh decouplingstrategy previously used in electromagnetic modeling. The global forward prob-lem is broken down into tiles each associated with nested Octree mesh. Secondly,141I leveraged open-sourced technologies for out-of-core storage of dense matrices.Used in concert, the two advancements allowed to run large forward problemswithout the need for compression.In Chapter 3 I reviewed some of the work I had previously investigated duringmy Master’s thesis regarding the implementation of a mixed `p-norm regulariza-tion. I identified a new scaling strategy based on the maximum partial derivativesof individual regularization functions such that multiple penalties could impact thesolution. The robust implementation of mixed norm assumptions is at the core ofthis research as it allows me to generate a suite of solutions and explore the modelspace in a consistent manner. This chapter resulted in the research paper Fournierand Oldenburg (2019a).Chapter 4 builds upon the magnetic vector inversion in spherical coordinatesfirst introduced by Lelie`vre (2009). The algorithm has received little attention dueto the difficulty in solving the non-linear inverse problem. With knowledge gainedin Chapter 3, I developed an iterative rescaling strategy based on the maximumpartial derivatives of the sensitivity function. The decoupling of the magnetizationstrength and orientation allowed me to apply sparsity assumptions for the recoveryof well-defined anomalies with coherent magnetization direction. It is an improve-ment over methods previously published as it is can deal with complex geologicalsettings comprising multiple anomalies with arbitrary shape.Chapter 5 is a generalization of the methodology introduced in Chapter 3 forthe recovery of oriented edges. The measure of model gradients along the Cartesianaxes was a choice of convenience but it poses a major limitation for the modelingof folded and dipping geological contacts. I introduced a 7-point gradient operatormeasuring the model gradients with diagonal cells. I have also shown how scarcestructural measurements can be interpolated and used to guide the inversion.Chapter 6 takes the model space inversion further by attempting to learn fromthe suite of solutions. I investigated ways to extract optimal inversion parametersbased on recurrence of dominant features. I utilized basic machine learning toolssuch as PCA, edge detection and image moment algorithms to identify patterns andbuild constraints. I elaborated a strategy to automate dip and strike estimation ofisolated geological units.All the technology brought forward in this thesis were put to test in Chapter 7142on the Kevitsa Ni-Cu-PGE deposit. I inverted ground gravity and airborne mag-netic data using rotated sparse norms. The resulting density and magnetizationmodel highlighted important features of the deposit. First from the density model,I confirmed the seismic interpretation of an extended ultra-mafic unit under cover.The magnetization model also corroborated previous studies in modeling the re-manent component of the central dunite. The inversion also indicated a possibleconnection with other units found at the base of the intrusion. Potentially the mostsignificant outcome of this case study is the orientation of magnetization recoveredin several folded units. While other studies have attempted to use magnetizationinversion to gain insight about the geological history of a deposit, it is likely thefirst time that this was done at such a large scale and over complex geology.8.1 Limitations and future workThe work presented in this thesis represents a significant gain in flexibility for in-verting gravity and magnetic data, but this comes at an added computational costneeded to a single inversion and also for generating a suite of models. This is exac-erbated by the cooling strategy of the threshold parameter introduced in Chapter 3such that sparsity assumptions are slowly phased in. It is important to note thateach model is a valid candidate with respect to the geophysical data. These modelscould potentially be used earlier in the learning process to avoid repetitive iterationsteps.Attempting to use magnetic vector inversion for large scale paleomagneticstudies is an appealing approach that warrants further investigation. Sparsity as-sumptions have helped in recovering coherent magnetization direction inside iso-lated anomalies. This process is less stable when performed over multiple magneticanomalies with overlapping signal. The use of magnetic gradient data may help inreducing the non-uniqueness and in better defining the boundaries of anomalies.Another aspect of magnetic vector inversion that I have not addressed in this thesisis to isolate the remanent component from the total magnetization. This could po-tentially be accomplished by jointly inverting electromagnetic and magnetic datafor the modeling of susceptibility, remanence and conductivity.The learning algorithms presented in this thesis have served their purpose in143identifying trends in the model space and extracting local parameters. A criticalcomponent of this process is the pattern recognition phase. This step still reliesheavily on the user to determine tuning parameters that are currently found by trialand error. Some level of quality control is also needed when defining rotation pa-rameters such that strike and dip angles are consistent with known dip directions.More advanced machine learning algorithm could potentially improve these pro-cedures. The learning process could be performed on multiple physical propertiesand serve as a link for joint and cooperative inversions.Research dedicated to potential field data is likely to continue and grow as thequality and quantity of surveys continues to increase. While I have made inroadsin reducing the computation cost of potential fields inversions, large (continental)scale inversion remains difficult. One of the main tuning parameters that I havemostly ignored in this project is the reference model. At a large scale, choosinga single reference value is undoubtedly a gross generalization that can adverselyaffect the solution. Future work should investigate the use of model decompo-sition techniques such that the inversion is performed on both the reference andanomalous properties. Large scale model decomposition inversion could be doneat various resolutions with the use of spherical Octree discretization.The methodology presented in this thesis will be extended to other geophysicaldata. Work is currently underway to test the model space inversion on electricalresistivity problems and airborne electromagnetics.144BibliographyAbedi, M., D. Fournier, S. G. Devriese, and D. W. Oldenburg, 2018a, Integratedinversion of airborne geophysics over a structural geological unit: A case studyfor delineation of a porphyry copper zone in Iran: Journal of AppliedGeophysics, 152, 188 – 202. → pages vi——–, 2018b, Potential field signatures along the Zagros collision zone in Iran:Tectonophysics, 722, 25 – 42. → pages viAjo-Franklin, J. B., B. J. Minsley, and T. M. Daley, 2007, Applying compactconstraints to differential traveltime tomography: Geophysics, 72, 67–75. →pages 8, 34, 46, 55Barber, C. B., D. P. Dobkin, and H. Huhdanpaa, 1996, The Quickhull algorithmfor convex hulls: ACM Transactions on Mathematical Software, 22, 469–483.→ pages 18Barbosa, V. C., and J. B. Silva, 2006, Interactive 2D magnetic inversion: A toolfor aiding forward modeling and testing geologic hypotheses: Geophysics, 71,L43–L50. → pages 84Barbosa, V. C. F., and J. B. C. Silva, 1994, Generalized compact gravity inversion:Geophysics, 59, 57–68. → pages 8, 34, 46, 55Blakely, R., 1996, Potential theory in gravity and magnetic applications:Cambridge University Press. → pages xii, 2Blaschek, R., A. Hordt, and A. Kemna, 2008, A new sensitivity-controlledfocusing regularization scheme for the inversion of induced polarization databased on the minimum gradient support: Geophysics, 73, 45–54. → pages 8, 46Bosch, M., and J. McGaughey, 2001, Joint inversion of gravity and magnetic dataunder lithologic constraints: The Leading Edge, 20, 877–881. → pages 7, 29Briggs, I. C., 1974, Machine contouring using minimum curvature: Geophysics,39, 39–48. → pages 96, 113Chartrand, R., 2007, Exact reconstruction of sparse signals via nonconvexminimization: IEEE Signal Processing Letters, 76, 167–188. → pages 8, 46Clark, D. A., 2014, Methods for determining remanent and total magnetisation ofmagnetic sources - a review: Exploration Geophysics, 45, 271–304. → pages14510Cordell, L., and V. J. S. Grauch, 2012, in Mapping basement magnetization zonesfrom aeromagnetic data in the San Juan Basin, New Mexico: 181–197. →pages 3Cox, L. H., G. A. Wilson, and M. S. Zhdanov, 2010, 3D inversion of entireairborne electromagnetic surveys using a moving footprint: ASEG ExtendedAbstracts, 1–10. → pages 22Dannemiller, N., and Y. Li, 2006, A new method for estimation of magnetizationdirection: Geophysics, 71, L69–L73. → pages 10Dask, D. T., 2016, Dask: Library for dynamic task scheduling:→ pages 22Daubechies, 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 8, 34, 42Davis, K., D. W. Oldenburg, and M. Hillier, 2012, Incorporating geologicalstructure into the inversion of magnetic data: Presented at the ASEG 22thConference and Exhibition. → pages 7, 84Domzalski, W., 1966, Importance of aeromagnetics in evaluation of structuralconstrol of mineralization: Geophysical Prospecting, 14, 273–291. → pages 1Ekblom, H., 1973, Calculation of linear best lp-approximation: BIT, 13, 292–300.→ pages 8, 34Engheta, N., W. D. Murphy, V. Rokhlin, and M. S. Vassiliou, 1992, The FastMultipole Method (FFM) for the electromagnetic scattering problem: IEEETransactions on Antenna and Propagation, 40, 634641. → pages 21Enkin, R., 2003, Paleomagnetic analysis of selected Ekati kimberlite.Unpublished report for BHP Billiton Diamonds. Geological Survey of Canada.→ pages 9Enkin, 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 9, 69Farquharson, C. G., and D. W. Oldenburg, 1998, Non-linear inversion usinggeneral measures of data misfit and model structure: Geophysical JournalInternational, 134, 213–227. → pages 8, 34Fedi, M., G. Florio, and A. Rapolla, 1994, A method to estimate the totalmagnetization direction from a distortion analysis of magnetic anomalies1:Geophysical Prospecting, 42, 261–274. → pages 10Foss, C., 2006, in Parametric inversion as an advanced technique for magneticsource depth estimation: 943–947. → pages 7Foss, C., and B. McKenzie, 2011, Inversion of anomalies due to remanentmagnetisation: an example from the Black Hill Norite of South Australia:146Australian Journal of Earth Sciences, 58, 391–405. → pages 10Fournier, D., 2015, A cooperative magnetic inversion method with lp-normregularization: Msc. thesis, The University of British Columbia. → pages 10,35, 48, 71Fournier, D., K. Davis, and D. W. Oldenburg, 2016, Robust and flexiblemixed-norm inversion: SEG Technical Program Expanded Abstracts,1542–1547. → pages vi, 11Fournier, D., and D. Oldenburg, 2019a, Sparse magnetic vector inversion inspherical coordinates: Application to the Kevitsa Ni-Cu-PGE magneticanomaly, Finland: Geophysics. (In revision). → pages 142Fournier, D., and D. W. Oldenburg, 2019b, Inversion using spatially variablemixed lp-norms: Geophysical Journal International, 218, 268–282. → pages viFullagar, P. K., and G. A. Pears, 2013, 3D magnetic modelling and inversionincorporating self-demagnetisation and interactions: ASEG ExtendedAbstracts, 2013, 1–4. → pages 10Fullagar, P. K., G. A. Pears, and B. McMonnies, 2008, Constrained inversion ofgeologic surfaces - pushing the boundaries: The Leading Edge, 27, 98–105. →pages 7, 29Gorodnitsky, 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 8, 34Grant, F., 1984, Aeromagnetic, geology and ore environments, I. Magnetite inigneous, sedimentary and metamorphic rocks: An overview: Geoexploration,23, 303–333. → pages 1Gregory, J., N. Journet, G. White, and M. Lappalainen, 2011, Kevitsa coppernickel project Finland: 43-101 report, First Quantum Mineral. → pages 117Haber, C., and C. Schwartzbach, 2014, Parallel inversion of large-scale airbornetime-domain electromagnetic data with multiple OcTree meshes: InverseProblems, 30, 1–28. → pages 22Haber, E., D. W. Oldenburg, T. Farnocombe, and A. Cellar, 1997, Directestimation of dynamic parameters in SPECT tomography: IEEE Transaction onNuclear Science, 44, 2425–2429. → pages 30Haber, E., and L. Tenorio, 2003, Learning regularization functionals: a supervisedtraining approach: Inverse Problems, 19. → pages 7Helbig, K., 1963, Some integrals of magnetic anomalies and their relation to theparameters of the disturbing body: Zeitschrift fr Geophysik, 29, 83–96. →pages 10Henkel, H., 1991, Petrophysical properties (density and magnetization) of rocksfrom the northern part of the baltic shield: Tectonophysics, 192, 1–19. → pages9, 69147Hestenes, M. R., and E. Stiefel, 1952, Method of Conjugate Gradient for solvinglinear systems: Journal of Research of the National Bureau of Standards, 49,409–436. → pages 29, 41Hotelling, H., 1933, Analysis of a complex of statistical variables into principalcomponents: Journal of Educational Psychology, 24, 417–441. → pages 104Hu, M.-K., 1962, Visual pattern recognition by moment invariants: IRETransactions on Information Theory, 8, 179–187. → pages 111Huber, P. J., 1964, Robust estimation of a location parameter: The Annals ofMathematical Statistics, 35, 73–101. → pages 34Jones, E., T. Oliphant, and P. Peterson, 2001, SciPy: Open source scientific toolsfor Python: → pages 42Kissel, C., and C. Laj, 1989, Paleomagnetic rotation and continental deformation:Kluwer Academic Publishers, volume 254 of Series C: Mathematical andPhysical Sciences. → pages 9Koivisto, E., A. Malehmir, N. Hellqvist, T. Voipio, and C. Wijns, 2015, Building a3D model of lithological contacts and near-mine structures in the Kevitsamining and exploration site, Northern Finland: constraints from 2D and 3Dreflection seismic data: Geophysical Prospecting, 63, 754–773. → pages xii,xxii, xxiii, 6, 7, 117, 118, 119, 120, 127, 128, 136Krahenbuhl, R., and Y. Li, 2015, Robust parametric inversion using adaptivequenched simulated annealing: SEG Technical Program Expanded Abstract,1526–1530. → pages 7Last, B. H., and K. Kubik, 1983, Compact gravity inversion: Geophysics, 48, no.8, 713–721. → pages 8, 34, 46, 55Lawson, C. L., 1961, Contribution to the theory of linear least maximumapproximation: PhD Thesis, University of California, Los Angeles. → pagesxiv, 8, 34, 35Le, C. V. A., B. D. Harris, and A. M. Pethick, 2016, Magnetotelluric inversion,carbonaceous phyllites and an ore zone: Kevitsa, finland: 25th GeophysicalConference and Exhibition, 525–529. → pages 120Lee, T.-C., R. Kashyap, and C.-N. Chu, 1994, Building skeleton models via 3-dmedial surface/axis thinning algorithms: Computer Vision, Graphics, andImage Processing, 56, 462–478. → pages 111Lelie`vre, P. G., 2009, Integrating geological and geophysical data throughadvanced constrained inversion: PhD Thesis, The University of BritishColumbia. → pages xix, 29, 69, 84, 88, 89, 90, 93, 142Lelievre, P. G., and D. W. Oldenburg, 2009, A 3D total magnetization inversionapplicable when significant, complicated remanence is present: Geophysics,74, L21–L30. → pages 10, 12, 70, 71, 74Lelie`vre, P. G., D. W. Oldenburg, and N. Williams, 2009, Integrating geologic and148geophysical data through advanced constrained inversions: Extended Abstracts,20th Geophysical conference and exhibition, ASEG, 1–6. → pages 29Li, Y., 1993, A globally convergent method for lp problems: SIAM J. ScI. Stat.Comput., 3, 609–629. → pages 8, 34Li, Y., A. Melo, C. Martinez, and J. Sun, 2019, Geology differentiation: A newfrontier in quantitative geophysical interpretation in mineral exploration: TheLeading Edge, 38, 60–66. → pages 6Li, Y., and D. W. Oldenburg, 1996, 3-D inversion of magnetic data: Geophysics,61, 394–408. → pages 5, 30, 31, 66——–, 1998, 3-D inversion of gravity data: Geophysics, 63, 109–119. → pages 5——–, 2000, Incorporating geological dip information into geophysicalinversions: Geophysics, 65, 148–157. → pages xix, 7, 12, 29, 84, 85, 87, 90——–, 2003, Fast inversion of large-scale magnetic data using wavelet transformsand a logarithmic barrier method: Geophysical Journal International, 152,251–265. → pages 21Lin, W., and M. Zhdanov, 2017, Joint inversion of seismic and gravitygradiometry data using gramian constraints: SEG Technical Program ExpandedAbstracts, 1734–1738. → pages 7Liu, 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 10, 71Liu, S., X. Hu, H. Zhang, M. Geng, and B. Zuo, 2017, 3D magnetization vectorinversion of magnetic data: improving and comparing methods: Pure andApplied Geophysics, 174, 4421–4444. → pages 74Lockhart, G., H. Grutter, and J. Carlson, 2004, Temporal, geomagnetic and relatedattributes of kimberlite magmatism at Ekati, Northwest Territories, Canada:Lithos, 77, 665–682. → pages 9Macleod, I. N., K. Jones, and T. F. Dai, 1993, 3-D Analytic Signal in theInterpretation of Total Magnetic Field Data at Low Magnetic Latitudes:Exploration Geophysics, 24. → pages 3Marjoribanks, R., 2010, Geological methods in mineral exploration and mining,second ed.: Springer. → pages 1McMillan, M., C. Schwarzbach, E. Haber, and D. Oldenburg, 2016, Multiple bodyparametric inversion of frequency- and time-domain airborne electromagnetics:SEG Technical Program Expanded Abstracts 2016, 846–851. → pages 7Miller, C., S. Kang, D. Fournier, and G. Hill, 2018, Distribution of and condensatein a hydrothermal system: Insights from selfpotential inversion at MountTongariro, New Zealand: Geophysical Research Letters, 45, 8190–8198. →pages viMiller, C. A., G. Williams-Jones, D. Fournier, and J. Witter, 2017, 3D gravity149inversion and thermodynamic modelling reveal properties of shallow silicicmagma reservoir beneath Laguna del Maule, Chile: Earth and PlanetaryScience Letters, 459, 14 – 27. → pages viMiller, H. G., and V. Singh, 1994, Potential field tilta new concept for location ofpotential field sources: Journal of Applied Geophysics, 32, 213–217. → pages3Montonen, M., 2012, Induced and remanent magnetization in two boreholes ofthe Kevitsa intrusion: Master’s thesis, University of Helsinki. → pages xi, 120,122, 124, 127, 128Mushayandebvu, M. F., P. van Driel, A. B. Reid, and J. D. Fairhead, 2001,Magnetic source parameters of two-dimensional structures using extendedEuler deconvolution: Geophysics, 66, 814–823. → pages 3Mutanen, T., 1997, Geology and ore petrology of the Akanvaara and Koitelainenmafic layered intrusions and the Keivitsa-Satovaara layered complex, northernFinland: Geological Survey of Finland, 395. → pages 117, 119Nabighian, M. N., M. E. Ander, V. J. S. Grauch, R. O. Hansen, T. R. LaFehr, Y.Li, W. C. Pearson, J. W. Peirce, J. D. Phillips, and M. E. Ruder, 2005,Historical development of the gravity method in exploration: Geophysics, 70,63ND–89ND. → pages 3Nagy, D., 1966, The gravitational attraction of a right rectangular prism:Geophysics, 31, 361–371. → pages 13NCEI, 2017, Emag2v3: Earth magnetic anomaly grid (2-arc-minute resolution):Online. ( → pages xii, 4Nocedal, J., and S. J. Wright, 1999, Numerical optimatization: Springer Science.→ pages 29, 41Norris, D. K., and R. F. Black, 1961, Application of paleomagnetism to thrustmechanics: Nature, 192, 933–935. → pages 9Osborne, M. R., 1985, Finit algorithms in optimization and data analysis: JohnWiley & Sons Lts. → pages 42Pearson, K., 1901, On lines and plane of closest fit to systems of points in space:Philosophical Magazine, Series 6, 2, 559–572. → pages 104Pedregosa, F., G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M.Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D.Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, 2011, Scikit-learn:Machine Learning in Python: Journal of Machine Learning Research, 12,2825–2830. → pages 105, 106Phillips, J. D., 2003, Can we estimate total magnetization directions fromaeromagnetic data using Helbig’s formulas: Abstracts Week B, IUGG 2003,B261. → pages 10——–, 2010, Estimating structural dip from gravity and magnetic profile data:150SEG Technical Program Expanded Abstracts 2010, 1202–1206. → pages 4Phillips, N. D., 1996, Geophysical inversion in an integrated exploration progra:examples from the San Nicolas deposit: Master’s thesis, The University ofBritish Columbia. → pages 7, 29Pilkington, M., 1997, 3-D magnetic imaging using conjugate gradients:Geophysics, 62, 1132–1142. → pages 21, 66Portniaguine, O., and M. S. Zhdanov, 2002, 3-D magnetic inversion with datacompression and image focusing: Geophysics, 67, 1532–1541. → pages 8, 46Portniaguine, 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 34Pratt, D. A., B. K. McKenzie, and T. S. White, 2014, Remote remanenceestimation (RRE): Exploration Geophysics, 45, 314–323. → pages 10Pueyo, E. L., F. Cifelli, A. J. Sussman, and B. Oliva-Urcia, 2016, Introduction:Palaeomagnetism in fold and Thrust belts: New perspective: GeologicalSociety, London, 425, 1–6. (Special Publications). → pages 10Queitsch, M., M. Schiffler, R. Stolz, C. Rolf, M. Meyer, and N. Kukowski, 2019,Investigation of three-dimensional magnetization of a dolerite intrusion usingairborne full tensor magnetic gradiometry (FTMG) data: Geophysical JournalInternational, 217, 1643–1655. → pages 10Ramon, M. J., E. L. Pueyo, J. L. Briz, A. Pocovi, and J. C. Ciria, 2012, Flexuralunfolding of horizon using paleomagnetic vectors: Journal of StructuralGeology, 35, 28–39. → pages 9Reid, A. B., J. M. Allsop, H. Granser, A. J. Millett, and I. W. Somerton, 1990,Magnetic interpretation in three dimensions using Euler deconvolution:Geophysics, 55, 80–91. → pages 3Reid, J. E., and A. Pfaffling, 2006, Airborne electromagnetic footprint in 1dearths: Geophysics, 71, 63–72. → pages 21Rocklin, M., 2015, Dask: Parallel computation with blocked algorithms and taskscheduling: Proceedings of the 14th Python in science conference, 126–132. →pages 22Roest, W. R., J. Verhoef, and M. Pilkington, 1992, Magnetic interpretation usingthe 3-D analytic signal: Geophysics, 57, 116–125. → pages 3Salem, A., S. Williams, J. D. Fairhead, D. Ravat, and R. Smith, 2007, Tilt-depthmethod : A simple depth estimation method using first-order magneticderivatives: The Leading Edge, 26, 1502. → pages 4Sanchez, M. G., M. M. Allan, C. J. R. Hart, and J. K. Mortensen, 2014, Extractingore-deposit-controlling structures from aeromagnetic, gravimetric, topographic,and regional geologic data in western Yukon and eastern Alaska: Interpretation,2, 75–102. → pages 3151Schermerhorn, W. D., B. Ritzinger, M. Anderson, J. Peacock, J. Witter, J. Glen,A. Steely, C. Forson, P. Stelling, and D. Fournier, 2017, Geophysicalinvestigation of the Mount Baker geothermal play: Presented at the GSAAnnual Meeting in Seattle. → pages viSharma, P. V., 1966, Rapid computation of magnetic anomalies anddemagnetization effects caused by bodies of arbitrary shape: Pure Appl.Geophys., 64, 89–109. → pages 15Shearer, S., 2005, Three-dimensional inversion of magnetic data in the presenceof remanent magnetization: Master’s thesis, Colorado School of Mines,Golden, Colorado. → pages 10Stocco, S., A. Godio, and L. Sambuelli, 2009, Modelling and compact inversionof magnetic data: A Matlab code: Computer and Geosciences, 35, 2111–2118.→ pages 8, 46, 55Sun, J., and Y. Li, 2013, A general framework for joint inversion withpetrophysical information as constraints: SEG Houston Annual Meeting,3093–3097. → pages 7——–, 2014, Adaptive lp inversion for simultaneous recovery of both blocky andsmooth feature in geophysical model: Geophysical Journal International, 197,882–899. → pages 8, 34, 35, 47, 56, 101, 108Thompson, D. T., 1982, EULDPH: A new technique for makingcomputer-assisted depth estimates from magnetic data: Geophysics, 47, 31–37.→ pages 3Thurston, J. B., and R. S. Smith, 1997, Automatic conversion of magnetic data todepth, dip, and susceptibility contrast using the SPI (TM) method: Geophysics,62, 807–813. → pages 4Verduzco, B., D. J. Fairhead, C. M. Green, and C. MacKenzie, 2004, New insightsinto magnetic derivatives for structural mapping: The Leading Edge, 23, 116.→ pages 3Villalain, J. J., A. M. Casas-Sainz, and R. Soto, 2015, Reconstruction of invertedsedimentary basins from syntectonic remagnetizations. a methodologicalproposal: Palaeomagnetism in Fold and Thrust Belts: New Perspectives. →pages 9Vine, F. J., and D. H. Matthews, 1963, Magnetic anomalies over the oceanicridges: Nature, 199, 947–949. → pages 9Vogel, C. R., 2002, Computational Methods for Inverse Problems (Frontiers inApplied Mathematics): Society for Industrial Mathematics. → pages 67Wijns, C., and P. Kowalczyk, 2007, Interactive geophysical inversion usingqualitative geological constraints: Exploration Geophysics, 38, 208–212. →pages 7Williams, N., 2008, Geologically constrained inversion: PhD Thesis, University152of British Columbia, Vancouver, British Columbia. → pages 7, 29Yang, D., and D. W. Oldenburg, 2014, 3D inversion of airborne electromagneticdata parallelized and accelerated by local mesh and adaptive soundings:Geophysical Journal International, 196, 1492–1507. → pages 22Zhdanov, M., and E. Tolstaya, 2004, Minimum support nonlinear parametrizationin the solution of a 3D magnetotelluric inverse problem: Inverse Problems, 20,937–952. → pages 56Zhu, Y., M. S. Zhdanov, and M. Cuma, 2015, Inversion of TMI data for themagnetization vector using Gramian constraints: Presented at the SEG AnnualMeeting. → pages 10, 71153AppendixThis section describes the partial volume calculation defined by the intersection oftwo upright prisms. Let us define the extent of a prism mi by the lower southwestand upper northeast nodal coordinates:mi = [nL nU ] =xL xUyL yUzL zU (A.1)To calculate the volume of intersection of two prisms mi and m j I first calculatethe length Lx intercepted along the x-axis asLx = max[xmax− xmin 0]xmin = min[nxU i nxU j]xmax = max[nxLi nxL j] (A.2)The same calculation is done along the y and z axes (Ly, Lz). The total volume isthan calculated asV = Lx ∗Ly ∗Lz (A.3)Intersecting length calculations in (A.2) can easily augmented for M pairs of cells.154


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items