UNIVERSIDAD  DE  JAÉN   ESCUELA  POLITÉCNICA  SUPERIOR   DE  LINARES     DEPARTAMENTO  DE  INGENIERÍA   DE  TELECOMUNICACIÓN                   TESIS  DOCTORAL       DEVELOPMENT  OF  SIGNAL   PROCESSING  METHODS  USING   GROUND-­PENETRATING  RADAR  TO   EVALUATE  THE  QUALITY  OF  STONE   MATERIALS       PRESENTADA  POR:   VIOLETA  MONTIEL  ZAFRA     DIRIGIDA  POR:   DR.  D.  NICOLÁS  RUIZ  REYES   DR.  D.  FRANCISCO  JESÚS  CAÑADAS  QUESADA     JAÉN,  16  DE  JUNIO  DE  2017     ISBN  978-­84-­9159-­120-­7   This work was carried out under the supervision of Dr. Nicolás Ruiz Reyes Telecommunication Engineering Department Higher Polytechnic School of Linares University of Jaén Dr. Francisco Jesús Cañadas Quesada Telecommunication Engineering Department Higher Polytechnic School of Linares University of Jaén A mi familia, en especial a mi hermana Azahara por su apoyo incondicional i Acknowledgements This Thesis has been supported by Andalusian Business, Science and Inno- vation Council under Project P11-TIC-7278, being Nicolás Ruiz Reyes the main researcher of the project. This dissertation would not have been possible without the guidance and help of several individuals who contributed and extended their valuable as- sistance. First of all, I would like to thank to Dr. Nicolás Ruiz Reyes for giving me the opportunity to take part in this project, his generosity and support and for sharing his knowledge and experience. I sincerely thank Dr. Francisco Jesús Cañadas Quesada for his technical guidance on signal processing and encouragement when required. I also thank him for his con- structive suggestions when reviewing the articles and this thesis. I wish to express my special gratitude to Dr. Pedro Vera Candeas for the patience, time and expertise contributed to the enhancement of this dissertation. Dur- ing these years, we held stimulating debates, which were a great source of inspiration to me. I am also grateful to the members of the team of the project, specially to Dr. Javier Rey Arrans and Dr. Julián Mart́ınez López for his encouragement and interest during my research. I would like to show my gratitude to the Telecommunication Engineering Department, for facilitating the access to its resources and infrastructures and for the support. Special thanks to Pedro Aguilar Aguilar for his kindness and his technical support when required. Also I thank Roćıo Pérez de Prado and Julio Carabias Orti for sharing with me my first steps in teaching. I would like to express my sincere thanks to my colleagues of the office Pablo, Francisco, Diego, Casto, Unai and Javi for making the everyday life enjoyable. I deeply thank Antonio for his continuous support and help during the final steps of this Thesis. iii Finally, I would like to express my gratitude to Dr. Luigi Zanzi from the Politécnico di Milano for offer me an opportunity to work with his research group. Thank you. Violeta Montiel Zafra June 2017 iv Agradecimientos Esta Tesis ha sido posible gracias a la beca asociada al Proyecto de Exce- lencia de la Junta de Andalućıa P11-TIC-7278, siendo el Dr. Nicolás Ruiz Reyes el investigador principal del proyecto. Este trabajo no habŕıa sido posible sin la gúıa y ayuda de varias personas que han contribuido a confeccionarlo y mejorarlo. Primero, quisiera agrade- cer al Dr. Nicolás Ruiz Reyes por darme la oportunidad de formar parte de este proyecto, por su generosidad y apoyo y por compartir su conocimiento y experiencia conmigo. Agradezco sinceramente al Dr. Francisco Jesús Cañadas Quesada por su conocimiento técnico en procesado de seńal y su apoyo siempre que lo he necesitado. Le agradezco también sus sugeren- cias constructivas revisando los art́ıculos y esta Tesis. Quisiera expresar mi agradecimiento más sincero al Dr. Pedro Vera Candeas por su paciencia, tiempo y experiencia que han contribuido más que significativamente a la mejora de este trabajo. Durante estos ańos hemos mantenido debates muy estimulantes que han sido una gran fuente de inspiración para mı́. Gracias también a los miembros del equipo del proyecto, especialmente al Dr. Javier Rey Arrans y al Dr. Julián Mart́ınez López por su especial apoyo e interés durante mi investigación. Quisiera mostrar mi agradecimiento al Departamento de Ingenieŕıa de Telecomunicación por facilitarme el acceso a sus recursos e infraestructuras y por el apoyo constante. Gracias a Pedro Aguilar Aguilar por su amabilidad y apoyo técnico cuando lo he necesitado. También quisiera agradecer a la Dra. Roćıo Pérez de Prado y al Dr. Julio Carabias Orti por compartir conmigo mis primeros pasos en la docencia. Quisiera expresar mis gracias más sinceras a mis compańeros de la sala Pablo, Francisco, Diego, Casto, Unai y Javi por hacer cada d́ıa de trabajo un d́ıa entretenido y divertido. Agradecer profundamente a Antonio por su apoyo continuo y su ayuda durante los últimos pasos de la consecución de v esta Tesis. Finalmente, quisiera expresar mi agradecimiento al Dr. Luigi Zanzi del Politécnico de Milán por ofrecerme la oportunidad de trabajar en su grupo de investigación. Gracias. Violeta Montiel Zafra Junio 2017 vi Abstract In the last decades, quality control in natural stone materials has become a critical task to evaluate correctly the density of the discontinuities, to facili- tate the exploitation, to minimize extraction and transport costs, to improve the cutting process, etc. Most of the anisotropies in stone materials come from their non-homogeneity, which can be classified as cavities, fractures, microfractures, karstification, etc. Non-destructive evaluation (NDE) is a novel procedure used in Science and Technology industry to test material properties without causing any damage. In particular, Ground-Penetrating Radar (GPR) is a high reso- lution geophysical method which is originally designed to investigate the subsurface of the earth. Under the convenient conditions, this technique can provide accurate information of the nature of buried objects. More- over, this technology has been employed to analyse the internal structure of ornamental rocks with promising results. In this thesis, a set of novel signal processing methods applied to GPR are proposed. These methods are mainly focused on providing a more use- ful GPR information for the user in order to improve the evaluation and characterization of the internal state of a stone material. Specifically, the global system developed has three main contributions: characterisation of stone materials, detection and classification of defects and removal of the typical noise active in GPR images. First of all, it is necessary to evaluate and characterize stone materials using GPR. Thus, a technique which highlights the defects inside the stone materials, both quarries and extracted blocks, is designed. For this purpose, a deconvolution method which has not yet been used in GPR images to our best knowledge is employed. Then, the defects are clearly emphasized and the probability of detection is increased. Besides, the stratification and mica vii schist contents of massifs can be located. The proposed method is evaluated using different GPR antennas applied to real data. It is proved that GPR is a suitable tool for the diagnosis of stone materials that has been widely used in the field of non-destructive evaluation and testing applied to the extraction or cutting process in the stone industry. Secondly, a novel GPR signal processing method is developed which au- tomatically detects defects. Besides, it establishes a classification of such defects according to its spatial orientation (horizontal, vertical or diagonal). Moreover, 3D maps have been designed for showing more clearly the defects. This algorithm employs a deconvolution technique trained with a synthetic database, which search for possible spatial orientation defects. This method has been evaluated using synthetic and real databases, with promising re- sults. Then, it can be used for determining the internal quality of a stone block. Besides, the classification of fractures can be useful for the cutting process of the extracted blocks, since during this process vertical and diago- nal fractures are critical for the breaking of the slabs and the cutting speed can be adjusted accordingly. Thirdly, we propose a novel, efficient and fast technique to remove back- ground noise present in GPR images. This noise reduces the resolution of these images and masks possible anisotropies inside the evaluated mate- rial. This method is based on exploiting the repetitive pattern shown by the horizontal noise in the direction of the movement of the antenna. It has been trained using synthetic GPR profiles, and it has been evaluated with synthetic and real data. The novel technique outperforms other classic background noise removal methods. Summarizing, the proposed methods can be considered useful and inter- esting to be applied in the stone industry because they can improve different processes such as, evaluation of the internal quality of stone blocks, reduc- tion of cutting and transport costs or minimization of the amount of waste products. Keywords: Ground-Penetrating Radar, GPR, Non-destructive Evaluation, NDE, Non-destructive Testing, NDT, Processing, Deconvolution, Migration, SI-PLCA, Background, Horizontal noise, non-neighboring. viii Resumen En las últimas décadas, la calidad de control en materiales pétreos se ha convertido en una tarea cŕıtica, de forma que se hace necesario evaluar cor- rectamente la densidad de las discontinuidades, para facilitar su explotación, para minimizar costes de extracción y de loǵıstica, para mejorar el proceso de corte, etc. La mayoŕıa de anisotroṕıas en piedra parten de su inhomogenei- dad, que pueden ser clasificadas como cavidades, fracturas, microfracturas, karstificación, etc. La Evaluación No Destructiva es un procedimiento novedoso en la indus- tria de la Ciencia y Tecnoloǵıa para testear las propiedades de un material sin causar daño en él. En particular, el Georradar es una técnica geof́ısica de alta resolución, originalmente diseñada para investigar la profundidad de la tierra. Bajo unas condiciones determinadas, esta técnica puede ofrecer información de la naturaleza de objetos enterrados. Además, esta tecnoloǵıa ha sido usada para analizar la estructura interna de piedras ornamentales con resultados prometedores. En esta Tesis se propone un conjunto de métodos novedosos de procesado de señal aplicados al georradar. Estos métodos están enfocados en ofrecer una información al usuario más útil, para mejorar la evaluación y caracteri- zación del estado interno de un material pétreo. Concretamente, el sistema total desarrollado presenta tres principales contribuciones: caracterización de materiales pétreos, detección y clasificación de defectos y eliminado del ruido normalmente activo en imágenes del georradar. Primero, es necesario evaluar y caracterizar materiales pétreos usando georradar. Entonces, se ha diseñado una técnica que resalta los defectos dentro de materiales pétreos, tanto canteras como bloques extráıdos. Para este propósito se usa un método de deconvolución, que no ha sido aún usado en imágenes de georradar para nuestro conocimiento. Entonces, los defectos ix se enfatizan claramente y la probabilidad de detección aumenta. Además, la estratificación y los micaesquistos pueden ser claramente localizados. El método propuesto se evalúa usando diferentes antenas aplicadas a datos reales. Se ha probado que el georradar es una herramienta realmente útil para la diagnosis de materiales pétreos que ha sido ampliamente usado en el campo de la evaluación no destructiva aplicada a la extracción o al proceso de corte en la industria de la piedra. Posteriormente, un método novedoso de procesado de señal ha sido diseñado, que automáticamente detecta defectos. Además, establece una clasificación de los defectos de acuerdo a su orientación espacial (horizon- tal, vertical o diagonal). También se han diseñado unos mapas 3D para mostrar más claramente los defectos. Este algoritmo usa una técnica de de- convolución que ha sido entrenada con una base de datos sintética, que busca defectos con las posibles orientaciones espaciales. Este método ha sido eval- uado usando bases de datos reales y sintéticas, con resultados prometedores. Entonces, este método se puede usar para determinar la calidad interna dde un bloque de piedra. Además, la clasificación de fracturas puede ser útil para el proceso de corte de bloques extráıdos, ya que durante dicho proceso las fracturas verticales y diagonales son cŕıticas en la rotura de las tablas extráıdas, aśı que la velocidad de corte puede ser ajustada convenientemente. Por otra parte, proponemos una técnica novedosa, eficiente y rápida para eliminar el ruido de fondo presente en las imágenes de georradar. Este ruido reduce la resolución en estas imágenes y enmascara las posi- bles anisotroṕıas dentro del material a evaluar. Este método está basado en la búsqueda de patrones repetitivos t́ıpicos en el ruido horizontal en la dirección del movimiento de la antena. Ha sido además entrenado usando perfiles sintéticos de georradar, y ha sido evaluado con datos sintéticos y reales. La técnica supera los resultados de otras técnicas clásicas de elimi- nado de ruido de fondo. Para concluir, los métodos propuestos pueden considerarse útiles e in- teresantes para poder ser aplicados a la industria de la piedra, ya que pueden mejorar diferentes procesos como la evaluación de la calidad interna de blo- ques de piedra, reducción de costes de corte y transporte y minimizado la cantidad de residuos. Keywords: Georradar, Evaluación No Destructiva, Procesado, Deconvolución, Migración, SI-PLCA, Fondo, Ruido horizontal. x List of Acronyms ACC Accuracy AGC Automatic Gain Control ATS Average Trace Subtraction BMS Background Matrix Subtraction CMP Common Mid-Point DTVM Directional Total Variation Minimisation EM Electromagnetic EUC Euclidean Distance FA False Alarm FDTD Finite-Difference Time-Domain FFT Fast Fourier Transform GPR Ground-Penetrating Radar HR Hit Rate HT Hilbert Transform ICA Independent Component Analysis NDE Non-destructive Evaluation NDT Non-destructive Technique PCA Principal Component Analysis PLCA Probabilistic Latent Component Analysis RMS Root Mean Square RRME Root Mean Square Error SI-PLCA Shift Invariant Probabilistic Latent Component Analysis SIMCA SIMulated Correlation Algorithm SVD Singular Value Decomposition SSI Subsurface Imaging TV Total Variation xi List of Symbols Symbol Units Definition Ē [V/m] Electric field strength vector H̄ [A/m] Magnetic field intensity vector D̄ [C/m2] Electric displacement vector B̄ [T ] Magnetic flux density vector J̄ [A/m2] Electric current density vector q [C/m3] Electric charge density ε [F/m] Permittivity µ [H/m] Permeability εr - Relative permittivity µr - Relative permeability σ [S/m] Conductivity v [m/s] Velocity of propagation ω [rad/s] Angular frequency β [rad/m] Phase coefficient α [Neper/m] Attenuation coefficient η [Ω] Intrinsic impedance r - Reflection coefficient φ [◦] Angle xiii Symbol Units Definition d [m] Distance to the target t [s] Two-way travel time to the target λ [m] Wavelength fc [Hz] Frequency R - Bandwidth to center frequency ratio B [Hz] Bandwidth ∆v [m] Vertical resolution ∆h [m] Horizontal resolution ∆x [m] Sampling interval vrms [m/s] RMS velocity k [m−1] Temporal wavenumber kx [m−1] Spatial x wavenumber kz [m−1] Spatial z wavenumber xiv List of Constants Symbol Value Definition ε0 8.85x10−12F/m Electric permittivity of free space µ0 1.26x10−6H/m Magnetic permeability of free space c 3x108m/s Speed of light xv List of Figures 2.1 Incident, reflected and refracted waves from a layer with con- trasting permittivities [4] . . . . . . . . . . . . . . . . . . . . 14 2.2 GPR basic scheme [6]. . . . . . . . . . . . . . . . . . . . . . . 15 2.3 GPR typical reflection schematic [22] . . . . . . . . . . . . . . 16 2.4 Scans description [22] . . . . . . . . . . . . . . . . . . . . . . 18 2.5 Resolution scheme [88]. . . . . . . . . . . . . . . . . . . . . . 21 2.6 Example of a radargram with clutter. . . . . . . . . . . . . . 22 2.7 Example of the clutter frequency concept. At 100 MHz the clutter is clearly visible, whereas at 50 MHz much of the clut- ter is suppressed [6]. . . . . . . . . . . . . . . . . . . . . . . . 23 2.8 The structure of a quarry and a 250 MHz radargram, where the layers are clearly highlighted [58]. . . . . . . . . . . . . . 26 2.9 A schematic presentation of two GPR radargram which show the propagation of the fractures in a quarry [60]. . . . . . . . 26 2.10 Cross sections of a stone block, where the horizontal and dip- ping fractures are marked by arrows [58]. . . . . . . . . . . . 28 2.11 3D map of a thin discontinuity present on a marble block [8]. 28 2.12 Example of the method proposed by Longoni et. al: raw radargram (top left), processed radargram with the disconti- nuities (top right), interpolated discontinuities (bottom left) and fitting planes (bottom right) [56]. . . . . . . . . . . . . . 30 2.13 An example of a radargram after applying each of the pre- processing methods above mentioned. . . . . . . . . . . . . . 32 2.14 An example of a radargram after applying background re- moval techniques [115]. . . . . . . . . . . . . . . . . . . . . . . 34 xvii 2.15 The harbor example shows that the gap acts as a Huygens’ secondary source (a) and the waves generated by Huygens’ secondary source have a hyperbolic traveltime trajectory (b) [126]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.16 Kirchhoff Migration. (a) Raw radargram, where a target ap- pears as the characteristic hyperbola. (b) Radargram after Kirchhoff [126]. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.17 Test examples of aperture widths in Kirchhoff migration [126]. 39 2.18 Test examples of velocity errors in Kirchhoff migration [126]. 40 2.19 F-K Migration. (a) A dipping reflector is represented by 0B in the (f, k) plane. (b) After migration, the radial line 0B has been mapped to 0B′ [126]. . . . . . . . . . . . . . . . . . 41 2.20 The diffraction hyperbola in the (x, t) plane is mapped onto a triangular area in the (f, k) plane [126]. . . . . . . . . . . . 42 2.21 Flow diagram of the proposed method by Al-Nuaimy [3]. . . . 43 2.22 Flow diagram of the proposed method by Sezgin [102]. . . . . 44 2.23 Flow diagram of the proposed method by Pasolli et. al [102]. 45 2.24 3D data volume which highlights the fractures at different depths ((a), (b) and (c)) [47]. . . . . . . . . . . . . . . . . . . 46 3.1 The four models that constitute the training database of pub- lication [P2]. In this case, white colour represents a permit- tivity value of 8 (the stone part of the block) whereas the black colour is related to the air, so the permittivity is 1. . . 53 3.2 The four radargrams related to the four models of the database of publication [P2] shown in Fig. 3.1. As can be observed, each resulting radargram comprises only non-negative values. 54 3.3 The set of models of a block of the database related to pub- lication [P3]. The selected relative permitivitties vary from 1 (black) to 8 (white). . . . . . . . . . . . . . . . . . . . . . . . 59 3.4 Types of noises of the database related to publication [P3]. . 62 xviii List of Tables 2.1 Typical range of EM parameters of various materials mea- sured at 100 MHz [23]. . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Relationship between EM parameters, operating parameters and the selected frequency [125]. . . . . . . . . . . . . . . . . 20 3.1 Setup of the real data acquisition related to publication [P1] . 48 3.2 Setup of the simulation related to publication [P2] . . . . . . 50 3.3 Setup of the real data acquisition related to publication [P2] . 51 3.4 Detailed information about length, width and orientation of each file of the database related to publication [P2]. . . . . . 55 3.5 Setup of the simulation related to publication [P3] . . . . . . 57 3.6 Setup of the real data acquisition related to publication [P3] . 58 3.7 Characteristics of the measured data of [P3]. . . . . . . . . . 60 xix Contents Abstract vii Resumen ix 1 Introduction 1 1.1 Context and motivation . . . . . . . . . . . . . . . . . . . . . 1 1.2 Scope of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Scientific Contributions . . . . . . . . . . . . . . . . . . . . . 5 1.4 Thesis structure . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 State of the art 9 2.1 Physical and mathematical background related to EM waves and GPR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.1 Propagation of electromagnetic waves in dielectric ma- terials . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.2 Ground-Penetrating Radar . . . . . . . . . . . . . . . 15 2.2 The use of GPR in the stone industry . . . . . . . . . . . . . 24 2.2.1 GPR applied to evaluate the quality related to quarries 24 2.2.2 GPR applied to evaluate the quality related to stone blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2.3 GPR applied to general stone assessment . . . . . . . 29 2.3 Signal Processing applied to GPR . . . . . . . . . . . . . . . 30 2.3.1 Pre-processing . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Background Removal . . . . . . . . . . . . . . . . . . . 33 2.3.3 Migration . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.4 Pattern recognition . . . . . . . . . . . . . . . . . . . . 42 2.4 3D modelling in GPR . . . . . . . . . . . . . . . . . . . . . . 45 xxi 3 Methodology 47 3.1 Experimental work description . . . . . . . . . . . . . . . . . 47 3.1.1 Simulation software . . . . . . . . . . . . . . . . . . . 47 3.1.2 Equipment for real acquisition data . . . . . . . . . . 48 3.2 Characterization of ornamental stone . . . . . . . . . . . . . . 48 3.2.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2.2 Material . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.3 Database . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Detection and Classification of fractures . . . . . . . . . . . . 49 3.3.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.2 Material . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3.3 Databases . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.3.4 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 GPR Background removal . . . . . . . . . . . . . . . . . . . . 57 3.4.1 Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.2 Material . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.3 Databases . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.4.4 Metrics . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4 Results, Conclusions and Future Work 63 4.1 Characterization of ornamental stone . . . . . . . . . . . . . . 63 4.2 Detection and Classification of fractures . . . . . . . . . . . . 64 4.3 GPR Background removal . . . . . . . . . . . . . . . . . . . . 66 Bibliography 69 Paper A: Ground-penetrating radar method used for the char- acterisation of ornamental stone quarries 83 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 2. Geological setting . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3. GPR data acquisition and processing . . . . . . . . . . . . . . . 89 4. Image enhancement . . . . . . . . . . . . . . . . . . . . . . . . . 90 5. Results and interpretation . . . . . . . . . . . . . . . . . . . . . 92 6. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 xxii Paper B: Detection and classification of internal defects in lime- stone blocks based on a deconvolution technique with SI- PLCA applied to GPR signals 105 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 2. Rock specimens . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3. Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4. Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5. Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6. Conclusions and Future Work . . . . . . . . . . . . . . . . . . . 134 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 Paper C: A novel method to remove GPR background noise based on the similarity of non-neighboring regions 143 1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 2. Proposed method . . . . . . . . . . . . . . . . . . . . . . . . . . 147 3. Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 4. Conclusions and Future Work . . . . . . . . . . . . . . . . . . . 175 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 xxiii List of Included Publications This thesis is a compound thesis consisting of the following publications, which are preceded by an introductory overview of the research field and a summary of the contributions. The publications will be referred in the text with [P1], [P2] and [P3]. • [P1] Rey J., Mart́ınez J., Vera P., Ruiz N., Cañadas F., Mon- tiel V., ”Ground-penetrating radar method used for the characteri- sation of ornamental stone quarries”, in Construction and Building Materials, vol. 77, pp. 439-447, 2015. • [P2] Montiel V., Cañadas F., Vera P., Ruiz N., Rey J., Mart́ınez J., ”Detection and classification of internal defects in limestone blocks based on a deconvolution technique with SI-PLCA applied to GPR signals”, Research in Nondestructive Evaluation. Status: under re- view. • [P3] Montiel V., Cañadas F., Vera P., Ruiz N., Rey J., Mart́ınez J., ”A novel method to remove GPR background noise based on the similarity of non-neighboring regions”, Journal of Applied Geophysics. Status: under review. xxv Chapter 1 Introduction 1.1 Context and motivation Subsurface imaging (SSI) is the imaging of an object buried below the surface of a medium [94]. Non-destructive Testing (NDT) or Non-Destructive Eval- uation (NDE) is the process of inspecting, testing, or evaluating materials in order to search for discontinuities, or differences in certain characteristics without destroying the functionality of the component under study. In this manner, the material has not suffered any damage when the test is com- pleted. Destructive testing implies damages suffered by the material. As a result, several samples are often necessary to inspect in order to determine its properties. Like other type of products, quality control in natural stone is a critical task, since evaluation of massive stone quality and homogeneity could lead to primary waste reduction and improving overall sustainability. [37]. This material can present structural discontinuities, such as fractures, cavities, etc., or studying its layers could be interesting. Determining the structure of a quarry could reduce extraction costs, since it could facilitate the plan- ning of the advance of its exploitation, minimizing the use of explosives. Besides, the size or position of the extracted standard-sized blocks could be adequately selected. In addition, the early detection of internal defects of stone blocks can reduce cutting and transport costs, amount of waste products and can optimize the production of slabs or smaller blocks. Specif- ically, the probability of occurrence of a defect inside an extracted block is really high in limestones. Based on information, a cutting protocol of the 1 extracted blocks can be established, or to consolidate the block by means of resins (or other components) or maybe to establish injection points, thus improving the performance of the entire process. Nowadays, NDT, as an upcoming technology, is used in several processes in order to assess the quality of product or to estimate the properties of stone. Several authors have proposed numerous NDT methods in order to evaluate stone products. Bodare [15] separated the NDT methods into two major categories, depending on the type of wave: seismic and electric methods. On the one hand, seismic methods can be divided, in turn, as ultrasonic and impact-echo methods. First, ultrasonic methods have been widely used to analyse the internal structure of stone products [17] [121] [71] [21] [93] [96] [14]. While it is true that the resulting resolution is really high and both big and small defects with excellent precision can be detected, ultrasonic tech- niques use sophisticated procedures (each transducer should be positioned along the material), and real measurements require rather time consuming [72]. Then, impact-echo methods have been widely used mainly to analyze concrete structures, but it has been also used to evaluate stone [92]. The main drawback of this method is that it can only provide information about the general quality of the material, without detailing the number, size, shape or spatial orientation of the possible defects. On the other hand, electric methods are divided into radar, resistive and electromagnetic methods. Radar techniques are undoubtedly the most used for these tests. In particular, Ground-Penetrating Radar (GPR), as a suitable diagnosis tool, is widely employed in the investigation of stone properties and fracturing [104] [40] [105] [57] [58] [82] [26] [34] [47] [8] [43] [116], due to the fact that the dielectric properties of stone materials allows this technique to perform well. GPR has become one of the most common used NDT methods for the assessment of natural stone. Under favorable conditions, this technique allows not only accurately identify targets, but characterize the environment, discriminate materials, detect inhomogeneities, etc., without any damage of the medium to test. In addition to its resolution (both lateral and vertical), GPR becomes a promising prospecting method due to its quick and easy operation. GPR antennas do not need to be in contact with the surface of the material to test, which enable a simple measure, and with an affordable cost. Finally, the reflection amplitudes and the travel times are compacted into 2 a high-resolution image, known as radargram, where the different objects are shown. In geological applications, this technique is used to analyse the structure of the stone at depths of 30-60 m. In other high resolution applications, high frequency antennas are able to offer a better noise control, as well as to provide a high quality image, with a depth of 5-10 m. However, the challenge of GPR processing is to get an accurate radar- gram analysis method, since this technique presents some limitations which need to be solved: the obtained radargrams are not a realistic picture of the medium and its interpretation is not intuitive. Then, the final user must interpret the image based on his own experience, so the analysis is subjec- tive. The beamwidth of the antenna, the wavelength of the radiation and the dimensions of the target will cause that the final image will not corre- spond to reality. In general, most stone materials present high attenuation, which makes the penetration of the waves more difficult. Moreover, the signal to clutter ratio of the radargram is the key to target detection [23], where clutter is a typical background noise present in almost all of these radargrams. In this context, the noise energy of the clutter masks the real targets existing inside the stone. Currently, there are several background removal methods as will be detailed in Chapter 2, in Section 2.3.2 whose main goal is to remove this type of noise. The main objective of any signal processing method applied to GPR is to provide a reconstructed radargram which can be easily interpreted by the user, removing any possible noise and also to automatically classify the target. Moreover, 3D GPR maps can contribute to a better understanding of the internal structure of a stone product, showing the depth, shape and the spatial localization of the detected fractures. There are different commercial programs (e.g. ReflexW [87]) which correlate the radargrams and offer a 3D vision of a material, but a classification of the defects cannot be done to our best knowledge. Several authors have also investigated in 3D GPR modelling as will be shown later (in Section 2.4 of the Chapter 2). This thesis applies GPR to evaluate the quality of stone materials. In particular, different signal processing methods are proposed to enhance the probability of detection of fractures, to automatically detect and classify the defects and to remove background noise. Then, the three methods compose a general system to remove noise, to detect and to classify fractures in stone materials. 3 For the first purpose, GPR is used to probe the texture of different types of stone materials, as well as to detect anisotropies. A deconvolu- tion technique is used for searching for activations of a GPR typical pulse in the radargram. Then, it enhances the energy related to the encoun- tered defects inside real GPR data, increasing the probability of detection of the anisotropies. The performance is evaluated over quarries and ex- tracted blocks, proving that GPR together with a deconvolution technique is useful for the evaluation of the general quality of stone materials as well as for the determination of their texture and composition. Subsequently, GPR is applied to automatically detect and classify frac- tures according to their spatial orientations. Then, a deconvolution tech- nique is used with a synthetic training database, which is composed of pos- sible defects with vertical, horizontal and diagonal orientations. Later, the proposed method is objectively evaluated using a synthetic testing database. Besides, real blocks data is used for the creation of 3D maps, one per orien- tation, which show their internal structure. The proposed method outper- forms the detection and classification performance compared to state-of-art methods. Besides, the results have shown that GPR can be used as a clas- sification tool for discarding blocks directly at the quarry or to optimize the cutting process. Finally, a background removal method is developed in order to remove the horizontal noise, which often masks the real anisotropies present in the GPR image. Since this noise is propagated in the direction of the move- ment of the GPR antenna, the method identifies the repeating regions to remove them. There is an optimization procedure using a synthetic training database. Then, experimental results over a test database, with synthetic and real data, prove that the proposed method obtains better results com- pared to the classic background removal techniques. The database is also composed with synthetic and real data and the results are compared with classic background methods, obtaining promising results. 1.2 Scope of the Thesis This thesis is focused on developing novel signal processing methods using GPR data in order to evaluate the quality of stone materials by means of 4 characterization, detection and classification of defects detected. In this Thesis, we explore two application scenarios: stone quarries and extracted blocks. The main objectives of the present thesis are: • To characterize stone quarries, by using a deconvolution technique in order to probe the texture and to enhance the anisotropies detection. • To create a method for detecting fractures and a 2D and 3D models for classification fractures according to their spatial orientation in stone blocks, by using a deconvolution technique. • To remove GPR background noise present in most radargrams, which is based on the non-local similarity of 2D regions over the distance in order to identify repetitive regions of the same depth. 1.3 Scientific Contributions Main scientific contributions of this thesis comprise: [P1] Ground-penetrating radar method used for the character- isation of ornamental stone quarries As mentioned before, evaluating the extension and the potential of a stone quarry can produce remarkable benefits in terms of reducing costs, risks and environmental impact as well as increasing the productivity. Then, this paper uses GPR to predict the texture and the anisotropies in stone for ornamental use. The proposed scheme consists of two stages: data acqui- sition and processing and radargram enhancement. First, data is acquired in a real limestone quarry and in commercial extracted blocks using differ- ent GPR antennas and a variety of pre-processing tools is applied. Then, the anisotropies are highlighted using a Shift-Invariant Probabilistic Latent Component Analysis (SI-PLCA). This algorithm separates the anisotropies from noise, by searching for activations of the typical GPR pulse across the two dimensions of the image. GPR is founded as a suitable tool for the diagnosis of stone materials and SI-PLCA is a good approach to highlight the relevant information. Besides, in this work it is observed that the GPR response is better for horizontal than for vertical fractures, so this result suggests that a classification of defects can be desirable. 5 [P2] Detection and classification of internal defects in limestone blocks based on a deconvolution technique with SI-PLCA applied to GPR signals It is well known that analysing the structure of stone blocks before per- forming the cutting process can improve the productivity. In this work, a novel approach to improve the detection and to classify defects accord- ing to their spatial orientation in damaged stone blocks is presented. For this purpose, SI-PLCA algorithm is used. First, the method is based on training fractures with a variety of orientations, such as vertical, horizontal an diagonals. Thus, SI-PLCA provides information about the activation of each pattern at each position of the block. Then, it is possible to know the number of defects in a block, as well as its importance and spatial ori- entation. Besides, a 3D map provides information about the fractures for each orientation pattern. Then, this type of classification is a novel result to our best knowledge in GPR signal processing. For the experiments, a synthetic and a real databases are employed and the results are objectively compared with state-of-art methods, with promising detection and classifi- cation results. Besides, the novel method obtains a suitable detection when the fracture is vertical, which is the most critical spatial orientation during the cutting process, as indicated in [P1] objectives. [P3] A novel method to remove GPR background noise based on the similarity of non-neighboring regions One of the major challenges in the research related to GPR image is the improvement of its quality. Due to interferences of the horizontal back- ground noise, the resolution quality is decreased and the true anisotropies cannot be detected. In order to remove this type of noise, this study pro- poses a novel, fast and robust background removal. It is based on observing the repetitive pattern shown by the noise in the direction of the movement of the antenna. Then, for a certain region, other repetitive ones at the same depth are searched. The method assumes that regions composed of horizontal background noise show higher similarity between them, whereas regions mainly composed of anisotropies present lower similarity. Besides, an optimization procedure for the selection of the parameters of the method needs to be applied. Experimental results show a promising performance compared to classic background removal techniques. 6 1.4 Thesis structure The rest of the thesis is organized as follows. Chapter 2 presents the state of art of the use of GPR for quality control purposes. Chapter 3 includes the methodology, where the stone specimens and the experimental work are described. Databases and metrics for evaluation are included here. Chapter 4 is focused on the evaluation of the proposed methods for characterisation of stone quarries, for detection and classification of defects in stone blocks, and finally for removal of background noise active in GPR data. Besides, future work is included. Finally, the articles published during this thesis are included in the appendix. 7 Chapter 2 State of the art The potential of GPR for analysing the internal structure of stone materials has been addressed. As explained in Chapter 1, GPR is widely used for prospecting stone materials due to its dielectric properties, both quarries and, less frequently, for extracted blocks. In order to obtain a fast evaluation of the material, it is often that classic activities use the raw GPR radargrams without any post-processing over quarries or blocks. In some studies, it is common to simply use pre- processing methods, which are present in most commercial software ([87] [84]). However, in other recent cases, several signal processing algorithms are applied to the GPR radargrams in order to remove noise and to recon- struct the real image. Besides, numerous specific signal processing methods have been created, in order to improve the detection of discontinuities or even to do a classification of the detected targets. Moreover, new algo- rithms have been developed for removing the background noise present in most radargrams. Finally, 3D visualization methods have been developed by several authors, since this type of reconstruction clearly shows the internal spatial structure of a material. In this Chapter, the state of the art related to GPR is presented. First, it is necessary to explain briefly the mathematical and physical background that supports GPR, including the propagation of the waves trough the ma- terial, the well-known Maxwell equations and the electromagnetic (EM) pa- rameters of the medium. The equations related to each operating parameter of the GPR system are presented. Besides, the EM properties of the stone material are detailed. Finally, a complete review of the most referenced and 9 recent published techniques are presented, separating the works related to the use of GPR over stone materials (both quarries and blocks), the gen- eral signal processing which most authors apply to their data and the 3D reconstruction methods. 2.1 Physical and mathematical background related to EM waves and GPR In this section, a brief review of the GPR fundamentals is presented. For this purpose, it is necessary to explain how the emitted waves are propagated trough the material. Moreover, as this Thesis is focused on stone materials, a review of the properties of these materials is depicted. To conclude this Chapter, an overview of the GPR principles is presented. 2.1.1 Propagation of electromagnetic waves in dielectric ma- terials The propagation of the EM waves along a medium was originally described by Maxwell [62], which summarized the work of prior researches in a more compacted form (equations 2.1 and 2.2). ∇̄ × Ē = −∂B̄ ∂t (2.1) ∇̄ × H̄ = J̄ + ∂D̄ ∂t (2.2) ∇̄ • D̄ = q (2.3) ∇̄ • B̄ = 0 (2.4) where Ē is the electric field strength vector, H̄ is the magnetic field intensity vector, D̄ is electric displacement vector, B̄ is the magnetic flux density vector, J̄ is the electric current density vector and finally q is the electric charge density. These four relations represent the properties and the dependence be- tween electric and magnetic fields. In other words, the Maxwell’s equations 10 establish that electric currents generate magnetic fields, and vice versa. Then, these relations introduce the property parameters (or EM parame- ters) of permittivity ε, magnetic permeability µ and conductivity σ. The material’s response to EM fields are described by the following constitutive relationships (equations 2.5 2.6 2.7). J̄ = σ̃Ē (2.5) D̄ = ε̃Ē (2.6) B̄ = µ̃H̄ (2.7) where σ̃, ε̃ and µ̃ are tensor quantities related to conductivity, permit- tivity and magnetic permeability respectively. Permittivity describes the capability of a material to store and release energy from an electric field. The higher the permittivity parameter, the more energy is absorbed by the medium or, in other words, the slower the wave travels through the medium. The relative permittivity is defined by εr = ε/ε0, where ε0 is the electric permittivity of free space. The wet mate- rials will slow down the EM wave, since the water presence causes that the permittivity increases. Magnetic permeability indicates the capability of a material to store and release energy from a magnetic field. It is related to the magnetic polariza- tion of a material. Then, when a material present a high permeability, it will interfere with the magnetic part of the wave. The relative permeability of a material is µr = µ/µ0, where µ0 is the magnetic permeability of free space. For non-ferromagnetic materials it is assumed that this value is equal to 1. Conductivity is the ability of a material to conduct electric current. This parameter affects to the penetration of the waves through the medium. Materials which present high conductivity conduct the electric part of the wave, and so the attenuation increases. Then, a low conductivity would be desirable. When the material is resistive (low conductivity), such as dry sand, ice or dry concrete, the wave stays intact longer and continues travelling through the material. However, if the material is conductive, such as salt water or wet concrete, the EM wave is absorbed. 11 Table 2.1 shows the conductivity and relative permittivity of several materials. As can be observed, wet materials present higher relative per- mittivity than dry mediums (and so the wave reduces its velocity). Besides, the conductivity is lower in materials such as dry concretes, limestones and sands, but it is higher in materials such as sea water. Material Conductivity Relative Permittivity Air 0 1 Asphalt dry 10−2 : 10−1 2–4 Asphalt wet 10−3 : 10−1 6–12 Clay dry 10−1 : 10−0 2–6 Clay wet 10−1 : 10−0 5–40 Coal dry 10−3 : 10−2 3.5 Coal wet 10−3 : 10−1 8 Concrete dry 10−3 : 10−2 4–10 Concrete wet 10−2 : 10−1 10–20 Freshwater 10−6 : 10−2 81 Freshwater ice 10−4 : 10−3 4 Granite dry 10−8 : 10−6 5 Granite wet 10−3 : 10−2 7 Limestone dry 10−8 : 10−6 7 Limestone wet 10−2 : 10−1 8 Permafrost 10−5 : 10−2 4–8 Rock salt dry 10−4 : 10−2 4–7 Sand dry 10−7 : 10−3 2–6 Sand wet 10−3 : 10−2 10–30 Sandstone dry 10−6 : 10−5 2–5 Sandstone wet 10−4 : 10−2 5–10 Sea water 102 81 Sea water ice 10−2 : 10−1 4–8 Shale dry 10−3 : 10−2 4–9 Shale saturated 10−3 : 10−1 9–16 Snow firm 10−6 : 10−5 6–12 Soil clay dry 10−2 : 10−1 4–10 Soil clay wet 10−3 : 10−0 10–30 Soil loamy dry 10−4 : 10−3 4–10 Soil loamy wet 10−2 : 10−1 10–30 Soil sandy dry 10−4 : 10−2 4–10 Soil sandy wet 10−2 : 10−1 10–30 Table 2.1: Typical range of EM parameters of various materials measured at 100 MHz [23]. In most common of applications, the key parameters are the permittivity and the conductivity, whereas variations of the permeability parameter are not significant [45]. Besides, as GPR operates best with low electrical loss materials, if conductivity is equal to 0, the working depth will be greater. The dielectric properties are critical GPR parameters, since they will 12 control the velocity of the EM waves through the material, and they will affect the horizontal and vertical resolutions. Then, a GPR survey can correctly be planned knowing these values previously mentioned, as well as the resulting radargrams will be properly interpreted. The velocity of propagation v of the EM wave is defined by equation 2.8. v = ω β (2.8) where ω is the angular frequency (ω = 2πf) and β is the phase parame- ter. The attenuation coefficient α is expressed in equation 2.9, whereas the phase parameter β is defined according to the equation 2.10. α = ω √ µε ( 1 2 [√ 1 + ( σ ωε )1/2] − 1 )1/2 (2.9) β = ω √ µε ( 1 2 [√ 1 + ( σ ωε )1/2] + 1 )1/2 (2.10) There is a simplifying assumption which is made for low-loss materials [45]. In this case, the attenuation is approximated as the equation 2.11 and the velocity of propagation can be approximated as the equation 2.12. α = σ 2 √ µ ε (2.11) v = c √ εr (2.12) where c is the velocity of light in free space, defined by equation 2.13. c = 1 √ µ0ε0 (2.13) Then, the conductivity has a large effect on attenuation [54] and the velocity is dependent on the permittivity of the material. Finally, the intrinsic impedance η is defined according to equation 2.14. η = √ µ ε (2.14) Besides, when an EM wave arrives at an interface between two materials 13 with different intrinsic impedance, part of the wave is reflected back to the source, whereas the other part is refracted [56]. The reflected signal is described by the reflection coefficient r (2.16). r = η2 − η1 η2 + η1 (2.15) where η1 and η2 are the impedances of the mediums 1 and 2 respectively. Then, when the material is low-lossy, the expression can be simplified as described in equation 2.16. r = √ εr2 − √ εr1√ εr2 + √ εr1 (2.16) Figure 2.1: Incident, reflected and refracted waves from a layer with con- trasting permittivities [4] A typical example can be observed in Fig. 2.1. Then, when εr2 > εr1, r has a positive value. This situation can be found in an air-filled void existing in a dielectric material [23]. Besides, the relationship between the velocities and the angles of the waves can be expressed in equation 2.17. v1 v2 = sin(φ1) sin(φ2) (2.17) These characteristics related to stone products have been investigated by several authors [98] [118] [70] [119] [101] [111]. However, in most practical situations the EM parameters will be unknown, so the determination of these properties remains largely experimental. Besides, stone is a complex 14 medium which is composed of many materials in varying proportion, so within different types of stone the dielectric parameters will vary [23]. Then, it can be concluded that, for the type of products this Thesis is focused on, the permittivity varies between 7 and 9. In real cases, it is also possible to estimate this value using processing software, by manual fitting an hyperbola. 2.1.2 Ground-Penetrating Radar GPR is a non-destructive evaluation technique which uses radio waves to inspect dielectric materials (Fig. 2.2). Figure 2.2: GPR basic scheme [6]. Originally, this method was applied to natural geologic materials. How- ever, today it is used extensively for a variety of applications in many fields, such as environmental [79], forensic [25], archaeological investigations [12], building construction [74] [99] [9], landmine detection [39], pipes or cable de- tection [80], sedimentology [67] and it is able to analyze numerous materials, such as concrete [42], wood [91], rock [124], soil [31], etc. 15 Basic Principles Figure 2.3: GPR typical reflection schematic [22] The most common use of GPR employs a display unit, that provides a display of the recorded data, a data collection unit and an emitting and a receiver antennas, which are moved along the surface of a material, in order to detect reflections of any type of target. Most commercial GPRs are shielded bistatic radar, that is, they have an emitting and a receiver 16 antennas. These bistatic antennas can separate its position for common mid-point (CMP) surveys, in order to determine the velocity of the wave through the material. On the other hand, shielded antennas are used to reduce interference from external sources which can produce EM radiation, being more robust than unshielded antennas. The fundamentals of GPR are based on the study of the propagation of polarized high-frequency EM waves (25 MHz - 4 GHz), while an antenna moves across the surface, emitting short EM pulses. In Fig. 2.3 a basic example of a typical reflection schematic is shown. In the top plot the emitting and the receiving antenna are placed at the top of a material, which is composed of a material and a buried object, and they present different dielectric properties. The EM wave is radiated by the emitting antenna and travels through the medium, with a velocity which is dependent on the permittivity of the material. When the wave hits a medium with different dielectric properties, it is reflected or scattered back to the surface, and then the receiving antenna records it for later interpretation. Besides, part of the wave energy continues to travel downward (refracted wave). The process is repeated until the EM wave is completely attenuated. Moreover, the greater the contrast between the two materials, the stronger the reflected signal. In the middle top of Fig. 2.3, the recorded trace is shown. This signal is known as A-scan or trace, where the only variable is the time, directly related to the depth of the EM wave through the medium. Then, as GPR antenna is moved along the distance, a set of traces constitute a B-scan, also referred as radargram, which is shown in the bottom-right top of Fig. 2.3. This typical image uses different colors to represent the amplitude of each point. In Fig. 2.4 a whole scans description is exposed, where in the top plot an A-scan is showed, in the middle plot a B-scan or radargram, as a collection of A-scans, is presented, and, finally, in the bottom plot a C-scan is displayed, where the spatial thee-dimensional data can be used to reconstruct a 3D image. Then, A-scan, B-scan and C-scan provide 1D, 2D and 3D data respectively. It is common to work with B-scans, and to apply processing to A-scans. 17 Figure 2.4: Scans description [22] GPR operating parameters A GPR survey is planed considering the expected subsurface material prop- erties, the nature and the depth of the targets and possible limitations [125]. Then, the performance of GPR will be dependent of the frequency and the selected operating parameters. One important factor is the time window, which is the period during which each A-scan is measured by the receiving antenna, and it needs to 18 be accurately adjusted to penetrate the desired depth, depending on the velocity of the wave through the material. Moreover, the relationship between the velocity of the wave and the material properties is the fundamental basis of GPR. Since the velocity of two mediums differs when they present different dielectric parameters, when a wave travels through these two materials, the travel times will be different. Knowing the velocity value, for an homogeneus medium, the depth or thickness can be measured from 2.18. d = vt 2 (2.18) where t denotes the two-way travel time between the surface and the target and d is the resulting distance to the target. Thus, time window can be approximately adjusted using this relation. However, estimating the velocity of the wave through a real material is not trivial. In most trial situations, the relative permittivity will be unknown and the velocity of propagation needs to be measured in situ [23]. On the other hand, wavelength λ is given by the relation between the velocity of the wave v and the selected frequency as indicated in equation 2.19. λ = v fc (2.19) where fc is the nominal frequency. There is a direct relationship between the frequency of the transmit- ted wave and the resolution. Besides, there is an inverse relationship be- tween frequency and penetration depth. Then, high frequency antennas are used for searching of small and/or shallow targets, whereas low frequency antennas are employed for detecting big and/or deep targets. Moreover, the medium and the target properties should be taken into account when selecting the transmitted frequency. Table 2.2 summarizes how the GPR parameters and the EM properties affect a GPR survey. Then, the choice of the nominal frequency fc is based on the trade-off be- tween desired depth and resolution [24] [89], as well as on the EM properties of the medium. The higher the nominal frequency, the less the wave is able to penetrate the material. Besides, the conductivity, as explained in Section 2.1.1, is proportional to the attenuation. Then, materials which present high 19 Permittivity Conductivity Frequency GPR Parameter low high low high low high Remark Velocity fast slow Velocity is high in materials such as dry sand, and slow in water- saturated materials. Attenuation low high Signal attenuation is influenced strongly by electrical conductivity at high frequencies. Penetration long short The lower the attenuation, the greater the penetration distance. Wavelength long short Short wavelengths are normally used for concrete structures; long wavelengths are applied to mapping geological layers. Resolution low high The shorter wavelength, the higher resolution of subsurface targets. Table 2.2: Relationship between EM parameters, operating parameters and the selected frequency [125]. conductivity suffer a high attenuation, so the penetration decreases. Moreover, GPR signals are characterized by the bandwidth to centre frequency ratio R (equation 2.20) [45]. R = B fc (2.20) The goal is to maximize the bandwidth B and minimize fc. Since R ≈ 1 in most commercial GPRs, they are usually characterized by the center frequency. In addition, resolution is a fundamental term which is related to the ability to distinguish two targets from one another [88], as can be seen in Fig. 2.5. The vertical resolution ∆v (in depth) is considered as defined in equation 2.21. ∆v = λ/4 (2.21) This resolution is independent on the distance between the antenna and the target in a real world. The lateral resolution ∆h (in distance) will be related to the properties of the radar wave, the EM parameters of the medium and the distance between the antenna and the target. Then, it can be calculated as defined in equation 2.22 [45]. ∆h = √ rλ/2 (2.22) 20 where r is the separation between the GPR antenna and the target. Some authors have investigated vertical and horizontal resolution in real cases [78]. Thus, the detection of any target by the GPR will become difficult, if it is very small compared to the wavelength, if it present dielectric properties very similar to the medium or if the distance between the target and the antenna is very small. Figure 2.5: Resolution scheme [88]. Besides, data acquisition system should satisfy the Nyquist sampling criteria, in order to avoid aliasing. For a given sinusoid frequency f , the sampling interval (or step size, in distance) ∆x must obey ∆x ≤ v/2f . It is common to use values that are half as large, so the sampling interval will be ∆x ≤ v/4f . For GPR signals with a bandwidth to nominal frequency ratio equal to unity, this is translated to ∆x ≤ λ/3 [45]. Then, it is recommended to select a value lower than satisfy ∆x ≤ λ/6. Besides, if the selected sampling interval is very small, it is necessary a slow speed of the movement of the antenna. Clutter Clutter is a noise which affects the GPR radargrams and it can be defined as signals that are unrelated to the target but are present in the same time window, with similar spectral properties. The reasons of this noise are diverse and they can be caused by reflections between the emitting antenna and the surface, interferences between the emitting and receiving antenna, which is called as cross-talk or by scattered signals of other objects [23]. This effect is usually presented as horizontal straps and it is common that the target energy of interest is masked by it. In general, clutter noise is 21 more significant at short range and decreases at longer times. A radargram which presents clutter can be seen in Fig. 2.6, where horizontal straps in the samples 40-60 in depth. This noise presents high energy and masks the real anomalies in the radargram (at a depth of 400 and 600 samples on the right). Moreover, the example presented in Fig. 2.7 clearly shows that the clutter influence increases with increasing frequency. Figure 2.6: Example of a radargram with clutter. 22 Figure 2.7: Example of the clutter frequency concept. At 100 MHz the clutter is clearly visible, whereas at 50 MHz much of the clutter is suppressed [6]. Methods which suppress the clutter effects are known as Background re- moval techniques. An exhaustive review of state-of-art background removal methods is presented in the Section 2.3.2 in this Chapter. Advantages and disadvantages of GPR The most important advantage of GPR is its non-destructive nature [16]. Moreover, it presents the capability of detecting targets from several cen- timeters to hundred of meters, with the appropriate antenna and configura- tion (according to the survey design). Besides, this technique is really fast of using and with low cost, allowing an on-site evaluation of the current situation. The antennas of a GPR do not need to be in contact with the surface of the earth [23]. Finally, it is able to detect many types of mate- rials beneath the subsurface, responding to both metallic and non-metallic objects [22]. However, one important limitation of GPR is that the antenna emits the signal in a 3D cone, so the receiving antenna can record data from reflections which can be originated from anywhere around this cone. Then, 23 when a reflection is recorded, the radargram is not a picture directly beneath the survey point and this can lead to problems with the interpretation of the depth and shape of the detected targets. Moreover, the performance of GPR is limited to the environment and with some types of material the emitted EM wave is rapidly attenuated. The selection of the operating parameters will be dependent on the material to study and the nature, depth and environment of the target to find. Finally, this technique is really sensitive to unwanted signals caused by the system, other EM transmissions or other factors [16]. 2.2 The use of GPR in the stone industry Determining the internal structure of a quarry can reduce extraction costs as well as improve its exploitation. Besides, once the stone blocks are extracted, the detection of possible defects can reduce transport and cutting costs and can improve the production of slabs. Then, the publications related to this paper [P1], [P2] and [P3] analyse stone quarries and blocks quality. Thus, this section briefly reviews the use of GPR to assess the quality of stone materials, specifically, quarries and stone blocks. Besides, a review of methods to evaluate stone materials is included. 2.2.1 GPR applied to evaluate the quality related to quarries GPR has been used in numerous works [104] [55] [113] [40] [63] [105] [29] [81] [73] [58] [7] [82] [47] [34] [60] [68] [48] in order to inspect the quality of stone quarries. Sigurdsson [104] conducted one of the first studies over quarries, per- forming lithographic characterisations in limestone quarries. Kong and By [55] discussed the advantages of using GPR with several applications, in- cluding the analysis of a granite quarry. This work proved the utility of this technique for detecting cracks in stone materials. Tillard and Dubois [113] used GPR for analysing granitic and limestone quarries in order to deter- mine the wave propagation velocity. Besides, the author outlined the utility of this technique for distinguish the depth or thickness of the multiple beds inside the same quarry, in order to detect two kinds of limestone. GPR was also applied by Grandjean and Gourry [40] for detection of fracture zones in a marble quarry in Greece. To carry out this study, the authors employed 24 frequencies of 900 and 300 MHz, obtaining suitable results and concluding that the resolution and penetration reached with these antennas were rea- sonably good. Meschede et. al [63] also used GPR with a 300 MHz antenna to investigate a limestone quarry, with promising results. Sigurdsson and Overgaard [105] conducted studies in limestone quarries, where zones and textural variations could be differentiated. Derovert and Abraham [29] com- bined GPR and seismic imaging to investigate a gypsum quarry using a 500 MHz antenna. Pipan et. al [81] employed GPR for analysing the bedding planes and fractures in limestone in Italy, using antennas of 50, 100, 200 and 250 MHz and applying pre-processing techniques. Moreover, GPR was used by Orlando [73] for doing a semi-quantitative evaluation of massive rock quality. Lualdi and Zanzi [58] developed a method to explore the po- tential of GPR applied to a limestone quarry in Italy, using antennas of 50 and 250 MHz. The results were compared with the geological structures. This work proved that GPR can be used to evaluate the extension and the potential of a quarry and to prevent extraction activities in non-profitable areas, since the limestone layers are clearly detected (Fig. 2.8). Apel and Dezelic [7] used a 1500 MHz antenna for mapping a quarries, differentiating the dolomite and limestone rocks. In [82] Porsani et. al introduced a method for detecting fractures, joints and massive blocks in an ornamental granite quarry in Brazil using antennas of 25, 50 and 100 MHz were employed. These frequencies are able to distinguish fractures and to localize massive blocks, information which can serve as a guide for planning the advance of the quarry. Kadioglu [47] used the technique for determining changes in layer thickness and fractures or air- or moisture-filled cavities in a marble quarry in Turkey. This study is performed using a 250 MHz antenna. More- over, Forte et. al [34] applied this technique to perform a general evaluation of limestone characteristics. Luodes [60] established a characterization of different natural stone assessment and their material properties using GPR. This technique is applied to granite, schist, soapstone and marble. The GPR measurements clearly show the fracturing of the rock (Fig. 2.9)). Results indicated the difficulty in detecting vertical and sub-vertical fractures due to their small reflections Nielsen et. al [68] developed a 3D vision of a lime- stone quarry to illuminate the architecture and growth patterns of carbonate mounds. Finally, Kana et. al [48] employed GPR for analysing a limestone quarry in order to characterize the filling material and the aperture of the 25 fractures. Figure 2.8: The structure of a quarry and a 250 MHz radargram, where the layers are clearly highlighted [58]. Figure 2.9: A schematic presentation of two GPR radargram which show the propagation of the fractures in a quarry [60]. Summarizing, the previous works suggest that GPR can be an effective and fast tool to evaluate the general quality of a stone quarry, detecting its layers and fractures. Besides, antennas of 100-250 MHz provide a good resolution with the necessary depth. 26 2.2.2 GPR applied to evaluate the quality related to stone blocks The use of GPR for analysing the quality of extracted stone blocks is not as extensive as of quarries studies. In the following, a short revision of GPR studies applied to blocks [58] [95] [8] [10] [11] is presented. Lualdi and Zanzi [58] developed a method to evaluate the volume of extracted stone blocks and to detect the areas affected by fractures using a 250 MHz antenna. The geometry of the detected fractures can be seen in Fig. 2.10. This work proves the utility of GPR for increasing the productivity and optimizing the cutting of the blocks. Sambuelli [95] used GPR for detecting thin fractures in marble blocks, in order to optimize their cutting, using a 2 GHz antenna. The author analyzed the amplitude and phase of the reflections to differentiate the fracture filling. Arosio et. al [8] established a method of quality control applied to stone blocks, in order to detect internal defects, using a 2 GHz GPR system and employing a simple data processing. Fig. 2.11 shows the utility of GPR for locating fractures inside stone blocks. Finally, Arosio et. al [10] [11] addressed a fracture characterization with GPR, using deterministic deconvolution, applied to small block samples. The fracture parameters which are characterized are related to the thickness and the filling material, by processing the amplitude and phase spectra of the thin bed reflection. The results suggest that GPR can be a fast and effective tool to analyse the fracture parameters. 27 Figure 2.10: Cross sections of a stone block, where the horizontal and dip- ping fractures are marked by arrows [58]. Figure 2.11: 3D map of a thin discontinuity present on a marble block [8]. This short review of works related to the use of GPR applied to block evaluation proves that GPR can be a suitable method to analyze the internal structure of stone blocks. Besides, the selected frequencies can be higher than the antennas employed in the quarry analysis, so the resolution can be improved. 28 2.2.3 GPR applied to general stone assessment Hereafter, a brief review of studies conducted by several authors to assess stone materials is presented here. Toshioka et. al [116] carried out measurements of vertical rock walls in order to design a mapping crack distribution. Roch et. al [90] conducted a study on the assessment of rock-fall hazards using a 100 MHz antenna, proving the potential of GPR for rock analysis. Theune et. al [112] used a 50 MHz GPR antenna for discovering highly fractured rocks. Deparis et. al [27] designed a fracture characterization method of unstable cliffs using a 100 MHz antenna, with a penetration of 25 m. Dorn et. al [32] performed a fracture reconstruction using 100 and 250 MHz GPR antennas applied to a granitic rock aquifer. Longoni et. al [56] designed a method to evaluate the characterization of a fractured rock mass using GPR. The proposed method integrates diverse techniques used to define the rock mass fracture pattern: joint orientation, spacing, and persistence. An example of its performance can be observed in Fig. 2.12. Finally, Forte et. al [35] used a 3D dataset to construct a tool to characterize sediments and rock masses. 29 Figure 2.12: Example of the method proposed by Longoni et. al: raw radar- gram (top left), processed radargram with the discontinuities (top right), interpolated discontinuities (bottom left) and fitting planes (bottom right) [56]. 2.3 Signal Processing applied to GPR The goal of applying signal processing methods to GPR is to achieve a more reliable radargram that provides an user’s opinion without ambiguity. As a result, it is compulsory to remove any noise and clutter existing in the radargram. Some processing methods are applied to each A-scan, since each trace can present its own properties that must be treated separately. However, there are other techniques which are used for B-scans or even for C-scans. 2.3.1 Pre-processing A pre-processing method is a routine which is applied over each trace of the radargram to correct anomalies from the recording. The most common used methods [16] [29] [36] [39] [40] [50] to pre-process a raw radargram (Fig. 2.13.(a)) are described as follows. 30 • DC Removal: this technique consists on deleting the constant com- ponent of the trace, in case there is one. Then, the method slides a window from the beginning to the end of the A-scan, subtracting the mean or the median within it from the centre value of such window. This component does not indicate any useful information, so it should be removed. An example is shown in Fig. 2.13.(b). • Background Removal: it is the subtraction of the mean determined in a selected window, in order to remove clutter, which usually blocks up the desired signal. After applying this technique, the real anomalies of the radargram can be detected, as can be observed in Fig. 2.13.(c). In this radargram, the clutter at 25 m of depth is almost completely removed. • Gain: this step is necessary in order to equalize the amplitude of the emitted wave, which suffers a significant attenuation along the medium. Automatic Gain Control (AGC) is usually selected to do this equalization, where, for each position of a running window, the average amplitude is calculated, and then the signal value in the selected point (which is usually the beginning of such window) is divided by this coefficient. Other types of gains include constant gain, exponential gain compensation. In all cases, the selection of the size of the window and the maximum gain allowed are crucial to get the best possible result. A clear example is shown in Fig. 2.13.(d), where it can be observed that the energy at the end of the radargram is now equalized. • Filtering: finally, the unwanted frequencies, due to system or human- induced noise, should be removed. Besides, there are many types of filters, from simple Band-Pass to sophisticated domain and transform filters. It is usually common to use Band-Pass filters, selecting the half part and the double of the nominal frequency as the upper and lower cut-off frequency points respectively. This final step is shown in Fig. 2.13.(e). 31 (a) Raw radargram (b) Radargram after applying DC removal (c) Radargram after applying Background re- moval (d) Radargram after applying Gain (e) Radargram after applying Filtering Figure 2.13: An example of a radargram after applying each of the pre- processing methods above mentioned. 32 2.3.2 Background Removal Background removal is the procedure of subtracting from the total field the part of clutter (or background noise), which has been introduced in Section 2.1.2. One of the most reference algorithm to remove the background noise present in the radargrams is the average trace subtraction (ATS) [108] [45] [69]. This method only takes the mean of several traces within a window and removes it from each trace. ATS assumes that this noise is constant. For this reason, it is not robust when evaluating rapid variations of the noise. Then, ATS could degrade the true anisotropies present in the radargram. Besides, other used estimator is the median trace. Moreover, by using a median filter, the algorithm can be made adaptive. Moreover, Singular Value Decomposition (SVD) is a commonly used background method [19] [16], known algorithm for diagonalisation of rectan- gular matrices. Using this technique decomposed the data into eigenimages and quantify how much a decomposed eigenimage correlates with the origi- nal dataset. If the data is severely contaminated by background noise, this noise can be regarded as a component which is the most correlated with the data. However, the components that can be extracted with this technique are always orthogonal, and they cannot be easily interpreted. A reference technique is the Principal Component Analysis (PCA) [122] [16] [2] [103]. This algorithm can be performed via several approaches, such as eigenvalue decomposition of covariance matrix or SVD. In this case, PCA decomposes measured data into two orthogonal matrices that provide information about the main components present in a dataset. Besides, it provides a third matrix where the diagonal elements indicate the amount of information of each principal component. This decomposition can be done in one dimension (1DPCA) or two dimensions (2DPCA). Problems related to orthogonality have an impact on PCA. 33 (a) Raw radargram (b) Radargram after applying ATS using mean estimator. (c) Radargram after applying ATS using me- dian estimator. (d) Radargram after applying median filter. (e) Radargram after applying PCA. Figure 2.14: An example of a radargram after applying background removal techniques [115]. Fig. 2.14 shows the most used background methods [115], where the comparison based on visual inspection suggests that in this example the more efficient method is PCA. Other approaches have been developed to improve the removal of the horizontal noise, such as Independent Component Analysis (ICA) [128], fil- tering [127] [59] [49] [83] [130], likelihood processing [65], deterministic de- convolution [123], wavelet domain [1], [44], parametric system identification [28] and eigenimage processing [51] [50]. Recently, some authors have published novel background removal meth- ods [85] [86]. For instance, a new method was developed by Rashed M. [85] known as Background Matrix Subtraction (BMS). This technique calculates the complete background matrix, with windowing, sample exclusion, weight- 34 ing and iteration. Specifically, at each depth of a vertical GPR section, it applies a 1D window and calculates an alpha-trimmed mean. Subsequently, it excludes any sample within this window which has a different sign from this mean. The samples are weighted according to their closeness to the alpha-trimmed mean and normalized. Then, this window is slid horizon- tally and vertically to conform a background noise matrix and the whole process is iterated until all the residuals are removed. Besides, a Directional Total Variation Minimisation (DTVM) was pro- posed by Rashed E. [86], which reduces the global variation of the signal while preserving a close match to the original form. This method is based on the concept of total variation (TV) minimisation, which is a commonly used approach for image restoration. Then, its main goal is to extract the clutter from the acquired data by forcing the TV regularisation term over the horizontal direction. Nevertheless, most algorithms have been designed to remove a purely horizontal noise, with no possible variation in its amplitude or sign. Thus, they have a suitable performance only with constant noise. 2.3.3 Migration It is well-known that the interpretation of a radargram is subjective and it is strictly dependent on the experience and training of the final user. Migration is the term used by geophysicists to describe the process of focusing the radargrams to more closely resemble the physical target dimensions [38]. It has been a basic tool for interpreters since at least the 1940s [97]. Due to the beamwidth of the antenna [120], targets appear as hyperbolas. Migration, based on spatial deconvolution, re-locates the apparent positions of the reflections (known as dipping reflectors) to their true positions [76] and collapses diffractions [126] resulting in a migrated radargram. This process increases spatial resolution, since a migrated image can offer an idea of the shape and the dimensions of the target [23], increasing target detectability and reducing the effect of the beamwidth. Migration requires the true medium velocity. If an incorrect velocity model is selected, the migrated section could be misleading. With a higher migration velocity than the real medium velocity, the diffraction curve is inverted more and more, and the image is overmigrated (taking the shape of a smile [129]). Rather, with lower velocities, the radargram is undermigrated 35 (the resulting shape can be distinguished as a frown). Although numerous methods have been proposed in the literature, the most reference approaches are based on Kirchoff anf F-K methods [126]. • Kirchhoff In Fig. 2.15.(a) the physical principles of migration are described us- ing the well-known harbor example [126]. There is a storm barrier at z3 from the beach and a gap in the barrier. Now, a calm breeze comes from the ocean as a plane incident wave, and the wavefront is parallel to the storm barrier. Then, walking along the beach, a different wave- front is seen. Thus, the gap acts as a secondary source and generates the semicircular wavefront along the beach. The experiment can be translated to a recorded time section (Fig. 2.15.(b)). The gap is known as a point aperture or Huygens’ secondary source, which responds to a plane incident wave and generates a semicircular wavefront in the (x, z) plane. Thus, the response in the (x, t) plane is the diffraction hyperbola. (a) (b) Figure 2.15: The harbor example shows that the gap acts as a Huygens’ secondary source (a) and the waves generated by Huygens’ secondary source have a hyperbolic traveltime trajectory (b) [126]. 36 Kirchhoff technique is based on the diffraction summation method, which is a summation of amplitudes along the hyperbolic trajectory. The curvature of the hyperbola is dependent on the velocity function. The process is based on the semicircle (the diffraction curve) superpo- sition, mapping the amplitude at a sample in the (x, t) plane onto a semicircle in the (x, z) plane. The migrated section will consist on the superposition of multiple semicircles. An example can be analysed in Fig. 2.16 [126], where a raw section is shown in Fig 2.16.(a). The apex of the hyperbola A corresponds to time t0. Then, the amplitude on point B is mapped onto apex A (Fig 2.16.(b)), following the hyperbolic traveltime equation 2.23. t2x = t20 + 4x2/v2rms (2.23) where vrms is the root mean square (RMS) velocity at A in t0. Figure 2.16: Kirchhoff Migration. (a) Raw radargram, where a target ap- pears as the characteristic hyperbola. (b) Radargram after Kirchhoff [126]. Besides, it is necessary to consider several factors which are associated 37 with the amplitude and phase of the waveform along the hyperbola. Firstly, considering the example in Fig. 2.15(a), the wave amplitude at location A on the z-axis is stronger than the amplitude at B, which is at an oblique angle from the z-axis. Then, this angle dependence of amplitudes should be corrected, which is known as the obliquity factor. Secondly, the wave amplitudes suffer a spherical divergence. The wave energy decays as 1/r2sw, where rsw is the distance from the source to the wavefront, while amplitudes decay as 1/rsw. Thus, all of the ampli- tudes should be scaled by this spherical spreading factor. Finally, the wavelet shaping factor, which is an inherent property of Huygens’ sec- ondary source waveform should restore the phase and the amplitude of the hyperbolic paths. When the diffraction summation method which has been explained incorporates the obliquity, the spherical spreading and the wavelet shaping factors, it is called Kirchhoff method. In practice, the choice of the aperture width and the maximum dip is critical to preserve dipping events and to decrease computational cost. In Fig. 2.17 a diffraction hyperbola and the resulting migrations using four aperture widths are shown. When the selected aperture is small, the migration is less capable in collapsing the hyperbola. Besides, velocity errors can occur as explained before, as can be observed in Fig. 2.18. When the selected velocity is lower than the medium velocity, undermigration appears (Fig. 2.18(a)), whereas if the velocity is higher than the real medium velocity it causes overmigration of the diffraction hyperbola (Fig. 2.18(b)). 38 Figure 2.17: Test examples of aperture widths in Kirchhoff migration [126]. 39 (a) Example of Kirchhoff migration using a lower velocity than the medium velocity (b) Example of Kirchhoff migration using a higher velocity than the medium velocity Figure 2.18: Test examples of velocity errors in Kirchhoff migration [126]. • F-K (Stolt) Stolt [109] [110] incorporated the Fourier Transform in migration. This method includes the frequency-wavenumber (ω−k) range. Originally, this method was based on a constant velocity assumption. 40 A F-K migration example [126] is represented in Fig. 2.19, where velocity is assumed to be 1. A pre-migrated target is shown in Fig. 2.19.(a) by the radial line 0B, where the vertical axis is the frequency axis ω. Then, F-K migration maps lines of constant frequency (AB) are changed to circles (AB′) in the new plane (kz, kx). Then, mapping is completed, and the dipping event 0B has been mapped to 0B′ after migration (Fig.2.19.(b)). Figure 2.19: F-K Migration. (a) A dipping reflector is represented by 0B in the (f, k) plane. (b) After migration, the radial line 0B has been mapped to 0B′ [126]. Now the diffraction hyperbola is examined, using the example shown in Fig. 2.20, which is collapsed to the apex after the F-K migration. This hyperbola is composed of a series of dipping segments, such as A, B, C, D and E, where A is the apex and E is the steepest dip along the asymptotes. Then, in the (f, k) space, the A segment is mapped along the frequency axis, the other segments are mapped along radial lines, increasingly away from the frequency axis. The last segment E is mapped along the radial line, which represents the boundary between 41 the propagation and the evanescent region (located at or greater than 90 degrees from the vertical). Figure 2.20: The diffraction hyperbola in the (x, t) plane is mapped onto a triangular area in the (f, k) plane [126]. If the medium velocity is constant, F-K migration can be expressed as a direct mapping from temporal frequency to vertical wavenumber, as observed in Fig. 2.19. However, Stolt extended his method to handle velocity variations. Due to the fact that it is Fourier based, it is the fastest known migra- tion method [13]. However, velocity errors can also cause under- or overmigrated radargrams (all the examples for the Kirchhoff migration presented in Fig. 2.18 can also be considered for the F-K method), together with problems related to the maximum dip selected. 2.3.4 Pattern recognition Therefore, more emphasis should be now placed on classification of pat- terns in radargrams, in order to discriminate different buried objects or anisotropies of a material. Then, some authors [3] [33] [114] [77] have de- signed pattern recognition algorithms applied to these GPR images. Al-Nuaimy et. al [3] proposed a system for automatically detect buried objects with GPR. The diagram of the processing steps is shown in Fig. 2.21. First, the radargram is pre-processed applying background removal, path loss compensation and low-pass filtering. Then, pattern recognition is used to discriminate the buried targets and the unwanted signals. Thus, a neural network classification is applied. 42 Figure 2.21: Flow diagram of the proposed method by Al-Nuaimy [3]. Falorni et. al [33] conducted a processing method to extract hyperbolic patterns inside the radargrams. Ting-jun et. al [114] also performed a model to extract hyperbolic signatures in GPR. Sezgin [102] developed a 2D template matching method for establishing a discrimination of buried objects, following the flow diagram of the Fig. 2.22. First, starting from the B-scan image, a background subtraction is applied and then, the radargram is normalized to positive values between 0 and 255. Then, as image should be binarized, the author applies a thresholding using the Minimum Error Thresholding [53]. Subsequently, a morphological pro- cess is used, applying an skelonization process. The whole process is applied to the templates and to the B-scan images. Finally, the author calculates the correlation between the templates and the radargrams using the Euclidean distance measure. The evaluation is performed using radargrams of buried metallic objects, with promising classification results. 43 Figure 2.22: Flow diagram of the proposed method by Sezgin [102]. In addition, Pasolli et. al [77] introduced a pattern recognition system to classify buried objects (Fig. 2.23). After a pre-processing step, the radar- gram is thresholded. Then, the system automatically detects the objects by means of a search of linear or hyperbolic patterns. These templates are formulated using a genetic optimization framework, where the apex and the curvature are modelled. Finally, a support vector machine approach classi- fies the object. The performance of the proposed method is evaluated using a synthetic database varying the number of buried objects (their position, size, shape and material). 44 Figure 2.23: Flow diagram of the proposed method by Pasolli et. al [102]. 2.4 3D modelling in GPR A 3D model of a stone product can offer a better understanding of its internal structure, offering a vision of the depth, shape and the spatial orientations of the detected fractures. However, the diffusion of 3D radar investigations is limited. One of the main problems is related to the acquisition technique: the grid of measurement data as well as the selected frequency will affect the accuracy of the model. Other relevant issue is the selected 3D processing software. In this section, a review of authors [40] [60] [47] [68] [100] [57] who have used GPR to create 3D models is presented. Grandjean and Gourry [40] developed a 3D fracture mapping. To ob- tain this type of representation, the author proposed the correlation of the fractures from each profile with the nearest one. Luodes [60] established a measurement grid on a quarry in two directions, offering an schematic 3D visualization, where the propagation of the fractures is showed in Fig. 2.24. Later, Kadioglu [47] presented a simple 3D visualization to determine changes in layer thickness and discontinuities. For this purpose, parallel 2D profiles are measured and displayed in a 3D grid. Nielsen et. al [68] collected a 3D GPR dataset over a limestone quarry using a 100 MHz antenna mount- ing the different profiles along the positions. Sengondan et. al [100] created a 3D reconstruction for the location of foundations in demolished buildings, comparing its results with the output of ReflexW software [87]. For this 45 purpose, the author used a simulated correlation algorithm (SIMCA) which is based on a comparison between the trace that would be returned by an ideal point reflector and the actual trace. Lualdi and Zanzi [57] carried out a GPR survey to create a 3D reconstruction for a cultural heritage application, improving the acquisition of the data to ensure the quality of the results. The authors stated that the 3D reconstruction is successful particularly in the detection of the hidden structure of timber floors or roofs. Besides, the system should ensure a correct antenna position and orientation, covering a regular grid of measuring points. Figure 2.24: 3D data volume which highlights the fractures at different depths ((a), (b) and (c)) [47]. 46 Chapter 3 Methodology This chapter describes the methodology used for each publication presented in this Thesis. First, the experimental work, with synthetic and real data, is detailed. Then, the setup, the involved materials, the databases and the metrics are described for each publication. 3.1 Experimental work description The acquisition system, for both synthetic and real data, is completely de- scribed for each of the associated publications. 3.1.1 Simulation software When a a signal processing method is proposed, a synthetic database is essential for training, as well as for evaluating its performance. Several simulations have been conducted using the modelling software ReflexW [87] [84] in [P2] and [P3]. This widely used commercial package [30] applies the Finite-Difference Time-Domain (FDTD) method to simulate EM waves. This interactive software is used to define layer boundaries, while the EM parameters are defined along such boundaries. These parameters are entered within a table. Besides, a random-layer option can be used in order to specify statistic parameters for a random perturbation of the physical parameters of the layers, using different spatial distributions. The objective is to set a configuration similar as a real GPR system, and with conditions similar to a real scenario. This simulated data has been used for designing the method and to test the results as argued before. 47 Moreover, as explained in Section 2.1 in Chapter 2, the EM parameters should be adjusted to the real properties of this type of stone. 3.1.2 Equipment for real acquisition data In this Section, the real acquisition system for [P1], [P2] and [P3] is de- scribed. A MALAGS Pro-Ex model RAMAC [61] GPR system has been employed to conduct every survey. The configuration has been modified in according to the selected antenna, the material and the desired depth to evaluate. For quarries studies, antennas of 100, 250 and 800 MHz were used. Re- sults indicated that the use of 250 MHz antenna is an optimal choice in order to analyze quarries due to its good resolution and accurate depth. Focusing on the analysis applied to extracted stone blocks, an antenna of frequency equal to 800 MHz is selected, due to its resolution and because the reached depth covers the size of the studied blocks. 3.2 Characterization of ornamental stone The following description is related to publication [P1]. This work is only focused to real data. 3.2.1 Setup Since real data of different depths need to be evaluated, three different an- tennas will be used, taking into account that for quarries, which present a higher depth, it was necessary to use low frequency antennas, whereas for blocks, which have a lower depth, the highest frequency antenna has been employed. The operating parameters were selected, as indicated in Table 3.1. Operating Parameter Value Antenna nominal frequencies 100, 250 and 800 MHz Sampling interval 0.03 m Sampling frequencies 1100, 2610 and 7964 MHz Time windows 410, 185 and 47 ns Table 3.1: Setup of the real data acquisition related to publication [P1] 48 For each antenna nominal frequency a particular sampling frequency (about ten times the nominal frequency) and a certain time window are selected. It is considered that the depth decreases as the nominal frequency increases. Besides, the velocity (about 100 mm/ns) has been experimentally estimated. Finally, the sampling interval has been calculated in order to satisfy the Nyquist criteria. 3.2.2 Material The study is focused on white-coloured marble quarries, known commercially as White Macael. Besides, two types of carbonated lithologies have been employed: a clear limestone known as Crema Marfil (White Ivory) and a highly porous red limestone, known as Red Travertino. The last lithology has been selected because the anisotropies cause problems in the exploitation of the quarry and in the cutting process. The White Macael quarry face presents two differentiated zones: a mar- ble unit and a mica schits unit. Besides, it has a sub-vertical cavity. The three blocks (Macael Marble, Crema Marfil and Red Travertino) present no defects. 3.2.3 Database On the one hand, several profiles have been obtained at the White Macael quarry faces using the 100, 250 and 800 MHz antennas. On the other hand, GPR has been applied to extracted blocks of the three types (White Macael, White Ivory and Red Travertino) using the 800 MHz antenna. Note that several signal pre-processing methods, e.g., filtering and gain, has been applied in the tests using the software RadExplorer [84]. 3.3 Detection and Classification of fractures The following description is related to publication [P2]. Synthetic and real data acquisition is presented. 49 3.3.1 Setup Simulated data These synthetic radargrams are used for training and for the evaluation of the proposed algorithm. Besides, two databases, one for synthetic data and other for real data, are created, in order to check the results. Since the blocks present a depth which is covered by the 800 MHz an- tenna, this frequency is selected for the analysis, since it provides the highest possible resolution with our equipment, and it is necessary to simulate with similar characteristics to the real measurements. In Table 3.2 the selected parameters are shown. Operating Parameter Value Antenna nominal frequency 800 MHz Sampling interval 0.005 m Sampling frequency 62.5 GHz Permittivity of stone 8 Permittivity of air 1 Time window 35 ns Table 3.2: Setup of the simulation related to publication [P2] The selected time window has been calculated according to the velocity (about 106 mm/ns) and the depth of the simulated profile. The permittivity of the stone is selected in accordance with the real material which is ana- lyzed. The sampling interval satisfies the Nyquist criteria. The sampling frequency is automatically selected by the ReflexW software. Finally, the time window is enough to cover the whole depth of the simulated block. Real data In this case, only a 800 MHz antenna has been employed, since it is the highest frequency of our equipment, covering the necessary depth with the highest resolution. The survey design is described in Table 3.3. Since the two measured blocks present different dimensions the time window has been adjusted to covering their corresponding depths. Besides, the sampling interval satisfies the Nyquist theorem. The sampling frequency is selected as ten times the nominal frequency. 50 Operating Parameter Value Antenna nominal frequency 800 MHz Sampling interval 0.026 m Sampling frequency 8 GHz Time windows 36 ns / 39 ns Table 3.3: Setup of the real data acquisition related to publication [P2] 3.3.2 Material Simulated data The objective of this simulation is to design profiles as closely to reality as possible with ReflexW. For this reason, the selected permittivity is 8. However, we can affirm, as mentioned in Section 2.1.1 in Chapter 2, that the value of permittivity of limestone varies between 7 and 9, with a normal distribution. Then, the part of stone of each profile is simulated as a layer with a random value between these limits. Then, the air-filled fracture, which presents different sizes, has been included, with a permittivity of 1. Real data With regard to real measurements, the studies that have been conducted in [P2] are focused on blocks of cream-coloured massive limestone, known as Crema Marfil (White Ivory), exploited in the Sierra de El Coto, next to the town of Pinoso (Alicante, Spain). Specifically, two blocks of different quality are tested. On one side, a top-quality block with dimensions of 160x300x140 cm is analysed. This block has an apparently good quality. The cutting process resulted in forty- six tables, and the report indicated that only four tables presented some minor problems (horizontal and diagonal filled fractures which did not cause breakage). On the other side, a second-quality block, with visible fractures, with dimensions of 120x260x160 cm is tested. It presented some visible fractures. After the cutting process, which resulted in fifty-three tables, several slabs presented major problems. Specifically, three tables suffered from an impor- tant vertical break, so they were completely broken in half. The intermediate tables had vertical fractures, but they did not cause breakage. Finally, the last table had a horizontal fracture, which caused the breakage of the slab.