Cross-Hole GPR Imaging: Traveltimeand Frequency-Domain Full-WaveformInversionbyDaryl Van VorstB.A.Sc., The University of British Columbia, 1998M.A.Sc., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2014c© Daryl Van Vorst 2014AbstractGround-penetrating radar (GPR) has the potential for high-resolution imaging of near-surfacematerial properties, including electrical conductivity and permittivity, which can be used forgeological interpretation of the near subsurface. This thesis presents ray-based traveltime inver-sion and frequency-domain full-waveform inversion (FWI) techniques for application to boreholeGPR surveys.Ray-based traveltime inversion is attractive for its speed, reliability, and ability to work in3D, but the ray approximation involved limits recoverable detail to greater than one wavelength.The traveltime method presented here uses an efficient and easily programmed fast-sweepingeikonal solver to compute traveltimes. The inversion method also incorporates the unknowntime offset between signal transmission and start of recording at the receiver as a model pa-rameter that is recovered simultaneously with the material slowness.The resolution of FWI approaches the diffraction limit of one half wavelength, but at asubstantial computational cost. The FWI inversion scheme presented here works in 2D and isunique in its simultaneous recovery of the source wavelet, conductivity, and permittivity. Itsfrequency-domain formulation allows for efficient factorization of the forward modeling operatorand its subsequent application to multiple right-hand sides in order to quickly construct theforward model Jacobian. Efficient calculation of the Jacobian allows the use of the Gauss-Newton technique rather than the gradient descent method that is common for other GPRFWI inversions.Measured data must be converted from 3D to 2D before use with this 2D FWI technique.I present a graphical derivation of the perpendicular ray Jacobian, which is an essential part of3D to 2D transformation. The graphical derivation provides the reader with an intuitive under-standing of the Jacobian that is difficult to obtain from traditional mathematical treatments.I also illustrate that 3D to 2D transfer functions previously derived for the acoustic case areiiAbstractapplicable to borehole GPR.iiiPrefacePublished WorksThis work has led to several submitted and published papers. Those that have been includedin this thesis are listed below along with the chapter number in which they can be found:• D. Van Vorst, M. Yedlin, P. Jeanne, Y. Guglielmi, F. Cappa, S. Gaffet, “Eikonal 3Dtraveltime inversion with strong robust GCV and its application to borehole GPR data,”submitted, December 2013. Chapter 2.• c©2012 Seismological Society of America. Reprinted with permission from D. Van Vorst,J. Virieux, and M. Yedlin, “A Graphical Derivation of the Perpendicular Ray Jacobian,”B. Seismol. Soc. Am., vol. 102, no. 6, pp. 2452–2457, 2012. Chapter 3.• c©2014 the authors. Reprinted with permission from D. Van Vorst, M. Yedlin, J. Virieux,and E. Krebes, “Three-dimensional to two-dimensional data conversion for electromagneticwave propagation using an acoustic transfer function: application to cross-hole GPRdata,” Geophys. J. Int., vol. 198, pp. 474–483, 2014. Chapter 4.• D. Van Vorst and M. Yedlin, “Full-waveform inversion of ground penetrating radar datain the frequency domain: simultaneous recovery of permittivity, conductivity, and sourcewavelet,” submitted, June 2014. Chapter 5.CollaborationIn the development of the above papers, I have collaborated with my advisor, Dr. MatthewYedlin, and his colleague Dr. Jean Virieux as well as several other researchers in the fieldof inversion and geophysics. For all of the papers, I performed the background research andivPrefacederivations, developed the computer code, ran the computer simulations, interpreted the results,wrote and submitted the papers. The only exception of note to this is that Dr. Jean Virieuxprovided the initial paraxial ray theory derivation of the perpendicular ray Jacobian in the2012 BSSA paper. I then performed the careful modifications that were necessary to make itquantitatively accurate. Dr. Edward Krebes provided two paragraphs in the 2014 GJI paperregarding the treatment of absorption in ray tracing. A few paragraphs were provided by Dr.Pierre Jeanne in the 3D traveltime tomography paper discussing the geological correlation withthe inversion results. Dr. Yedlin and Dr. Virieux (and to a lesser degree the other co-authorslisted above) provided insight and feedback for my work as a whole. Dr. Yedlin provideddirection for my research, helped relate my work to previous work in the field, and helped editthe papers.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction and Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 The Inverse Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.1 Gradient Descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2.2 Gauss-Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.3 Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.3 Ground-Penetrating Radar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.3.1 Electrical Properties of Rock . . . . . . . . . . . . . . . . . . . . . . . . . 141.4 The Ray Approximation and Traveltime Inversion . . . . . . . . . . . . . . . . . 151.5 Full-Waveform Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.6 Review of Related Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191.6.1 Traveltime Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191.6.2 3D to 2D Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 20viTable of Contents1.6.3 Full-Waveform Inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . 211.7 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251.8 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262 Eikonal 3D Traveltime Inversion With Strong Robust GCV and its Applicationto Borehole GPR Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.2 Inversion Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.2.1 T0 Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3 Eikonal Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.4 Jacobian Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.5 Regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.6 Synthetic Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.7 Inversion of GPR Data from LSBB . . . . . . . . . . . . . . . . . . . . . . . . . 462.7.1 Sources of Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.7.2 Measurement Inversion Results . . . . . . . . . . . . . . . . . . . . . . . . 472.8 Geological Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 482.9 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502.10 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513 A Graphical Derivation of the Perpendicular Ray Jacobian . . . . . . . . . . 533.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.2 Graphical Derivation of J⊥ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.2.1 Jacobian Factorization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.2.2 Derivation of J⊥ for a Straight Ray . . . . . . . . . . . . . . . . . . . . . 573.2.3 Extension to Curved Rays . . . . . . . . . . . . . . . . . . . . . . . . . . 593.3 Paraxial Computation of J⊥ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.4 Simplified Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.6 Data and Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.7 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66viiTable of Contents4 Three-Dimensional to Two-Dimensional Data Conversion for ElectromagneticWave Propagation Using an Acoustic Transfer Function: Application tocross-hole GPR data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2 Borehole GPR as an Acoustic Problem . . . . . . . . . . . . . . . . . . . . . . . 694.3 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3.1 Gradient of Permittivity . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3.2 Dielectric Lens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.3.3 Gradient of Permittivity with Loss . . . . . . . . . . . . . . . . . . . . . . 814.4 Full-wave Inversion Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.6 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 875 Full-Waveform Inversion of Ground Penetrating Radar Data in the FrequencyDomain: Simultaneous Recovery of Permittivity, Conductivity, and SourceWavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 885.2 Forward Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905.2.1 FDFD Helmholtz Formulation . . . . . . . . . . . . . . . . . . . . . . . . 915.2.2 Perfectly Matched Layer . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.2.3 Antenna Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.3 Inversion Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.4 Jacobian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1035.5 3D to 2D Conversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1045.6 Synthetic Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065.7 Inversion of Radar Data from LSBB . . . . . . . . . . . . . . . . . . . . . . . . . 1105.7.1 Data Fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.7.2 Source Wavelet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1205.7.3 Water Content . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.8 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123viiiTable of Contents5.9 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1236 Summary and Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . 1256.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131AppendixA Full-Waveform Jacobian Derivation . . . . . . . . . . . . . . . . . . . . . . . . . 144A.1 Jacobian for Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146ixList of Tables1.1 Summary of cross-hole full-waveform GPR inversion research. . . . . . . . . . . . 235.1 Experimental results comparing an antenna inside a borehole to the assumedcurrent model without borehole. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99xList of Figures1.1 Impulse borehole radar system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111.2 Simplified block diagram of a typical impulse GPR transmitter. . . . . . . . . . . 121.3 Simplified block diagram of a typical impulse GPR receiver. . . . . . . . . . . . . 142.1 Ray tracing algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.2 The true model for the synthetic case. . . . . . . . . . . . . . . . . . . . . . . . . 412.3 GCV, RGCV, R1GCV curves for the synthetic experiment. . . . . . . . . . . . . 422.4 L-curve for the synthetic experiment. . . . . . . . . . . . . . . . . . . . . . . . . . 422.5 True model and synthetic results shown as a slice at x = 0.5 m. . . . . . . . . . . 442.6 Misfit vs. regularization parameter for the synthetic experiments. . . . . . . . . . 452.7 Synthetic result from the FDTD data set that best fits the true model, shown asa slice at x = 0.5 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.8 Plane view of the borehole positions. . . . . . . . . . . . . . . . . . . . . . . . . . 462.9 GCV, RGCV, and R1GCV curves for the inversion of measured data from LSBB. 482.10 L-curve for the measured data from LSBB. . . . . . . . . . . . . . . . . . . . . . 482.11 Inversion results for the measured data found with R1GCV shown as a solidvolume and with isosurfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 492.12 Vertical profiles of (a) relative water content, (b) permeability, (c) storativity, (d)porosity, and (e) P-wave velocity measured along the borehole B3. (f) Boreholeviewer image illustrating the lithology contrast along B3. (g) 3D view of thegeological structure [53, 52, 51]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.1 A ray tube cut off at a particular arc length, s. . . . . . . . . . . . . . . . . . . . 55xiList of Figures3.2 Calculation of the angle, dψ0, between the tilted ray segment, with take-off anglesθ0 and dφ, and its in-plane projection with length s0. The angle dφ is assumedsmall. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.3 A ray normally incident on a layered structure. . . . . . . . . . . . . . . . . . . . 583.4 Out-of-plane spreading, hi, along a ray that is initially straight and normallyincident to a layered structure with constant velocities, vi. . . . . . . . . . . . . . 584.1 Schematic representation of the physical problem and action of the transfer func-tion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.2 Computed rays and 2-D fields at 94.8 MHz for a 10 m × 10 m area with varyingslowness: ǫr = 10− 0.8z. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3 Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions (in the medium corresponding to Fig. 4.2), measured at thereceiving antenna as a function of depth (f = 94.8 MHz). . . . . . . . . . . . . . 774.4 Computed rays and 2-D fields at 94.8 MHz for a 10 m × 10 m area with a varyingslowness and small lens: ǫr =[1√10−0.8z − 0.2e−((x−3)2+(z−6)2)]−2. . . . . . . . . . 794.5 Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions (in the medium corresponding to Fig. 4.4), measured at thereceiving antenna as a function of depth (f = 94.8 MHz). . . . . . . . . . . . . . 804.6 Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions for a permittivity gradient with σ = 0.005 S m−1, measuredat the receiving antenna as a function of depth (f = 94.8 MHz). . . . . . . . . . 824.7 True model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 844.8 Recovered permittivity using 2-D ray inversion. . . . . . . . . . . . . . . . . . . . 854.9 Recovered model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 865.1 The computation grid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 935.2 Electric field streamlines computed in COMSOL for a thin cylindrical antennainside of a cylindrical borehole. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.3 Radial electric field lines around an antenna inside a borehole. . . . . . . . . . . 985.4 Source wavelet for the synthetic inversion experiments. . . . . . . . . . . . . . . . 107xiiList of Figures5.5 Synthetic experiment true model and inversion results. . . . . . . . . . . . . . . . 1085.6 L-curves for the (a) 2D and (b) 3D data sets. . . . . . . . . . . . . . . . . . . . . 1095.7 Recovered source wavelet coefficients for the 2D data set. . . . . . . . . . . . . . 1105.8 Recovered source wavelet coefficients for the 3D data set. . . . . . . . . . . . . . 1115.9 Borehole locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.10 Ray traveltime inversion results using the LSBB radar data set. . . . . . . . . . . 1135.11 Full-wave inversion results using the LSBB radar data set. . . . . . . . . . . . . . 1135.12 Full-wave inversion results for opposing hole pairs using the LSBB radar data set.1145.13 Recovered permittivity and conductivity for the 3-2 hole combination. . . . . . . 1165.14 Power content in the received GPR signal at three frequencies with 7.3 MHzbandwidths plotted against the power content at 157.4 MHz. . . . . . . . . . . . 1175.15 L-curves for the full-wave inversions of LSBB data. . . . . . . . . . . . . . . . . . 1185.16 To test for dispersion, the argument of the ratio of two received signals is plottedagainst frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.17 Received and predicted signal amplitudes and phases at 121 MHz. . . . . . . . . 1195.18 Recovered source wavelet coefficients for (a) the 5-2 hole combination and (b)the 5-3 hole combination. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1215.19 Recovered source wavelet coefficients for all hole combinations. . . . . . . . . . . 1225.20 Water content for the 5-2 and 3-4 hole combinations estimated using the Topprelation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122xiiiList of AbbreviationsADC analog-to-digital converterBJT bipolar junction transistorFFT fast Fourier transformFMCW frequency modulated continuous waveFWI full-waveform inversionGCV generalized cross-validationGPR ground penetrating radarLNA low-noise amplifierLSBB Laboratoire Souterrain a` Bas Bruit (Underground Low Noise Laboratory)PRBS pseudo-random binary sequenceRGCV robust generalized cross-validationR1GCV strong robust generalized cross-validationSFCW stepped frequency continuous wavexivAcknowledgementsI am grateful for the incredible level of support, mentoring, and encouragement that I have re-ceived from my advisor Matthew J. Yedlin. Without our numerous discussions this work wouldnot have been possible. I am also grateful to Douglas Oldenburg for teaching me the importantconcepts of inversion and Eldad Haber for providing continued insights and suggestions for myresearch. Jean Virieux has provided a great deal of insight and analytical expertise which Iam extremely lucky and grateful to have received. External reviewer, Niklas Linde, providednumerous helpful suggestions that improved the quality of this thesis.My lab mates, Peter Pawliuk and Jonathan Fraser, have indulged me in countless (andsometimes heated) discussions about many topics, mostly related to our research. These dis-cussions have led to valuable advances in my understanding and research. Thank you Jon andPeter.Funding for this research has come from numerous sources. Most substantially it has beenfrom the Natural Sciences and Engineering Research Council of Canada. Funding has beenprovided by ANR “Captage de CO2” through the HPPP-CO2 project, and by PACA countythrough the PETRO-PRO project. LSBB and the re´gion PACA have provided the infrastruc-ture. The ANR MAXWELL project provided direct funds for travel and per diem.I am thankful for the University of Pau and Pays de l’Adour’s generosity in providing acomplete borehole radar system for my measurements at LSBB. I am also thankful for MALA˚Geoscience’s generosity in lending me a set of 250 MHz borehole radar antennas at very shortnotice in order to complete my measurements when a problem was discovered with the systemI was using.xvChapter 1Introduction and BackgroundInversion is the process of determining the physical makeup of a system from indirect measure-ments. It is the inverse of solving the forward problem. Simulating a system, a task familiarto every engineer, is solving the forward problem. Solving the forward problem provides theanswer to the question: “Assuming I have this system, what will I measure when I perform myexperiment?” An example of this is calculating how a lump of Jell-O with grapes in it wobbleswhen tapped with a spoon (assuming all of the mechanical properties of the Jell-O are known).The inverse problem answers the (usually) far more difficult question: “Given that I have theseresults from my experiment, what is the system that produced them?” This is analogous totapping the Jell-O with a spoon and measuring how it wobbles and from this determining theJell-O’s mechanical properties along with that of the grapes and their locations.In this thesis, I present my research on solving the cross-hole GPR inverse problem. Themeasurement setup for the cross-hole GPR problem consists of a radar transmitter and receiverpair whose positions are moved inside boreholes in the earth while measurements are made.The desired outcome is a map, or image, of the material parameters (electrical permittivity andconductivity) of the material between the boreholes.The research started with 2D and then 3D traveltime inversions that use a recent advance-ment in solving the eikonal equation by Zhao [124], which is efficient and easily implemented.From there, my work progressed toward full-waveform inversion in 2D in the frequency domainusing the Gauss-Newton method. Performing 2D inversion requires adjusting measured dataso that it matches synthetic 2D data (as closely as possible). This led to additional researchinto 3D to 2D data conversion in the uniform ray approximation.I start off here with the motivation for my research, followed by brief backgrounds of theinverse problem, ground-penetrating radar, the ray approximation, and full-waveform inversion.11.1. MotivationInterested readers are encouraged to explore the references provided for a much more in-depthdiscussion. Following this is a summary of research related to this thesis, an overview of mywork that connects it to the papers presented in later chapters and to the motivation, and aconcise summary of my contributions to the field.1.1 MotivationThe primary motivation at the start of this work was simply to improve my own understandingof inversion, forward modeling methods, electromagnetics (and acoustics/seismics), and theassociated mathematical tools. The motivation for tackling the specific issues in inversionaddressed by this thesis are outlined here.Simplicity and Reliability: By simplicity I mean being easy to understand, explain, andimplement [117]. By reliability I mean the probability of a system completing its expectedfunction [116]. Simplicity and reliability often (but not always) go hand in hand. The fewerthings there are to go wrong usually means that fewer things will go wrong. Although notspecifically necessary for application, simplicity also leads to clearer and more complete under-standing and is of key importance to most people starting out in a new area (save for the fewgeniuses among us). Reliability is of great practical importance for inversion, both in the sensethat it is unacceptable if an inversion fails outright and in that the inversion should producean acceptable result. Both the forward modeling operation and the type of inversion algorithmstrongly affect reliability. I have attempted to address both of these issues in all of the publishedand submitted papers that appear in the following chapters of this thesis.Speed: Speed is of critical importance in inversions, particularly where the forward mod-eling operation is of concern. The forward modeling operation is often executed hundreds orthousands (or even more) times through the course of an inversion, and therefore its speed hasa dramatic effect on the time an inversion takes. The type of algorithm used to achieve theminimization of an objective function, and thereby yield the sought-after inversion result, canalso dramatically affect the number of iterations required and consequently the overall time aswell. Although all practitioners of inversion need be relatively patient people, there is a limit.Speed has therefore had considerable influence on the inversion approaches taken in this work.21.1. MotivationThe methodology used in the following research produces inversion results in a timely fashion.In the case of 2D traveltime inversion this can be seconds, but for full-waveform inversion it ishours or days.Resolution: By resolution I mean the ability to resolve detail [115]. Detailed resolutionanalysis for GPR has been presented by Oberro¨hrmann et al. [70]. The clearer a subterraneanimage is, the easier it is to determine what an anomaly is and where it is located. There isgenerally a trade-off between resolution and processing time. Traveltime inversion algorithmscan provide low resolution but extremely fast inversion results. Such results may be adequateon their own, or they may raise questions that higher-resolution approaches, such as full-waveform inversion, can address. Resolution is ultimately limited by the information contentand accuracy of the measurements on which the inversion is based, which are in turn limitedby the governing physics (such as the diffraction limit) and by the number and arrangementof the measurements. Increasing resolution while minimizing the impact on processing time iswhat full-waveform inversion algorithms attempt to address.Precision & Accuracy: Accuracy is the degree of closeness measurements are to the quan-tity’s actual value [113]. Precision is the degree to which measurements of the same quantityunder unchanged conditions show the same results [113]. A limiting factor that directly affectsthe accuracy of an inversion is the degree to which the physical world is accurately representedby the forward modeling operation.Understanding: I have attempted to write papers, and this thesis, in a way that both incre-mentally advances the field of GPR inversion and renders the theory of it easily understandable(within reasonable limits). Academic works should always at least attempt to address the lat-ter, although the word “easily” is perhaps optional. A related issue is implementability, whichis of particular concern in applied science. I have also attempted to provide enough detail thatall of the described methods can be implemented by practitioners in the field without greatdifficulty.31.2. The Inverse Problem1.2 The Inverse ProblemMuch of the material in this section is adapted from Doug Oldenburg’s graduate course on geo-physical inversion [71]. There is a vast repository of literature describing the inverse problem,including [69, 93, 110, 77, 72]. The inverse problem seeks to find a model that suitably fits themeasured data. The problem is generally ill-posed due to non-uniqueness of models that fit thedata, inconsistencies in the data due to noise, forward-modeling approximations and discretiza-tion, and any other errors. Additionally, the forward model in most geophysical problems has asmoothing nature, and so its inverse must necessarily roughen, making it particularly sensitiveto data inconsistencies. In order to achieve reasonable inversion results, some preconceivednotion of what the model should be like must be incorporated into the solution.There are many ways to incorporate some notion of what the model should be like into theinverse problem. Examples include statistical approaches such as those pioneered by Tarantola[95, 94], the Backus-Gilbert method [11, 42], and Tikhonov regularization (see, e.g., [97, 2,44, 36]). Tikhonov regularization is the most commonly used method for regularizing ill-posedinverse problems. It directly penalizes some measure of the model that is undesirable. Unde-sirable qualities of the model that are usually penalized are deviation from an initial guess androughness.The Tikhonov method sets the inversion up as a least-squares minimization of a combinationof data misfit, φd(m), and a model norm, φm(m). A typical objective function that is minimizedisφ(m) = φd(m) + βφm(m)= 12 ‖Wd(F[m]− dobs)‖2 + β 12Nreg∑i=1‖Wi(m−mref )‖2 , (1.1)where m is a vector containing the discretized unknown model parameters (e.g., velocities).The forward modeling operation, F[m], generates the predicted data, dpred:dpred = F[m] (1.2)41.2. The Inverse ProblemThe measured (observed) data are dobs, and Wd weights the data misfit according estimatederror magnitude. The notions about the model are encapsulated in the reference model, mref ,and the matrices Wi that, for example, affect derivative operations if model roughness is to bepenalized. Nreg is the number of regularization measurements (such as three if roughness inboth x and z and the difference between m and mref are to be penalized). The regularizationparameter, β, sets the trade-off between data fitting and model smoothness and/or deviationfrom mref . The 1/2 factor is optional and is only there to simplify derivatives later on.In writing equation (1.1) we are making the assumption that the data misfit, F[m]− dobs,obeys Gaussian statistics [93] (which may not be true). The equation also makes assumptionsabout the structure of the subsurface through the matrices Wi and the reference model mrefthat may not match the true model statistics. Despite the potential difference in statistics, thistype of formulation is frequently used due to its simplicity of implementation.During an inversion, the regularization parameter is often set to a very high value at thestart and then lowered as the inversion progresses. This achieves a ‘cooling’ effect where therecovered model starts out very smooth and, as the regularization parameter is reduced, morestructure ‘crystallizes’ out of the model. Reducing β too far produces a model that is excessivelyrough as the inversion starts to fit the noise. Some appropriate value of β must be chosen so thata good model is found. The simplest method to do this is simply by eye. That is, the inversionis allowed to progress to a very small β and the user views the recovered images at variouslevels of regularization and chooses one that they feel has the most realistic structure. If thestatistics of the noise are known accurately, then a target data misfit, φd = φ∗d, can be sought.If the statistics are Gaussian and uncorrelated, and the entries of the diagonal matrix Wd areequal to the inverse of the standard deviation, then the target misfit is φ∗d = N/2 (where N isthe number of data points). This finds the simplest model, or the model closest to the referencemodel (depending on the choice of Wi matrices) under the constraint that the data are fittedto the presumed noise level. Alternative automated approaches include the l-curve method [45,36] and generalized cross validation (GCV) [44, 36]. The l-curve method is a graphical methodthat involves analyzing a plot of φd vs. φm to find the corner of a characteristic ‘L’ shape.However, sometimes the corner is not clearly defined. GCV is a statistical approach in whichβ is chosen such that the recovered model best approximates missing data. The GCV method51.2. The Inverse Probleminvolves finding the minimum of a potential function that often has a very flat bottom, andfinding the minimum can therefore be difficult. GCV is also sensitive to correlated data [63].Other methods exist that are related to GCV, called robust GCV (RGCV) and strong robustGCV (R1GCV), which have potential functions with more clearly defined minimums and areless sensitive to correlated data [62, 63].The forward modeling operator can be either linear or non-linear depending on the geophys-ical problem and the degree of approximation used. For linear problems, such as straight-raytomography or microgravity inversions [42], it is possible to solve equation 1.1 directly. For non-linear inverse problems, such as traveltime tomography (with curved rays) and full-waveforminversion, equation (1.1) is minimized iteratively.There are many methods available for solving least-squares minimization problems. Forgeophysical inverse problems, there are three related methods that are most commonly used:gradient descent, Gauss-Newton, and full Newton. An excellent review for the seismic case isgiven by Virieux and Operto [110]. Other methods exist, such as singular value decomposi-tion and direct search methods, but these see less use because of difficulty in scaling to largeproblems. Of the three related methods, gradient descent is extremely common because of theease with which the gradient can be calculated. Gauss-Newton is also quite common and offersbetter convergence than gradient descent, but requires more memory and computation time.These three methods are briefly introduced in what follows.1.2.1 Gradient DescentGradient descent is the most intuitively simple approach for solving non-linear inverse problems,and is used extensively in geophysical inversion [54, 35, 33, 59, 67, 15, 57, 58, 81, 61, 125]. Whenapplying gradient descent, the objective function is usually simplified so that it just containsthe data misfit,φ(m) = 12 ‖Wd(F[m]− dobs)‖2 . (1.3)The gradient descent method iteratively minimizes equation (1.3) by perturbing the currentmodel in the negative direction of the gradient a distance that minimizes φ(m). That is, at61.2. The Inverse Problemeach step m is updated withmk+1 = mk − γ∇mφ(mk), (1.4)where k is the iteration number and γ is a scalar parameter, referred to as the step length, thatis found at each iteration. A line search is often used to find γ such that φ(mk+1) is minimized,and requires at least a couple of forward calculations to do so. A more efficient method thatrequires one additional forward modeling step is suggested by Pica et al. [79] and used by Meleset al. [67]. The gradient is often smoothed to prevent excessive structure from entering themodel. This is similar to using a regularized objective function, such as equation (1.1).Convergence is reached, in principle, when the gradient is zero. This will not usually occurin a numerical optimization method, so the decision to stop the iterative process must be madesome other way. Checking that the data misfit and/or the gradient are below acceptable levelsis one option. Stopping the process after some number of iterations is another.Gradient descent is an attractive method of solving inverse problems because of its very lowprocessing and memory requirements. For waveform inversion problems, only one additionalforward modeling operation is required to calculate the gradient [81]. However, gradient descentcan converge very slowly and require many iterations. For some pathological problems this canbe prohibitively expensive. The Gauss-Newton and full Newton methods described below havebetter convergence properties, and are less likely to stagnate in this way.1.2.2 Gauss-NewtonThe Gauss-Newton method finds the minimum of the objective function by first linearizing theforward model, then finding the exact minimum of the linearized objective function at eachstep. The Gauss-Newton method has also received much attention in geophysical inversionpublications (see, e.g., [123, 101, 103, 82] and [40, p. 128]). The linearization of the forwardmodel is achieved by a vector-valued Taylor series expansion about the current model,F[m+ δm] ≈ F[m] + Jδm, (1.5)71.2. The Inverse Problemwhere δm is a small model perturbation and the Jacobian, J, is a matrix that contains thepartial derivatives of each data point, di, with respect to each model parameter, mj :Jij =∂di∂mj. (1.6)The linearized objective function is obtained by inserting equation (1.5) into equation (1.1)to obtainφ(m+ δm) ≈ 12 ‖Wd(F[m] + Jδm− dobs)‖2 + β 12Nreg∑i=1‖Wi(m+ δm−mref )‖2 . (1.7)To find an equation for the model update at each iteration, the gradient of equation (1.7) withrespect to δm is set to zero. This results inδm =JTWdTWdJ+ βNreg∑i=1WiTWi−1 JTWdTWd(dobs − F[m]) + βNreg∑i=1WiTWi(mref −m) .(1.8)The model is then updated at each iteration withmk+1 = mk + γδm, (1.9)where the constant γ is usually chosen such that 0 < γ ≤ 1. Choosing γ less than unity tendsto improve the stability of the procedure but may slow down convergence. Iterations are againstopped using a check of the gradient size, stagnation, or that an excessive number of iterationshave been performed.This method is related to the gradient descent method. The term JTWdTWd(F[m]−dobs)appearing in equation (1.8) is the gradient of equation (1.3). The inverted matrix that appearsbefore it in equation (1.8) both smooths the gradient and transforms it with an approximationof the curvature of the objective function. The transformation that is applied to the gradientby the inverse matrix in equation (1.8) provides a better update direction that performs welland approaches quadratic convergence.The Gauss-Newton method requires significantly more computation time than gradient81.2. The Inverse Problemdescent, and so may not always be feasible for very large inversion problems. Computing theJacobian occupies most of the additional computation time. However, with the steady increasein computing speed and memory size that we have seen thus far, it is becoming possible to applythis method to larger problems. I apply this method to 3D traveltime inversion in Chapter 2and to the 2D frequency-domain full-waveform inversion problem in Chapter 5.1.2.3 NewtonThe Newton method has seen somewhat less application in geophysical inverse problems thanthe gradient descent and Gauss-Newton methods (see, e.g., [82] and [40, pp. 126–128]). TheNewton method uses an expansion of the objective function about the current model up to thesecond order:φ(m+ δm) ≈ φ(m) + [∇mφ(m)]T δm+12δmTHδm, (1.10)where the Hessian matrix, H, is symmetric and contains all of the second partial derivatives ofthe objective function with respect to the model parameters:Hij =∂φ(m)∂mi∂mj. (1.11)Setting the gradient of equation (1.10) with respect to δm to zero results in the updateexpressionδm = −H−1∇mφ(m). (1.12)When applying this to the regularized objective function of equation (1.1), the gradient termin equation (1.12) is∇mφ(m) = JTWdTWd(F[m]− dobs) + βNreg∑i=1WiTWi(m−mref ), (1.13)which appears in the Gauss-Newton method in equation (1.8). The Hessian isH = ∇m(JT )WdTWd(F[m]− dobs) + JTWdTWdJT + βNreg∑i=1WiTWi. (1.14)91.3. Ground-Penetrating RadarThe last two terms of the Hessian also appear in the Gauss-Newton method in equation (1.8).The Gauss-Newton method is then seen to be the Newton method but with the first term ofequation (1.14) neglected. The first term of equation (1.14) can be difficult to compute andit depends on the data misfit, which should shrink as we approach a solution. So this term isoften neglected, and the resulting penalty will likely be slower convergence early in the iterativeprocess.1.3 Ground-Penetrating RadarThere are a wide variety of ground-penetrating radar systems in use (see, e.g., [28]), withapplications in geology, civil engineering, archaeology, unexploded ordinance detection, andexploration. Perhaps the most common GPR units run along the ground surface and haveeither a single antenna or essentially co-located receiving and transmitting antennas. Theysend signals down into the ground and record the signals reflected off of material discontinuities.Next most common, and what I focus on in this thesis, are borehole radar units. Borehole radarunits usually have separate receive and transmit antennas so that the transmitter and receivercan be placed in different holes, thus allowing signal transmission through the intervening rockmedium to be recorded. Borehole radar units can also be configured so that both antennas arein the same hole and then operate in a similar fashion to the down-looking surface radar, butwith the view extending out radially from the borehole.The 250 MHz MALA˚ Geoscience borehole radar system that was used to make the mea-surements for this research consists of two thin battery-powered radar antennas controlled viafiber-optic cables from a control unit. The antennas are suspended in the boreholes by thefiber-optic cables, which run over ticker wheels that measure the antennas’ depths. The tickerwheel associated with the receiving antenna is also connected to the control unit and triggersthe sequential acquisition of a trace at fixed intervals of travel down the hole. A schematic ofthe particular system that is used for this research is shown in Fig. 1.1. Flexible plastic central-izers are taped to the antennas to keep them centered within the boreholes. This reduces errorsthat could be introduced into the recorded waveform due to lateral movement of the antennaswithin the borehole (14 cm diameter for this survey).101.3. Ground-Penetrating RadarRockBoreholesRX AntennaTX AntennaControlUnitTickerWheelsPlasticCentralizerFiber-Optic CableFigure 1.1: Impulse borehole radar system.111.3. Ground-Penetrating RadarThere are several GPR modulation techniques, with the most common being impulseradar. Other methods include stepped frequency continuous wave (SFCW) (see, e.g., [90, 74]),frequency modulated continuous wave (FMCW), pseudo-random binary sequence (PRBS), andsingle frequency systems. A thorough summary of these techniques and a few others is givenby Daniels [28, pp. 185–246]. The design of impulse radar systems is more straightforwardthan other systems, but due to the limitations of analog-to-digital converters (ADCs) (thathave moderate power consumption and cost) the signal must be re-transmitted many times inorder for it to be sampled completely. This means that a significant portion of the transmittedsignal power is lost. The other methods largely mitigate this problem with more complex down-conversion receiver architectures and low-frequency ADCs, but they often need high-speed gateson both the transmitter and receiver in order to prevent receiver saturation from direct couplingfrom the transmitter.In impulse radar, the most common method of generating impulses is by the avalanchebreakdown of a bipolar junction transistor (BJT) [28, pp. 201–202]. The avalanche breakdownof BJTs typically occurs in the range of hundreds of picoseconds to one nanosecond [17, 12].The breakdown of the transistor delivers the charge stored on a small capacitor to the antenna.The capacitor both limits the duration of the applied pulse and limits heating of the transistorto prevent damage. The differentiating action of the antenna converts the rapidly rising edgeof this pulse into an impulse. A simplified block diagram of a GPR transmitter is shown inFig. 1.2.AntennaTrigger fromControl UnitMatchingand BalunHV SupplyFigure 1.2: Simplified block diagram of a typical impulse GPR transmitter. A trigger pulsefrom the control unit triggers avalanche breakdown in the BJT, which delivers the charge storedon the capacitor to the antenna.121.3. Ground-Penetrating RadarThe differentiating action of a transmitting antenna is apparent upon examination of thegoverning equation for the far electric field of a short dipole antenna [13, pp. 280,318],Eθ ≈jωηIl sin(θ)c4πr e−jkr, (1.15)where Eθ is the electric field (in the θˆ spherical coordinate direction), I is the current, l is thelength of the current element, η is the impedance of free space, c is the speed of light, andk is the phase constant. The jω in the numerator indicates a time derivative is taking place.Feynman [39, p. 28-2] gives an expression for the electric field due to a point charge, whichclearly illustrates the differentiating action of a transmitting antenna,E = −q4πǫ0[er′r′2 +r′cddt(er′r′2)+ 1c2d2dt2er′]. (1.16)The last term in equation (1.16) indicates that accelerating charge is responsible for the farelectric field. As current is a steady movement of charge, accelerated charge is the change ofcurrent with time and this produces the radiated electric field. The pulse produced by the BJTis a voltage pulse, but if one assumes that the antenna terminals are resistive within its usefulbandwidth, then the antenna current is proportional to the driving voltage. In equation (1.16),er′ is the unit vector that points from the position where the field is measured to the charge.The prime indicates retarded quantities.An impulse radar receiver is essentially just an antenna, low-noise amplifier (LNA), andADC. The primary design challenge with such a receiver is the ADC, as it must sample at a highrate. For example, the ADC in the 250 MHz units that were used to make the measurements forthis research ran at approximately 2.5 Gsps. Recent advances in ADC technology are makingthis less of an issue, but typically the receiver structure includes a fast sample-and-hold circuitbefore the ADC. This allows the use of a low-cost and low-power ADC that operates at a ratesubstantially below the effective sample rate. A downside to this approach is that the radarpulse must be transmitted many times in order to build up a complete trace. A simplified blockdiagram of such a radar receiver is shown in Fig. 1.3.131.3. Ground-Penetrating RadarLNA ADCSample andHoldAntennaControlCircuitsTrigger fromControl UnitData toControl UnitFigure 1.3: Simplified block diagram of a typical impulse GPR receiver. A trigger from thecontrol unit initiates the sampling process.1.3.1 Electrical Properties of RockGPR signals are affected by the electrical properties of the rock through which they travel, andtherefore it is these properties are most directly determined during an inversion. It is of interestto find some relationship between the electrical properties and mechanical or other propertiesof rock to aid in identifying it based on inversion results. A complete analysis of the electricalproperties of rock is complicated and beyond the scope of this research. However, I review heresome of the common relationships that are used to relate electrical properties to the physicalproperties of rock.Archie’s law [114] is a low-frequency relationship between the effective electrical conductivityof a brine saturated rock, σef to the conductivity of the brine, σw, the porosity of the rock, φ,the cementation index, m, the tortuosity of the rock, a, and to the saturation (expressed bySnw):σef =1aσwφmSnw (1.17)For higher frequencies the complex refractive index method (CRIM) should produce betterresults, and is given by [32],(ǫef )1/2 =N∑i=1ci(ǫi)1/2. (1.18)In equation (1.18) the effective complex permittivity is ǫef , N is the number of components inthe composite, ci is the volume fraction of the ith component, and ǫi is the complex permittivity141.4. The Ray Approximation and Traveltime Inversionof the ith component.Another expression of interest is the Topp relation [99]. The Topp relation is an empiricallyderived expression for materials that range from sandy loams to clay. It relates water contentto the effective relative electric permittivity of the material, and is given by:ǫr,ef = 3.03 + 9.3θv + 146.0θ2v − 76.7θ3v (1.19)In equation (1.19) ǫr,ef is the effective relative permittivity and θv is the volumetric watercontent.1.4 The Ray Approximation and Traveltime InversionThe ray approximation can be used to simplify either vector (electromagnetic or seismic) waveequations, or the scalar (acoustic) wave equation. For simplicity, I introduce the ray approx-imation using the scalar wave equation, but the results are also applicable to vector waveequations. The ray approximation reduces the computational complexity of solving the wavescalar equation,∇2P (x) + w2v2(x)P (x) = δ(x), (1.20)by using the Wentzel-Kramers-Brillouin (WKB) approximation [98, 100], where the solutionis assumed to be in the form of a series in powers of ω−n (with n = 0 . . .∞) [121]. In equa-tion (1.20), ω = 2πf is the angular frequency, v is velocity, P is pressure, and δ(x) is the Diracdelta function. The approximation is valid as long as∣∣∣∣δvv∣∣∣∣≪ 2πλmin, (1.21)where λmin = 2πvmin/ω and vmin is the lowest velocity. Equation (1.21) states that the materialproperties must vary slowly in space relative to the wavelength. This fundamentally limits theresolution that can be accurately obtained when using the ray approximation to generate imagesusing inversion.151.4. The Ray Approximation and Traveltime InversionThe classic asymptotic series [56] that is assumed for the solution of equation (1.20) isP (x) = e−iωT (x)∞∑n=0An(x)(iω)n , (1.22)where T (x) is the traveltime function and An are amplitude coefficients. A drawback of usingequation (1.22) is that the source singularity is not captured within it, which leads to problemscalculating amplitudes in the vicinity of the source. A modification to this [121, 8, 7, 10, 9]that directly includes the source singularity isP (x) = − e−iωT (x)4πv(x)T (x)∞∑n=0An(x)(iω)n . (1.23)Substituting either of equations (1.22) or (1.23) into equation (1.20) leads to a hierarchy ofequations. Setting the coefficient of ω2 to zero results in the eikonal equation,|∇T (x)|2 = 1v2(x) = s2(x), (1.24)where s(x) is slowness. The remaining equations of the hierarchy are referred to as transportequations and allow the determination of the coefficients An. Usually only the first transportequation is used and allows calculation of the amplitude along the ray. The higher termsattempt to take into account diffraction.The eikonal equation (equation (1.24)) forms the basis of ray tracing and from it boththe traveltime from source to receiver and the ray path itself can be determined. (See, e.g.,[64, pp. 82–87] and [22, pp. 99–109].) Traveltime inversion matches the solution of the eikonalequation to measured traveltimes in order to estimate the velocity profile of the material betweentransmitter and receiver locations. In this way, the forward modeling operation in equation (1.1)is the eikonal equation, and the observed data are traveltimes. Due to the limit set out inequation (1.21), ray-based traveltime inversions are only able to resolve anomalies that arelarge in comparison to the wavelength in the background medium. Further details on traveltimeinversion are provided in Chapter 2.Ray theory can also be used to recover conductivity (or whatever dissipative process mayexist) by using ray amplitudes that are computed from the transport equation (the equation161.5. Full-Waveform Inversionthat governs A0) (see, e.g., [49] and [22, chap. 5]). However, as Holliger et al. [49] point out,uncertainty in the radiation pattern of the antennas leads to difficulty in obtaining accurateamplitude inversion results.Analysis of the behaviour of the ray amplitude for 2D and 3D cases leads to a simpleand effective approximate method for translating 3D radar data into equivalent 2D data [121,100, 104]. This allows 2D full-waveform inversion methods to be applied to measured data.When simulating the electromagnetic fields from a radar antenna in two dimensions, the sourceextends to infinity in the third dimension. This fundamentally changes the character of thetransmitted signal. Rather than decaying as 1/r in 3D, the 2D signal decays as 1/√r. Further,the transmitted field from an impulsive source has a tail in 2D that it does not in 3D. The tailexists as a result of progressively weaker and delayed signals arriving from the far extents ofthe source, which is stretched out to infinity. The decay of the far-field signal as a function ofdistance from the source can be handled using ray theory, and the tail can be accounted for byapplying a π/4 phase change and a weighting proportional to 1/√ω to the data. Details of thetransformation are discussed in [121] and in Chapters 3 and 4.1.5 Full-Waveform InversionFull-waveform inversion attempts to recover a more precise estimate of the unknown model bydirectly simulating the measured waveform rather than some aspect of it, such as traveltime oramplitude. Although a full waveform can be simulated in the ray approximation, full-waveforminversion usually implies solving the wave equation (Maxwell’s equations in the case of GPR)without the ray approximation during the forward modeling step. This extends the resolutionof the inversion beyond the requirement of a slowly varying model, as set out in equation 1.21.In principle, full-waveform inversion should be capable of resolving detail up to the diffractionlimit [83, 78]d = λ2 , (1.25)where d is the size of the smallest observable detail, and λ is the wavelength of the signal inthe background medium.The solution of Maxwell’s equations requires significantly more memory and processing time171.5. Full-Waveform Inversionthan computing traveltimes by ray tracing or by solving the eikonal equation directly. As theforward modeling step must be performed many times during the inversion, this can lead to adramatic increase in the processing time required for completion. This is particularly true of theGauss-Newton and full Newton methods where computing the Jacobian and Hessian requiresmany forward-modeling operations. As a result, full-waveform inversion is often limited togradient descent.An estimate of the transmitted signal (the source wavelet) must be made during the inver-sion, and is a key part of producing a realistic predicted waveform. In traveltime inversion, allthat must be estimated is the point in time when the pulse is transmitted. In full-waveforminversion, this is estimated in conjunction with the wave shape and amplitude. The sourcewavelet (or the time of transmission) represents more unknowns that must be determined dur-ing the inversion, but their effect on the data is different from that of the material properties.As a result, these unknowns are sometimes found at the start of the inversion and are keptunchanged as the recovered model evolves. The methods demonstrated in Chapters 2 and 5treat these unknowns in a similar way to the permittivity (or velocity) and conductivity modelparameters and update them simultaneously as the inversion proceeds.The antenna radiation pattern plays an important role in full-waveform inversion, partic-ularly when finding conductivity (or any other absorption mechanism). With the assumptionthat material loss is relatively low so that the rock behaves as a good dielectric for wave prop-agation, the phase of the received signal strongly affects the recovered permittivity and theamplitude of the signal strongly affects recovered conductivity. Complete numerical modelingof an antenna system is a difficult task and requires a very fine numerical grid in the vicinity ofthe antenna due to the size scale of the antenna structure. This can mean a significant increasein problem size (memory requirements and processing time), so some approximation must bemade when performing inversions. The proprietary nature of commercial GPR systems alsopresents an impediment to a complete simulation of the antennas. A simple, and not entirelyunrealistic, approach is to assume that the antennas are infinitesimal dipoles. In Chapter 5, Ipresent a simple approach for achieving an approximation of real antenna behaviour for boththe receive and transmit antennas.As the above discussion shows, the full-waveform inversion problem is not only more compu-181.6. Review of Related Researchtationally expensive than the traveltime problem, the required detail in modeling the physicalsystem is also higher. In order to help ensure that a full-waveform inversion finds a reasonablemodel, traveltime inversion is often used to produce a starting model. This allows good esti-mates of the source wavelets to be generated at the start of the inversion and greatly improvesthe likelihood of avoiding a local minima during the inversion. A good starting model, especiallyfor FWI problems, is therefore a key requirement of achieving an acceptable inversion result.1.6 Review of Related ResearchIn this section, I provide a review of research that is related to the work presented in Chapters 2through 5.1.6.1 Traveltime InversionTraveltime inversion has received a great deal of attention because of its relatively low processingrequirements. In a traveltime inversion scheme the forward problem is computing the traveltime,and must be made as computationally efficient as possible to produce a timely inversion result.Ray tracing has been used extensively for this purpose in seismics (see, e.g., [109, 22, 123]). Ithas the advantage over other methods that multiple arrivals can be found by following differentray paths. A chief difficulty with this method is finding the ray paths. This generally involvesiteratively perturbing an initial guess of a ray take-off angle and computing new rays until oneis found that lands on the receiver (see, e.g., Virieux and Farra [109]).An alternative approach avoids the use of rays to find traveltimes and instead solves theeikonal equation directly. One method for this is presented by Vidale [107, 108] (for anapplication see, e.g., Aldridge and Oldenburg [2]). Vidale’s method has been updated byHole and Zelt [47] to handle sharp velocity contrasts. The algorithm involves computing thetraveltime on rings in 2D (or surfaces in 3D) at increasing distance from the source point andrequires a sorting operation. Vidale’s scheme is referred to as a fast-marching approach [86]. Asimple and efficient approach that also handles sharp velocity contrasts is given by Zhao [124],and used here in Chapter 2. Zhao’s method is a fast-sweeping approach in which the entire gridis swept in alternate directions and updates are computed at each node based on surrounding191.6. Review of Related Researchtraveltimes. The alternating sweep directions ensure that rays traveling in any direction arecaptured and accurate first-arrival times are computed.An advantage of the fast-sweeping and fast-marching methods is that they are guaranteedto find a traveltime from the source to a receiver at any location. This contrasts with theray method where something in the model structure, such as a lens, could prevent a ray fromreaching a particular receiver even though there is a measured traveltime at that receiver. Sucha breakdown of the forward modeler poses a problem for the inversion algorithm because theobjective function can no longer be evaluated without modification.The model update during the inversion can be computed using any of the methods outlinedin section 1.2. Recently, an adjoint method for calculating the gradient that is related toZhao’s fast-sweeping algorithm was presented by Leung and Qian [61]. Alternatively, whenthe fast-marching or fast-sweeping methods are used, an approximate sparse Jacobian can begenerated using ray paths [2, 101]. The sparse nature of the Jacobian allows quick calculation ofa Gauss-Newton update direction that converges much faster than the gradient descent method.1.6.2 3D to 2D TransformationPerforming realistic full-waveform inversions in 3D is a task that is still not achievable onmodern desktop machines due to memory and processing time requirements, and is somethingthat must be relegated to large compute clusters or supercomputers. To mitigate this problem,full-waveform inversions are often performed in two dimensions. Simulating the world in twodimensions is equivalent to simulating a three-dimensional world in which nothing changesin the direction of the missing dimension. In other words, two things happen: point sourcesbecome infinite line sources and an assumption that material properties do not vary along themissing dimension is imposed. The former of these two things can be addressed mathematicallyby transforming measured 3D data into equivalent 2D data. In the following, this 3D to 2Dtransformation is referred to as a transfer function. The deviation from reality that is imposedby the second item is unfortunate, but it is generally simply accepted.Excellent reviews of the 3D to 2D conversion problem are given by Auer et al. [6] and byWilliamson and Pratt [118]. A transfer function based on acoustic diffraction theory is presentedby Deregowski and Brown [30]. A couple of years later, Vidale et al. [106] present a transfer201.6. Review of Related Researchfunction derived in the time domain for the seismic case. Later, Bleistein [19] formulates atransfer function using source matching methods for the acoustic case. A somewhat simplertransfer function for GPR that uses a uniform material approximation is employed by Ernstet al. [33] . Yedlin et al. [121] directly take the ratio of frequency domain asymptotic Green’sfunctions for the acoustic case and arrive at a transfer function that incorporates both near-and far-field wave behaviour. Their result is valid for inhomogeneous media and is equivalent toBleistein’s result in the far field. In Chapter 4, I show that this transfer function is applicableto the borehole GPR case.The methods of Yedlin and Bleistein require calculation of the perpendicular portion of theray Jacobian (i.e., the amount a planar ray would move out of the plane if its takeoff angle outof the plane were perturbed from zero). Cˇerveny´ [22] presents a thorough derivation of the rayJacobian, but it is complex and quite difficult to follow. His treatment encompasses much morethan is necessary to address 3D to 2D transformation. I present two straightforward derivationsof the perpendicular portion of the ray Jacobian in Chapter 3.1.6.3 Full-Waveform InversionFull-waveform inversion began in the 1980s with Tarantola’s application of an acoustic methodto seismic data [94, 95]. Since then, desktop computing power has increased dramatically witha nearly four-hundred-thousand-fold increase in memory and a twenty-thousand-fold increasein memory bandwidth (comparing my Apple IIe to a Core i7 machine similar to what is cur-rently on my desk). This has opened up possibilities for the solution of full-waveform inversionproblems that were beyond reach in the 1980s.Research into full-waveform inversion of GPR data has received a lot of interest in thepast few years. This is likely due in large part to the continued increase in performance ofdesktop computers, which can now be used to perform realistic 2D full-waveform inversions.Additionally, once boreholes are drilled, cross-hole GPR measurements are relatively easy toundertake.Distributed computing clusters are attractive and reasonably widely available and havethe potential to extend FWI inversion beyond 2D or allow the solution of larger problems.However, they are expensive and require more sophisticated (distributed) programming than211.6. Review of Related Researchdesktop machines do. Therefore, the ability to perform FWI on desktop machines increasesthe accessibility in environments where compute resources might be limited, such as in anengineering office.Table 1.1 summarizes research into GPR full-waveform inversion, and is ordered approxi-mately chronologically. For clarity, I begin with the time domain approaches and summarizethe frequency domain approaches next. The first application of a full-waveform method tocross-hole radar measurements appears to be quite recent, by Ernst et al. [33] in 2007. Theyleverage previous work on finite-difference time-domain modeling [34, 48] that accurately mod-els a transmitting antenna situated in a borehole. However, the receiving antenna is modeledas an infinitesimal dipole. The inversion uses a gradient descent method where permittivityand conductivity are updated one after the other in a sequential fashion. The source waveletis determined at the start of the inversion, with the possibility of update during the inver-sion if required. Prior to the work of Ernst, in 2002 Jia et al. [54] present a gradient descenttime-domain inversion method and apply it to synthetic data. Later, Kuroda et al. [59] (2007)present another time-domain gradient descent scheme that recovers permittivity and apply itto synthetic data. A few years later, Meles et al. [67] (2010) present a time-domain gradientdescent method for updating permittivity and conductivity in a quasi-simultaneous fashion.Meles et al. calculate the gradient through the electromagnetic vector wave equation ratherthan a scalar approximation such as that used by Ernst et al. Belina et al. [15, 16] (2012)explore a deconvolution approach for finding the source wavelet and apply the inversion schemeof Ernst et al. to synthetic and measured data. Belina et al. also study the impact of sourcewavelet variability on an inversion. Cordua et al. [27] (2012) present a time-domain MonteCarlo approach to GPR inversion. Klotzsche et al. [58, 57] (2010, 2013) apply the method ofMeles et al. to field data.In 2011 Ellefsen et al. [31] apply a gradient descent frequency domain method to laboratorydata. In their approach phases and amplitudes are inverted separately for permittivity andconductivity. In 2012 Yang et al. [119] present a frequency-domain Gauss-Newton method forthe simultaneous recovery of permittivity and conductivity, and apply it to a synthetic data set.The source wavelet does not appear to be recovered during the inversion. Lastly, in 2014 Lavoue´et al. [60] apply a frequency-domain quasi-Newton method to synthetic multi-offset surface GPR221.6. Review of Related Researchdata. Their method simultaneously recovers permittivity and conductivity (the source waveletis assumed known). They iteratively approximate the action of the inverse Hessian operatorto achieve a Newton-like descent method. The polarization of the 2D electromagnetic problemthat is solved for surface GPR in their work is the opposite of what is required for boreholeGPR.Table 1.1: Summary of cross-hole full-waveform GPR inversion research.Author Domain Quantities Type† YearJia et al. [54]a Time Quasi-Simultaneous:bPermittivity, ConductivityGD 2002Kuroda et al. [59]a Time Permittivity, Conductivityc GD 2007Ernst et al. [33, 35] Time Sequential:Permittivity, ConductivityInitial:Source WaveleteGD 2007Meles et al. [67]a Time Quasi-Simultaneous:bPermittivity, ConductivityGD 2010Ellefsen et al. [31] Freq.f Sequential:Permittivity, ConductivityGD 2011Cordua et al. [27]a Time Permittivity MC 2012Yang et al. [119]a Freq. Simultaneous:Permittivity, ConductivityGN 2012Belina et al. [15, 16] Time Quasi-Simultaneous:bPermittivity, Conductivity,Source Waveletd,eGD 2012Klotzsche et al. [58, 57] Time Quasi-Simultaneous:bPermittivity, ConductivityInitial:Source WaveleteGD 2010–2013Lavoue´ et al. [60]a,g Freq. Simultaneous:Permittivity, ConductivityQN 2014Van Vorst et al. [102] Freq. Simultaneous:Permittivity, Conductivity,Source WaveletGN 2014† Inversion types: Gradient descent (GD), Gauss-Newton (GN), Quasi-Newton (QN), and Monte Carlo (MC).a Simulation only.b Quasi-simultaneous means that the quantities are updated at each iteration, but in an independent fashion.c Conductivity was not recovered in the numerical experiment.d Discusses both initially determining the source wavelet and including it as an unknown in the inversion.e When finding the source wavelet as an initial step, the outlined inversion algorithm provides a mechanismfor refinement of the source wavelet during the inversion process if necessary.f Phase and amplitude are inverted separately.g Surface radar with opposite polarization to cross-hole.Aside from the inversion formulation, the solution of the forward problem (Maxwell’s equa-231.6. Review of Related Researchtions) has received extensive research attention. Unlike the solution of the eikonal equationthrough ray tracing or one of the sweeping methods, where points outside of the computationdomain are simply ignored, the boundary of the computational domain when solving Maxwell’sequations requires special treatment. In some engineering problems the domain boundary canbe treated as a perfect conductor, and simply setting the fields to zero outside the domain isall that is required. For unbounded problems, where waves are expected to radiate outwardfrom the domain and not return, some form of absorbing boundary condition must be used.Problems like this in the frequency domain tend to result in a poorly conditioned coefficientmatrix that must be inverted [87]. The poor conditioning of the matrix leads to difficulty withiterative solutions.For GPR inverse problems, time-domain methods are, so far, the norm (see, e.g., [33, 35,67]). The finite-difference time-domain (FDTD) approach (Taflove and Hagness [91] have anexcellent book on the topic) avoids the matrix inversion problem by solving the discretized setof Maxwell’s equations directly for an explicit expression for the fields at each time step. Thedownside to this approach is that the time steps must be short and many of them are requiredto propagate the transmitted wavelet through the entire domain. However, once completed,the results can be used to find the impulse response of the system for a very broad range offrequencies. This long simulation time leads to the use of gradient descent methods insteadof the Newton or Gauss-Newton methods because of the large number of forward modelingoperations that would be required.For 2D problems, a finite-difference frequency-domain (FDFD) method such as is describedby Rappaport et al. [85] leads to a coefficient matrix that is small enough to be factored onmodern desktop machines. The factorization can then be used to find the solution to theforward problem very efficiently for a large number of sources. This allows a Jacobian for theforward problem at several frequencies to be built up in an acceptable amount of time so thatthe Gauss-Newton method can be used. This is the approach taken in Chapter 5.241.7. Overview1.7 OverviewChapters 2 through 5 are reproductions of my published or submitted works. They representmy contributions to the field of GPR inversion.In Chapter 2, I present my work on 3D ray tomography. It summarizes and extends mypreviously published work in ray tomography [101, 103]. It differs from previous work by otherauthors through the use of Zhao’s fast-sweeping eikonal solver and the incorporation of findingthe time offset between transmission and start of reception (signified by T0). There is littlediscussion in the literature about incorporating this quantity as an unknown in the inversion.The paper also extends two methods that are related to generalized cross validation (GCV)for finding the regularization parameter for linear problems to nonlinear inversion. These arereferred to as robust GCV (RGCV) [62] and strong robust GCV (R1GCV) [63]. They allowthe GCV method to be applied in situations where data errors are correlated and also producea much sharper minimum that is easier to pinpoint than the standard GCV function.In Chapter 3, I present both graphical and analytical derivations of the perpendicularportion of the ray Jacobian. These are far easier to follow than the treatment provided byCˇerveny´ [22] and provide a degree of insight into the nature of the ray Jacobian.In Chapter 4, I show that the 3D to 2D frequency domain data transfer function derived byYedlin et al. [121] for the acoustic case is directly applicable to the borehole GPR problem. Thisis accomplished analytically by showing that under the WKB approximation of equation (1.21)the ratio of the solutions of Maxwell’s equations in 2D and 3D is the same in the far field asthe acoustic problem when the same assumptions about the model are applied.In Chapter 5, I present a frequency-domain Gauss-Newton inversion formulation for cross-hole GPR. This appears to be the first frequency-domain Gauss-Newton method for cross-holeGPR data to simultaneously recover permittivity, conductivity, and the source wavelet. Iapply the method to both synthetic and field data. I also present a reciprocity-based antennamodel for both receiver and transmitter that is easy to implement and more accurate than theinfinitesimal dipole antenna model that is typically used for FWI.At the start of this chapter, I give a list of motivations for this research. I feel that I haveaddressed them as follows:251.8. ContributionsSimplicity: Simplicity has been addressed in a number of ways. The ray tomography problemis simplified through the use of the fast-sweeping eikonal solver by Zhao along with the directinclusion of T0 in the formulation. The graphical derivation of the perpendicular ray Jacobianis simple and direct. For the full-waveform inversion problem, the matrix formulation reducesthe complexity of much of the math associated with computing the Jacobian.Reliability: Zhao’s fast-sweeping solver improves the reliability of finding traveltimes com-pared to ray tracing. It will always find a traveltime even when the ray method might breakdown. Both inversion formulations use the more robust and faster converging Gauss-Newtonmethod instead of gradient descent.Speed: The fast-sweeping eikonal solver is indeed fast. The overall Gauss-Newton formu-lation of the ray tomography inversion with a sparse Jacobian is extremely fast, taking mereminutes to render an inversion result. The LU or LDLT decomposition of the forward modeloperator for the full-waveform inversion problem allows efficient calculation of the Jacobian.This, in turn, provides convergence in fewer iterations using the Gauss-Newton method insteadof gradient descent.Resolution: The full-waveform inversion formulation approaches the diffraction limit inachievable resolution.Accuracy & Precision: The more complete physical model that is simulated for the full-waveform inversion increases the overall accuracy of the inversion results. The use of a reciprocity-based antenna model for both receiver and transmitter also improves the accuracy of the re-covered conductivity profile.Understanding: The arguments above for simplicity all also apply to understanding. I haveattempted to present all of my research in a way that should be easily understood by myexpected readers.1.8 ContributionsThe contributions that I have discussed throughout the preceding text are summarized here:• Developing 2D and 3D traveltime Gauss-Newton inversion algorithms that use a fast-sweeping solver and incorporate finding the transmitter time offset, T0.261.8. Contributions• Extending robust GCV and strong robust GCV to nonlinear inversion.• Providing both a graphical and analytical derivation of the perpendicular ray Jacobian.The basis of the analytical derivation was provided by Jean Virieux. I further refined itto be quantitatively correct.• Presenting an analysis that shows that the acoustic 3D to 2D transfer function is validfor cross-hole GPR.• Formulating a frequency-domain full-wave inversion that simultaneously finds permittiv-ity, conductivity, and the source wavelet using the Gauss-Newton method.27Chapter 2Eikonal 3D Traveltime InversionWith Strong Robust GCV and itsApplication to Borehole GPR Data2.1 IntroductionIn this work we present a 3D traveltime inversion algorithm for cross-well measurements that isan extension to our preliminary 2D work [101] and some very preliminary 3D work [103]. Thealgorithm uses a robust and easily implemented fast-sweeping solver to compute traveltimesbased on a model estimate [124], an alternative to Vidale’s method [107, 108] and the methodof Podvin and Lecomte [80]. The inversion is formulated as a Gauss-Newton least-squaresminimization where the Jacobian of the forward model is estimated using rays [2, 107, 80]. Theapproximate Jacobian is very quickly computed and sparse. Its use thereby greatly reducesthe memory and processing requirements of the inversion when compared to the use of a fullJacobian [see, for example, 69, p. 128]. The inversion algorithm also solves for the potentiallyunknown time offset, referred to as T0, between signal transmission and start of sampling atthe receiver.The inversion algorithm is applied to a synthetic model and to borehole GPR measure-ments taken at the Vaucluse karst aquifer at LSBB1. A discussion of the geology of the area isincluded and provides ground-truth for the inversion results. The regularization parameter ischosen using a modified form of generalized cross-validation (GCV) called strong robust gener-alized cross-validation, R1GCV, [63] that has been extended, for the first time (to the authors’1(http://lsbb.oca.eu)282.2. Inversion Formulationknowledge), for use in this non-linear inversion context.In what follows the inversion formulation is presented, followed by the 3D eikonal solver, Ja-cobian estimation, regularization and GCV, synthetic and experimental results, and a geologicaldiscussion.2.2 Inversion FormulationThe inversion is formulated in a typical way using Tikhonov regularization [see e.g., 2, 44, 36],where the following objective function is minimized:φ(m, β) = φd(m) + βφm(m)=∥∥∥Wd(F[m]− dobs)∥∥∥2+ βNreg∑i=1αi ‖Wi(m−mref )‖2 (2.1)The objective function, φ, consists of a data norm, φd, that measures misfit between measuredand predicted data, and a model norm, φm, that, in this case, measures smoothness of themodel. The regularization parameter, β, sets the trade-off between fitting the measured dataand model smoothness.In equation (2.1) dobs is the measured data and F[m] = dpred is the predicted data (theforward model). These data are traveltimes from sources to receivers. The model vector, m,is the slowness discretized over a volume and reorganized as a vector for ease of computation.The matrix Wd is a diagonal matrix that contains estimates of the reciprocal of the standarddeviation of each measurement. The matrices Wi are discrete forms of each component of thegradient operator, and can also include other measures such an identity matrix to encouragethe recovered model to be close to the reference model, mref . The parameters αi allow eachregularization term to be weighted independently.The traveltimes are estimated by the eikonal equation, which is a non-linear function of mand is discussed in section 2.3. The objective function is minimized iteratively by linearizingthe forward model about the current model,F[mk + δm] ≈ F[mk] + Jkδm, (2.2)292.2. Inversion Formulationwhere k is the iteration number and J is the Jacobian. Setting the gradient of the objectivefunction (with respect to δm) to zero results in [2]WdJk√βα1W1...√βαNregWNregδm =Wdδd√βα1W1(mref −mk)...√βαNregWNreg(mref −mk), (2.3)where δd = dobs − F[mk]. Equation (2.3) is large and sparse (thanks to a sparse approximateJacobian discussed in section 2.4) and can be solved using the LSQR algorithm [75, 2]. Themodel is then updated according tomk+1 = mk + ζδm. (2.4)In equation (2.4) the parameter 0 < ζ ≤ 1 improves convergence, particularly when δmis large and the linearization in equation (2.2) is poor [44]. We also verify that φk+1 < φkand reduce ζ if necessary [44]. This involves another forward modeling operation, but it isnot particularly wasteful because the results can be used for the calculation of the next modelupdate as long as it is not necessary to reduce ζ (thereby requiring another forward computationto check that φk+1 < φk).The process continues until m converges to the one that minimizes equation (2.1). At theminimum point the gradient of the objective function is zero, but due to scaling and precisionissues, we can not realistically achieve ∇mφ = 0. Instead, we use a relative gradient test thatrelates the relative change in φ to the relative change in a model element, mi, that is also usedby Haber and Oldenburg [44] [and derived in 29, pp. 159–161],max[ |(∇mφk)i(mk)i|φ(mk); i = 1..M]< tol, (2.5)where M is the number of model parameters. As β is reduced, it becomes more difficult to findthe minimum of the objective function. This is likely due to the effects of approximate Jacobianbecoming more apparent as smoothing is reduced. At some point the scaled solution ζδm fails302.2. Inversion Formulationto reduce φ, even for very small ζ. Additionally, the convergence test in equation (2.5) reliesinternally on the gradient being accurate, which is computed using the approximate Jacobian.To prevent infinite looping and to allow the inversion to proceed despite the lack of precision,the relative change of the model is also used as a stopping criterion (from [44]):‖mk+1 −mk‖max(‖mk+1‖ , ‖mk‖)< δ (2.6)If either ζ has been reduced to the point that the model is not changing significantly, or themodel has become stationary because a local minimum has been reached, the iterations arestopped. Finally, the process is stopped if the number of iterations becomes excessive.2.2.1 T0 EstimationIn our case, an additional unknown that must be included in m is an additive constant, T0, thatis part of all of the measured traveltimes. This variable represents various unknown parametersin the measurement setup that combine to produce an unknown offset between generation ofthe transmitted impulse and start of sampling at the receiver. These unknowns are assumedconstant and because only a single pair of radar antennas were used for all measurements, weonly need to find a single offset that applies to all of the measured data.T0 is included in the inverse problem by appending it to the model vector, m, so that it isthe last element. When computing the predicted data, this value is added to each traveltime(see Section 2.3). This way, the Jacobian matrix and the Wi matrices are easily augmented toinclude T0. The Jacobian is padded with a column of ones,J =∂F1∂m1 · · ·∂F1∂mM−1 1... ... ...∂FN∂m1 · · ·∂FN∂mM−1 1, (2.7)where the subscripts represent indexes into the column vectors, N is the number of measure-ments and M is the number of model parameters including T0.Similarly, a column of zeros is added to the regularization matrices, Wi, so as not toinclude T0. In our preliminary paper [101] there was an unnecessary row of zeros added to312.3. Eikonal Solverthese matrices.2.3 Eikonal SolverThe forward modeling operation must be computed each time equation (2.3) is solved, so it isimportant that it be as efficient as possible. We estimate the travel times from a given sourceto an array of receivers with the eikonal equation,|∇T (x)| = 1v(x) = s(x), (2.8)where T (x) is the travel time, v(x) is the velocity and s(x) is the slowness and all are functionsof position. This is solved with the boundary condition that T (x0) = 0, where x0 is the sourcelocation. The offset, T0, is added after finding T (x) with this boundary condition.The eikonal equation has been used extensively as a basis for traveltime inversion (forexample: [2, 101, 107, 61]; and [22, p. 55]). Here we use an easily-implemented fast sweepingmethod of solving the equation, due to Zhao [124]. Leung and Qian [61] also use this methodand use an adjoint formulation to find the gradient of the objective function. Here, instead,we use a simple ray-tracing technique to estimate the sensitivities, which results in a sparsesensitivity matrix that we use for a non-linear Gauss-Newton solver.The fast sweeping method of Zhao [124] is much easier to implement than a fast marchingmethod (such as [107, 47]) and will always produce a solution. Ray paths are used to computethe sensitivity matrix, J, and are found by following the traveltime gradient from receiver tosource (discussed in section 2.4). This method has the substantial benefit over two-point raytracing that a ray path can always be found (unless the model is excessively rough) and it isnot an iterative process, but it is restricted to first-arrivals.Here we outline the application of Zhao’s fast sweeping upwind-difference scheme [124] tothe solution of the 3D eikonal equation. For details about the derivation and performance ofthe method, the reader is referred to Zhao’s paper.In 3D, Zhao’s algorithm works on a cubic grid with cell spacing h where the traveltime isstored at each grid point. At the start, each grid point is initialized to a large value (larger322.3. Eikonal Solverthan any expected traveltime in the final solution) except for the starting point, which is setto zero. In our case the starting point is not generally exactly on a grid point, so we computethe traveltime to nearby points assuming a straight ray and a weighted average of the slownessat the nearby points. Then the entire grid is swept in alternating directions and updatedtraveltimes are computed at each grid point. The alternating sweep directions ensure thatthe traveltime is updated along characteristics (rays) in any direction. For the 3D case, thereare a total of eight sweep directions to accommodate all combinations of positive or negativesweeping on each coordinate axis. This results in eight sets of triply-nested loops. The sweepingis continued until the changes in traveltimes are below the desired tolerance. The calculationthat is performed at each node (at the inside-most level of these loops at position (x, y, z)) isoutlined in the following.The minimum traveltime of nodes that are adjacent to the current node in each of thecoordinate directions is found and stored in the variables a1 through a3:a1 = Tx min, a2 = Ty min, a3 = Tz min (2.9)In equation (2.9), at position (x, y, z), Tx min = min(Tx+h,y,z, Tx−h,y,z) and similarly for theother coordinate directions. Tx,y,z is the traveltime value that is currently stored at (x, y, z).At the edges of the domain the node outside the domain is ignored (i.e. the minimum is replacedwith the value of the adjacent node inside the domain).Next, these values are sorted such thata1 ≤ a2 ≤ a3. (2.10)If the traveltime at the node is already less than the smallest neighboring value, then we havealready found the minimum time at the current node: If a1 ≥ Tx,y,z then move on to the nextpoint determined by the sweep direction.If Tx,y,z > a1 then we must compute a new time. Up to three calculations are performedper cell, depending on the traveltime values obtained when performing the calculations. Thesecalculations (equations (2.11), (2.12) and (2.13)) are solutions of Zhao’s upwind-difference ap-332.4. Jacobian Estimationproximation of the eikonal equation, and are arrived at by following his method. The firstcomputation finds a tentative traveltimeTt = a1 + sh, (2.11)where s is the slowness at the current node. If Tt ≤ a2 then we are finished with this node andset Tx,y,z = min(Tt, Tx,y,z).Otherwise, we computeTt =a1 + a2 +√2s2h2 − (a1 − a2)22 (2.12)and then test to see if Tt ≤ a3. If so, then we are finished and set Tx,y,z = min(Tt, Tx,y,z).If not, we perform the last computation,Tt =a1 + a2 + a3 +√3s2h2 + 2(a1a2 + a1a3 + a2a3 − a21 − a22 − a23)3 , (2.13)and set Tx,y,z = min(Tt, Tx,y,z).The grid size used for the fast sweeping solver is set at an integer division of the model grid.This allows a greater accuracy when computing travel times without an associated increasein the matrix system that must be solved at each iteration of the inversion. The slownessand traveltimes are scaled to be near order one by setting s(x) =√ǫr(x) (ǫr is the relativepermittivity).2.4 Jacobian EstimationThe length of the ray path in each cell is an approximation of the sensitivity [2, 101]:Jij =∂Ti∂mj= ∂∂mj∑mjlij (2.14)≈ lij342.5. RegularizationIn equation (2.14) the source-receiver combinations are indexed by subscript i and the modelcells are indexed by subscript j. The approximation is due to the fact that the length lij is notindependent of m (in that cell or even nearby cells).In a similar manner to Aldridge and Oldenburg [2] (and [107] and [80]) the rays are tracedfrom receiver to source by following the steepest descent through the travel time field. Thegradient in each cell is assumed constant.The tracing algorithm is described by the flow chart in Figure 2.1. Starting at the receiver,a straight segment is drawn from the receiver to a cell face in the direction of the negativetraveltime gradient (−∇T ). From there until the source cell is reached, straight segments aredrawn from face to face in each cell along that cell’s gradient. When the cell containing thesource location is reached, a straight segment is drawn directly from where the ray entered thecell to the source, regardless of the gradient direction within that cell.Occasionally, the gradient in an adjacent cell directs the ray back into the cell from whereit came, thus trapping the ray along a face or an edge of a cell. This can happen when thegradient is nearly parallel to a face between cells and −∇T in either cell points slightly towardthe face between them. Simply following −∇T from cell to cell would result in an infinitenumber traversals across that cell boundary. To deal with this situation, a list of cells that theray has passed through is maintained and the ray is not allowed to enter a cell twice. Whenthe trajectory of a ray takes it into a cell where it has been, the trajectory is modified bysetting the component of the gradient in the direction of the previously entered cell to zero andthen recomputing the trajectory. This is referred to as finding the next best direction in theray tracing flow chart shown in Figure 2.1. This causes the ray to follow a face or an edge inbetween cells when the gradient in adjacent cells pushes the ray toward the face or edge in eachcell.2.5 RegularizationThe inversion process proceeds starting with a large β and then reducing it until a suitable levelof regularization is reached. Thus the model is, in a sense, cooled and structure crystallizesout of it. At what point this cooling process should be stopped or, alternatively, what value of352.5. RegularizationStart:Position = Receiver locationDraw segment along negativegradient to cell boundaryLeavingDomain?CrossingPath?Change segmentto next bestdirectionInsourcecell?Compute and savefinal segmentDoneYesNoYesNoYesNoSave ray segmentRay Tracing Flow ChartFigure 2.1: Ray tracing algorithm.β provides the best model, is an important question. Unfortunately the answer is not alwaysclearly defined, and it is often easier to let the cooling process continue past the optimal valueand then use one’s judgment to determine the best value simply by looking at the inversionresults. This potentially leads to a biased selection of β. Automated methods of finding β areof interest because they can reduce this bias. In this section we briefly discuss two commonmethods for finding the regularization parameter: the l-curve [45, 46, 36] and generalized cross-validation (GCV) [44, 36]. We also discuss some robust extensions to GCV [62, 63]. If the dataerrors are well known, it is also possible to use an expected value of φd to stop the cooling.However, the true magnitude of the errors is rarely known accurately.2 Plots and further2We were somewhat mistaken in stating in the preliminary work [101] that the magnitude of these errors waswell known (and also assumed normally distributed), leading to the use of an expected value of φd for selectingthe regularization parameter.362.5. Regularizationdiscussion are found in sections 2.6 and 2.7.The l-curve is a plot of φd(β) vs φm(β), usually on a log-log scale, and generally results ina characteristic “L” shape, often with a clearly defined corner. The value of β at the corner(point of maximum convex curvature) of the plot is considered to be a good trade-off betweenmodel smoothness and data fitting [36, 46]. If one were to use this method, the curvature canbe estimated during the cooling process to adjust the selection of β to home in on the optimalvalue [36]. Usually the l-curve method works very well, but in some situations the log-log plotof φd and φm is not “L”-shaped and other methods should be looked at.The GCV method seeks to find a model that best predicts missing data, but can resultin an under-smoothed model and is sensitive to correlated data [63]. Modifications to GCVsuch as robust GCV (RGCV) [62] and strong robust GCV (R1GCV) [63] help resolve thoseissues. Both of these modifications to GCV add a measure of the influence of each data pointon the regularized solution as a penalty function, but differ in the norm used to compute theinfluence. The degree to which the penalty function affects the RGCV and R1GCV functionsis determined by a manually chosen parameter, γ, with the range 0 < γ ≤ 1. γ = 1 representsno influence from the penalty function and both RGCV and R1GCV reduce to the usual GCVfunction.Farquharson and Oldenburg [36] extend GCV to a non-linear inversion problem with avery similar objective function to ours. The extension is achieved by linearizing the objectivefunction about the current model using equation (2.2) in the same way that the model update,δm, is found in equation (2.3). Their result is repeated here for reference:V (β) =∥∥Wdδd−WdJM−1(JTWdTWdδd+ r)∥∥2[tr(I−A)]2 , (2.15)where V (β) is the GCV function to be minimized and A and r are given by equations (2.21)and (2.23). The value of β that minimizes equation (2.15) is considered to be the best value.When this function and those in equations (2.16) and (2.17) are evaluated, all of the values arethe ones computed during the most recent iteration (except for β).Using the extension of GCV to a non-linear inversion problem by Farquharson and Olden-burg [36] as a guide, we obtain RGCV, VR(β, γ), and R1GCV, VR1(β, γ), which are modified372.5. Regularizationforms of the functions of Lukas [62, 63]:VR(β, γ) =γ + (1− γ)µ2(1− µ1)2N−1∥∥Wdδd−WdJM−1(JTWdTWdδd+ r)∥∥2 (2.16)andVR1(β, γ) =γ + (1− γ)µ12(1− µ1)2N−1∥∥Wdδd−WdJM−1(JTWdTWdδd+ r)∥∥2 , (2.17)whereµ1(β) = N−1 tr(A), (2.18)µ2(β) = N−1 tr(ATA), (2.19)µ12(β) = N(µ1 − µ2)/β, (2.20)A(β) = WdJM−1JTWdT , (2.21)M(β) = JTWdTWdJ+ βNreg∑i=1WiTWi, (2.22)andr(β) = βNreg∑i=1WiTWi(mref −m). (2.23)In equations (2.15) through (2.23), all quantities are those found during the most recent iter-ation, except for β. In equation (2.20) the parameter λ in the original work [62, 63] has beenreplaced with β/N due to a different scaling of the regularization parameter in the objectivefunction. In equations (2.16) through (2.20), N is the number of data points. The influencematrix (equation (2.21)) is defined byAij =∂(Wddpred)i∂(Wddobs)j, (2.24)with similarity to Wahba [112, pp. 51–52, 116] [111].When the number of measurements is large, the matrix A cannot easily be formed explicitly,but the result of applying it to a vector can be computed (ex. iteratively or with an LU382.6. Synthetic Experimentsdecomposition). To compute the trace we use the stochastic estimator of Hutchinson [50].To produce the plots found in this paper, several estimations of the trace are computed andaveraged.The β found by minimizing any of the GCV functions is an estimate of the optimal β andcan be used during the inversion process to control the reduction of β and to halt the inversionprocess once the appropriate value is reached. Farquharson and Oldenburg [36] and Haberand Oldenburg [44] show that at early iterations the estimate of β provided by GCV may beincorrect, but as the inversion process proceeds the value found by GCV converges. It shouldalso be noted that an approximate Jacobian is used here, and may influence the GCV results.2.6 Synthetic ExperimentsA synthetic experiment is performed using essentially the same borehole layout and transmitterand receiver placement as the data collected at LSBB (see Section 2.7). The borehole arrange-ment is visible in Figs. 2.2 and 2.8. The boreholes are approximately 20 m deep. The cluster offour holes (B2-B5) are used to generate the synthetic data. For the holes directly across fromeach other, five transmitters are used at equally spaced positions in each hole. For the holesdiagonally adjacent to each other, three equally spaced transmitter positions are used. Thereceived signal is sampled at 5 cm intervals in each hole. The total number of recorded tracesis 17709.Two sets of synthetic data are generated. The first set of traveltime data is computed usingthe fast sweeping eikonal solver discussed in section 2.3. Random noise with zero mean and astandard deviation of 10 % of the mean travel time was added to this data set before inversion.This relatively large error is intended to demonstrate the capability of the inversion to functionin the presence of substantial noise. Of the total number of measurements, a total of 12398 datapoints are used to match the selection of traces that are picked successfully for the measureddata discussed in Section 2.7.As a more realistic test of the inversion algorithm, a second set of traveltime data is com-puted using the MEEP FDTD solver [73]. In order to simplify the use of MEEP, straight andperfectly vertical boreholes are used in the simulation. The FDTD simulation produces time-392.6. Synthetic Experimentsdomain waveforms at the receiver locations that are essentially the same as those producedby a real borehole impulse radar system (such as the MALA˚ units used to make the actualmeasurements at LSBB). These time-domain waveforms were picked using a correlation-basedpicker, discussed in [101]. Due to difficulty in accurately determining the first-arrival timewhen the first-arrival is not the strongest arrival, synthetic traveltime data is not available atall receiver locations. The fast-sweeping eikonal solver discussed here computes first arrivalsonly, so picking results where the picker failed to find the first arrival are discarded. This isdone by manually checking the picked data. The total number of data points for this data setis 15553. No noise was added to these data. The picked FDTD data naturally have substantialcorrelated errors if the solution of the eikonal equation is taken to be the true traveltime. It isthe effect of these correlated errors on the inversion that is being investigated with this dataset. Further, the picking process adds a small amount of error to the synthetic traveltime.The starting model used for the inversion of both data sets is a uniform material with arelative permittivity of 11. The true model is shown in Figure 2.2. The background mediumconsists of two layers of differing permittivity. There is a truncated ball-shaped inclusionin the top portion of the model of a higher permittivity and a rod-shaped inclusion of lowerpermittivity that runs at an angle through the model. The isosurfaces in Figure 2.2 are clusteredat the material interfaces because of the infinite gradients there.The inversion process is allowed to run freely until the solver breaks down and is unable toreduce φ for a reduced β. During the inversion, the regularization parameter, β, is reduced bya factor of 0.7 at each iteration of the outer-most loop. This provides us with enough inversionresults to clearly see how the various criteria for selecting β in section 2.5 relate to each otheron the l-curve.Figure 2.3 shows the GCV, RGCV, and R1GCV curves computed using an intermediateinversion result (i.e. somewhat before the final β) for both synthetic data sets. The curves varyonly slightly during most of the inversion process. Farquharson and Oldenburg [36] and Haberand Oldenburg [44] offer some examples of the convergence of GCV in a non-linear setting.The parameter γ used with the RGCV and R1GCV functions is chosen to be γ = 0.8, therebygiving a relatively low level of robustness (i.e. closer to the standard GCV). With this selectionof γ, the GCV and RGCV curves are nearly identical, with only a slight difference at small β.402.6. Synthetic Experimentsx (m)y (m)-20-15z (m)-10-5123450 0910111213Relative Permittivity1 21415True ModelFigure 2.2: The true model for the synthetic case. The volume is split vertically with ǫr = 10above a region with ǫr = 13. There is a ball-shaped inclusion with ǫr = 15 visible in the upperpart of the model (and cut in half by the x = 0 plane). A rod-shaped inclusion with ǫr = 8 runsdiagonally through the model. The pane on the right shows various isosurfaces of the relativepermittivity. The boreholes are shown as vertical black bars in the isosurface plot.A reduction of γ would increase this difference and also shift the RGCV minimum to the right.The R1GCV function curves sharply upward as β is reduced, thus making the minimum veryclearly defined. Figure 2.4 shows where the three cross-validation results lie on the l-curves foreach data set. The R1GCV result happens to land very close to the corner of the l-curve forthe eikonal data set, but the FDTD data set inversion fails before a corner is reached. TheGCV and RGCV results for the eikonal data set are slightly to the right of the corner (a lowerdegree of regularization), but predict extremely small regularization parameters for the FDTDdata set (the inversion algorithm breaks down before they are reached). For comparison, thetrue model has a much higher model norm (φm ≈ 439) than the recovered models, indicatingthat the inversion process is finding smooth models that fit the data.There are a few potential reasons for the failure of the FDTD data inversion before thel-curve corner. Most fundamental, the eikonal equation (solved for first arrivals) does not pre-412.6. Synthetic Experiments10010110-6 10-4 10-2 100 102 104 106 108Eikonal Synthetic Data GCV CurvesβVGCVRGCVR1GCV(a)10-310-210-110010110210310-10 10-5 100 105 1010FDTD Synthetic Data GCV CurvesβVGCVRGCVR1GCV(b)Figure 2.3: GCV, RGCV, R1GCV curves for the synthetic experiment using data generated by(a) the fast-sweeping eikonal solver and (b) picked FDTD waveforms from MEEP. The functionshave been scaled such that their value at large β is unity, and γ = 0.8 has been used for RGCVand R1GCV. The circle, cross, and triangle indicate the minimums of the GCV, RGCV, andR1GCV curves respectively.cisely model all of the finite-frequency physics that the FDTD simulation models. Additionally,there may be very weak first arrivals that are missed during the picking process where a laterand much stronger arrival is erroneously chosen. Both of these issues can cause correlatederrors, which are not accounted for in the inversion formulation.10510-3 10-2 10-1 100 101 102Eikonal Synthetic Data L-curveφmφ dGCVRGCVR1GCV(a)10310410-2 10-1 100 101FDTD Synthetic Data L-curveφmφ dR1GCV(b)Figure 2.4: L-curve for the synthetic experiment using data generated by (a) the fast-sweepingeikonal solver and (b) picked FDTD waveforms from MEEP. The points on the curve found byGCV, RGCV, and R1GCV are indicated. For the FDTD data set, GCV and RGCV producepoints on the l-curve that were not reached during the inversion process.422.6. Synthetic ExperimentsSlices at x = 0.5 m through the recovered models for the two data sets are shown alongsidethe true model in Figure 2.5. Both data sets are created such that T0 ≈ 0. The value of T0 thatis recovered by the inversion ranges from −0.63 ns to −0.86 ns for the eikonal data set, and is−2.2 ns for the FDTD data set. These values are much lower than the minimum traveltime ofabout 20 ns.Although the model recovered from the eikonal data set using GCV (or RGCV) shown inFigure 2.5b appears somewhat noisy compared to the true model, Figure 2.5a, and also tothat recovered using R1GCV, Figure 2.5c, it is a better fit to the true model. That this istrue is clearly visible by computing the norm of the difference between the true model and therecovered models. A plot of this misfit for l1 and l2 norms is shown in Figure 2.6a. Althoughthe R1GCV choice of β appears closer to the corner of the l-curve in Figure 2.4a, it produces aslightly over-smoothed recovered model and does not resolve the anomalies quite as clearly asthe model found at the GCV value of β. This is only true with the eikonal data set where thesynthetic observed data is unrealistic. The R1GCV method finds a reasonable β for the FDTDdata set when both GCV and RGCV do not.Although R1GCV finds a good model for the FDTD data set, the model is not quite asclose to the true model as it is for other values of β, as indicated in Figure 2.6b. There is acloser l2 fit with a sightly larger β. The recovered model that best fits the true model is shownin Figure 2.7. In this case R1GCV slightly under-estimates β rather than over-estimates it aswith the eikonal data set. R1GCV could be made to produce a larger value of β by decreasingγ thereby selecting a greater degree of robustness. Regardless, the value of β produced withγ = 0.8 is very good and higher values of β only slightly improve the inversion results, asindicated with the l2 fit in Figure 2.6b.432.6. Synthetic Experimentsz (m)y (m)-20-15-10-50 1 2 3 4 589101112131415Trueǫ r(a)z (m)y (m)-20-15-10-50 1 2 3 4 589101112131415GCV, RGCVǫ r(b)z (m)y (m)-20-15-10-50 1 2 3 4 589101112131415R1GCVǫ r(c)z (m)y (m)-20-15-10-50 1 2 3 4 589101112131415ǫ rFDTD R1GCV(d)Figure 2.5: True model and synthetic results shown as a slice at x = 0.5 m: (a) true model,(b) recovered model using GCV or RGCV and the eikonal data set, (c) recovered model usingR1GCV and the eikonal data set, (d) recovered model using R1GCV and the FDTD data set.442.6. Synthetic Experiments20002500300035004000450010-1 100 101 102 103 104 105 106Eikonal Synthetic Data Inversion Misfitp = 1p = 2GCVRGCVR1GCVβ‖ǫr,true−ǫ r,rec‖p(a)1500200025003000350040004500500010-1 100 101 102 103 104 105 106FDTD Synthetic Data Inversion Misfitp = 1p = 2R1GCVβ‖ǫr,true−ǫ r,rec‖p(b)Figure 2.6: Misfit vs. regularization parameter for the synthetic experiments (a) eikonal and(b) FDTD. The p = 2 norm is scaled for comparison with the p = 1 norm. The associatedpoints on each curve for GCV, RGCV and R1GCV are also shown. The GCV and RGCV pointsare not shown for the FDTD results because they are at smaller β than was reached during theinversion. The true model is ǫr,true and the recovered model is ǫr,rec.z (m)y (m)L2 Minimum-20-15-10-50 1 2 3 4 589101112131415ǫ rFigure 2.7: Synthetic result from the FDTD data set that best fits the true model, shown as aslice at x = 0.5 m.452.7. Inversion of GPR Data from LSBB2.7 Inversion of GPR Data from LSBB250 MHz GPR data were collected at LSBB in a total of five boreholes (see Figure 2.8), four ofwhich are in a diamond shaped cluster [101]. The boreholes are approximately 20 m deep. Thedata from these four holes is used in this work. For the pairs of holes directly across from eachother, five equally spaced transmitter locations (from top to bottom of the holes) were usedwhile data were collected in the opposite holes at intervals of 5 cm. For the holes diagonallyadjacent to each other, three equally spaced transmitter locations were used. The same boreholearrangement is used for the synthetic model and is visible in Figure 2.2. Borehole deviationsfrom vertical were measured and are used to adjust the transmitter and receiver locations usedin the inversion.B1B5B4B2B3Tunnel WallScale: 5mBorehole LocationsFigure 2.8: Plane view of the borehole positions.The impulse radar data is picked for first arrivals with a simple correlation procedure [101].Measurements are correlated with a synthetic waveform with similar characteristics to theactual data. The radar trace with the strongest correlation is used as the starting point, andfurther traces are picked by again correlating the synthetic impulse with the measured data,but the chosen correlation peak is limited to a small window in time around the picked timeon the adjacent trace. This way, the picking procedure follows a single crest of the impulsefrom one receiver to the next. If the correlation peak does not meet a certain threshold, thenthe measurement is discarded. The results of the picking process are checked manually to462.7. Inversion of GPR Data from LSBBensure that the correct crest of the impulse is picked. Of the 17709 traces collected, 12398 aresuccessfully picked and used for this inversion.2.7.1 Sources of ErrorThe measurement data has a number of additional sources of error in comparison to the syn-thetic case. The exact physical location of the boreholes at depth, and therefore the radarantenna positions, is not known precisely. There are errors incurred during the traveltime pick-ing process. Further, the eikonal equation and solver only provides an approximation of thetrue nature of the electromagnetic system consisting of the earth, radar units, boreholes, andother objects.A simple straight-ray relationship between antenna location and traveltime in a uniformvelocity model is used to relate error in antenna depth and antenna separation to error inmeasured travel time. Rough estimates of the error associated with the vertical and horizontalpositions of the antennas are then used to compute estimates of the standard deviation ofthe traveltime due to antenna position uncertainty. This is combined with the uncertainty ofthe traveltime picker (assumed to have a standard deviation equal to the sampling interval ofthe GPR system) to obtain an overall estimate of the standard deviation of each traveltimemeasurement. The reciprocal of the estimated standard deviation of each measurement is placedat the associated location on the diagonal of Wd. The matrix Wd is a diagonal matrix.2.7.2 Measurement Inversion ResultsAs in the synthetic case, the inversion was allowed to run until the solver broke down (stagnated)to allow plotting as much of the l-curve as possible. The same starting model (uniform withǫr = 11) was used. The GCV, RGCV, and R1GCV curves are plotted in Figure 2.9. Theminimums of the GCV and RGCV curves are poorly defined and have a β value so low thatthe inversion results are excessively under-smoothed. The R1GCV curve is the only one thatproduces a reasonable regularization parameter. This poor behavior is also seen with thesynthetic FDTD data set.The resulting l-curve is shown in Figure 2.10. A beginning of a corner of the “L” is startingto form but this is where the solver breaks down and is no longer able to reduce φ.472.8. Geological Correlation10-210-110010110-6 10-4 10-2 100 102 104 106 108Measured Data GCV CurvesβVGCVRGCVR1GCVFigure 2.9: GCV, RGCV, and R1GCV curves for the inversion of measured data from LSBB.The functions have been scaled such that their value at large β is unity, and γ = 0.8 is used forRGCV and R1GCV. The circle, cross, and triangle indicate the minimums of the GCV, RGCV,and R1GCV curves respectively.10310410-3 10-2 10-1 100 101 102Measured Data L-curveφmφ dGCVRGCVR1GCVFigure 2.10: L-curve for the measured data from LSBB. The point on the curve found byR1GCV is indicated. The points for GCV and RGCV represent the nearest β encounteredduring the inversion process.The recovered electrical permittivity is shown in Figure 2.11 for the level of regularizationgiven by R1GCV. At this level of regularization there is structure that looks suspiciously likeinversion artifacts, visible as either high or low permittivity regions along the boreholes nearthe transmitter locations, so further reduction of β is probably not warranted.2.8 Geological CorrelationThe study area is located in the Underground Low Noise Laboratory (LSBB) gallery (http://lsbb.oca.eu)in southern France. The LSBB is nested in a karst aquifer at depths of 200 to 500 m. At 250-m482.8. Geological Correlationx (m)y (m)-20-15z (m) -10-5123450 08101214Relative Permittivity1 2R1GCVz (m)y (m)R1GCV-20-15-10-50 1 2 3 4 5ER891011121314Figure 2.11: Inversion results for the measured data found with R1GCV shown as a solid volumeand as a slice taken at x = 1 m.depth, within a N170-trending gallery, five vertical cored boreholes of 20 to 22 m deep weredrilled in the gallery floor in order to explore the geological structure and to develop fluid injec-tion experiments (Figure 2.8 [43]). Previous studies have characterized in detail the geological,petrophysical and hydromechanical properties of the area where the GPR data were collected.In Figure 2.12, we summarize the main results obtained by: (1) field observation [53]; (2) lab-oratory measurement on samples from boreholes B3: VP and porosity (φ) [52]; (3) boreholelogging in B3 (borehole viewer images and water content variations relative to a reference valuetaken in the boreholes close to the surface); and (4) injection tests performed between twoinflatable packers at different depths in B3 and analyzed with the hydraulic numerical codeTOUGH2 [84] in order to estimate the hydraulic permeability (k) and the storage coefficient(S) [51].The boreholes intersect a sedimentary series made of fractured limestone (lower cretaceouscarbonates) oriented N80◦E and dipping 10 to 20◦ to the South. The fracture network mainlyconsists of sub parallel N030 and few N120 joints. This series can be divided into three me-492.9. Concluding Remarkschanical units: one porous layer (φ = 20 %) 5 to 6 m thick with few fractures (Unit 2, Fig-ure 2.12), intercalated between two fractured and low porosity layers (φ = 3 %) (Units 1 and 3,Figure 2.12). These variations of the petrophysical properties have a strong impact on the me-chanical properties distribution inside the sedimentary series. The low-porosity layers presenthigh VP values (7300 ms−1) and the porous layer has low VP values (3600 ms−1). In the sameway, the evolution of the permeability and the specific storage coefficient are strongly influencedby the petrophysical properties. The low porosity layers have a low specific storage coefficient(up to 6.9 · 10−3 m−1), and permeabilities ranging between 1.8 · 10−12 m2 to 1.75 · 10−15 m2depending on the fracture density. Inversely, the porous and low fractured layer has a higherspecific storage (2.0 · 10−2 m−1) and intermediate permeability values (around 3.0 · 10−13 m2).These variations of the hydraulic properties are well highlighted by the borehole logging. Thereis higher water content in the porous layer and lower water content in the low-porosity layers.These geological experiments were performed to study the architectural and hydromechani-cal properties of a small fault zone cutting the gallery near the borehole B1. It has been shownthat the porosity inside the porous layers decreased up to 50 % in B1. We can suppose thatthis reduction of the porosity is linear and can start near the boreholes B2, which is located atthe boundary between the host rock and the damage zone.The foregoing geological measurements and analyses provide the ground-truth for our in-version. The high water content in the middle highly porous layer of the rock structure isconsistent with the high relative electrical permittivity visible in the center of Figure 2.11.Water has a high electric permittivity and its presence in rock strongly affects the observedelectric permittivity of the rock. This can be seen through the complex refractive index method(CRIM) [32] or Topp [99] relationships.2.9 Concluding RemarksA straightforward traveltime inversion algorithm has been presented, along with its applicationto synthetic and measured data. The extension and use of a modified form of GCV calledstrong robust GCV, R1GCV, to the non-linear Gauss-Newton inversion problem has also beenpresented for the first time (to the authors’ knowledge).502.10. Acknowledgements-5-1015-2054321 0 0 12y (m)x (m)z (m)-20B3B2B4B5B3-510 20 3000 7000Vp (m.s )-1Vp and Porosity measurements on sample taken on core Porosity (%)00-5-10-15-18-10-150-5-10-15-180-5-10-15-18Borehole logging Hydraulic properties estimated from TOUGH20 1 S (10 m )2 -23 -10K (10 m )2-122 1 0-5-10-15-18-5 5 0 0-5-10-15-18Water content (%)a) b) c) d) e)f) g)NFigure 2.12: Vertical profiles of (a) relative water content, (b) permeability, (c) storativity, (d)porosity, and (e) P-wave velocity measured along the borehole B3. (f) Borehole viewer imageillustrating the lithology contrast along B3. (g) 3D view of the geological structure [53, 52, 51].A synthetic experiment using FDTD data produces similar behaviour with regards to thel-curve and GCV, RGCV, and R1GCV results when inverted as the measured data set. Thisis likely due to the implicit approximation that using the eikonal equation to model finite-frequency traveltimes entails. Despite the poor l-curve and GCV/RGCV results, the syntheticFDTD experiment succeeds in finding a reasonable approximation of the true model. The in-verted measured data also corresponds well with the geological studies that have been performedin the area.2.10 AcknowledgementsThe authors would like to thank the providers of funding and support for this project. Thiswork is financed by the ANR Captage de CO2 through the HPPP-CO2 Project, by the PACAcounty through the PETRO-PRO Project. Golder Inc. provided in-kind donation for the fieldwork. MALA˚ Geoscience provided two 250 MHz borehole antennas as in-kind support. LSBBand the re´gion PACA provided the infrastructure. Direct funds for travel and perdiem wereprovided by the project ANR MAXWELL.512.10. AcknowledgementsThe Natural Sciences and Engineering Research Council of Canada (NSERC) providedsubstantial support through scholarship and research grants.Calculations were performed using MATLAB, GNU Octave, and MEEP.52Chapter 3A Graphical Derivation of thePerpendicular Ray Jacobian3.1 IntroductionIn the following we present a graphical derivation of the perpendicular component of the 3Dray Jacobian for planar rays. This part of the ray Jacobian describes the out-of-plane spreadingassociated with a ray and is the only difference between the 2D and 2.5D cases in ray theory[22, p. 394]. The difference between the two cases for full-wave simulations is more complex,but knowledge of the perpendicular component of the Jacobian allows one to preprocess real-world acoustic and ground-penetrating radar data into a form that is usable for 2D full-waveinversion, provided the geology is sufficiently 2.5D.Planar rays imply that the derivative of velocity in the direction normal to the plane iszero. Here, we satisfy this by assuming that the material properties do not vary normal tothe plane, but the arguments are also applicable as long as the normal derivative vanishes. Inthat case, the velocity does not vary in the direction normal to the plane if one stays in aninfinitesimal region near the plane. Previous related work focused on mapping 3D acousticdata to 2D includes that of Cˇerveny´ and Ravindra [24, pp. 74-93], Cˇerveny´ and Hron [23], andBleistein [19].The 3D to 2D data transfer function for the acoustic-wave problem derived in [121] (see2This article has been accepted for publication in the Bulletin of the Seismological Society of America: D.Van Vorst, J. Virieux, and M. Yedlin, “A Graphical Derivation of the Perpendicular Ray Jacobian”, vol. 102,no. 6, pp. 2452–2457, 2012, c©Seismological Society of America.533.2. Graphical Derivation of J⊥[120]) isTF = u2D(w)u3D(w)= −iπH(2)0 [ωT (s)]√v(0)J3D(s)T (s)sin(θ0)J2D(s)eiωT (s). (3.1)In equation (3.1), H(2)0 is a Hankel function, v(0) is the velocity at the source, s is the arc lengthof the ray from source to receiver, T (s) is the travel time from source to receiver, J3D(s) is the3D Jacobian at the receiver, J2D(s) is the 2D Jacobian at the receiver, and θ0 is the take-offinclination. The transfer function transforms the measured field, u3D(w), into an equivalent2D field, u2D(w). Here, and in the rest of the paper, it is assumed that the complex timedependency is eiωt.Ray-based inversions usually precede full-wave inversions and provide much of the infor-mation that is required to compute equation (3.1): ray paths and an approximate velocitydistribution (travel time is already known from picking first arrivals). In the following we show,through graphical and paraxial methods, that aside from travel time, the remaining quantitiesin the radical in equation (3.1) amount to an integral of velocity along a ray.For both the graphical and paraxial derivations, we assume that the coordinate system hasbeen selected such that the plane where the 2D inversion is to be performed is y = 0. Whenworking in 3D, spherical coordinates are used where θ is the inclination, and φ is the azimuth.We are interested in how rays behave very close to the y = 0 plane, so the azimuth angle isvery small and is represented as dφ, as indicated in Figures 3.1 and 3.2. It is also assumed thatthe material parameters do not vary as a function of y.3.2 Graphical Derivation of J⊥We begin with the graphical derivation of the ratio of ray Jacobians. First, we show that the3D ray Jacobian can be factored into a parallel part, J‖ = J2D, and a perpendicular part J⊥.Next, a simple expression to compute J⊥ is derived with the assumption that the in-plane rayis straight. Finally, it is shown that the expression is valid for curved rays.543.2. Graphical Derivation of J⊥3.2.1 Jacobian FactorizationThe 3D ray Jacobian, J3D(s, θ0, φ0), relates the cross sectional area of the ray tube, dS, to theassociated infinitesimal changes in the take-off angles dθ and dφ:dS = J3Ddθdφ. (3.2)An illustration of the 3D ray Jacobian for a planar ray is shown in Fig. 3.1.xyzθ0dφA: (θ0, 0)B: (θ0, dφ)C: (θ0 + dθ, dφ)D: (θ0 + dθ, 0)dSFigure 3.1: A ray tube cut off at a particular arc length, s. The shaded end with vertices A, B,C and D corresponds to the area, dS = J(s, θ0)|φ0=0 dθdφ, covered by the ray as the take-offangles (θ0 and φ = 0) are perturbed by dθ and dφ respectively. The arc from the origin topoint A is the original unperturbed ray.xyzθ0 s0dφdψ0s0 sin(θ0)s0 sin(θ0)dφOriginal In-plane RayTilted RayFigure 3.2: Calculation of the angle, dψ0, between the tilted ray segment, with take-off anglesθ0 and dφ, and its in-plane projection with length s0. The angle dφ is assumed small.Similarly, the 2D ray Jacobian relates the length of the arc spanned by the end of the rayto a perturbation of the take-off angle, dθ:dS2D = J2Ddθ. (3.3)In Figure 3.1, this corresponds to the length of the line segment AD.553.2. Graphical Derivation of J⊥To find the relationship between the 2D and 3D Jacobians, it is useful to determine theshape of the surface element dS. This can be accomplished by applying Snell’s law to themodel. If we approximate the true continuum of velocity with a fine grid of rectangular prisms,each with constant velocity, then we can apply Snell’s law directly. The prisms finely discretizethe model in the x and z directions (they break up the y = 0 plane into small rectangles) andare infinitely long in the y direction, thus maintaining a 2.5D representation of the model. Thesurface normals of these prisms have no y component because there is no material variation inthe y direction. We can apply the vector form of Snell’s Law,tˆ1 × nˆv1= tˆ2 × nˆv2, (3.4)to any interface in this discretization. In equation (3.4) the velocities on either side of theinterface are v1 and v2, the unit tangent vectors to the ray are tˆ1 and tˆ2, and the surfacenormal is nˆ. The components of tˆ1 or tˆ2 are related to the spherical coordinates θ and φ bytx = sin(θ) cos(φ), (3.5)ty = sin(θ) sin(φ), (3.6)andtz = cos(θ). (3.7)Equating the y components of the cross products for small azimuthal angles dφ1 and dφ2, wehavenx cos(θ1)− nz sin(θ1) cos(dφ1)v1= nx cos(θ2)− nz sin(θ2) cos(dφ2)v2. (3.8)The angles dφ1 and dφ2 are infinitesimal because the ray is either in-plane or only perturbedslightly to be out of plane, thus allowing us to simplify equation (3.8):nx cos(θ1)− nz sin(θ1)v1= nx cos(θ2)− nz sin(θ2)v2+O(dφ21) +O(dφ22). (3.9)563.2. Graphical Derivation of J⊥In the limit as dφ approaches zero, the error terms in equation (3.9) vanish. For small dφ, theprojection of the ray onto the plane is governed by equation (3.9). The projection is independentof dφ because tx and tz only depend on θ. As dφ is increased from zero, the ray must thereforelift vertically out of the plane, as indicated by ty ≈ sin(θ)dφ. So the surface element dS inFigure 3.1 is indeed a rectangle, and the element area is equal to the product of the lengths ofthe sides. If there had been a dependency on dφ in equation (3.9), then a change in dφ wouldcause the end of the ray to move at an angle out of the plane, and dS would no longer berectangular.The rectangular nature of dS allows us to factor the 3D Jacobian and writeJ3D = J‖J⊥, (3.10)andJ2D = J‖, (3.11)where J‖ is the in-plane (parallel) portion of the 3D Jacobian, and J⊥ is the out-of-plane(perpendicular) portion of the 3D Jacobian. Graphically, these components are (referring toFig. 3.1),J‖ = Length ADdθ =Length BCdθ , (3.12)andJ⊥ = Length ABdφ =Length CDdφ . (3.13)Thus, the ratio of Jacobians in equation (3.1) simplifies to J⊥.3.2.2 Derivation of J⊥ for a Straight RayThe perpendicular portion of the ray Jacobian, J⊥, can be determined geometrically by relatingthe length of the segment AB to dφ (Fig. 3.1). For simplicity, we depart from the prismaticdiscretization that we used to demonstrate the factorizability of J3D (equation (3.10)) andcompute J⊥ for a straight ray (at an arbitrary inclination θ0). To ensure that the unperturbedin-plane ray is straight, the material must consist of slabs with surface normals parallel to theray, as shown in Figure 3.3. In the “Extension to Curved Rays” section, we show that this573.2. Graphical Derivation of J⊥result also applies to curved rays through an arbitrary 2.5D medium.xzRaySlabss0s1s2s3v0v1 v2 v3sFigure 3.3: A ray normally incident on a layered structure. Each layer consists of a constantvelocity vi and is of thickness si. The total length of the ray is s.A profile of the ray as it travels through layers of varying velocity is illustrated in Figure 3.4.The in-plane ray is initially straight (the horizontal axis in the figure) and at normal incidenceto each layer. It is then tilted up out of the plane by a small angle, dψ0. The gray section inFigure 3.2 corresponds to the far left layer, v0, of Figure 3.4.v0 v1 v2OutofPlaneMovementArc Length Along Raydψ0dψ1dψ2ss0s1s2h0h1h2Figure 3.4: Out-of-plane spreading, hi, along a ray that is initially straight and normallyincident to a layered structure with constant velocities, vi. The ray is tilted up by a smallangle, dψ0. Each layer has a thickness, si.From the definition of spherical coordinates, perturbing a point (s0, θ0, 0) by dφ to (s0, θ0, dφ)results in a displacement dy = s0 sin(θ0)dφ out of the plane. The displacement is related toan effective angle, dψ0, by dy = s0dψ0. The first ray segment with inclination θ0 and azimuthperturbation dφ makes an angle dψ0 with the y = 0 plane (refer to Fig. 3.2):dψ0 = sin(θ0)dφ. (3.14)This is recognizable as part of the canonical solid-angle dΩ = sin(θ)dθdφ and is the root583.2. Graphical Derivation of J⊥of the difference between the solid-angle Jacobian, JΩ3D, and the Jacobian, J3D, related toperturbations of the spherical coordinates θ and φ.As the ray travels through the layers of material with velocity vi, the angle, dψi, that eachsegment, i, makes with the y = 0 plane changes according to Snell’s law,sin(dψi)vi= sin(dψi+1)vi+1, (3.15)ordψivi≈ dψi+1vi+1, (3.16)for small angles. This is illustrated in Figure 3.4.The elevation of the ray off of the plane increases by an amount hi as it travels througheach layer. For each segment, this elevation is related to the angle dψi byhi = si tan(dψi) ≈ sidψi, (3.17)where si is the length of the original ray segment.Combining equations (3.17) and (3.16) and summing gives the total elevation, h, of the ray(segment AB in Fig. 3.1):h =N∑i=0hi ≈dψ0v0N∑i=0sivi ≈sin(θ0)dφv0N∑i=0sivi. (3.18)Taking the limit as the segment lengths and the azimuth angle approach zero gives the result(see equation (3.13))J⊥ = ∂y∂φ∣∣∣∣φ=0= sin(θ0)v0∫ s0v(s)ds, (3.19)where the integral is along the ray from the source location to the receiver location.3.2.3 Extension to Curved RaysThe derivation of J⊥ in the previous section that resulted in equation (3.19) assumes a straightray at normal incidence to a layered material structure. To extend the result to curved rays,it must be shown that equation (3.16) is still valid if the x and z components of the surface593.3. Paraxial Computation of J⊥normals of the velocity layers are allowed to vary arbitrarily.Equating the x or z components of equation (3.4) givessin(θ1) sin(dφ1)v1= sin(θ2) sin(dφ2)v2. (3.20)For small dφ, this simplifies tosin(θ1)dφ1v1= sin(θ2)dφ2v2. (3.21)Substituting equation (3.14) into equation (3.21) again produces equation (3.16). Therefore,the angle (in the plane) at which a ray hits a velocity interface does not affect how far the rayrises out of the plane (for small dφ). Thus, the derivation that resulted in equation (3.19) alsoapplies to curved rays in a 2.5D medium.3.3 Paraxial Computation of J⊥In this section we use the paraxial ray method to analytically derive J⊥. The approach followsthat of Virieux and Farra [109], but extra attention is payed to the perturbations of the raytake-off parameters so that the ray Jacobian can be computed precisely. The notation largelymirrors that of Virieux and Farra [109] and Farra and Madariaga [38].We start by defining a Hamiltonian,H(x,p) = 12[p · p− u2(x)], (3.22)that is constructed directly from the eikonal equation,|∇T |2 = 1v2 = u2, (3.23)where u is the slowness of the medium (v is velocity) and p = ∇T . The vector p is also referredto as slowness. Governing differential equations for rays can be derived from equation (3.22)using the method of characteristics and can be found in [109]. The rays are parameterizedcurves in six-dimensional phase space described by the vector yT0 (σ) = [x0(σ),p0(σ)], where603.3. Paraxial Computation of J⊥x0 is the parameterized position and the subscript indicates that this is a reference ray. Thereare a several common options for parameterizing the position on the ray [22, p. 126]. Here wechoose σ (τ in [109]), which is defined bydT = p · dx = u2(x)dσ. (3.24)Paraxial rays are perturbations of the reference ray and are defined by x(σ) = x0(σ)+δx(σ)and p(σ) = p0(σ) + δp(σ). The perturbations along the central ray, δyT = (δx, δp), aregoverned by another set of differential equations [109],dδydσ = Aδy, (3.25)whereA =∇x∇pH ∇p∇pH−∇x∇xH −∇p∇xH. (3.26)Additionally, the paraxial rays must satisfyδH = ∇pH · δp+∇xH · δx = 0, (3.27)otherwise the computed trajectories are not actually rays. If this condition is satisfied anywhereon a paraxial ray, then it is satisfied everywhere on that ray.To compute the ray Jacobian, we start by perturbing the slowness, p, at the start of thecentral ray to find two new paraxial rays. These two new rays define the cross-sectional area ofthe ray tube. In order to find a Jacobian per solid angle, we must ensure that our perturbationstrace out a known portion of a solid angle.To satisfy equation (3.27) at the start of the paraxial ray, we must ensure that the pertur-bations of initial slowness for the two paraxial rays, δpi and δp′i, are perpendicular to the initialslowness, pi. We are not perturbing the starting position of the ray, only the take-off angle,so δxi = 0. Further, if we ensure that each perturbation is of length ǫ|pi|, where ǫ is a smallpositive number and the perturbations are perpendicular to each other, then they trace out anarea of ǫ2|pi|2 on a sphere of radius |pi|. This corresponds to a solid angle of ǫ2 steradians.613.3. Paraxial Computation of J⊥Requiring that the perturbations are perpendicular to each other and are of the same length isnot essential, but doing so simplifies the computation of the solid angle spanned by the rays.If we define the position perturbations associated with the slowness perturbations δpi and δp′ito be δx and δx′, respectively, then we can express the ray Jacobian per solid angle JΩ3D as,ǫ2JΩ3D = pˆ · (δx× δx′), (3.28)where we have assumed that δpi × δp′i is in the direction of pi rather than −pi. This ensuresthat the sign of the Jacobian is compatible with the geometrical derivation (i.e., positive for auniform medium).A set of initial perturbations that satisfies the above requirements isδpi = ǫ|pi|√p2zi + p2yi(0,−pzi, pyi) (3.29)andδp′i = ǫ|pi|√(p2yi + p2zi)2 + p2xi(p2yi + p2zi)(−p2yi − p2zi, pxipyi, pxipzi). (3.30)Following Virieux and Farra [109], we choose three reference paraxial trajectories with unitperturbations of initial slowness in each coordinate direction:δp1 = (1, 0, 0), (3.31)δp2 = (0, 1, 0), (3.32)andδp3 = (0, 0, 1). (3.33)The corresponding perturbations in position along each trajectory are δq1, δq2, and δq3. Ingeneral, the trajectories are not rays because equation (3.27) is not satisfied. However, theperturbations in equation (3.28) for the two paraxial rays, δx and δx′, can be constructed usinglinear combinations of these three trajectories. This allows equation (3.28) to be expressed as623.3. Paraxial Computation of J⊥a determinant [109]:JΩ3D =−|pi||p|∣∣∣∣∣∣∣∣∣∣∣∣∣px δq1x δq2x δq3xpy δq1y δq2y δq3ypz δq1z δq2z δq3z0 pxi pyi pzi∣∣∣∣∣∣∣∣∣∣∣∣∣. (3.34)Using analogous reasoning, we can express the 2D ray Jacobian asJ2D =1|p|∣∣∣∣∣∣∣∣∣∣px δq1x δq3xpz δq1z δq3z0 pxi pzi∣∣∣∣∣∣∣∣∣∣, (3.35)where p has no y-component.In the 2.5D case the material properties do not vary in the y direction. The paraxial raytracing equation (3.25) can then be simplified and expressed asddσδqxδqyδqzδpxδpyδpz=0 0 0 1 0 00 0 0 0 1 00 0 0 0 0 112∂2u2∂x2 012∂2u2∂x∂z 0 0 00 0 0 0 0 012∂2u2∂x∂z 0 12 ∂2u2∂z2 0 0 0δqxδqyδqzδpxδpyδpz. (3.36)Inspection of equation (3.36) indicates that the y component of the slowness vector is aconstant along the ray:dδpydσ = 0. (3.37)Additionally, from equation (3.36), we can relate the y component of the position pertur-bation to the perturbation of initial slowness by solvingdδqydσ = δpy (3.38)633.3. Paraxial Computation of J⊥to obtainδqy = δpyσ. (3.39)From equation (3.39), we see that δq2y = σ, and δq1y = δq3y = 0, thus allowing us to simplify thesecond row of equation (3.34). Similarly, δq2x = δq2z = 0. The 3D Jacobian (equation (3.34))can now be simplified to:JΩ3D =−|pi||p|∣∣∣∣∣∣∣∣∣∣∣∣∣px δq1x 0 δq3xpy 0 σ 0pz δq1z 0 δq3z0 pxi pyi pzi∣∣∣∣∣∣∣∣∣∣∣∣∣. (3.40)This can be manipulated and expressed in a form that contains the determinant of the 2D rayJacobian (equation (3.35)),JΩ3D = σ|pi||p| Det2D − p2yi|pi||p|∣∣∣∣∣∣∣δq1x δq3xδq1z δq3z∣∣∣∣∣∣∣, (3.41)where Det2D is the determinant in equation (3.35), and we have made use of equation (3.37).The vector p in equation (3.35) does not have a y component, so it would not be correct toinclude J2D in equation (3.41). However, if we set pyi = 0, as is the case for planar rays, wecan writeJΩ3D = σ|pi|J2D. (3.42)The 3D Jacobian can be factored into a parallel and perpendicular part. The parallel part issimply the 2D Jacobian, and the perpendicular part isJΩ⊥ = σ|pi| =1v0∫ s0v(s)ds, (3.43)where s is arc-length along the ray and v0 is the velocity at the start of the ray. The definitionof the ray parameter (equation (3.24)) has been used to express σ as an integral, and the eikonalequation (3.23) gives us |pi| = 1/v0.We can use the solid-angle relationship dΩ = sin(θ)dθdφ to convert the solid-angle Jacobian643.4. Simplified Transfer FunctionJΩ3D into J3D that operates on perturbations of the spherical coordinates dθ and dφ,J3D = sin(θ0)JΩ3D, (3.44)orJ⊥ = sin(θ0)v0∫ s0v(s)ds, (3.45)which is identical to equation (3.19), thus verifying the geometrical derivation.3.4 Simplified Transfer FunctionIncorporating equations (3.10), (3.11), and (3.19) or (3.43) into equation (3.1) results in thesimplified expressionTF = −iπH(2)0 [ωT (s)] eiωT (s)√T (s)∫ s0v(s)ds. (3.46)If we replace the Hankel function in equation (3.46) with its large-argument expansion [1, p.364], we obtainTF ≈ e−ipi4√1f∫ s0v(s)ds. (3.47)Although equation (3.46) can be used, equation (3.47) provides accurate results when thesource and receiver are sufficiently far apart.3.5 ConclusionThe perpendicular portion, J⊥, of the ray Jacobian, J , was derived for the 2.5D planar ray casethrough a graphical approach and confirmed by an analytical approach. Both approaches arereasonably straightforward and produce the same result, but the graphical approach requiresonly a knowledge of Snell’s law. It also provides a more physically intuitive view of ray spreadingout of the plane and the factorizability of the ray Jacobian.653.6. Data and Resources3.6 Data and ResourcesNo data were collected for or used in this article. The article was prepared in LATEX, and thefigures were created with Inkscape. Maple was used to verify some of the equations.3.7 AcknowledgementsThe authors would like to thank the National Science and Engineering Research Council ofCanada for its direct financial support. This work has also been partially supported by theFrench Agence Nationale de la Recherche under grant ANR-2011-BS56-017.66Chapter 4Three-Dimensional toTwo-Dimensional Data Conversionfor Electromagnetic WavePropagation Using an AcousticTransfer Function: Application tocross-hole GPR data4.1 IntroductionFull-wave inversion techniques are becoming increasingly popular for acoustic [76, 125, 68],seismic [82, 81, 88, 110], and electromagnetic [33, 35, 67, 66, 57] imaging problems. However,the memory and computational requirements often restrict its use to a 2-D domain. One of twoprocedures can be used to overcome the fundamental difference between real-world (3-D) dataand restriction of the computation to two dimensions. The first is to compute 3-D accuratefields during the inversion. The second is to first modify the real-world data so that it matcheswhat would be produced by a 2-D simulation.In this work, we investigate the second of the two options and focus our efforts on anacoustic 3-D to 2-D transformation derived by Yedlin et al. [121]. Bleistein [19] arrives at the2This article has been accepted for publication in Geophysical Journal International c©: 2014 The Authors.Published by Oxford University Press on behalf of The Royal Astronomical Society. All rights reserved.674.1. IntroductionxzyBoreholesTransmitter(infinitesimal dipole)Receiver Locations(Only a few are shown)Material Boundary(No variation in y)xzDepth (-z)Depth (-z)V/mV/m3D2DReceived Signal3Dto2DSurfaceDepth of SurveyTransmitter becomesa ribbon of current withinfinite y-dimension andinfinitesimal z-dimensionTransfer Function:ℜ{Ez}Figure 4.1: Schematic representation of the physical problem and action of the transfer function.The transfer function converts measured 3-D GPR data into what would be generated in anequivalent 2-D problem. The material is allowed to vary arbitrarily in x and z, and is assumednot to vary in y.same acoustic transfer function but using standard asymptotic ray theory and source matchingtechniques rather than directly taking the ratio of 3-D and 2-D asymptotic Green’s functions.A thorough discussion on various 3-D to 2-D transformation techniques is presented by Aueret al. [6] and by Williamson and Pratt [118]. The utility of the method that we focus onin this work, but for inversion of acoustic and seismic data, is demonstrated for a syntheticcase by Auer et al. [6]. Similar methods have been used successfully for real cross-well groundpenetrating radar (GPR) surveys by Klotzsche et al. [57] and Ernst et al. [33] (reading thelatter prompted our investigation into these transformations). Although an acoustic transferfunction has been successfully used for vector electromagnetic GPR applications as mentionedabove, we show analytically and numerically that this is indeed justified for a cross-well setup.The analysis assumes that the material does not vary in the y-direction (a 2.5-D structure) andthat the material is lossless and isotropic. A schematic of the physical problem and the actionof the transfer function are illustrated in Fig. 4.1.684.2. Borehole GPR as an Acoustic ProblemThese transfer functions are normally evaluated using first arrivals. A drawback of usinga transfer function that is based on first-arrivals only is that when there are multiple arrivalsthe approximation may be poor. It is possible to apply this type of transformation to multiplearrivals as long as the corresponding rays can be found, but the arrivals must be separatedin time and that may not always be the case [6]. Despite this limitation, this type of trans-fer function (evaluated using only first arrivals) is used successfully in seismic inversion andcross-hole GPR inversion. We provide additional synthetic inversion results in this paper thatdemonstrate the effectiveness of the transfer function for frequency-domain cross-hole GPR in-version where permittivity and conductivity are simultaneously recovered. When the transferfunction is not used, the accuracy of the recovered conductivity is adversely affected while therecovered permittivity remains essentially unaffected (at least if the conductivity is small). Thisis because there is a strong link between conductivity and amplitude of a received signal andthe transfer function addresses the change in amplitude in 3-D versus 2-D worlds.4.2 Borehole GPR as an Acoustic ProblemIn this section, we show that the 2-D and 3-D vector electromagnetic GPR problems can beapproximated as acoustic problems. Once this is accomplished, we show how the acoustictransfer function of Yedlin et al. [121] can be used for the electromagnetic case.For both the acoustic and GPR cases it is assumed that∣∣∣∣∇vv∣∣∣∣≪ 2πλmin, (4.1)so that the Wentzel-Kramers-Brillouin (WKB) approximation holds [98]. In eq. (4.1), λmin =2πvmin/ω, and vmin is the slowest velocity found in the model. Eq. (4.1) states that thelength scale of the wavelength is much smaller than the length scale of material propertyvariation. When eq. (4.1) is satisfied the waves behave locally like plane waves and we can usean asymptotic series in powers of ω−n (with n = 0 . . .∞) to describe the solution [25, pp. 135,278]. The use of such an expansion allows us to simplify the electromagnetic equations (at theexpense of accuracy) and put them into a form compatible with an easily computed transfer694.2. Borehole GPR as an Acoustic Problemfunction.The governing equations for the borehole GPR problem are those of Maxwell,∇×E = −jωµ0H, (4.2)∇×H = jωǫE+ Js, (4.3)∇ · µ0H = 0, and (4.4)∇ · ǫE = ρ, (4.5)where Js is a driving current, E is the electric field, and H is the magnetic field. The magneticpermeability, µ0, is assumed to be constant and equal to the free space value. The electricpermittivity, ǫ, is allowed to vary with position. The free charge density, ρ, is assumed to bezero everywhere except at the source (in order to support the existence of the source current).It is also assumed that the material is lossless and isotropic. Although this derivation assumesno loss, a numerical experiment shows that the results are also applicable to the low-loss case.A synthetic low-to-moderate-loss example is given in Section 4.3.3, and a synthetic inversionexample with low loss is given in Section 4.4.The borehole GPR experiment consists of a pair of transmit and receive antennas. Withthe assumption that the current distribution on the receiving antenna in transmit mode isknown, reciprocity allows one to easily compute the voltage at the terminals of the antennawhen receiving in terms of the incident electric field [55, p. 353]. In the case of a shortreceiving antenna, the output voltage is the vertical component of the electric field incident onthe antenna, scaled appropriately. To simplify the work that follows, electrically short dipoleantennas are assumed. Here, we assume that the antennas are oriented vertically (in the zˆdirection), and the effects of the borehole itself are ignored. Due to the orientation of theantennas relative to the 2-D plane between them, computing E is necessarily a vector problembecause it will generally have non-zero Ex and Ez. So we start our investigation by findingequations for the magnetic field for both the 3-D and the 2-D cases.If we apply the curl operation to Ampere’s law, eq. (4.3), and expand ∇× ǫE, we obtain∇(∇ ·H)−∇2H = jω∇ǫ×E+ ω2ǫµ0H+∇× Js, (4.6)704.2. Borehole GPR as an Acoustic Problemwhere eq. (4.2) has been substituted for ∇×E. Assuming that eq. (4.1) holds allows us to usean asymptotic series for the field quantities. Doing so, and remembering that we are workingat a high, but finite, frequency allows us to neglect the term involving ∇ǫ to get a zeroth-orderapproximation. The divergence term vanishes due to eq. (4.4). The result is three independentscalar equations, one for each component of H,∇2H+ k2H ≈ −∇× Js, (4.7)where k = ω√µ0ǫ is the wavenumber.A transmitting antenna in a borehole GPR setup can be roughly approximated as an in-finitesimal vertical current element,Js = zˆδ(x− x0) [A/m2 or A/m], (4.8)where x0 = (x0, 0, z0). Here, we assume that the 2-D plane of interest is where y = 0 and thatthe transmit antenna is oriented in the zˆ direction. Note that in 2-D Js is a ribbon of current.The forcing term for eq. (4.7) is the curl of the source current. In three dimensions this is∇× Js = xˆδ(x− x0)δ′(y)δ(z − z0)− yˆδ′(x− x0)δ(y)δ(z − z0) [A/m3]. (4.9)The 3-D magnetic field is thus described by∇2Hx + k2Hx ≈ −δ(x− x0)δ′(y)δ(z − z0), (4.10)∇2Hy + k2Hy ≈ δ′(x− x0)δ(y)δ(z − z0), and (4.11)Hz ≈ 0. (4.12)In two dimensions the forcing term is∇× Js = −yˆδ′(x− x0)δ(z − z0) [A/m2]. (4.13)714.2. Borehole GPR as an Acoustic ProblemThis source term allows us to find a scalar equation for Hy (with Hx = Hz = 0) from eq. (4.7):∇2Hy + k2Hy ≈ δ′(x− x0)δ(z − z0) (4.14)Above, we have found approximate scalar equations for each component of the magneticfield from a GPR transmitting antenna, but what we are really interested in is how the 3-Delectric field relates to the 2-D electric field. The electric field is proportional to the curl of Hfrom eq. (4.3),jωǫE = ∇×H = −xˆ∂Hy∂z + yˆ∂Hx∂z + zˆ(∂Hy∂x −∂Hx∂y). (4.15)With a little work we can relate these derivatives of H to the acoustic problem discussed in[121],∇2P + k2P = δ(x− x0)δ(y)δ(z − z0), (4.16)as follows. Eq. (4.16) differs from eqs. (4.10) and (4.11) only in the source term. If we differen-tiate both sides of eq. (4.16), we obtain∇2 ∂∂xP +( ∂∂xk2)P + k2 ∂∂xP = δ′(x− x0)δ(y)δ(z − z0), (4.17)which has a comparable source term to eq. (4.11) and a similar structure to eq. (4.6), but herewe wish to relate solutions of the electromagnetic problem directly to solutions of the acousticproblem of eq. (4.16) rather than solve a new differential equation. We can again use the WKBapproximation and neglect the derivative of k2, allowing us to write∇2 ∂∂xP + k2 ∂∂xP ≈ δ′(x− x0)δ(y)δ(z − z0). (4.18)By inspection it is clear that solutions to eq. (4.11) are derivatives of the solutions to eq. (4.16),so ∂P/∂x ≈ Hy as long as the velocity varies slowly enough in space.Following a similar approach by differentiating with respect to y, we obtain −∂P/∂y ≈ Hx.In this case, no approximation needs to be made when differentiating eq. (4.16) because k doesnot depend on y, but overall there is an approximation due to eq. (4.7).Due to the symmetry of the material and of the doublet source in eq. (4.10), the x-component724.2. Borehole GPR as an Acoustic Problemof H is zero in the y = 0 plane. Clearly, if it is zero everywhere in the plane then ∂Hx/∂zmust also be zero within the plane, and from eq. (4.15) we see that Ey = 0. The ∂Hx/∂yterm is more complicated. It is equivalent to differentiating the solution, P , of eq. (4.16) twicewith respect to y. We can show that this term can be neglected by differentiating the uniformasymptotic ansatz discussed in [121],P (x) = − e−iωT (x)4πv(x)T (x)∞∑n=0An(x)(iω)n . (4.19)In eq. (4.19), the velocity, v(x), and An(x) are smooth functions independent of ω. Thederivative of traveltime, T (x), with respect to y in the plane is zero due to the planar raysand the symmetry of the velocity profile. As a result, the ω2 term that appears when applying∂2/∂y2 to the exponential in eq. (4.19) disappears because ∂T (x)/∂y = 0. Therefore, we canneglect ∂Hx/∂y in eq. (4.15). The ω2 terms that appear when computing the derivatives ofHy in eq. (4.15) in a similar manner do not disappear because, in general, ∂T (x)/∂x 6= 0 and∂T (x)/∂z 6= 0. What we are left with is,∇×H = −xˆ∂Hy∂z + zˆ∂Hy∂x , (4.20)which only involves derivatives of Hy tangent to the plane. Knowledge of the behaviour of Hyoutside of the plane is, therefore, not required to compute E. Combining eqs. (4.20), (4.18),and (4.3), away from the source, we find a differential operator that when applied to solutionsof eq. (4.16) will yield the electric field:E ≈ 1jωǫ(zˆ∂2∂x2 − xˆ∂2∂z∂x)P. (4.21)The first order asymptotic solutions to eq. (4.16) for the 2-D and 3-D cases derived in [121]areP2D(s) ≈14√2v(s)πωJ2D(s)ej 3pi4 e−jωT (s), (4.22)734.3. Numerical Experimentswhere the Hankel function has been replaced by its asymptotic form, andP3D(s) ≈ −14π√sin(θ0)v(s)J3D(s)v(0)e−jωT (s). (4.23)In eqs. (4.22) and (4.23) the parameter s is the arc length along a given ray. The Jacobians,J2D and J3D, are with respect to take-off angles, v(s) is the velocity along the ray, v(0) is thevelocity at the source and θ0 is the take-off inclination. The sin(θ0) term appears because J3Dis defined for the spherical coordinate system azimuth and inclination rather than for a solidangle.When we apply the operator in eq. (4.21) to eq. (4.22) or eq. (4.23), we can again applyeq. (4.1) to drop all terms except the ones involving derivatives of the exponentials becauseof their dependence on ω. The two exponential terms are identical, so when the ratio of the3-D to the 2-D fields is taken, the derivatives cancel out and we are left with the same transferfunction as for the acoustic case:TF = E(2D)zE(3D)z≈ E(2D)xE(3D)x≈ e−j pi4√2πJ3D(s)v(s0)ωJ2D(s) sin(θ0). (4.24)This is only valid in the far field, and many approximations were made to arrive at it, but it isadequate justification for applying the acoustic transfer function to the electromagnetic case.The ratio of Jacobians in eq. (4.24) can be computed by integrating velocity along the ray path[100], resulting inTF ≈ e−ipi4√1f∫ st0v(s)ds = e−ipi4√1f∫ st0ds√ǫ(s)µ, (4.25)where st is the total arc-length of the ray.4.3 Numerical ExperimentsIn this section, we demonstrate the effectiveness of the transfer function that was arrived at ineq. (4.25) for the electromagnetic case of cross-borehole GPR by applying it to three synthetictest models. For comparison, we also apply the method of Ernst et al. [33] (we use the complex744.3. Numerical Experimentsconjugate of their eq. A-1 to match our time convention), repeated here for reference:TF =√2πTiωǫmeanµ = e−ipi4√Tfǫmeanµ, (4.26)where T is the traveltime from source to receiver and ǫmean is the average permittivity in thedomain. In the experiments below we refer to this as ‘quasi-straight’ because some effect ofray curvature is captured through the traveltime, T , even though a constant permittivity isused. Eq. (4.25) is referred to as ‘curved’ because of the direct inclusion of ray bending. Thefast-sweeping eikonal solver of Zhao [124] is used to compute traveltimes. The ray paths foreq. (4.25) are found by following the traveltime gradient from receiver to source.In the tests, the computed fields from a full-wave 3-D finite-difference time-domain (FDTD)electromagnetic simulation using MEEP [73] are converted to their 2-D equivalent using thetransfer functions. In this section, the transfer functions are evaluated on the true test mod-els and not approximate ones arrived at through a ray inversion. The converted fields arecompared with 2-D full-wave electromagnetic simulations that are generated using a customfinite-difference frequency-domain (FDFD) code. The computational domain is surrounded bya perfectly matched layer (PML) in both cases.The 2-D material variation in the xz plane is effectively extruded in the y-direction forthe 3-D simulations, resulting in a 2.5-D model. The experiments are performed using a shortvertical dipole current element as the source on the left side of the model. The received signalis taken to be the vertical component of the electric field at the right-hand side of the model, asa function of depth. The material in each model is non-magnetic, isotropic and lossless (exceptfor the last case).The FDFD grid used for the 2-D simulations is set to 60 points per shortest wavelength of thebackground medium [see eq. (4.27)]. The grid size for the 3-D simulations is set to 40 pointsper shortest wavelength of the background medium. The frequency used for all simulationsis 94.8 MHz (this was arrived at to match wavelengths with previous acoustic experiments[121] given the chosen electric permittivities). The 3-D results were verified using the finite-element method software package COMSOL, but at a substantially coarser mesh due to memoryconstraints and were found to closely match.754.3. Numerical Experiments02468100 2 4 6 8 101.61.822.22.42.62.83RelativeSlownessx (m)Depth(m)Depth(m)Receiver LocationsTransmitterX(a)02468100 2 4 6 8 10-15-10-5051015mV/mx (m)Depth(m)Receiver LocationsTransmitterX(b)Figure 4.2: Computed rays and 2-D fields at 94.8 MHz for a 10 m × 10 m area with varyingslowness: ǫr = 10− 0.8z: (a) A selection of the computed rays plotted on the relative materialslowness, √ǫr; (b) Real part of the electric field. The transmitter and receiver locations arealso shown.4.3.1 Gradient of PermittivityThe first model has a relative permittivity that depends on depth asǫr(z) = 10− 0.8z. (4.27)The 2-D rays and electric field throughout the domain are shown in Fig. 4.2.Fig. 4.3 compares the transformed 3-D data (using the transfer functions) to the 2-D data.The data are collected along the right-hand side of the model as if from a borehole antennamoved from top to bottom of the hole. See the illustrations in Figs. 4.1 and 4.2 for clarificationof experimental setup and data representation. It can be seen that excellent agreement isachieved with the curved-ray transfer function of eq. (4.25), with an overall relative error of2.0 per cent3. The quasi-straight-ray transfer function of eq. (4.26) also provides good results,but not quite as accurate with a relative error of 14.8 per cent. A small phase error for bothtransfer functions (both functions have a constant −45◦ phase angle) and amplitude discrepancyfor the curved-ray transfer function, visible in Fig. 4.3c, is likely due to the relatively strong3The relative error is the rms error of the transfer function,√1N ‖TF − TFtrue‖2, divided by the rms valueof the true transfer function,√1N ‖TFtrue‖2. Here, the transfer functions are taken to be vectors with N datapoints (one for each receiver location). TF is given by eq. (4.25) and TFtrue is the ratio between the computed3-D and 2-D fields. Phase errors increase the rms error and so are included, at least to some degree, in thiscalculation.764.3. Numerical Experiments-0.01-0.005 0 0.005 0 1 2 3 4 5 6 7 8 9 10ℜ{2D FDFD}Curved ℜ{3D → 2D}|2D FDFD|Curved |3D → 2D|Quasi-straight ℜ{3D → 2D}Quasi-straight |3D → 2D|Depth (m)ElectricField(V/m)(a)-200 0 200 400 600 800 1000 1200 0 1 2 3 4 5 6 7 8 9 106 2D FDFD6 3D → 2DDepth (m)Phase(Deg)(b)-46-45.5-45-44.5-44-43.5 0 1 2 3 4 5 6 7 8 9 10 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 0 1 2 3 4 5 6 7 8 9 106 2D/3D6 TF|2D/3D|Curved |TF|Quasi-straight |TF|Depth (m)Depth (m)Phase(Deg)Magnitude(c)Figure 4.3: Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions (in the medium corresponding to Fig. 4.2), measured at the receiving antennaas a function of depth (f = 94.8 MHz): (a) Electric field magnitude and real part; (b) Electricfield phase (unwrapped for clarity); (c) Phase and magnitude of the transfer functions comparedwith phase and magnitude of the ratio of the computed 2-D and 3-D fields. The phase of thetwo transfer functions is fixed at −45◦, so only one curve is shown for phase. The two transferfunctions are indicated by the words ‘Curved’ for eq. (4.25) and ‘Quasi-straight’ for eq. (4.26).The thick grey curves represent the target 2-D FDFD and transfer function values.774.3. Numerical Experimentspermittivity gradient that does not well satisfy eq. (4.1). For the very deepest receivers, around9 m and below, a greater discrepancy is clear in the magnitude plot of Fig. 4.3c. This is likelydue to a combination of increased velocity gradient at depth, the artificial confinement of therays within the domain that is visible at the bottom of Fig. 4.2a, and possibly interaction withthe PML layers where the permittivity gradient stops in the 2-D case (due to implementation)and a non-physical attenuating action exists. The somewhat larger amplitude discrepancy forthe quasi-straight-ray transfer function that is visible in Fig. 4.3c likely arises from the largevariation of permittivity from the top to the bottom of the domain and the use of the averagepermittivity instead of a permittivity that is more localized to the rays.4.3.2 Dielectric LensThe second test is performed using the same background permittivity as the first case, but withan added lens:ǫr(x, z) =[ 1√10− 0.8z− 0.2e−((x−3)2+(z−6)2)]−2. (4.28)Plots of the computed rays, electric field and intensity are shown in Fig. 4.4. The plot ofintensity, Fig. 4.4c, is included to illustrate the effect of the lens. The focusing action of thelens indicates that a ray caustic is present.Fig. 4.5 again compares the transformed 3-D data to the 2-D data. The agreement is not asclose as in the first test case, but is still good with a relative error of 12.2 per cent for the curved-ray transfer function and 13.1 per cent for the quasi-straight-ray one. There is a discontinuityin the magnitude of the curved-ray transfer function (Fig. 4.5c) near a depth of 5 m. Thisdiscontinuity is a result of the rays taking a greatly different path around the anomaly, as isclearly visible in Fig. 4.4a. The abrupt change in path results in an abrupt change in the resultof the integration in eq. (4.25). The phase of the transfer functions, on the other hand, is quiteaccurate except where the signal is very weak. This region is at a depth near 3 m and is markedin Fig. 4.5 with vertical dotted lines.784.3. Numerical Experiments02468100 2 4 6 8 101.522.533.54RelativeSlownessx (m)Depth(m)Depth(m)(a)02468100 2 4 6 8 10-15-10-5051015mV/mx (m)Depth(m)(b)02468100 2 4 6 8 1001e-062e-063e-064e-065e-06W/m2x (m)Depth(m)(c)Figure 4.4: Computed rays and 2-D fields at 94.8 MHz for a 10 m × 10 m area with a varyingslowness and small lens: ǫr =[1√10−0.8z − 0.2e−((x−3)2+(z−6)2)]−2: (a) A selection of the com-puted rays plotted on the relative material slowness, √ǫr; (b) Real part of the electric field; (c)Magnitude of the Poynting vector.794.3. Numerical Experiments-0.015-0.01-0.005 0 0.005 0.01 0 1 2 3 4 5 6 7 8 9 10ℜ{2D FDFD}Curved ℜ{3D → 2D}|2D FDFD|Curved |3D → 2D|Quasi-straight |3D → 2D|Quasi-straight ℜ{3D → 2D}Depth (m)ElectricField(V/m)(a)-200 0 200 400 600 800 1000 1200 0 1 2 3 4 5 6 7 8 9 106 2D FDFD6 3D → 2DDepth (m)Phase(Deg)(b)-60-55-50-45-40-35-30 0 1 2 3 4 5 6 7 8 9 10 3 3.5 4 4.5 5 5.5 6 0 1 2 3 4 5 6 7 8 9 106 2D/3D6 TF|2D/3D|Curved |TF|Quasi-straight |TF|Depth (m)Depth (m)Phase(Deg)Magnitude(c)Figure 4.5: Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions (in the medium corresponding to Fig. 4.4), measured at the receiving antennaas a function of depth (f = 94.8 MHz): (a) Electric field magnitude and real part; (b) Electricfield phase (unwrapped for clarity); (c) Phase and magnitude of the transfer functions comparedwith phase and magnitude of the ratio of the computed 2-D and 3-D fields. The vertical dottedlines indicate a region of low signal strength and large transfer function error. The phase of thetwo transfer functions is fixed at −45◦, so only one curve is shown for phase. The two transferfunctions are indicated by the words ‘Curved’ for eq. (4.25) and ‘Quasi-straight’ for eq. (4.26).The thick grey curves represent the target 2-D FDFD and transfer function values.804.3. Numerical Experiments4.3.3 Gradient of Permittivity with LossThe last simulation is of a permittivity model equal to the first case [see eq. (4.27)] but with lossfrom a constant conductivity of σ = 0.005 S m−1. For this case the real part of the complexpermittivity, ǫ = ǫrǫ0 − iσ/ω, is used to compute the traveltime for eq. (4.26) and the raypaths that are necessary to evaluate the integral in eq. (4.25). The velocity in the integrand iscomplex and computed directly from the complex permittivity by v = 1/√ǫµ0. Similarly, ǫmeanin eq. (4.26) is the mean complex permittivity. The 3-D and 2-D electromagnetic fields thatare compared are again computed using FDTD and an FDFD code with the specified constantconductivity.By doing this we have neglected some of the effect of the loss by keeping the rays completelyreal. This is partly justifiable because the conductivity is small enough that at the frequency ofinterest the material behaves as a moderate-to-good dielectric with phase velocity nearly equalto the lossless case. A good dielectric is defined as having (ωǫ/σ)2 ≫ 1 [13, p. 150], whichis satisfied except for where the permittivity is very low in this model. Including the complexvelocity in the computation of the transfer function at least partly takes into account the loss,and is further justified through this numerical experiment. The comparison of transformed 3-Ddata to 2-D data in Fig. 4.6 shows that there is excellent agreement, but with an increased phaseerror compared to the lossless case and a relative error of 3.6 per cent for the curved-ray transferfunction and 13.4 per cent for the quasi-straight-ray one. Including the complex slowness in theintegrand in eq. (4.25) gives a transfer function phase that is different from the constant −45◦that is seen in the previous test cases. For both transfer functions, the result is a much closerestimate of the true phase difference between 2-D and 3-D for this low-to-moderate-loss case.A more careful computation to include absorption losses would entail tracing complex rays.Studies involving complex rays have been carried out for both electromagnetic and seismic wavepropagation in dissipative media (see, e.g. [126, 5], and the references within them).However, difficulties in the physical interpretation of complex rays have led to studies involv-ing real rays in absorbing media (see, e.g. [105]). Discussions of ray tracing in absorbing mediacan also be found in various texts, for example, Carcione [21] and Cˇerveny´ [22]. These studiesshow that if the absorption is weak, then real rays combined with complex frequency-dependent814.3. Numerical Experiments-0.0002-0.00015-0.0001-5e-05 0 5e-05 0.0001 0.00015 0 1 2 3 4 5 6 7 8 9 10ℜ{2D FDFD}Curved ℜ{3D → 2D}Quasi-straight ℜ{3D → 2D}|2D FDFD|Curved |3D → 2D|Quasi-straight |3D → 2D|Depth (m)ElectricField(V/m)(a)-200 0 200 400 600 800 1000 1200 0 1 2 3 4 5 6 7 8 9 106 2D FDFDCurved 6 3D → 2DQuasi-straight 6 3D → 2DDepth (m)Phase(Deg)(b)-43-42-41-40-39-38-37-36 0 1 2 3 4 5 6 7 8 9 10 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 0 1 2 3 4 5 6 7 8 9 106 2D/3DCurved 6 TFQuasi-straight 6 TF|2D/3D|Curved |TF|Quasi-straight |TF|Depth (m)Depth (m)Phase(Deg)Magnitude(c)Figure 4.6: Comparison of the 2-D electric field to the 3-D electric field transformed by thetransfer functions for a permittivity gradient with σ = 0.005 S m−1, measured at the receivingantenna as a function of depth (f = 94.8 MHz): (a) Electric field magnitude and real part; (b)Electric field phase (unwrapped for clarity); (c) Phase and magnitude of the transfer functionscompared with phase and magnitude of the ratio of the computed 2-D and 3-D fields. The twotransfer functions are indicated by the words ‘Curved’ for eq. (4.25) and ‘Quasi-straight’ foreq. (4.26). The thick grey curves represent the target 2-D FDFD and transfer function values.824.4. Full-wave Inversion Experimentmedium parameters (e.g. wave speeds) generally produce wave amplitudes and phases that aresufficiently accurate for most purposes.4.4 Full-wave Inversion ExperimentTo demonstrate the effectiveness of this technique to the inversion of borehole radar data, weprovide in this section a synthetic full-wave inversion example. The inversion algorithm uses a2-D FDFD solver to compute the forward model. The inversion is formulated as a regularizedminimization problem using the Gauss-Newton method [see, e.g.: 97, 82] and simultaneouslyfinds permittivity, conductivity and the source wavelet. The FDFD forward model is solvedusing an LU factorization, and this factorization is also used for the direct calculation of theJacobian matrix that relates perturbations in permittivity, conductivity and source wavelet tochanges in predicted data [see, e.g., 82]. The objective function that is minimized isφ(m) = φd + βφm=∥∥∥Wd(F[m]− dobs)∥∥∥2+ βNreg∑i=1‖Wi(m−mref )‖2 ,(4.29)where the model vector, m, contains the permittivity, conductivity and source wavelet. Theregularization operators, Wi, measure roughness and deviation from the reference model, mref ,of the permittivity and conductivity. The regularization operators ignore the source waveletcoefficients in m and mref . In eq. (4.29), the forward modeling operator F represents thepredicted data found using Maxwell’s equations (approximated with FDFD), and the observeddata are dobs. The data are weighted by Wd, β is the regularization parameter, and Nreg is thenumber of regularization matrices (three in this case: two for spatial derivatives and anotherthat weights the difference between m and mref ). The regularization is weighted such that theconductivity is allowed to vary somewhat more (in a relative sense) than permittivity.Synthetic test data are generated using the MEEP FDTD solver [73] in three dimensionson a model that varies only in two dimensions. Gaussian pulse source excitation is used with acenter frequency of 150 MHz and a full-width half-maximum of 10 ns. A total of ten infinitesimaldipole transmitter locations are used (5 on each side of the model) with finely spaced receivers834.4. Full-wave Inversion Experiment02468100 1 2 3 4 57891011121314RelativePermittivityx (m)Depth(m)Transmitter/receiverlocations(a)02468100 1 2 3 4 51234567Conductivity(mS/m)x (m)Depth(m)(b)Figure 4.7: True model: (a) Relative permittivity and (b) conductivity. The transmitter andreceiver locations are also shown.on the opposite side of the model from a given transmitter. Gaussian random noise is added tothe synthetic data such that the signal-to-noise ratio approximates that of actual radar takenfrom a karst aquifer with similar borehole spacing [101]. The FDTD grid is set to 20 points pershortest wavelength at 150 MHz. The receivers measure the z-component of the electric field.The true model is shown in Fig. 4.7.To provide a starting model for the full-wave inversion, a 2-D ray inversion is performed[101]. The recovered permittivity from the ray inversion is shown in Fig. 4.8. The regulariza-tion parameter for the starting model is chosen manually here, but could be automated withgeneralized cross-validation or a related method. The ray inversion is only able to recover thepermittivity on a very coarse scale. It is able to resolve the background permittivity half-spacesand partially resolve the two larger cylindrical anomalies. Most of the anomalies that are vis-ible in the true model are not visible in the ray inversion. The conductivity is assumed to beconstant at the start of the full-wave inversion, and is found by fitting received power to auniform material model (of permittivity equal to the average of the ray inversion results) anda short dipole antenna radiation pattern. This method produces a conductivity of 2.9 mS m−1for this model.Six frequencies evenly spaced from 75 MHz to 200 MHz are used. The forward modelingoperation is performed at 10 points per shortest wavelength assuming a relative permittivityof 20. Three inversions are performed, one for each of the two transfer functions and anotherwithout the use of a transfer function. The results of the full-wave inversion are shown in844.4. Full-wave Inversion Experiment02468100 1 2 3 4 57891011121314RelativePermittivityx (m)Depth(m)Figure 4.8: Recovered permittivity using 2-D ray inversion. The background half-spaces arerecovered and the two large cylindrical anomalies in the true model are partially resolved. Thesmaller anomalies are not resolved.Fig. 4.9. The regularization parameter is chosen manually using the l-curve method [45]. Somevariation in the results that are obtained with each method is introduced at this step andchoosing a value slightly above or below the corner of the l-curve would change the resultsslightly. In all cases the inversion is able to recover the permittivity very well (all but thesmallest anomaly are clearly visible). The relative error on the recovered permittivity with thecurved-ray transfer function is 2.7 per cent, with the quasi-straight-ray one it is 2.8 per cent,and with no transfer function at all it is 2.9 per cent. This small variation may be related to theselection of the regularization parameter mentioned above. The structure of the conductivity isrecovered moderately well. The numerical value of the conductivity is found accurately (witha relative error of 10.8 per cent for the curved-ray transfer function and 9.3 per cent for thequasi-straight-ray one) when the transfer function is used, but the conductivity is stronglyover-estimated when the transfer function is not used (with a relative error of 53 per cent).It is interesting to note that quasi-straight-ray transfer function produces slightly betterresults with this example than the curved-ray one, even though the opposite is true with thegradient tests in Section 4.3. The model used for the inversion test has a narrower range ofpermittivity values, particularly in the background medium, than the gradient tests. This likelyincreases the accuracy of the quasi-straight-ray transfer function with regard to the use of themean permittivity, as compared to the gradient cases. The transfer functions are evaluatedon the ray-based inversion result, which adds another degree of uncertainty into the mix. Theinversion also uses a somewhat coarse grid that likely introduces some error into both the854.4. Full-wave Inversion Experiment024680 1 2 3 4 57891011121314RelativePermittivityx (m)Depth(m)(a)024680 1 2 3 4 51234567Conductivity(mS/m)x (m)Depth(m)(b)024680 1 2 3 4 57891011121314RelativePermittivityx (m)Depth(m)(c)024680 1 2 3 4 51234567Conductivity(mS/m)x (m)Depth(m)(d)024680 1 2 3 4 57891011121314RelativePermittivityx (m)Depth(m)(e)024680 1 2 3 4 51234567Conductivity(mS/m)x (m)Depth(m)(f)Figure 4.9: Recovered model: The recovered permittivity and conductivity when using thecurved-ray transfer function of eq. (4.25) are shown in (a) and (b), and with the quasi-straight-ray transfer function of eq. (4.26) in (c) and (d). The recovered permittivity and conductivitywhen no transfer function is used are shown in (e) and (f). Without the transfer function, thepermittivity is still recovered and so is the structure of the conductivity, but the values are overestimated.864.5. ConclusionFDFD calculation and the evaluation of the transfer functions. The inclusion of conductivityin the quasi-straight-ray transfer function acts as a constant complex multiplier and would beabsorbed in a linear way into the source estimation that is employed here and so has no effecton this inversion (and so wasn’t included in this test). Ignoring conductivity when evaluatingthe curved-ray transfer function also has negligible effect on the inversion results in this case,but this is likely related to the low conductivity having only a minor effect on phase.4.5 ConclusionWe have shown through analytical and numerical means that the acoustic transfer functionderived by Yedlin et al. [121] and Van Vorst et al. [100], and also by Bleistein [19], basedon curved rays, is applicable to the borehole GPR problem and can be extended to includeconductivity. We have compared its performance to a similar transfer function by Ernst etal. [33] that is based on quasi-straight rays. The analysis assumed lossless slowly varyingisotropic material, but numerical experiments show that for low loss and for situations thatviolate eq. (4.1) to some degree, the methods are still useful.The intended application of the transfer functions is for full-wave inversion of GPR data inmaterial that is known to have relatively low loss. Before performing the full-wave inversion itis assumed that a ray-based inversion was previously performed, which supplies the ray pathsnecessary for the evaluation of the integral in eq. (4.25) [or the traveltimes in eq. (4.26)]. Thesetransfer functions provide the necessary adjustment to 3-D data that is required for recoveringconductivity in a 2-D full-wave inversion.4.6 AcknowledgementsWe wish to thank the reviewers, Giovanni A. Meles and two anonymous, for their very helpfulcriticisms and for referring us to recent works in this area.The Natural Sciences and Engineering Research Council of Canada (NSERC) providedsubstantial support through scholarship and research grants.Calculations were performed using MATLAB, GNU Octave, COMSOL, and MEEP.87Chapter 5Full-Waveform Inversion of GroundPenetrating Radar Data in theFrequency Domain: SimultaneousRecovery of Permittivity,Conductivity, and Source Wavelet5.1 IntroductionFull-waveform inversion (FWI) offers an imaging resolution higher than corresponding ray-based traveltime inversions, but at a significant computational expense. The improvementresults from an increase in the degree of accuracy in the physical model beyond what is possiblewith ray-based simulations. Furthermore, full-waveform inversion directly employs more of theinformation that is contained in the received waveform than traveltime inversion.Despite the added cost, full-waveform inversion has received significant attention in recentyears, primarily for acoustic and seismic inversions. Virieux and Operto [110] present an excel-lent review of FWI, with a focus on acoustic and seismic techniques.The first research into full-waveform inversion of cross-hole GPR data appears to be byJia et al. [54], where they leveraged their previous work on a related electromagnetic imagingproblem [92]. It is a time-domain gradient-descent method where the reconstruction is carriedout in a small region within the larger domain between the boreholes. Later, Kuroda et al. [59]present a method for imaging the full domain. The first application of a full-wave technique885.1. Introductionto cross-hole radar measurements appears to be by Ernst et al. [33] using a similar method toKuroda et al. Kuroda et al. and Ernst et al. use a time-domain gradient-descent method to findpermittivity and conductivity in a sequential manner [35]. The source wavelet is determinedusing a deconvolution algorithm applied to the starting model, and an iterative approach isgiven by Belina et al. [15]. Meles et al. [67] present a time-domain gradient descent approachfor the quasi-simultaneous inversion of permittivity and conductivity where step lengths forthe two material parameter updates are determined separately. Klotzsche et al. [57] applyMeles’s algorithm to measured data in an array of boreholes. Cordua et al. [27] approach theGPR problem from a statistical point of view and present an algorithm that provides useful aposteriori statistical information.There have been very few publications of frequency-domain GPR inversion techniques. Ellef-sen et al. [31] present a gradient-descent frequency-domain approach and verify it using a lab-oratory model. They split the inversion into two steps where the data phases are fit first, andthen the amplitudes, and it is unclear how the source wavelet is estimated. Yang et al. [119]present a Gauss-Newton frequency-domain approach that simultaneously finds permittivityand conductivity and apply it to synthetic data. In their conference publications, no mention ismade about source wavelet estimation and very few technical details of the inversion are given.Lavoue´ et al. [60] present a frequency-domain inversion method for application to surface radar(as opposed to cross-hole GPR in the above papers). They apply a quasi-Newton approachwhere the Hessian is approximated using previous gradients and updated models. The forwardproblem that is solved has the opposite polarization of what is required for 2D borehole radarsimulations, and determination of the source wavelet is not investigated.On modern desktop machines it is now feasible to solve frequency-domain forward-modelingproblems using direct matrix solvers, and the resulting factorization can be re-used for manyright-hand sides with relatively little additional cost. This leads to an advantage of a frequencydomain approach over time-domain in that the Jacobian matrix can be calculated in a rea-sonably short period of time and thus a Gauss-Newton algorithm can be employed instead ofgradient descent.In this work we present a frequency-domain full-waveform inversion algorithm for cross-hole GPR measurements that simultaneously finds the relative permittivity, conductivity, and895.2. Forward Modelsource wavelet associated with the GPR transmitter. We employ the Gauss-Newton methodto minimize an objective function that includes a regularization term, thus allowing a coolingprocess to be applied during the inversion where structure is slowly added by reducing the levelof regularization as the inversion progresses (see, e.g., [37]). An analysis of the interaction of aGPR antenna with an air-filled borehole is also presented and it results in a simplified antennamodel that can easily be included in finite-different frequency-domain (FDFD) simulations.The algorithm is successfully applied to a synthetic test model where data is generatedusing an FDTD solver and inverted using the FDFD algorithm presented here. The novelsimultaneous inversion of permittivity, conductivity, and source wavelet is shown to effectivelyrecover both the model and the source wavelet. The algorithm is also applied to data collectedfrom a karst aquifer at the Laboratoire Souterrain a` Bas Bruit (LSBB) (http://lsbb.oca.eu).This is a preliminary test of the algorithm and errors that are introduced into the field data,such as uncertainty in position and some significant distortion of the received signal, cause somedifficulty in the inversion. In particular, the conductivity results are corrupted by the distortedsignal but the permittivity results show a significant increase in resolution. Analysis of theimpact of various measurement uncertainties on the inversion (and handling them) is a topicfor future research.Unless explicitly stated otherwise, anywhere that vectors are raised to a power, or a constant,such as e, is raised to the power of a vector, element-by-element operations are implied. Forexample eǫ˜ =[eǫ˜1 eǫ˜2 · · · eǫ˜n]T or k−2 =[k−21 k−22 · · · k−2n]T .5.2 Forward ModelIn any inversion formulation, the forward modeling operation produces simulated data by mod-eling the physical processes that produce the measured data. For full-waveform GPR inversion,the forward modeling operation attempts to capture all of the electromagnetic aspects of theGPR experiment, and is accomplished using a finite-difference frequency-domain (FDFD) for-mulation. The outcome of the forward modeling operation, the predicted data, is the voltageat the terminals of the receiving antenna because an amplified version of this is recorded by theGPR unit.905.2. Forward ModelSome aspects are too complex to realistically model, such as the electronics inside of the GPRantenna units. The interaction between the borehole and the GPR units, although possible tomodel with FDFD, would require a very fine grid in the vicinity of the antenna. Further, if onewere to model both the receive and transmit antennas, then a separate forward problem wouldneed to be solved for each source-receiver combination because the receive antenna becomespart of the model. To simplify the FDFD formulation a regular square grid and simplifiedantenna model are used. To reduce memory and computational resources during the inversion,the material parameters are stored on a coarser grid than the FDFD grid and interpolated whenneeded.The FDFD implementation is discussed in Section 5.2.1 and the associated perfectly matchedlayer (PML) in Section 5.2.2. The simplified antenna model is discussed in Section 5.2.3. Inthis work we use the eiωt time convention.5.2.1 FDFD Helmholtz FormulationMaxwell’s equations govern the electromagnetic fields associated with the GPR problem. Maxwell’scurl equations, in the frequency domain, are repeated here for reference. The first is theMaxwell-Ampe´re’s equation,∇×H(x, ω) = J(x, ω) + iωD(x, ω), (5.1)where H is the magnetic field, D is the electric flux density, and J is the current density. Thesecond equation is the Maxwell-Faraday equation,∇×E(x, ω) = −iωB(x, ω), (5.2)where E is the electric field and B is the magnetic flux density. Additionally we have threeconstitutive relations that describe the material properties,D(x, ω) = ǫ′(x)E(x, ω), (5.3)J(x, ω) = σ(x)E(x, ω) + Js(x, ω), (5.4)915.2. Forward ModelandB(x, ω) = µ0H(x, ω). (5.5)In equations (5.3) through (5.5), ǫ′ is the electric permittivity, σ is the conductivity, µ0 isthe magnetic permeability of free space, and Js is the source current. For the remainder ofthis derivation, we drop the explicit spatial and frequency dependence. In this inversion for-mulation the material properties are assumed frequency independent, but a frequency-domainimplementation such as this allows frequency-dependent material parameters to be included ina straightforward manner. The prime on ǫ′ indicates that it is the real part of the complexpermittivity that is commonly used in electrical engineering,ǫ = ǫ′ − iσω . (5.6)Combining equations (5.1), (5.3), (5.4), and (5.6) allows us to write∇×H = Js + iωǫE. (5.7)We aim to solve these equations in two dimensions on the xz-plane. To do this, we startwith the 2.5D assumption that the material parameters and source (and therefore the fields)do not vary as a function of y. A further simplification that we make is that the material isnon-magnetic and that the permittivity and conductivity are frequency independent. Theselast two points are indicated in equations (5.3) through (5.5).With a source current excitation in the z-direction (such as that produced by a GPR an-tenna), the symmetry of the 2.5D problem makes Hx and Hz zero in the xz-plane that containsthe source. For simplicity, we make this the y = 0 plane. By dividing equation (5.7) throughby ωǫ, then taking the curl of both sides and substituting equation (5.2) and equation (5.5)and finally dividing through by ωµ0 allows us to write∇× 1k2∇×H−H = ∇×Jsk2 , (5.8)925.2. Forward ModelxzHyEx, k−2Ez, k−2, Jsz∆x∆zFigure 5.1: The computation grid. The y-component of the magnetic field is on the referencegrid, while the electric field components are staggered off of it by half of a grid spacing. Thelocation of the source current is also indicated. k−2 is also staggered by half a grid spacing toline up with the electric field. This is not the coarse grid where k−2 is stored, but rather wherethe interpolation operations, MxnMc and MznMc, place it in equation (5.12).wherek2 = ω2ǫµ0. (5.9)In deriving equation (5.8) we have used the assumption that the material is nonmagnetic andµ0 is independent of position. Expanding equation (5.8) and taking into account the symmetryresults in∂∂x1k2∂∂xHy +∂∂z1k2∂∂zHy +Hy =∂∂xJszk2 . (5.10)We can discretize this equation onto a Yee grid [122] and define the Hy positions as thereference grid with coordinates (x,z). Although not explicitly part of the equation, the first-derivatives of Hy correspond to the electric field, and their locations are shifted off of thereference grid by half of a grid spacing. The locations of the associated wave numbers inequation (5.10) are also shifted. The z-component of the electric field, Ez, has values at (x −dx/2,z) and Ex has values at (x,z − dz/2). Similarly, Jsz also has values at (x− dx/2,z). Thegrid is illustrated in Fig. 5.1.In matrix form we can write this asAhy = q, (5.11)where hy is a column vector containing Hy throughout the domain, q is a column vector935.2. Forward Modelcontaining the discretized values of the right-hand side of equation (5.10), and A is a matrixcontaining the discrete form of the differential operator on the left-hand side of equation (5.10).The discrete operator A can be expressed asA =Dxp diag(MxnMck−2)Dxn+Dzp diag(MznMck−2)Dzn + I.(5.12)The subscripts on the matrix operators describe the coordinate axis and direction in which theoperation acts (p for positive and n for negative). All three terms of A operate on the referencegrid and produce results on that grid. The D matrices are central differences that shift theresult forward or backward half a grid space in the specified direction. For example, Dxp takesthe x-derivative such that the resulting vector has values that are accurate half a grid space tothe right of whatever it operated on. The combination DxpDxn is therefore a second derivativewith no coordinate shifting. The vector k−2 contains 1/(ω2µǫ) evaluated on a coarse grid in theinterior region, and Mc interpolates it onto the reference grid and expands this interpolationinto the PML regions (see Section 5.2.2). The averaging operators Mxn and Mzn transfer k−2onto the Ez and Ex grids respectively.The source vector, q, is related to the source current density, js, byq = Dxp diag(MxnMck−2)js. (5.13)The elements of js are the discretized values of Jsz. This source current represents the transmitantenna, and is discussed in Section 5.2.3.From the Ampe`re equation (equation (5.1)) the electric field can be found by taking thecurl of hy. This results inez =1iω diag(MxnMc~ǫ−1)Dxnhy, (5.14)where ez is the z-component of the electric field. The vector ~ǫ−1 contains ǫ−1 evaluated on thecoarse grid.Finally, the predicted data are computed using the electric field by means of a measurement945.2. Forward Modeloperator, Q:dpred = Qez (5.15)The measurement operator approximates the output voltage of a receiving antenna by integrat-ing the incident electric field along the length of the antenna. The approximation that is usedis discussed in Section 5.2.3.The data that are computed in equation (5.15) is a subset of the complete data set associatedwith the current transmitter location and the current frequency. As matrix Q is associated withreceiver positions and the receiver positions may change depending on the source position (forexample when the transmitter moves from one borehole to the other), a different measurementoperator may be required per source location. Additionally, Q is a function of frequency becausethe antenna current in equation (5.19) is a function of frequency.5.2.2 Perfectly Matched LayerThe perfectly matched layer (PML) is implemented using complex stretched coordinates fol-lowing Chew and Weedon [26], and Rappaport et al. [85]. This type of PML is related toBerenger’s [18] and can also be implemented as an anisotropic material [96].This type of absorbing boundary condition is very effective and is easy to implement inthe FDFD method. Essentially, all derivatives in the forward operator A are prefixed with acomplex scaling factor. For example, the x-derivative becomes∂∂x →1sx∂∂x, (5.16)where sx is 1 + iσx/(wǫ′) (the real part can also be greater than one to cause real coordinatestretching within the PML region). This scaling maintains perfect matching (i.e. it is reflection-less) for waves incident on the PML at any angle while it also causes them to decay within thePML. In this work, we use the average permittivity in the domain when computing sx, sy, andsz. Although this may not be ideal, it simplifies the construction of the derivative matrices andmakes them essentially independent of the material parameters, which simplifies the Jacobiancalculation. The parameter σx (along with the associated one for z) is a function of position955.2. Forward Modeland follows an exponential grading from zero at the boundary between PML and the interiorregion to a large value at the edge of the computational grid. Within the central region of thedomain (i.e. outside of the PML) σx = σz = 0 and the derivatives behave as usual.The PML scaling, as described in the previous paragraph, makes the forward modelingmatrix A asymmetric. Although that is not a major difficulty, it is more convenient if it issymmetric. It is found that premultiplication by sxsz makes it symmetric. This premultipli-cation can be included within the derivative matrices and as a diagonal matrix in place of theidentity matrix that makes up the last term of A in equation (5.12).For additional details, the reader is referred to the references. Taflove and Hagness [91, pp.305–308] give some guidelines for computing appropriate PML coefficient values.5.2.3 Antenna ModelTo fully simulate the antenna structure and borehole for all source and receiver positions ateach iteration of an inversion would require a great deal of computation time. Doing so isnot realistic with current desktop machines. In this section we introduce a simplified antennamodel, and the reasoning behind it, that aims to be computationally inexpensive and easy toimplement while maintaining an acceptable level of accuracy. The model, in its current form, isonly valid for air-filled boreholes surrounded by rock. Water filled (or partially filled) boreholesand antenna contact with the formation require further research.The structure that is formed when an antenna is centered inside of a borehole is much likethat of a coaxial transmission line. If the dielectric constant of the rock outside of the hole islarge then the walls of the borehole affect the fields inside the borehole much like a conductorwould. An important difference is that energy leaks out of this structure where it does not in thecase of an ideal coaxial line. In both cases, where the center conductor (the antenna) terminatesand the outer cylinder (borehole) continue, the dimension of the resulting circular waveguideis often too small to allow propagation of typical GPR frequencies (this case is assumed here).COMSOL experiments justify this line of reasoning and show that the electric field alongthe length of the antenna is nearly radial (see Fig. 5.2), just as it is inside of a coaxial cable.The field continues to be so as it crosses the borehole boundary (unlike in the coaxial cablewhere it abruptly falls to zero in the outer conductor). This suggests that the borehole-antenna965.2. Forward Model−0.3−0.2−0.1 0 0.1 0.2 0.3−0.3 −0.2 −0.1 0 0.1 0.2 0.3x (m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)z(m)ǫr = 10ǫ0Antenna in boreholeFigure 5.2: Electric field streamlines computed in COMSOL for a thin cylindrical antenna insideof a cylindrical borehole. The surrounding rock has ǫr = 10 and the frequency is 150 MHz. Thelines are nearly horizontal and continue to be so as they cross the borehole wall.interaction for the transmit antenna can be simplified to an assumed current distribution (cal-culated using transmission line theory) at the antenna’s location, but without the antenna andborehole structure present in the actual FDFD computation.This reasoning can be partially justified by looking at Ampe`re’s equation for H integratedalong a loop that is just slightly larger than the borehole, as illustrated in Fig. 5.3. ApplyingStokes theorem around the loop gives∮H · dl =∫∫(Js + iωǫE) · dA. (5.17)If E is radial (that is, perpendicular to the normal vector of the disc that spans the integrationcircle), then the contribution of the displacement current density, iωǫE, to the surface integralin equation 5.17 is zero. The electric field is not exactly radial, but iωǫE ·dA is small comparedto Js ·dA. Removing the displacement current contribution from equation (5.17) leaves us withan expression for the magnetic field due to a current element in the low-frequency limit (orvery close to the current element). This tells us that from outside the borehole, as long as theborehole is small enough, the borehole and antenna combination behaves similarly to a currentelement with the same current that is on the antenna (but now situated within the medium,975.2. Forward ModelBorehole WallIntegration PathRockAirAntennaElectric FieldLinesJsEFigure 5.3: Radial electric field lines around an antenna inside a borehole. An integration pathfor Ampe`re’s law is shown as a dotted line just outside of the borehole boundary.rather than within a borehole).The current distribution along the antenna, assuming an ideal transmission line model, isgiven byI(z) = a sin[k0( lant2 − |z|)], |z| ≤ lant2 , (5.18)where a is a complex constant (independent of z) that sets the magnitude and phase of thedriving current and k0 = w√µ0ǫ0. This current distribution is also often used to model thindipole antennas in free space even though it is not the exact current distribution [14, p. 170].The sinusoidal weighting that appears in equation (5.18) is the standing wave pattern thatarises in an open-ended transmission line of length lant/2.Although removing the borehole and replacing the antenna with the assumed transmissionline current results in a different field distribution in the vicinity of the antenna, numericalexperiments show that the far-field pattern is very similar to that produced by the antennainside of a borehole. The pattern differs only slightly in amplitude and phase, and this can bewell accommodated by adjusting the parameter a. This factor is of little consequence for theinversion because it becomes part of the source wavelet, which we are inverting for. However,variation in this factor due to variation in material parameters would be an issue. Numericalexperiments show that changing the permittivity of the rock has a small effect on the requiredadjustment to phase and amplitude of the current on the antenna.To verify the accuracy of this model, we compare it to COMSOL simulation results wherewe have a 14 cm diameter borehole with a 0.5 m long transmitting antenna at 150 MHz andan otherwise uniform permittivity. To determine the sensitivity of our model to the surround-985.2. Forward ModelTable 5.1: Experimental results comparing an antenna inside a borehole to the assumed cur-rent model without borehole. The experiment is performed at 150 MHz. The RMS error ofthe electric field is calculated by taking the norm of the difference between the COMSOL sim-ulation results (with borehole) and the idealized dipole case, and dividing by the norm of theCOMSOL results. For the simplified model, the coefficient a corresponds to equation (5.18).For comparison, the results for a short dipole of length 0.01 m with a current of a Amps arealso included.ǫr Coef. a (A) 6 Coef. a |Coef. a| (A) Rx E-field RMS Error (%)Simplified 0.5 m dipole model:10 0.117 + 0.600i 79.0◦ 0.611 6.1614 0.057 + 0.573i 84.3◦ 0.576 5.18Short dipole (lant = 0.01 m) with I = a A:10 2.004 + 10.419i 79.1◦ 10.609 10.1014 0.945 + 9.662i 84.4◦ 9.709 12.73ing rock, we perform the experiment twice with different permittivities (ǫr = 10 and 14) forcomparison. The transmitting antenna is centered within a borehole and centered verticallyin a domain with a height of 8 m. The received signal is taken to be the vertical electric fieldmeasured along a vertical line from top to bottom of the domain and 3.5 m away from thetransmitting antenna. Using GNU Octave, we do a least-squares fit between the COMSOLdata and fields of an ideal current element in a uniform medium with the same permittivity.The current on the element is described by equation (5.18) and the least-squares fit finds thebest magnitude, a (i.e. the drive current for an assumed current without the borehole). Forcomparison, we also calculate how well a short dipole performs. The results of this experimentare summarized in Table 5.1. The 40 % change in permittivity leads to a 5.3◦ phase and a6 % amplitude adjustment in the drive current on the assumed current distribution withoutthe borehole. The overall error with the best fit drive current is 5 or 6 %. The error associatedwith the short dipole model is approximately twice that of the simplified model. These resultsshow that large changes in the electric permittivity of the rock surrounding the borehole havea relatively small effect on the equivalent current. This allows us to use the same source ap-proximation at each position in each borehole and reduces the number of unknowns we mustinvert for.The principle of reciprocity tells us that not only is the receiving antenna just as importantas the transmit antenna, but also that if we know the radiation properties of the transmit995.3. Inversion Formulationantenna then we can determine the behavior of the receive antenna. We can write the outputvoltage of a receiving antenna in terms of the incident electric field (see, e.g., [55, pp. 345–353]),Vo = −1Io∫ lant/2−lant/2Ez(z)I(z)dz, (5.19)where I(z) is the current on the antenna in transmit mode, Io is the driving current at thefeed point in transmit mode, and Vo is the open-circuit voltage at the terminals of the receiveantenna given an incident field Ez. As long as the amplifier that the antenna is connected to islinear then its sampled output is related to the antenna terminal voltage by a complex constant(at a given frequency). The loading that the amplifier puts on the antenna will cause the actualterminal voltage to deviate from Vo, but as long as the antenna’s impedance does not change asit is moved up and down inside of a borehole, then this can also be accounted for by a complexconstant. The experiments performed for the transmit antenna show that the current on theantenna is affected only slightly by moderate variation in the surrounding rock permittivity,thus suggesting that the antenna’s circuital behavior also remains relatively constant.The measurement operator Q implements the integral in equation (5.19) and interpolationof ez to allow receiver positions that are not exactly on the grid. Similarly, interpolation is usedto obtain the transmit current density, js, from equation (5.18) to allow transmitter positionsthat are offset from the grid.5.3 Inversion FormulationThe inversion is cast as a regularized minimization problem using the Gauss-Newton method[see, e.g.: 97, 82]. The objective function that we minimize is:φ(m) = φd + β2φm= 12 ‖Wd(F[m]− dobs)‖2 + β2 12Nreg∑i=1‖Wi(m−mref )‖2(5.20)In equation (5.20), F is the forward modeling operator, dobs is the observed data, Wi areregularization matrices (of number Nreg), Wd weights the data residual, δd = F[m]− dobs,1005.3. Inversion Formulationm is the model vector, mref is a reference model, and β is the regularization parameter. Theweighting in Wd can be used to weight the data according to measurement uncertainty and/orother parameters such as frequency or source-receiver offset to emphasize certain data overother data. Here we weight it with the reciprocal of an estimate of the standard deviation ofthe measurement error. The objective function can be broken into two pieces: the data norm,φd, and the model norm, φm. The model vector contains log-scaled values (indicated with atilde over the variables) of the relative permittivity, ~ǫr = eǫ˜, and conductivity, ~σ = eσ˜σ0 (withσ0 = 1 Sm−1), along with the source wavelet coefficients (the complex frequency components):m =[ǫ˜T σ˜T sT]T(5.21)The regularization matrices only act on the part of m that contains the permittivity andconductivity, and ignore the wavelet coefficients. Note that ~σ, σ˜, ~ǫr, and ǫ˜ are vector quantitiesas they are discretized over the area of interest for the inversion.The data, dobs, contains the discrete Fourier transform coefficients of the measured wave-form at all frequencies of interest for the inversion. The inversion acts on all of these frequenciessimultaneously. The synthetic experiment in Section 5.6 illustrates a case where six frequen-cies are sufficient to produce an inversion result that is a vast improvement over traveltimetomography.Relative permittivity and conductivity often differ by orders of magnitude, leading to adifference in scale in their log-scaled values. Therefore some relative scaling is required in theregularization matrices so that regularization is applied at a controlled level to both. In thiswork, the regularization seeks to find both a smooth model and one that is close to the referencemodel. For smoothness, the regularization matrices take the x and z first derivatives of themodel. The portion of these matrices that operate on ~σ are scaled by ǫ˜av/σ˜avζ. The ratioof the average values of ǫ˜ and σ˜ scales the derivatives of ǫ˜ and σ˜ such that they are similarwhen either parameter varies by the same factor. The additional scaling parameter, ζ, allowsthis relationship to be changed if one of the parameters is expected to have a greater relativevariation than the other. The same scaling can be used for closeness to the reference modelso that similar percent differences between recovered and reference models for each parameter1015.3. Inversion Formulationaffect the model norm to a similar degree. This is the only place in this method where anykind of manual adjustment of the weighting of permittivity and conductivity is required (i.e.the selection of ζ). No scaling is applied to the Jacobian.The inversion proceeds by minimizing equation (5.20), starting with a high value of theregularization parameter and reducing it until an appropriate level of regularization is found.The starting value is chosen such that φm dominates the objective function so that a verysmooth model is found. This ensures that the inversion begins by adding only the simplest (lowspatial frequency) structure first to help guide the inversion toward the global minimum. Theminimization of equation (5.20) is an iterative process where a model update, δm, is found andapplied to the model until convergence is achieved. The model update is given byδm = −[ℜ{J†W†dWdJ}+ β2Wr]−1·[ℜ{J†W†dWdδd}+ β2Wr(m−mref )](5.22)where the superscript † represents complex conjugate transpose, J is the Jacobian (see Sec-tion 5.4) andWr =Nreg∑i=1WTi Wi. (5.23)The ǫ˜av/σ˜avζ scaling discussed previously is contained within Wr. For the 2D problems dis-cussed here, equation (5.22) is small enough to be solved directly on a current desktop machinewith adequate memory. The update specified in equation (5.22) is a Gauss-Newton updatewhere an approximate Hessian matrix (the matrix being inverted) is used in place of the fullHessian found in the Newton method. This avoids calculation of second derivatives and usuallyhas only a small effect on convergence because the neglected portion of the Hessian depends onthe data residual, which decreases as the inversion progresses.Convergence is checked using the relative gradient method [44]. Iterations are also stoppedif the process stagnates. The regularization parameter can be determined by means of thel-curve [45].1025.4. Jacobian5.4 JacobianThe Jacobian matrix contains all first partial derivatives of data with respect to model param-eters. The model parameters for this inversion formulation are the log-scaled permittivity andconductivity, and the source wavelet coefficients. As the inversion operates on all frequencies ofinterest simultaneously, the data contains all of these frequencies and the Jacobian must alsoact on them.With the assumption that the source signature does not change as a function of position,we can extend the model vector to contain a single source wavelet,m =[ǫ˜T σ˜T s(1)R · · · s(Nf )R s(1)I · · · s(Nf )I]T, (5.24)where sR and sI are the real and imaginary parts of the complex source amplitude (the discreteFourier transform coefficients of the source wavelet) at each frequency and the superscripts arethe frequency index. These amplitudes correspond to the coefficient, a, in equation (5.18).Assuming linear material parameters, we can writed(f)pred = d♭(f)preds(f) = d♭(f)pred(s(f)R + is(f)I), (5.25)where d♭(f)pred are the predicted data calculated with a unit driving current and s(f) is the complexsource amplitude at the given frequency. The superscript f indicates that the expression refersto all data at a particular frequency. So, we can write∂d(f,i)pred∂s(f)R= d♭(f,i)pred =d(f,i)preds(f) (5.26)and∂d(f,i)pred∂s(f)I= id♭(f,i)pred = id(f,i)preds(f) , (5.27)where the index i refers to a particular source/receiver combination.1035.5. 3D to 2D ConversionThis allows us to write the complete Jacobian asJ =J(1)ǫ˜ J(1)σ˜d(1)preds(1) 0 · · · 0 id(1)preds(1) 0 · · · 0J(2)ǫ˜ J(2)σ˜ 0d(2)preds2 · · · 0 0 id(2)preds2 · · · 0... ... ... ... . . . ... ... ... . . . ...J(Nf )ǫ˜ J(Nf )σ˜ 0 0 · · ·d(Nf )preds(Nf )0 0 · · · id(Nf )preds(Nf ), (5.28)where each of the internal Jacobians, J(f)ǫ˜ and J(f)σ˜ , include all source/receiver combinationsat the specified frequency (see equation (A.7)). Expressions for J(f)ǫ˜ and J(f)σ˜ are derived inAppendix A.5.5 3D to 2D ConversionThe measured 3D data are converted to 2D equivalent data using the following transfer function[104]4:TF = E2DE3D= a0e−ipi4√2πω∫ stot0v(s)ds, (5.29)where a0 = 1 m−1. The integration in equation 5.29 is carried along rays traced from sourcesto receivers, and s is the arc length along the ray. The total arc length from source to receiveris stot. Rays are traced by following the traveltime gradient backward from receiver to sourceon a traveltime field that is computed using the eikonal equation. This process is described in[101] and a similar one in [2].In order to approximately account for attenuation in the medium, the velocity that isused in equation (5.29) is complex and computed directly from the complex permittivity ofequation (5.6), v = 1/√µǫ. Our ray inversion results that are used as a starting point forthe full-wave inversion do not include an estimate of the conductivity. For the purposes ofperforming the 3D to 2D conversion, and as a starting point for the full-wave inversion, anestimate of a constant background conductivity is found that best fits a simple dipole-to-dipoletransmission model.4A mistake in the units used for the sources was discovered after publication. The factor a0 is added here tocorrect for this. This function assumes numerically equal source currents, but for the 3D case the current is aninfinitesimal dipole and for the 2D case it is a current filament.1045.5. 3D to 2D ConversionAssuming that the radiation from the transmitter travels in a straight line to the receiver,that the antennas can be approximated as infinitesimal vertically-oriented dipoles, and that thevariation of material parameters across the domain is small, then we can estimate the receivedpower, Prx (in Watts), to bePrx = Ptx( λ4πr)2e−2αr[32 sin2(θ)]2, (5.30)where Ptx is the transmit power, r is the distance from transmitter to receiver, λ is the wave-length, and α is the attenuation constant. This is a modified form of the Friis transmissionequation [41] that includes path loss and the dipole radiation pattern. The 3/2 sin2(θ) term isthe gain of a small dipole. For simplicity, we find a best fit of equation (5.30) to the receiveddata at the dominant frequency rather than at all frequencies. A scaled version of the receivedpower at the dominant frequency is easily calculated by squaring the magnitude of this com-plex frequency component of the received signal. Such a scaling as no effect on the computedattenuation constant.After taking the log of equation (5.30), it can be written as a linear system for each source-receiver combination at a single frequency:−2r1 1... ...−2rN 1αln(C)=ln(Prx1) + 2 ln(r1)− ln[sin4(θ1)]...ln(PrxN ) + 2 ln(rN )− ln[sin4(θN )], (5.31)where the subscript indicates source/receiver combination andC = Ptx(3λ8π)2. (5.32)The least-squares solution to equation (5.31) produces an estimate of the attenuation coef-ficient α. From α, the conductivity can be computed asσ ≈ 2α√ ǫµ0, (5.33)where the permittivity, ǫ, is taken to be the average value throughout the domain computed1055.6. Synthetic Experimentsusing the ray inversion results. Equation (5.33) is valid at frequencies where (σ/(ωǫ))2 ≪ 1,which is the case for our real and synthetic data sets at the dominant frequency. For lowerfrequencies or higher conductivities, a more exact expression may need to be used (see, e.g.,[13, p. 150]).The numerical experiments in section 5.6 show that this conductivity estimate provides agood starting point for the inversion. However, it has been observed for some data sets thatequation (5.31) produces a negative conductivity and a different estimate should be used.5.6 Synthetic ExperimentsAn experiment is performed where synthetic data is produced by both a 2D and a 3D FDTDsimulation using MEEP [73] and both data sets are inverted for permittivity, conductivity andthe source wavelet using our method (the transfer function, equation (5.29), is applied to the3D data before inversion). The true 3D model has no variation in the y-direction, and sois 2.5D in nature. Inverting both 2D and 3D data sets allows us to illustrate the successfulapplication of the 3D to 2D transfer function. The transmitting antennas are represented bynearly infinitesimal current sources and the received signal is taken to be the vertical componentof the electric field. A Gaussian-modulated sine wave with a center frequency of 150 MHz isused as the source, and is shown in Fig. 5.4. A total of ten source locations (five at each side ofthe model) are used. Receivers are spaced at 25 cm intervals at the opposed side of the modelfrom a given transmitter.The time domain output of the FDTD solver is converted using FFT to a set of six frequen-cies equally spaced from 75 MHz to 200 MHz. Frequency sampling strategies for acoustic andseismic FWI have been explored by a number of authors (see, e.g.:[65, 68, 88]). The radar caseis more narrow-band than the acoustic case, and for simplicity we choose an equally spaced setof frequencies over a bandwidth as wide as possible before the signal becomes weak and ap-proaches the noise floor. Further research into frequency selection for cross-hole GPR FWI mayyield different results. These six frequencies are inverted simultaneously, rather than invertingfirst low frequencies and then higher ones as suggested by Pratt [81] and Bunks et al. [20]. Theuse of cooling (reduction of the regularization parameter as the inversion proceeds) helps avoid1065.6. Synthetic Experiments−1−0.500.515 10 15 20 25 30CurrentDensity(A/m2 )Time (ns)Figure 5.4: Source wavelet for the synthetic inversion experiments.local minima.The true model and the results of the two inversions are shown in Fig. 5.5. The regularizationparameter is chosen by means of the l-curve, shown in Fig. 5.6 for each data set. For thesynthetic data, a distinct knee is visible in the l-curves, making the choice of β clear. Theresolution of the FWI permittivity results are drastically improved compared to the first-arrivalray-based inversion results. The ray inversion is only capable of recovering the largest cylindricalanomalies and the material slabs, but the FWI permittivity results clearly show all but thesmallest anomaly. Even the smallest anomaly is visible when the 2D synthetic data are inverted,but it is essentially indistinguishable from noise. The smallest anomaly has a diameter of 0.25 mand is slightly less than half of the dominant wavelength of 0.63 m in the background medium(ǫr = 10). The next largest anomaly has a diameter of 0.5 m, indicating that the resolutionlimit of this FWI scheme and data acquisition geometry is about one wavelength. A denseracquisition geometry may bring the resolution of the inversion closer to the diffraction limit ofλ/2. The resolution of the FWI conductivity results are somewhat less than for permittivity.Other researchers have found similar limitations and difficulties when inverting for conductivity(see, e.g.: [57]). For reference, the starting value of the conductivity for the inversion is 2 mS/m.Only the largest high conductivity anomalies are clearly recovered. The more resistive smalleranomalies are not clearly resolved. This could be related to either the size of the anomalies orthat they are more resistive or both.The transfer function of equation (5.29) applies a constant phase shift with a computed1075.6. Synthetic Experiments0 1 2 3 4 5 6 7 8 9 100 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 7 8 9 1011121314True Model Ray Inversion FWI (2D Data) FWI (3D Data)RelativePermittivityx (m)x (m)x (m)x (m)z(m)(a)0 1 2 3 4 5 6 7 8 9 100 1 2 3 4 5 0 1 2 3 4 5 0 1 2 3 4 5 1.01.52.02.53.03.54.04.55.0True Model FWI (2D Data) FWI (3D Data)Conductivity(mS/m)x (m)x (m)x (m)z(m)(b)Figure 5.5: Synthetic experiment true model and inversion results: (a) permittivity and (b)conductivity. The left-most panels show the true permittivity and conductivity. The recoveredrelative permittivity from a ray-based inversion is shown in (a). The FWI results for the 2Dand 3D synthetic data sets are shown.1085.6. Synthetic Experiments101001000100000.0001 0.001 0.01 0.1 1 10 100φmφ dφ d(a)101001000100000.0001 0.001 0.01 0.1 1 10 100φmφ dφ d(b)Figure 5.6: L-curves for the (a) 2D and (b) 3D data sets. The chosen value (closest β to thecorner of the curve) is indicated with a circle.amplitude adjustment to the measured data. The phase shift is inconsequential to this inver-sion scheme because any multiplicative complex constant applied to all of the data at a givenfrequency will get absorbed into the source wavelet and will not affect the inversion. The am-plitude adjustment is critical for accurately finding the conductivity. Without it, the inversionsignificantly over-estimates conductivity [104].The source wavelet is found simultaneously with permittivity and conductivity during theinversion, and so it changes as β is reduced. The evolution of the source wavelet is shown forthe 2D data set in Fig. 5.7 and for the 3D data set in Fig. 5.8. In Figs. 5.7 and 5.8, the leftplot shows unscaled wavelet coefficients and the right plot shows the coefficients scaled suchthat the true values are unity. The initial value of the source wavelet is calculated using themethod suggested by Pratt [81] and is indicated on the figures as a circle. As the inversionprogresses and β is decreased, the source wavelet coefficients follow the paths to their chosenvalues (the best β determined by means of the l-curve) marked by crosses. The paths extendpast the chosen values as β is further reduced. A direct comparison of the time-domain sourcewavelets is impractical because of the limited accuracy of an inverse Fourier transform with sofew frequencies.For both data sets, the recovered source wavelet is close to the original, as indicated by theproximity to unity in Figs. 5.7b and 5.8b. Indeed, their amplitudes are found quite accurately1095.7. Inversion of Radar Data from LSBB−0.06−0.04−0.020.000.020.040.06−0.06−0.04−0.02 0.00 0.02 0.04 0.0674.98 MHz99.98 MHz149.97 MHz199.95 MHz99.98 MHz199.95 MHz124.97 MHz174.96 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsR(a)−1.50−1.00−0.500.000.501.001.50−1.50 −1.00 −0.50 0.00 0.50 1.00 1.50124.97 MHz149.97 MHz174.96 MHz199.95 MHz99.98 MHz 74.98 MHz74.98 MHz199.95 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsR(b)Figure 5.7: Recovered source wavelet coefficients (real, sR, and imaginary, sI , parts) for the2D data set: (a) raw values; (b) normalized values (the black dot at unity represents thenormalized value of the true source wavelet coefficients). Circles indicated the initial sourcewavelet estimate and the paths indicated the progression to the chosen values, marked withcrosses.with a frequency-dependent phase shift. This suggests that the recovered permittivity may beoffset slightly from the true value, but these errors are also approaching the level of accuracythat is attainable with the discretization that is used for the FDTD and FDFD models. Thepaths traced out as β is reduced exhibit distinct sharp bends that correspond closely to thecorner of the l-curve (Fig. 5.6) and therefore the chosen values of the wavelet coefficients.5.7 Inversion of Radar Data from LSBBGPR data were collected from the Vaucluse karst aquifer at LSBB near Rustrel, France. Studiesof the architectural and petrophysical properties of a fracture zone within the aquifer [53, 52]and of CO2 injection [51] have been carried out. The GPR data were collected to allow imagingof structure and water content within the fracture zone. These data are inverted using thealgorithm discussed in section 5.3. However, this is a preliminary application of the methodologyand it does not take into account some of the issues associated with the real data collected inthis survey. In particular, antenna position uncertainty and distortion of the received signalare not specifically addressed. The distortion appears to lead to corruption of the recoveredconductivity.A total of five boreholes (Fig. 5.9), with four of them in a diamond-shaped cluster, were1105.7. Inversion of Radar Data from LSBB−0.30−0.20−0.100.000.100.200.30−0.30 −0.20 −0.10 0.00 0.10 0.20 0.3074.98 MHz99.98 MHz149.97 MHz199.95 MHz99.98 MHz199.95 MHz124.97 MHz174.96 MHzs I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000s I·1000sR · 1000(a)−1.50−1.00−0.500.000.501.001.50−1.50 −1.00 −0.50 0.00 0.50 1.00 1.50124.97 MHz149.97 MHz174.96 MHz199.95 MHz99.98 MHz74.98 MHz74.98 MHz199.95 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsR(b)Figure 5.8: Recovered source wavelet coefficients for the 3D data set: (a) raw values (scaled forplotting); (b) normalized values (the black dot at unity represents the normalized value of thetrue source wavelet coefficients). Circles indicated the initial source wavelet estimate and thepaths indicated the progression to the chosen values, marked with crosses.drilled in the vicinity of a fracture zone. We apply the inversion algorithm to measurementstaken in the four holes in the diamond-shaped cluster. The GPR measurements were takenusing 250 MHz MALA˚ Geoscience borehole radar units. For the boreholes directly opposite eachother, the transmitter was held at five equally spaced positions from top to bottom of each holewhile the receiver was moved from top to bottom of the opposite hole (taking measurementsevery 5 cm). For the diagonally adjacent holes, the same procedure was followed but withonly three transmitter locations. To reduce processing time receiver measurements at 25 cmspacing are used for the full-wave inversions. As in the synthetic case, six equally spacedfrequencies are used, but over a slightly narrower range of 100 MHz to 200 MHz. As anexample, with this receiver spacing, the inversion takes approximately 18 hours to sweep througha very broad range of regularization parameter values on a quad-core i7 machine and usesabout 9 GB of RAM as currently implemented. The time could be reduced by selecting anarrower range of regularization values.5 The 5 cm spacing is very small and including allof the data points is unlikely to enhance the inversion results significantly. This is the samespacing that is used successfully for the synthetic tests. However, the transmitter spacing islarger. Reducing the number of data points reduces the number of rows in the Jacobian, thus it5For comparison, the method of Ernst et al. [33] took less than 12 hours using N + 1 processors, where N isthe number of transmitters.1115.7. Inversion of Radar Data from LSBBB B B BH 2BH 3Tunnel Wall Scale: 5mFault CoreFault ZoneFigure 5.9: Borehole locations. The hole arrangement is drawn to scale, but the tunnel wallsand fault core are not. The boreholes are approximately 20 m deep.strongly affects memory usage and processing time. With five transmitters per hole we shouldexpect an improvement over the traveltime inversion results, but not quite to the same levelas the synthetic case due to the sparser ray coverage. For the diagonally adjacent holes, whereonly three transmitters are used per hole, the full-wave inversion algorithm may not significantlyimprove upon the traveltime inversion results.2D traveltime inversions [101] are performed for each pair of boreholes to find initial esti-mates of the relative permittivity. The results of the traveltime inversions are shown in Fig. 5.10.An initial estimate of the conductivity is calculated using radar data collected between holes 5and 2 using equation (5.31) and is found to be σstart = 2.0 mS/m.Fig. 5.11 shows the full-wave inversion results for permittivity and conductivity, plottedin the same arrangement. For clarity, the two sections between opposing holes are shownseparately in Fig. 5.12. For these hole pairs, where five transmitters per hole are used, asubstantial improvement in resolution is clearly visible when comparing the full-wave inversionresults to the traveltime inversion results.5.7.1 Data FitThe conductivity results for the inversions involving three transmitters per hole vs five trans-mitters per hole show inconsistent regions (horizontal streaks) of high loss in the vicinity of thetransmitters. Close inspection of Fig. 5.12 and Fig. 5.13 show this. The inconsistency in theinversion results indicates that the inversion algorithm is compensating for an inconsistencybetween our numerical model and the measured data. There are a number of areas where sucha discrepancy could arise, including the antenna model or the 3D to 2D transformation. In this1125.7. Inversion of Radar Data from LSBB 0 1 2 0 1 2 3 4 5-20-15-10-5 0 8 10 12 14 16 18 20x(m)y(m)z(m)RelativePermittivityFigure 5.10: Ray traveltime inversion results using the LSBB radar data set. Borehole locationsare at the four corners of the diamond. The four slices that run from corner to corner eachinvolve six transmitter locations. The two slices that run across the diamond involve tentransmitter locations. 0 1 2 0 1 2 3 4 5-20-15-10-5 0 8 10 12 14 16 18 20x(m)y(m)z(m)RelativePermittivity(a) 0 1 2 0 1 2 3 4 5-20-15-10-5 0 0 0.001 0.002 0.003 0.004 0.005x(m)y(m)z(m)Conductivity(S/m)(b)Figure 5.11: Full-wave inversion results using the LSBB radar data set: (a) relative permittivityand (b) conductivity. Borehole locations are at the four corners of the diamond. The four slicesthat run from corner to corner each involve six transmitter locations. The two slices that runacross the diamond involve ten transmitter locations.1135.7. Inversion of Radar Data from LSBB−21−20−19−18−17−16−15−14−13−12−11−10−9−8−7−6−5−4−3−2−10 0 1 2 3 4 5 0 1 2 3 4 5 8 1012141618200 1 2 3 4 5 0 1 2 3 4 5 x(m)x(m)x(m)z(m)RelativePermittivityRay ǫr FWI ǫr FWI σConductivity(mS/m)(a)−21−20−19−18−17−16−15−14−13−12−11−10−9−8−7−6−5−4−3−2−10 0 1 2 0 1 2 8 1012141618200 1 2 0 1 2 3 4 5 x(m)x(m)x(m)z(m)RelativePermittivityRay ǫr FWI ǫr FWI σConductivity(mS/m)(b)Figure 5.12: Full-wave inversion results for opposing hole pairs using the LSBB radar data set:(a) holes 5 and 2, and (b) holes 3 and 4.1145.7. Inversion of Radar Data from LSBBcase the discrepancy appears to be primarily due to signal compression and distortion in theradar receiver when the transmit and receive antennas are relatively close (i.e. directly acrossfrom each other). Signal compression results in a lower recorded received signal level than isexpected with a linear receiver and the inversion attempts to compensate for this by addinghorizontal streaks of high loss in the vicinity of the transmit antennas.A GPR receiver is essentially an antenna and amplifier followed by a digital to analogconverter. Nonlinearity in an amplifier can be modeled in the simplest case with a polynomial,vo(vi) = avi + bv2i + cv3i + · · · , (5.34)where vo is the output voltage of the amplifier and vi is the input voltage. Equation (5.34)ignores any frequency dependence that may be present. The presence of the square and cubicterms leads to intermodulation distortion and signal compression. To characterize the non-linear behavior of a receiver or amplifier, a pair of test sinusoids closely spaced in frequencyis usually applied to the input (see, e.g., [89, pp.93–99]). The strength of the intermodulationproducts caused by the square and cubic terms is then measured, and the associated distortionmetrics (e.g. second- and third-order intercept points) are calculated. In this case we nolonger have access to the receivers so we cannot carry out a typical receiver analysis on them.However, we can analyze the sampled output to determine whether or not there is significantharmonic distortion. Fig. 5.14 is a plot of the power at three frequencies (f1 = 157.4 MHz,f2 = 314.9 MHz, and f3 = 17.0 MHz). Normally, such a plot would be made with the inputpower level of one of those frequencies on the horizontal axis, but we do not have access to theinput power level. So in this case, we put the output power at f1 on the horizontal axis. Doingso, the plot of f1 is a straight line with a slope of unity, as expected. The power content at theother two frequencies appears to rise twice as fast (at roughly 2 dB/dB). The second harmonicof f1 is f2, and a 2 dB/dB rise in this signal is consistent with it being generated by thesquare term in equation (5.34). The lower frequency, f3, is where second-order intermodulationproducts of strong signals near f1 could be found. The 2 dB/dB increase is also consistent withthis. However, it must be remembered that the data for this plot were obtained by movingthe receiving antenna, and therefore there are other potential explanations for the behavior1155.7. Inversion of Radar Data from LSBB−21−20−19−18−17−16−15−14−13−12−11−10−9−8−7−6−5−4−3−2−10 0 1 2 8 1012141618200 1 2 0 1 2 3 4 5 x(m)x(m)z(m)RelativePermittivityFWI ǫr FWI σConductivity(mS/m)Figure 5.13: Recovered permittivity and conductivity for the 3-2 hole combination.seen in the plot. The plot could also be explained by frequency-dependent attenuation andnear-field effects (particularly for f3). Also, second order distortion does not lead directly tosignal compression (for a single frequency input), but third order distortion can. The presenceof such strong second order distortion suggests that third order distortion is also likely present.This, coupled with the streaks seen in the recovered conductivity images, strongly suggests thatsignificant signal compression is occurring.Despite the streaks in the recovered permittivity images, there is a region of higher con-ductivity (about 3 mS/m) that is visible and coincides with the higher permittivity region atmid-depth (visible in Fig. 5.12b between −16 m and −9 m). This is likely not an inversionartifact.The FWI l-curves for the measured data collected at LSBB do not have well defined kneeslike those visible in Fig. 5.6 for the synthetic data sets. A few example l-curves from themeasured data inversions are shown in Fig. 5.15. The curvature of the log-log plots of φd vs φm(l-curves) are calculated and a point is selected where there is a high positive curvature. Theregion of the curve where a suitable result may be found is located by looking at the inversionresults and narrowing the area of interest to where the results are not overly smooth or overly1165.7. Inversion of Radar Data from LSBB 50 60 70 80 90 100 110 120 90 95 100 105 110 115Output Power at f1 = 157.4 MHz (dB)OutputPower(dB)f1 = 157.4 MHzf2 = 314.9 MHzf3 = 17.0 MHzFigure 5.14: Power content in the received GPR signal at three frequencies with 7.3 MHzbandwidths plotted against the power content at 157.4 MHz. This is taken from the 3-4 holecombination with the transmitter at a depth of 15 m in hole 4 while the receiver is swept from13 m to 16 m in hole 3.rough, and then the curvature is used to pick β. At the selected points there is only a very slightdent in the curves. These l-curve plots also show that as β is reduced the inversion algorithm isable to continue finding better and better data fits with rougher and rougher models, unlike thesynthetic case where the l-curve abruptly levels off and little improvement in data fit is foundas the model is allowed to become under-regularized.This lack of clear corner in the l-curve suggests that there is a disconnect between thecomputer model and the measured data, and is consistent with the non-linear behavior of thereceiver discussed previously. Such a disconnect could also arise from the lack of a true 3Dmodel, or, rather that the rock structure is not actually 2.5D in nature. Another source of errorcould be the antenna model or errors in transmitter or receiver positions. Yet another could bethe assumption that the permittivity and conductivity are not frequency independent withinthe radar band being used. A simple test for dispersion involves dividing (in the frequencydomain) the signal at two receiver positions and then plotting the argument of that value withrespect to frequency. A linear plot indicates a lack of dispersion, and this is visible in Fig. 5.16.This plot verifies that permittivity is not frequency dependent and that the conductivity issmall enough that it does not significantly affect group velocity, but the conductivity itself maystill be frequency dependent. Anisotropy is another potential source of error.Fig. 5.17 compares a subset of the observed and predicted data at a single frequency for1175.7. Inversion of Radar Data from LSBB1001000100001000000.00010.001 0.01 0.1 1 10 100 1000φmφ dφ dHoles 5 and 2(a)1001000100001000000.00010.001 0.01 0.1 1 10 100 1000φmφ dφ dHoles 3 and 4(b)1001000100001000000.00010.001 0.01 0.1 1 10 100 1000φmφ dφ dHoles 3 and 2(c)1001000100001000000.00010.001 0.01 0.1 1 10 100 1000φmφ dφ dHoles 5 and 4(d)Figure 5.15: L-curves for the full-wave inversions of LSBB data: (a) holes 5 and 2; (b) holes 3and 4; (c) holes 3 and 2; (d) holes 5 and 4. The chosen β is indicated with a circle. The toprow of figures corresponds to hole combinations where five transmitter positions were used forthe data collection, and the bottom two rows are for two of the hole combinations where onlythree transmitter locations were used. 150 200 250 300 350 400 450 500 550 60 80 100 120 140 160 180 200 220Frequency (MHz)UnwrappedPhase(Deg)LSBB DataLinear FitFigure 5.16: To test for dispersion, the argument of the ratio of two received signals is plottedagainst frequency. The linear plot indicates negligible dispersion.1185.7. Inversion of Radar Data from LSBB00.20.40.60.811.20 5 10 15 20Receiver Depth (m)SignalAmplitudeSignalAmplitudeSignalAmplitudeSignalAmplitudeSignalAmplitudeHoles 5 and 2ObservedPredicted(a)-5000-4000-3000-2000-1000010000 5 10 15 20Receiver Depth (m)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)Holes 5 and 2ObservedPredicted(b)00.20.40.60.811.20 5 10 15 20Receiver Depth (m)SignalAmplitudeSignalAmplitudeSignalAmplitudeSignalAmplitudeSignalAmplitudeHoles 3 and 4ObservedPredicted(c)-5000-4000-3000-2000-1000010000 5 10 15 20Receiver Depth (m)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)SignalPhase(Deg)Holes 3 and 4ObservedPredicted(d)Figure 5.17: Received and predicted signal amplitudes and phases at 121 MHz for: (a) amplitudeand (b) phase for holes 5 and 2; (c) amplitude and (d) phase for holes 3 and 4. The phases areunwrapped and shifted such that the maximum phase is near zero degrees.the 5-2 and 3-4 hole combinations. For holes 3 and 4, which are close together, the amplitudeplots show a distinct depression in the signal strength where the transmitter and receiver aredirectly opposite from each other (i.e. where the signal is strong). This is likely an artifact ofthe signal compression discussed above. This behavior is much less severe for holes 5 and 2where the transmitter and receiver are further apart. The phase fit for both hole combinationsis extremely good where the signal is reasonably strong. These results are consistent with abetter recovery of the permittivity than the conductivity.1195.7. Inversion of Radar Data from LSBB5.7.2 Source WaveletIn order to improve the inversion results for the diagonally adjacent holes, where there areonly three transmitter locations in each hole, the recovered source wavelet from the 5-2 holecombination inversion is used as a starting guess. During the inversion, deviation from thisstarting guess is penalized through the regularization matrices, Wi. The results of the 5-2 holecombination inversion are used as a starting point for the three-transmitter inversions becauseof the large hole spacing that reduces the received signal strength and therefore the distortionin the receiver, and best approximates the far-field behavior that is assumed in the 3D to 2Dtransformation.Fig. 5.18a shows the progression of the source wavelet as the regularization parameter isdecreased for the 5-2 hole combination. The progression for the 5-3 hole combination is shownin Fig. 5.18b (in this case the chosen values from the 5-2 inversion are used as a starting point).In both cases the progression of the source wavelet is different than what was found in thesynthetic experiments (see Figs. 5.7 and 5.8). In the 5-2 case we do see the general increase inamplitude (distance from the origin) as β is decreased and the data misfit improves. For the5-3 case, the starting point is already very close to the final value, so little movement is seenbefore the model becomes under-regularized.Fig. 5.19 is a scatter plot showing the recovered source wavelet coefficients for each of thesix inversions. The close proximity of the recovered wavelet coefficients in the scatter plotsupport our use of a single source wavelet when inverting the data collected from each pairof boreholes. The data for this set were collected in the period of one day, thus reducing thelikelihood of strong variation between wavelets. A study by Belina et al. [16] also illustratesthat even if there is variation in source wavelets with different transmitter locations there isvery little negative impact on the inversion results, unless there is a very strong variation in thesource wavelet. Further, they suggest that attempting to find source wavelets separately foreach transmitter location may degrade the inversion results due to greater uncertainty (lowerredundancy) when determining the source wavelet. They also suggest that a strong materialcontrast near a transmitter location is likely to negatively influence results by strongly biasingthe source wavelet estimation for that transmitter, further supporting the use of a single source1205.7. Inversion of Radar Data from LSBB−10−5 0 5 10−4 −2 0 2 499.30 MHz121.10 MHz140.48 MHz159.86 MHz179.23 MHz201.03 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsR(a)−10−5 0 5 10−4 −2 0 2 499.30 MHz121.10 MHz140.48 MHz159.86 MHz179.23 MHz201.03 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsR(b)Figure 5.18: Recovered source wavelet coefficients for (a) the 5-2 hole combination and (b) the5-3 hole combination. The circles indicate the initial estimate of the source wavelet and thecrosses indicate the chosen values. The starting values for the 5-3 hole combination are thechosen values from the 5-2 hole combination.wavelet.5.7.3 Water ContentLastly, the Topp relation [99] is used to determine the approximate water content of the rock.The water content for the hole combinations 5-2 and 3-4 is shown in Fig. 5.20. Measurementsin the host rock [53, 52] indicate a porous layer (15–20 % porosity) intercalated between twonon-porous (< 5 % porosity) but highly fractured layers. The water content estimated with theTopp relation is approximately consistent with the rock being saturated with water, but thewater content may be slightly over estimated. The general conditions in the tunnel were wetand muddy when the measurements were taken, suggesting that we can expect a high watercontent. The non-porous layers are indicated with an estimated 15-20 % water content (aboveabout −9 m and below −15 m), which is higher than the porosity. This could be due to waterwithin the fractures. The estimated water content of the porous layer (from −15 m to about−9 m) is also higher than the porosity at 20-30 % water content. These high water estimatescould be due to an over-estimated electrical permittivity, or due to the approximate nature ofTopp relation when applied outside of its calibrated range (sandy loam to clay), or both.1215.7. Inversion of Radar Data from LSBB−10−5 0 5 10−4 −2 0 2 499.30 MHz121.10 MHz140.48 MHz159.86 MHz179.23 MHz201.03 MHzs Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is Is IsRFigure 5.19: Recovered source wavelet coefficients for all hole combinations. The circles indicatethe coefficients for each of the six frequencies.5−2−21−20−19−18−17−16−15−14−13−12−11−10−9−8−7−6−5−4−3−2−10 0 1 2 3 4 5 3−40 1 2 10.0015.0020.0025.0030.00x(m)x(m)z(m)Watercontent(%)Figure 5.20: Water content for the 5-2 and 3-4 hole combinations estimated using the Topprelation.1225.8. Concluding Remarks5.8 Concluding RemarksWe have presented and demonstrated the effective use of a frequency-domain full-wave inver-sion algorithm that simultaneously recovers permittivity, conductivity, and the source wavelet.The algorithm was applied with very promising results to both synthetic and measured datasets. The measured data set exhibits significant non-linear distortion that produces streaks inthe recovered conductivity and is also likely responsible for poor l-curve plots. Despite the dis-tortion, the permittivity is recovered well and water content estimated with the Topp relationcorrelates well with the rock porosity measured in the boreholes. The permittivity is stronglyconnected to the phase of the received signal, and is therefore much less sensitive to distortionin the radar receiver.A simple reciprocity-based antenna model that is easily implemented in an FDFD codewas also presented and verified numerically. The antenna model uses a current distributionthat was arrived at through analogy with a coaxial transmission line. Simulations indicatedthat the current distribution provides realistic radiation characteristics for an antenna in anair-filled borehole, and through reciprocity the current distribution can also be used for modelthe reception characteristics.Calculating the full Jacobian is a computationally intensive task. Performing the calcula-tions in the frequency domain allows re-use of the factorization of the forward model to greatlyreduce the required computation time. Finding the Jacobian allows the use of a Gauss-Newtoniterative solution in which gradual cooling (reduction of the regularization parameter) can beused to help guide the solution to first find a coarse (low spatial frequency) model and thenprogressively include higher spatial frequency components.5.9 AcknowledgementsFunding for this work comes from a number of sources, and the authors would like to thankall of them. ANR “Captage de CO2” has provided funding through the HPPP-CO2 project,and the PACA county through the PETRO-PRO project. An in-kind donation of field workwas provided by Golder Inc. MALA˚ Geoscience provided two 250 MHz borehole antennas asin-kind support. LSBB and the re´gion PACA provided the infrastructure. The project ANR1235.9. AcknowledgementsMAXWELL provided direct funds for travel and per diem.Substantial support through scholarship and research grants has been provided by TheNatural Sciences and Engineering Research Council of Canada (NSERC).Calculations were performed using MATLAB, GNU Octave, MUMPS [3, 4], and MEEP[73].124Chapter 6Summary and Concluding RemarksIn this thesis I have presented my research on the inversion of cross-hole GPR data, which hasculminated in an effective 2D frequency-domain full-waveform inversion algorithm that works ona single modern desktop computer. This is possible largely due to the continued advancementof computer hardware and software, particularly as they relate to the solution of large linearalgebra problems. I have contributed to the field of GPR inversion in a number of importantareas. In the following paragraphs, I summarize these contributions.In Chapter 2, I presented a 3D traveltime inversion method for the recovery of electricalpermittivity (derived from slowness) from cross-hole GPR data. The methodology could alsobe applied to acoustic measurements. It is based on my previous 2D traveltime inversionresearch [101]. These traveltime inversion algorithms are based on the eikonal equation, forwhich Zhao has recently published a fast, efficient, accurate, and easily implemented solver[124]. My contribution related to this is the presentation and implementation of an efficientand particularly easy to implement 3D traveltime inversion algorithm.The eikonal equation is an infinite-frequency approximation of Maxwell’s equations, and assuch leads to a disconnect between measured traveltimes, which are taken at finite frequency,and the predicated traveltimes. This disconnect applies to any eikonal (or ray-based) inversionalgorithm and causes problems with determining the regularization parameter using the l-curveor GCV, which I illustrate with a synthetic example in section 2.6. As a solution to thisproblem, I have extended the work of Lukas on robust GCV and strong robust GCV [62, 63] tothe non-linear inversion problem. These functions, particularly strong robust GCV, appear tobe less sensitive to the correlated data errors that result from the eikonal approximation. Theyproduce potential functions with well-defined minima.When considering the application of a 2D full-waveform inversion algorithm to measured125Chapter 6. Summary and Concluding RemarksGPR data, one must reconcile the differences between 2D wave propagation in the computermodel and the real-world 3D wave propagation. In uniform media in the frequency domainin 2D the waves can be described with a Hankel function, and in 3D with a simple complexexponential with r−1 amplitude decay due to spreading. When we are interested in performinga full-waveform inversion, the media is most likely non-uniform and we have, presumably, theresults of a traveltime inversion at hand. Even though the traveltime results are arrived atthrough a ray approximation, they an be used effectively to increase the accuracy of a 3D to2D transformation beyond the uniform media assumption. This improves the accuracy of thefull-waveform inversion that follows.When transforming frequency domain data from 3D to 2D, it is customary to use the rayapproximation (see, e.g., [19, 6, 118]). One of the terms that appears in the 3D to 2D transferfunction is the ratio of 3D to 2D ray Jacobians and I refer to this as the perpendicular rayJacobian. The ray Jacobian describes the spreading of the ray due to perturbations of thetakeoff angles, and its derivation is generally quite complex (see, e.g., [22]). In Chapter 3,I provide a graphical derivation of the perpendicular ray Jacobian. The derivation clearlyand directly demonstrates that the 3D Jacobian can be factorized into a parallel (2D) and aperpendicular portion and finds an expression for the perpendicular portion. The derivationprovides the reader with an intuitive understanding of the ray Jacobian that is extremelydifficult to obtain through the traditional derivation methods.In Chapter 3, I also present a paraxial derivation of the perpendicular ray Jacobian that,although not as intuitive as the graphical derivation, provides an analytic verification of thegraphical results. The paraxial derivation is focused in particular on the perpendicular rayJacobian and bypasses a great deal of the complexity that is found in the work of Cˇerveny´ [22].This too greatly aids the reader in understanding the ray Jacobian.Although 3D to 2D transformations based on acoustic models have been used successfully forthe cross-hole GPR problem, a thorough derivation of such a transfer function for that case hasnot been performed (to my knowledge). In Chapter 4, I demonstrated that the acoustic transferfunction arrived at by Bleistein [19], and that we arrived at using a more direct approach [121],is directly applicable to the GPR problem. These derivations are based on the WKB (ray)approximation and thus are only strictly valid when the length scale of material variation is126Chapter 6. Summary and Concluding Remarkslarge compared to a wavelength. In Chapter 4, I also presented synthetic FWI simulation resultsthat illustrate the effectiveness of two 3D to 2D transfer functions (the one I derive and one byErnst et al. [33]), and it is clear from the results that even though the true model violates theWKB approximation and the traveltime inversion results with which the transfer functions arecomputed do not resolve the fine detail of the model, the FWI results show much of the detailthat is seen in the true model (particularly for electrical permittivity) and accurately find theconductivity. This illustrates that the transfer function is highly effective for FWI even thoughits derivation is in the ray regime.The transfer functions that I discuss in Chapter 4 apply a phase shift of very close to 45◦ fora low-loss medium, and a strongly varying amplitude adjustment that is a function of the raypath from source to receiver. Any constant complex multiplier that is applied to the frequencydomain data gets lumped into the source wavelet in the FWI algorithm that I used to producethe synthetic test results in Chapter 4 (and derived in Chapter 5). This implies that thetransfer function is of particular importance for recovering conductivity, which I demonstrateby inverting a synthetic 3D data set without applying a transfer function. The results showthat the recovered permittivity is hardly affected, while recovered conductivity is significantlyoverestimated.In Chapter 5, I presented a frequency-domain FWI algorithm that simultaneously recov-ers conductivity, permittivity, and the source wavelet. It is a Gauss-Newton formulation thatinverts multiple frequencies simultaneously and employs a cooling strategy to help steer theinversion result closer to the global minimum of the objective function. With both syntheticand measured data, the inversion algorithm is able to provide a substantial improvement inrecovered detail over what is obtained with traveltime inversion. The field data were unfortu-nately corrupted by signal compression within the receiver, leading to streaks in the recoveredconductivity images.The FWI algorithm ties together the research that is presented in the previous chapters. Ituses a 2D version of the traveltime inversion algorithm that is presented in Chapter 2 to generatethe starting model. The expression for the perpendicular ray Jacobian that is developed inChapter 3 and the transfer function derived in Chapter 4 are used to convert both syntheticand measured 3D data into 2D equivalent data.1276.1. Future WorkEarly in Chapter 5, I developed an FDFD solution to Maxwell’s equations with PML absorb-ing boundary conditions and gave a matrix representation for the resulting coefficient matrix,A. An LDLT factorization is used to dramatically speed up the many applications of A−1that must be performed during the inversion, both during the forward modeling operation andwhen calculating the Jacobian.The matrix formulation leads to a straightforward (albeit somewhat lengthy) formulation forfinding the Jacobian matrix. The Jacobian calculation involves back-propagation of a modifiedform of the measurement operator and this result can be re-used where the same measurementoperator is used for portions of the data, thus leading to a significant savings in processingtime.Another contribution of Chapter 5 is a simplified antenna model for borehole GPR FWIproblems. This new model is an improvement over the infinitesimal dipole model that is com-monly used for GPR FWI. The model approximates the radiation characteristics of a trueantenna, and because conductivity inversion relies on accurate signal amplitudes this modelshould improve conductivity inversion results over those obtained with the simple dipole model.Lastly, the inclusion of the source wavelet in the objective function and Jacobian of the FWIformulation is new to GPR FWI and has proven, so far, to be a reliable method of determiningthe source wavelet. This type of formulation allows the estimated source wavelet to evolve alongwith the model to match the measured data set.The research that has been presented in Chapters 2 through 5 amount to a complete cross-hole GPR imaging solution that starts with an extremely fast but low-resolution traveltimeinversion and ends with a high-resolution full-waveform inversion. The algorithms that areoutlined have been tested with both synthetic and measured data and give results in line withwhat is expected from these kinds of imaging modalities.6.1 Future WorkThere are a number of avenues for future research related to the research presented in this thesis.Of most interest is that related to the FWI algorithm, largely because of its high resolutionand the continued increase in performance of personal computers that allows faster and faster1286.1. Future Worksolution of large problems. Extension of the FWI algorithm to 3D is an attractive direction forfuture research. The mathematical formulation itself could easily be extended to 3D, includingthe simplified antenna model. The difficulty comes with the size of the problem. In 3D, evensolving the forward problem, as currently presented, presents a challenge.The current implementation of the 2D FWI algorithm uses MUMPS to factorize and solvethe forward problem, as well as perform part of the Jacobian calculation. MUMPS has thecapability to efficiently use a cluster to carry out the factorization of large sparse matriceslike the one used for the forward model, so this may suffice for solving a 3D version of theforward problem. The use of an iterative solver is most likely out of the question because ofthe large number of solutions that must be found when calculating the Jacobian. The forwardmodeling operation is also poorly conditioned and an appropriate conditioner would need to beimplemented in order for it to be solved iteratively at all. It may be necessary to explore otherforward modeling options, such as the finite element method.The calculation of the Jacobian presents a slightly different problem. As currently imple-mented, the Jacobian calculation is localized to a single processor. When a set of vectors ispassed to MUMPS to find a set of solutions, each vector is distributed across the working nodesfor solution and the results are sent back to the host processor. When calculating the Jacobian,this is a substantial amount of network traffic and an experiment has shown that adding anothermachine (connected with a single 1 Gbps ethernet link) actually slows down the solution of theinverse problem as compared to using a single machine. To remedy this, the entire Jacobiancalculation would need to be distributed across the compute cluster, rather than performingit on a single host. Most of the matrices involved in the Jacobian calculation are sparse, soduplicating them on each working processor (which might be required) should not occupy anexcessive amount of memory.The solution of the Gauss-Newton model update in equation (1.8) is another potentialproblem for moving to 3D. In 2D, I have noticed that iterative solvers have difficulty calculatingthe correct model update. This could be related to the inclusion of the source wavelet and/orthe relative scaling of the permittivity, conductivity, and source wavelet components of theJacobian and model vector. An iterative solution might be desirable because then the JTJproduct would not need to be formed explicitly. Additionally, JTJ is not sparse, so using1296.1. Future WorkMUMPS to solve for the model update would likely be unhelpful.Other possible directions for future research include:• Full automation of the selection of the regularization parameter for the FWI algorithm.The application of the R1GCV method to the FWI problem might be a good way of han-dling this, but as it is currently implemented the calculation time is quite long. Perhapsapplying the algorithm to a subset of the data is an option.• Decompression of the LSBB data. The LSBB data that I collected are compressed (dis-torted) where the received signal is strong and this led to a corrupt recovered conductivityin the FWI. It might be possible to preprocess the measured data to remove some of thedistortion and obtain a better conductivity image. This would entail somehow charac-terising the GPR receiver nonlinearity using the measured data that I currently have (orgaining access to the same receiver unit for further characterization).• Exploration of other techniques that may be used to speed up the FWI algorithm. Matrixprobing techniques to estimate the Hessian, or the L-BFGS-B algorithm that Lavoue´ etal. [60] use are possibilities. The goal here would be to speed up the inversion withoutadversely affecting the reliability of the Gauss-Newton method.• Adaptation of the FWI technique to other areas, such as medical imaging, through-the-wall imaging, or the imaging of building structures such as the internals of walls withoutcutting holes.There are a great many directions for future work that I or someone else could pursueusing the research presented here and the many excellent references given as a basis. Althoughresearch into FWI started in the 1980s, continued advancement in computer technology nowmakes FWI possible without the need for a large compute cluster. And the potential with largecompute clusters seems almost unlimited, as long as the problem can be efficiently distributedand solved in a distributed manner. There is potential for FWI-type algorithms in many areas,and we are just beginning to scratch the surface.130Bibliography[1] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions with Formu-las, Graphs, and Mathematical Tables. 10th printing December 1972. US GovernmentPrinting Office, 1964.[2] D. F. Aldridge and D. W. Oldenburg. “Two-dimensional tomographic inversion withfinite-difference traveltimes”. In: Journal of Seismic Exploration 2.3 (1993), pp. 257–274.[3] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster. “A fully asynchronous mul-tifrontal solver using distributed dynamic scheduling”. In: SIAM Journal on MatrixAnalysis and Applications 23.1 (2001), pp. 15–41.[4] P. R. Amestoy, A. Guermouche, J.-Y. LExcellent, and S. Pralet. “Hybrid scheduling forthe parallel solution of linear systems”. In: Parallel Computing 32.2 (2006), pp. 136–156.[5] D. Amodei, H. Keers, D. Vasco, and L. Johnson. “Computation of uniform wave formsusing complex rays”. In: Physical Review E 73.3 (2006), pp. 036704–1–14.[6] L. Auer, A. M. Nuber, S. A. Greenhalgh, H. Maurer, and S. Marelli. “A critical appraisalof asymptotic 3D-to-2D data transformation in full-waveform seismic crosshole tomog-raphy”. In: Geophysics 78.6 (Oct. 2013), R235–R247. issn: 0016-8033, 1942-2156. doi:10.1190/geo2012-0382.1.[7] G. S. S. Avila and J. B. Keller. “The high-frequency asymptotic field of a point sourcein an inhomogeneous medium”. In: Communications on Pure and Applied Mathematics16.4 (1963), pp. 363–381.[8] V. M. Babich. “On asymptotic solutions of the problems of moving sources of diffusionand oscillations”. In: Journal of Soviet Mathematics 57.3 (Nov. 1991), pp. 3067–3071.doi: 10.1007/BF01098969.131BIBLIOGRAPHY[9] V. M. Babich. “The higher-dimensional WKB method or ray method”. In: Encyclopediaof Mathematical Sciences 34 (1997), pp. 91–131.[10] V. M. Babich. “The short wave asymptotic form of the solution for the problem of a pointsource in an inhomogeneous medium (english)”. In: USSR Compututational Mathematicsand Mathematical Physics 5.5 (1965), pp. 247–251.[11] G. E. Backus and J. F. Gilbert. “Numerical applications of a formalism for geophysicalinverse problems”. In: Geophysical Journal International 13.1-3 (1967), pp. 247–276.[12] R. J. Baker. “High voltage pulse generation using current mode second breakdown in abipolar junction transistor”. In: Review of Scientific Instruments 62.4 (1991), pp. 1031–1036.[13] C. A. Balanis. Advanced Engineering Electromagnetics. Wiley, 1989. isbn: 9780471621942.[14] C. A. Balanis. Antenna Theory: Analysis and Design. third. John Wiley & Sons, 2005.[15] F. Belina, J. Irving, J. Ernst, and K. Holliger. “Analysis of an iterative deconvolution ap-proach for estimating the source wavelet during waveform inversion of crosshole georadardata”. In: Journal of Applied Geophysics 78 (2012), pp. 20–30.[16] F. A. Belina, J. Irving, J. R. Ernst, and K. Holliger. “Waveform Inversion of CrossholeGeoradar Data: Influence of Source Wavelet Variability and the Suitability of a SingleWavelet Assumption”. In: IEEE Transactions on Geoscience and Remote Sensing 50.11(2012), pp. 4610–4625.[17] D. M. Benzel and M. D. Pocha. “1000-V, 300-ps pulse-generation circuit using siliconavalanche devices”. In: Review of Scientific Instruments 56.7 (1985), pp. 1456–1458.[18] J.-P. Berenger. “A perfectly matched layer for the absorption of electromagnetic waves”.In: Journal of Computational Physics 114.2 (1994), pp. 185–200.[19] N. Bleistein. “TWO-AND-ONE-HALF DIMENSIONAL IN-PLANE WAVE PROPA-GATION”. In: Geophysical Prospecting 34.5 (1986), pp. 686–703.[20] C. Bunks, F. Saleck, S. Zaleski, and G. Chavent. “Multiscale seismic waveform inversion”.In: Geophysics 60.5 (Sept. 1995), pp. 1457–1473. issn: 0016-8033. doi: 10.1190/1.1443880.132BIBLIOGRAPHY[21] J. J. M. Carcione. Wave Fields in Eeal Media: Wave Propagation in Anisotropic, Anelas-tic and Porous Media. Vol. 31. Pergamon, Amsterdam, 2001.[22] V. Cˇerveny´. Seismic Ray Theory. Cambridge University Press, Cambridge, 2001.[23] V. Cˇerveny´ and F. Hron. “The ray series method and dynamic ray tracing system forthree-dimensional inhomogeneous media”. In: Bulletin of the Seismological Society ofAmerica 70.1 (1980), pp. 47–77.[24] V. Cˇerveny´ and R. Ravindra. Theory of Seismic Head Waves. University of TorontoPress, 1971.[25] C. H. Chapman. Fundamentals of Seismic Wave Propagation. Cambridge UniversityPress, Cambridge, 2004.[26] W. C Chew and W. H Weedon. “A 3-D perfectly matched medium from modifiedMaxwell’s equations with stretched coordinates”. In: Microwave and Optical Technol-ogy Letters 7.13 (1994), pp. 599–604.[27] K. S. Cordua, T. M. Hansen, and K. Mosegaard. “Monte Carlo full-waveform inver-sion of crosshole GPR data using multiple-point geostatistical a priori information”. In:Geophysics 77.2 (2012), H19–H31.[28] D. J. Daniels, ed. Ground Penetrating Radar. 2nd. Vol. 15. IET Radar, Sonar, Navigationand Avionics Series. The Institute of Engineering and Technology, 2007.[29] J. J. E. Dennis and R. B. Schnabel. Numerical Methods for Unconstrained Optimiza-tion and Nonlinear Equations. Classics in Applied Mathematics. SIAM, 1983. isbn:9781611971200.[30] S. M. Deregowski and S. M. Brown. “A theory of acoustic diffractors applied to 2-Dmodels”. In: Geophysical Prospecting 31.2 (1983), pp. 293–333.[31] K. J. Ellefsen, A. T. Mazzella, R. J. Horton, and J. R. McKenna. “Phase and amplitudeinversion of crosswell radar data”. In: Geophysics 76.3 (2011), J1–J12.[32] A. L. Endres and R. Knight. “A theoretical treatment of the effect of microscopic fluiddistribution on the dielectric properties of partially saturated rocks”. In: GeophysicalProspecting 40.3 (1992), pp. 307–324.133BIBLIOGRAPHY[33] J. R. Ernst, A. G. Green, H. Maurer, and K. Holliger. “Application of a new 2D time-domain full-waveform inversion scheme to crosshole radar data”. In: Geophysics 72.5(2007), J53–J64.[34] J. R. Ernst, K. Holliger, and H. Maurer. “FDTD modeling of borehole georadar data”.In: Proceedings of the Tenth International Conference on Ground Penetrating Radar,2004. GPR 2004. Vol. 1. June 2004, pp. 229–232.[35] J. R. Ernst, H. Maurer, A. G. Green, and K. Holliger. “Full-Waveform Inversion of Cross-hole Radar Data Based on 2-D Finite-Difference Time-Domain Solutions of Maxwell’sEquations”. In: IEEE Transactions on Geoscience and Remote Sensing 45.9 (2007),pp. 2807–2828.[36] C. G. Farquharson and D. W. Oldenburg. “A comparison of automatic techniques forestimating the regularization parameter in non-linear inverse problems”. In: GeophysicalJournal International 156.3 (2004), pp. 411–425.[37] C. G. Farquharson and D. W. Oldenburg. “Inversion of time-domain electromagneticdata for a horizontally layered Earth”. In: Geophysical Journal International 114.3(1993), pp. 433–442.[38] V. Farra and R. Madariaga. “Seismic waveform modeling in heterogeneous media by rayperturbation theory”. In: Journal of Geophysical Research 92.B3 (1987), pp. 2697–2712.[39] R. P. Feynman, R. B. Leighton, and M. Sands. The Feynman Lectures on Physics,The Definitive Edition Volume 1 (Feynman Lectures on Physics (Hardcover)). AddisonWesley, 2006.[40] A. Fichtner, F. Bleibinhaus, and Y. Capdeville. Full Seismic Waveform Modelling andInversion. Springer, 2011.[41] H. T. Friis. “A note on a simple transmission formula”. In: Proceedings of the IRE 34.5(1946), pp. 254–256.[42] W. Green. “Inversion of gravity profiles by use of a Backus-Gilbert approach”. In: Geo-physics 40.5 (Oct. 1975), pp. 763–772.134BIBLIOGRAPHY[43] Y. Guglielmi, F. Cappa, H. Lanc¸on, J. B. Janowczyk, J. Rutqvist, C. F. Tsang, andJ. S. Y. Wang. “ISRM Suggested Method for Step-Rate Injection Method for FractureIn-Situ Properties (SIMFIP): Using a 3-Components Borehole Deformation Sensor”.In: The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring:2007–2014. Ed. by R. Ulusay. Springer International Publishing, 2015, pp. 179–186.[44] E. Haber and D. Oldenburg. “A GCV based method for nonlinear ill-posed problems”.In: Computational Geosciences 4.1 (2000), pp. 41–63.[45] P. C. Hansen. “Analysis of discrete ill-posed problems by means of the L-curve”. In:SIAM Review 34.4 (1992), pp. 561–580.[46] P. C. Hansen. The L-curve and its use in the Numerical Treatment of Inverse Problems.IMM, Department of Mathematical Modelling, Technical University of Denmark, 1999.[47] J. A. Hole and B. C. Zelt. “3-D finite-difference reflection travel times”. In: GeophysicalJournal International 121.2 (1995), pp. 427–434.[48] K. Holliger and T. Bergmann. “Numerical modeling of borehole georadar data”. In:Geophysics 67.4 (July 2002), pp. 1249–1257. issn: 0016-8033, 1942-2156. doi: 10.1190/1.1500387.[49] K. Holliger, M. Musil, and H. Maurer. “Ray-based amplitude tomography for crossholegeoradar data: a numerical assessment”. In: Journal of Applied Geophysics 47.3 (2001),pp. 285–298.[50] M. F. Hutchinson. “A stochastic estimator of the trace of the influence matrix for Lapla-cian smoothing splines”. In: Communications in Statistics-Simulation and Computation19.2 (1990), pp. 433–450.[51] P. Jeanne, Y. Guglielmi, and F. Cappa. “Dissimilar properties within a carbonate-reservoir’s small fault zone, and their impact on the pressurization and leakage asso-ciated with CO2 injection”. In: Journal of Structural Geology 47 (2013), pp. 25–35. issn:0191-8141. doi: http://dx.doi.org/10.1016/j.jsg.2012.10.010.135BIBLIOGRAPHY[52] P. Jeanne, Y. Guglielmi, and F. Cappa. “Multiscale seismic signature of a small faultzone in a carbonate reservoir: relationships between VP imaging, fault zone architectureand cohesion”. In: Tectonophysics 554-557 (2012), pp. 185–201. issn: 0040-1951. doi:http://dx.doi.org/10.1016/j.tecto.2012.05.012.[53] P. Jeanne, Y. Guglielmi, J. Lamarche, F. Cappa, and L. Marie´. “Architectural charac-teristics and petrophysical properties evolution of a strike-slip fault zone in a fracturedporous carbonate reservoir”. In: Journal of Structural Geology 44 (2012), pp. 93–109.issn: 0191-8141. doi: http://dx.doi.org/10.1016/j.jsg.2012.08.016.[54] H. Jia, T. Takenaka, and T. Tanaka. “Time-domain inverse scattering method for cross-borehole radar imaging”. In: IEEE Transactions on Geoscience and Remote Sensing40.7 (2002), pp. 1640–1647.[55] E. C. Jordan and K. G. Balmain. Electromagnetic Waves and Radiating Systems. Second.Prentice-Hall, 1968.[56] F. C. Karal Jr and J. B. Keller. “Elastic wave propagation in homogeneous and inho-mogeneous media”. In: The Journal of the Acoustical Society of America 31.6 (1959),pp. 694–705.[57] A. Klotzsche, J. van der Kruk, N. Linde, J. Doetsch, and H. Vereecken. “3-D characteriza-tion of high-permeability zones in a gravel aquifer using 2-D crosshole GPR full-waveforminversion and waveguide detection”. In: Geophysical Journal International 195.2 (2013),pp. 932–944.[58] A. Klotzsche, J. van der Kruk, G. A. Meles, J. Doetsch, H. Maurer, and N. Linde. “Full-waveform inversion of cross-hole ground-penetrating radar data to characterize a gravelaquifer close to the Thur River, Switzerland”. In: Near Surface Geophysics 8.6 (2010),pp. 635–649.[59] S. Kuroda, M. Takeuchi, and H. J. Kim. “Full-waveform inversion algorithm for in-terpreting crosshole radar data: a theoretical approach”. In: Geosciences Journal 11.3(2007), pp. 211–217.136BIBLIOGRAPHY[60] F. Lavoue´, R. Brossier, L. Me´tivier, S. Garambois, and J. Virieux. “Two-dimensional per-mittivity and conductivity imaging by full waveform inversion of multioffset GPR data: afrequency-domain quasi-Newton approach”. In: Geophysical Journal International 197.1(Apr. 2014), pp. 248–268. issn: 0956-540X, 1365-246X. doi: 10.1093/gji/ggt528.[61] S. Leung and J. Qian. “An adjoint state method for three-dimensional transmission trav-eltime tomography using first-arrivals”. In: Communications in Mathematical Sciences4.1 (2006), pp. 249–266.[62] M. A. Lukas. “Robust generalized cross-validation for choosing the regularization pa-rameter”. In: Inverse problems 22.5 (2006), pp. 1883–1902.[63] M. A. Lukas. “Strong robust generalized cross-validation for choosing the regularizationparameter”. In: Inverse Problems 24.3 (2008), p. 034006.[64] D. Marcuse. Light Transmission Optics. Electrical/Computer Science and EngineeringSeries. Van Nostrand Reinhold, 1982. isbn: 0-442-26309-0.[65] H. Maurer, S. Greenhalgh, and S. Latzel. “Frequency and spatial sampling strategies forcrosshole seismic waveform spectral inversion experiments”. In: Geophysics 74.6 (2009),WCC79–WCC89.[66] G. Meles, S. Greenhalgh, J. Van der Kruk, A. Green, and H. Maurer. “Taming the non-linearity problem in GPR full-waveform inversion for high contrast media”. In: Journalof Applied Geophysics 73.2 (2011), pp. 174–186.[67] G. A. Meles, J. Van der Kruk, S. A. Greenhalgh, J. R. Ernst, H. Maurer, and A. G.Green. “A new vector waveform inversion algorithm for simultaneous updating of con-ductivity and permittivity parameters from combination crosshole/borehole-to-surfaceGPR data”. In: Geoscience and Remote Sensing, IEEE Transactions on 48.9 (2010),pp. 3391–3407.[68] W. A. Mulder and R. E. Plessix. “Exploring some issues in acoustic full waveform in-version”. In: Geophysical Prospecting 56.6 (2008), pp. 827–841.[69] G. Nolet. A Breviary of Seismic Tomography. Cambridge University Press Cambridge,UK, 2008.137BIBLIOGRAPHY[70] M. Oberro¨hrmann, A. Klotzsche, H. Vereecken, and J. van der Kruk. “Optimization ofacquisition setup for cross-hole GPR full-waveform inversion using checkerboard analy-sis”. In: Near Surface Geophysics 11.2 (2013), pp. 197–209.[71] D. Oldenburg. “EOSC 550 Linear Inverse Theory”. A graduate course offered at theUniversity of British Columbia. 2009.[72] D. W. Oldenburg. “An introduction to linear inverse theory”. In: IEEE Transactions onGeoscience and Remote Sensing GE-22.6 (Nov. 1984), pp. 665–674.[73] A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. John-son. “MEEP: A flexible free-software package for electromagnetic simulations by theFDTD method”. In: Computer Physics Communications 181.3 (2010), pp. 687–702.[74] M. J. Oyan, S.-E. Hamran, L. Hanssen, T. Berger, and D. Plettemeier. “Ultrawidebandgated step frequency ground-penetrating radar”. In: IEEE Transactions on Geoscienceand Remote Sensing 50.1 (2012), pp. 212–220.[75] C. C. Paige and M. A. Saunders. “LSQR: An Algorithm for Sparse Linear Equationsand Sparse Least Squares”. In: ACM Transactions on Mathematical Software 8.1 (1982),pp. 43–71.[76] G. S. Pan, R. A. Phinney, and R. I. Odom. “Full-waveform inversion of plane-wave seis-mograms in stratified acoustic media: theory and feasibility”. In: Geophysics 53 (1988),p. 21.[77] R. L. Parker. Geophysical Inverse Theory. Princeton University Press, 1994.[78] P. C. Pawliuk. “Evanescent field interactions in the scattering from cylinders with ap-plications in super-resolution”. PhD thesis. The University of British Columbia, 2013.[79] A. Pica, J. P. Diet, and A. Tarantola. “Nonlinear inversion of seismic reflection data ina laterally invariant medium”. In: Geophysics 55.3 (1990), pp. 284–292.[80] P. Podvin and I. Lecomte. “Finite difference computation of traveltimes in very con-trasted velocity models: a massively parallel approach and its associated tools”. In:Geophysical Journal International 105.1 (1991), pp. 271–284.138BIBLIOGRAPHY[81] R. G. Pratt. “Seismic waveform inversion in the frequency domain, Part 1: theory andverification in a physical scale model”. In: Geophysics 64.3 (1999), pp. 888–901.[82] R. G. Pratt, C. Shin, and G. J. Hicks. “Gauss-Newton and full Newton methods infrequency-space seismic waveform inversion”. In:Geophysical Journal International 133.2(1998), pp. 341–362.[83] R. G. Pratt and R. M. Shipp. “Seismic waveform inversion in the frequency domain,Part 2: fault delineation in sediments using crosshole data”. In: Geophysics 64.3 (1999),pp. 902–914.[84] K. Pruess, C. Oldenburg, and G. Moridis. TOUGH2 users guide, version 2.0. Tech. rep.LBNL-43134. Berkeley, Calif.: Lawrence Berkeley National Laboratory, 1999.[85] C. M. Rappaport, M. Kilmer, and E. Miller. “Accuracy considerations in using thePML ABC with FDFD Helmholtz equation computation”. In: International Journal ofNumerical Modelling: Electronic Networks, Devices and Fields 13.5 (2000), pp. 471–482.[86] J. A. Sethian. Level Set Methods and Fast Marching Methods: Evolving Interfaces inComputational Geometry, Fluid Mechanics, Computer Vision, and Materials Science.Vol. 3. Cambridge University Press, 1999.[87] W. Shin and S. Fan. “Choice of the perfectly matched layer boundary condition forfrequency-domain Maxwells equations solvers”. In: Journal of Computational Physics231.8 (2012), pp. 3406–3431.[88] L. Sirgue and R. G. Pratt. “Efficient waveform inversion and imaging: a strategy forselecting temporal frequencies”. In: Geophysics 69.1 (2004), pp. 231–248.[89] J. Smith. Modern Communication Circuits. Second. McGraw-Hill, Inc., 1998.[90] G. F. Stickley, D. A. Noon, M. Cherniakov, and I. D. Longstaff. “Gated stepped-frequency ground penetrating radar”. In: Journal of Applied Geophysics 43.2 (2000),pp. 259–269.[91] A. Taflove and S. C. Hagness. Computational Electrodynamics: The Finite-DifferenceTime-Domain Method. 2nd. Artech House, Incorporated, Jan. 2000. isbn: 9781580530767.139BIBLIOGRAPHY[92] T. Takenaka, H. Jia, and T. Tanaka. “Microwave imaging of electrical property distri-butions by a forward-backward time-stepping method”. In: Journal of ElectromagneticWaves and Applications 14.12 (2000), pp. 1609–1626.[93] A. Tarantola. Inverse Problem Theory and Methods for Model Parameter Estimation.SIAM, 2005.[94] A. Tarantola. “Inversion of seismic reflection data in the acoustic approximation”. In:Geophysics 49.8 (1984), 12591266.[95] A. Tarantola and B. Valette. “Generalized nonlinear inverse problems solved usingthe least squares criterion”. In: Reviews of Geophysics and Space Physics 20.2 (1982),pp. 219–232.[96] F. L. Teixeira and W. C. Chew. “General closed-form PML constitutive tensors to matcharbitrary bianisotropic and dispersive linear media”. In: IEEE Microwave and GuidedWave Letters 8.6 (1998), pp. 223–225.[97] A. N. Tikhonov, A. Goncharsky, V. V. Stepanov, and A. G. Yagola. Numerical Meth-ods for the Solution of Ill-Posed Problems. Vol. 328. Mathematics and its Applications.Springer, 1995.[98] I. Tolstoy. “The WKB approximation, turning points, and the measurement of phasevelocities”. In: The Journal of the Acoustical Society of America 52 (1972), p. 356.[99] G. C. Topp, J. L. Davis, and A. P. Annan. “Electromagnetic determination of soil watercontent: Measurements in coaxial transmission lines”. In: Water Resources Research 16.3(1980), pp. 574–582. issn: 0043-1397.[100] D. Van Vorst, J. Virieux, and M. J. Yedlin. “A Graphical Derivation of the Perpendic-ular Ray Jacobian”. In: Bulletin of the Seismological Society of America 102.6 (2012),pp. 2452–2457.[101] D. Van Vorst, M. Yedlin, Y. Guglielmi, F. Cappa, S. Gaffet, M. Maxwell, and D. Olden-burg. “GPR imaging of a fracture zone in the Vaucluse Karst Aquifer using 2D eikonalinversion”. In: EDP Sciences, Mar. 2011, 01002–p.1–01002–p.8. isbn: 978-2-7598-0630-0.doi: 10.1051/idust/201101002.140BIBLIOGRAPHY[102] D. Van Vorst and M. J. Yedlin. “Full-Waveform Inversion of Ground Penetrating RadarData in the Frequency Domain: Simultaneous Recovery of Permittivity, Conductivity,and Source Wavelet”. In: Geophysical Journal International (Submitted 2014).[103] D. Van Vorst, M. J. Yedlin, Y. Guglielmi, F. Cappa, S. Gaffet, M. Maxwell, and D.Oldenburg. “GPR Imaging of a Fracture Zone in the Vaucluse Karst Aquifer Using 2Dand 3D Eikonal Inversion”. In: i-DUST Conference. Extended abstract. 2010. url: http://www.i-dust.org/doc_journal/doc/idust/2010/iDUST2010_A11_abstract.pdf.[104] D. Van Vorst, M. J. Yedlin, J. Virieux, and E. S. Krebes. “Three-dimensional to two-dimensional data conversion for electromagnetic wave propagation using an acoustictransfer function: application to cross-hole GPR data”. In: Geophysical Journal Inter-national 198.1 (July 2014), pp. 474–483.[105] V. Vavrycˇuk. “Real ray tracing in anisotropic viscoelastic media”. In: Geophysical Jour-nal International 175.2 (2008), pp. 617–626.[106] J. Vidale, D. V. Helmberger, and R. W. Clayton. “Finite-difference seismograms for SHwaves”. In: Bulletin of the Seismological Society of America 75.6 (1985), pp. 1765–1782.[107] J. E. Vidale. “Finite-difference calculation of travel times”. In: Bulletin of the Seismo-logical Society of America 78.6 (1988), pp. 2062–2076.[108] J. E. Vidale. “Finite-difference calculation of traveltimes in three dimensions”. In: Geo-physics 55.5 (May 1990), pp. 521–526.[109] J. Virieux and V. Farra. “Ray tracing in 3-D complex isotropic media: an analysis of theproblem”. In: Geophysics 56.12 (1991), pp. 2057–2069.[110] J. Virieux and S. Operto. “An overview of full-waveform inversion in exploration geo-physics”. In: Geophysics 74.6 (2009), WCC1–WCC26.[111] G. Wahba. “Constrained Regularization for Ill Posed Linear Operator Equations, withApplications in Meteorology and Medicine”. In: Statistical Decision Theory and RelatedTopics, III. Ed. by S. S. Gupta and J. O. Berger. New York: Academic Press, 1982,pp. 383–418.141BIBLIOGRAPHY[112] G. Wahba. Spline Models for Observational Data. Vol. 59. Society for Industrial andApplied Mathematics, 1990.[113] Wikipedia. Accuracy and precision — Wikipedia, The Free Encyclopedia. [Online; ac-cessed 18-December-2014]. 2014. url: http://en.wikipedia.org/w/index.php?title=Accuracy_and_precision&oldid=638426604.[114] Wikipedia. Archie’s law — Wikipedia, The Free Encyclopedia. [Online; accessed 18-December-2014]. 2014. url: http://en.wikipedia.org/w/index.php?title=Archie\%27s_law&oldid=635633886.[115] Wikipedia. Optical resolution — Wikipedia, The Free Encyclopedia. [Online; accessed18-December-2014]. 2014. url: http://en.wikipedia.org/w/index.php?title=Optical_resolution&oldid=635962905.[116] Wikipedia. Reliability theory — Wikipedia, The Free Encyclopedia. [Online; accessed18-December-2014]. 2014. url: http://en.wikipedia.org/w/index.php?title=Reliability_theory&oldid=592914358.[117] Wikipedia. Simplicity — Wikipedia, The Free Encyclopedia. [Online; accessed 18-December-2014]. 2014. url: http://en.wikipedia.org/w/index.php?title=Simplicity&oldid=634382795.[118] P. R. Williamson and R. G. Pratt. “A critical review of acoustic wave modeling proce-dures in 2.5 dimensions”. In: Geophysics 60.2 (1995), pp. 591–595.[119] X. Yang, J. van der Kruk, J. Bikowski, P. Kumbhar, and G. A. Meles. “Frequency-domain Full-waveform Inversion of GPR Data”. English. In: Near-surface Geophysicsand Environment Protection. 2012, pp. 344–348. url: http://www.research.ed.ac.uk/portal/en/publications/frequencydomain-fullwaveform-inversion-of-gpr-data(9d9fcfbd-0ae3-423d-bf0b-6e89f66cd152).html (visited on 09/11/2014).[120] M. J. Yedlin. “Uniform asymptotic solution for the Green’s function for the two-dimensionalacoustic equation”. In: Journal of the Acoustical Society of America 81 (1987), pp. 238–243.142[121] M. J. Yedlin, D. Van Vorst, and J. Virieux. “Uniform asymptotic conversion of Helmholtzdata from 3D to 2D”. In: Journal of Applied Geophysics 78 (2012), pp. 2–8.[122] K. S. Yee. “Numerical solution of initial boundary value problems involving Maxwell’sequations”. In: IEEE Transactions on Antennas and Propagation 14.3 (1966), pp. 302–307.[123] C. A. Zelt and R. B. Smith. “Seismic traveltime inversion for 2-D crustal velocity struc-ture”. In: Geophysical Journal International 108.1 (1992), pp. 16–34.[124] H. Zhao. “A fast sweeping method for eikonal equations”. In: Mathematics of Compu-tation 74.250 (2005), pp. 603–627.[125] C. Zhou, W. Cai, Y. Luo, G. Schuster, and S. Hassanzadeh. “Acoustic wave-equationtraveltime and waveform inversion of crosshole seismic data”. In: Geophysics 60.3 (1995),pp. 765–773.[126] T. Zhu and K.-Y. Chun. “Complex rays in elastic and anelastic media”. In: GeophysicalJournal International 119.1 (1994), pp. 269–276.143Appendix AFull-Waveform Jacobian DerivationThe Jacobian is the matrix of all combinations of first partial derivatives of predicted datawith respect to model parameters, including the unknown source wavelet. The parts of theJacobian associated directly with permittivity and conductivity are arrived at by followingPratt’s approach [82, 81].We start by finding the part of the Jacobian that relates changes in permittivity to changesin predicted data, which is incorporated into the full Jacobian of equation (5.28). The corre-sponding part of the Jacobian for conductivity is derived in a similar manner and the resultsare shown in section A.1. Taking the derivative of equation (5.15) (using equation (5.14)) wehaveJǫ˜(f,s) =∂d(f,s)pred∂ǫ˜=Q(f,s)iω[diag(Dxnh(f,s)y)MxnMc diag(−ǫ0eǫ˜ ◦ ~ǫ−2)+diag(MxnMc~ǫ−1)Dxn∂h(f,s)y∂ǫ˜],(A.1)where we have explicitly indicated that this Jacobian is for a particular frequency and sourceposition with the superscripts f and s. The ◦ operator represents the Hadamard (element byelement) product. The complex permittivity, ~ǫ, and the log scaled relative permittivity, ǫ˜, arevectors containing the discretized material properties over the domain. The magnetic field, hy,and the measurement operator, Q, are those associated with that source and frequency. Inequation (A.1), eǫ˜ =[eǫ˜1 eǫ˜2 · · · eǫ˜n]T and ~ǫ−2 =[ǫ−21 ǫ−22 · · · ǫ−2n]T .Equation (A.1) indicates two contributing components of the Jacobian. The first termof the sum in equation (A.1) is related to perturbations in the permittivity that multipliesE in equation 5.1, and this term is localized to that point in space. The second term is144Appendix A. Full-Waveform Jacobian Derivationrelated to perturbations in the field itself due to non-local changes in permittivity that affectwave propagation. To find ∂hy/∂ǫ˜, we take the derivative of equation (5.11) with respect topermittivity,∂∂ǫ˜(A(f)h(f,s)y,fixed)+A(f)∂h(f,s)y∂ǫ˜ =∂q(f,s)∂ǫ˜ , (A.2)where the hy,fixed term is hy held constant. Written this way, the left-most derivative is appliedto a vector rather than to a matrix. Isolating ∂hy/∂ǫ˜ results in∂h(f,s)y∂ǫ˜ = A(f)−1F(f,s)ǫ˜ , (A.3)withF(f,s)ǫ˜ =[∂q(f,s)∂ǫ˜ −∂∂ǫ˜(Ah(f,s)y,fixed)]. (A.4)The first term in equation (A.4) is∂q(f,s)∂ǫ˜ =−1w2µDxp diag(j(f,s)s)MxnMc diag(ǫ0eǫ˜ ◦ ~ǫ−2), (A.5)where we have neglected changes in the source current, js, due to changes in permittivity. Thesecond term in equation (A.4) is found through differentiation of equation (5.12) to be∂Ah(f,s)y,fixed∂ǫ˜ = −[Dxp diag(Dxnh(f,s)y)MxnMc+Dzp diag(Dznh(f,s)y)MznMc]· diag(k(f,s)−2 ◦ ~ǫ−1ǫ0 ◦ eǫ˜).(A.6)Following the same approach but taking the derivatives with respect to σ˜ allows the calcu-lation of Jσ˜ (see section A.1).Evaluating the second term of equation (A.1) as it is written involves applying A−1 inequation (A.3) to Fǫ˜, which is a matrix with the number of columns equal to the numberof unknowns on the coarse model grid. An alternative is to find the transpose of the secondterm of equation (A.1) (and of the Jacobian) so that A−T in equation (A.3) is applied toDxnT diag(MxnMc~ǫ−1)QT , which is a matrix with the number of columns equal to the number145A.1. Jacobian for Conductivityof data points for the given transmitter and frequency (the size of QT ). Even if the totalnumber of data points is greater than the number of unknowns, the number of data points pertransmitter per frequency is likely less, thus making Jǫ˜T faster to compute than Jǫ˜. Additionally,several source locations may share the same set of receiver positions (i.e. they have the sameQ) allowing the results of the application of A−T to be reused for other source locations. In thepresent work, this is the approach taken. The model matrix A is symmetric, so A−1 = A−T ,thus allowing the algorithm that is used to apply A−1 when computing the forward model tobe used without modification to compute the Jacobian.Working in the frequency domain allows us to factor the forward model operator and thenuse this factorization many times to compute the Jacobian efficiently. This is critically im-portant because the forward model must be computed once for every transmitter position andevery receiver position (and at each frequency) in order to build the Jacobian. In this work, anLDLT factorization of A is found and applied to multiple right-hand sides using MUMPS [3,4].For a given frequency, the Jacobians for each source position are combined into a singlematrix by vertically stacking them,J(f)ǫ˜ =[JT (f,1)ǫ˜ JT (f,2)ǫ˜ · · · JT (f,Ns)ǫ˜]T, (A.7)where the numbered superscripts represent the source positions and f is the frequency index.The Jacobian expressed in equation (A.7) is for a single frequency, and several of these areagain stacked (see equation (5.28)) to obtain the complete Jacobian for all frequencies and allsource positions.A.1 Jacobian for ConductivityThe Jacobian for conductivity can be found by following the same approach as that used inSection 5.4 for permittivity. We start by taking the derivative of equation 5.15 with respect toconductivity:146A.1. Jacobian for ConductivityJ(f,s)σ˜ =∂d(f,s)pred∂σ˜=Q(f,s)iω[diag(Dxnh(f,s)y )MxnMc diag( iωσ0eσ˜ ◦ ~ǫ−2)+diag(MxnMcǫ−1)Dxn∂h(f,s)y∂σ˜].(A.8)To find ∂hy/∂σ˜ we take the derivative of equation (5.11),∂∂σ˜(A(f,s)h(f,s)y,fixed)+A(f,s)∂h(f,s)y∂σ˜ =∂q(f,s)∂σ˜ , (A.9)which can be re-arranged as∂h(f,s)y∂σ˜ = A(f,s)−1F(f,s)σ˜ , (A.10)andF(f,s)σ˜ =[∂q(f,s)∂σ˜ −∂∂σ˜(A(f,s)h(f,s)y,fixed)]. (A.11)The first term of equation A.11 is∂q(f,s)∂σ˜ =1w2µDxp diag(j(f,s)s)MxnMc diag( iωσ0eσ˜ ◦ ~ǫ−2), (A.12)where we have again neglected changes in the source current. The second term of equa-tion (A.11) is again found through differentiation of equation (5.12) to be∂A(f,s)h(f,s)y,fixed∂σ˜ =[Dxp diag(Dxnh(f,s)y)MxnMc+Dzp diag(Dznh(f,s)y)MznMc]· diag( iωk(f,s)−2 ◦ ~ǫ−1 ◦ σ0eσ˜).(A.13)147