Threshold Estimation Using Wavelets and Curvelets on Ground Penetrating Radar Data for Noise and Clutter Suppression by Dus t in Harr ison B . S . , South Dako ta School of Mines and Technology, 2000 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F Master of Applied Science in T H E F A C U L T Y O F G R A D U A T E S T U D I E S (Electr ical and Computer Engineering) The University of British Columbia M a y 2005 © Dus t in Harr ison, 2005 Abstract T h i s thesis provides an evaluation of the Redundant Discrete Wavelet Trans-form wi th applicat ion to the removal of additive white or colored Gaussian noise on a synthetic G P R signal. Special attention is given to the parameter that controls the number of decomposition levels. Eva lua t ion is performed using a level-dependent threshold to estimate and remove noise. Results are presented using noisy synthetic G r o u n d Penetrat ing Radar pulses to compare Wiener filtering and thresholding the Redundant and Non-redundant Discrete Wavelet transform. A d d i t i o n a l results are presented on the effects of choosing a number of decomposition levels using signal-to-noise ratio measurements, which suggest the importance of choosing this parameter. Recommendations are made and supported which determine the order of thresholding before or after the practice of trace averaging. Us ing G P R images, an application of a novel 2D threshold model in the newly discovered curvelet domain is compared to average trace subtraction. P romis ing results are presented on both synthetic and actual landmine data, which shows thresholding as a viable method of clutter suppression. 11 Contents Abstract ii Contents iii List of Tables vi List of Figures vii List of Abbreviations ix Acknowledgements x Dedication xi 1 Introduction 1 1.1 Research Objectives 4 1.2 Thesis Overview 5 2 Background 7 2.1 Ground Penetrating Radar 7 2.1.1 Principles of Operation 9 2.1.2 Data formatting and visualization 12 iii 2.1.3 Sources of Noise 13 2.2 Selected Review of G P R Signal Processing 17 2.2.1 Review of A-scan Processing Methods 18 2.2.2 Review of B-scan Processing Methods 20 2.3 Review of Wavelet Transforms 22 2.3.1 Discrete Wavelet Transform 23 2.3.2 Redundant Discrete Wavelet Transform 26 2.3.3 Curvelet Transform 28 2.4 Noise Removal using Threshold Est imat ions 31 3 Removing Noise from a GPR A-scan 35 3.1 Determining a Basis for A - S c a n Approx ima t ion 37 3.2 Determining the Number of Levels of Decomposi t ion 39 3.3 Es t imat ion of a Threshold from Noisy D a t a 41 3.4 G P R A - S c a n M o d e l 44 3.4.1 An tenna M o d e l 46 3.4.2 Noise M o d e l .' ' 47 3.5 Summary 49 4 Clutter Suppression in GPR B-Scans 50 4.1 B-scan M o d e l using Synthetic G P R Pulse 53 4.1.1 Clu t te r M o d e l 53 4.1.2 Target M o d e l 55 4.2 Proposed M e t h o d of Clu t te r Suppression 57 4.3 M e t h o d of Performance Evalua t ion 60 iv 5 Exper imenta l Results and Discussion 62 5.1 Experimental Results of Noise Removal in A-scans 63 5.1.1 Performance of Gaussian White Noise Removal 65 5.1.2 Performance of Gaussian Colored Noise Removal 67 5.1.3 Dependencies of Number of Decomposition Levels 67 5.1.4 Thresholding Multiple A-scans before Stacking 69 5.2 Discussion of Noise Removal Results 72 5.3 Simulation of Clutter Suppression in B-scans 74 5.3.1 Synthetic B-scan Experiments 75 5.3.2 Landmine B-scan Experiments : 78 5.4 Discussion of Clutter Suppression Results 81 6 Conclusions 83 6.1 Future Work 85 Bibl iography 87 A p p e n d i x A Image Plot t ing Funct ion 94 v List of Tables 3.1 Ent ropy measurement of different basis functions 38 5.1 Effect of noise coloring on SNRimp using the R D W T and D W T . . . 63 5.2 P S N R of proposed method and average trace subtract ion 75 v i List of Figures 1.1 Example Landmines in situ 3 2.1 Pictorial representation of F W H M 10 2.2 Examples of A , B and C scans 14 2.3 Pictorial representation of the D W T 25 2.4 Pictorial representation of the R D W T 27 2.5 Example of a Curvelet 31 3.1 Example showing level-constant thresholding 43 3.2 Estimating an example antenna response 47 3.3 Example of a typical received pulse in a Ultra-Wideband G P R system 48 4.1 Demonstration of curvelets versus wavelets 52 4.2 Example clutter '. 54 4.3 Example synthetic targets 56 5.1 Comparison of D W T and R D W T in white and colored noise 64 5.2 R D W T , D W T and Wiener filtering with A G W N 66 5.3 Risk decay of R D W T , D W T and Wiener A G C N removal 66 5.4 SNRimp of R D W T level-constant thresholding for A G W N 68 vii 5.5 SNRimp of R D W T level-constant thresholding for A G C N 69 5.6 Effect of thresholding before and after trace averaging 70 5.7 Results of stacking with median operator after thresholding 71 5.8 Understanding the choice of noise models for clutter suppression. . . 76 5.9 Comparison of residuals between proposed method and'ATS for point target 77 5.10 Comparison of residuals between proposed method and ATS for large target 79 5.11 Proposed method demonstrated on P M N landmine data. . . . . . . . 80 vm List of Abbreviations A G C N Additive Gaussian Colored Noise A G W N Additive Gaussian White Noise ATS Average Trace Subtraction DSP - D i g i t a l Signal Processing-DWT Discrete Wavelet Transform F W H M Full Width Half Maximum G C N Gaussian Colored Noise GPR Ground Penetrating Radar G W N Gaussian White Noise HumDem Humanitarian Demining JRC Joint Research Centre PSNR Peak Signal-to-Noise Ratio RDWT Redundant Discrete Wavelet Transform SNR Signal-to-Noise Ratio UWB Ultra-wideband ix Acknowledgements • D r . M . R . Ito provided me wi th funding and academic advisement which facilitated my completion. • D r . Fel ix Her rmann teaches an inspir ing course on wavelets, which motivated this research and his great knowledge of the topic gave me the techniques and insight to finish. • T h e Univers i ty of B r i t i s h C o l u m b i a also saw enough potential in my determi-nation and skills to provide addit ional funding for this project. • Emmanue l Candes, Laurent Demanet, D a v i d Donoho, Lex ing Y i n g deserve gracious thanks for making CurveLab available for research. • Joint Research Centre for having the bri l l iant foresight to facilitate human-i tar ian landmine detection research by the creation of the J R C Landmines Signature Database [11]. W i t h o u t this generous support I would not have been able to complete such an arduous task. D U S T I N H A R R I S O N The University of British Columbia May 2005 x To my wife, D o r a L i for her patience and love and my sister, A l i s sa for showing me I can do it xi Chapter 1 Introduction There are mil l ions of landmines worldwide, which are difficult to find and remove. The concept of humanitarian demining, which is sometimes referred to as H u m D e m , is to provide simple, sustainable and affordable method of detection and removal of landmines. T h e groups that implement these programs are often non-governmental organizations who d id not place the mines and generally have l i t t le information on mine type or placement locations. T h e most common and highly visible aspect of demining are the activities that in poor countries devastated by war. T h e objective of H u m D e m groups in these countries are to develop programs wi th in that country that allows present and future landmines to be found and destroyed even after the H u m D e m group has left. The sustainability objective of many H u m D e m programs is cr i t ical to imple-menting a method of landmine clearance that can be performed by locally trained demining groups. Previous estimations state that 100 mi l l ion or more landmines are currently laid in the world, however that number has been reduced to 60 mi l l ion or less [61]. T h e number of landmines, danger of detection and the removal proce-dures increase the cost and t ime requirements to nearly impossible levels for poor 1 countries. The cost of detection and removal is also related to the amount of t ime a landmine is in the ground. Delays in landmine clearance causes vegetation coverage, which increases the challenge to detect minefields and mine locations. Due to the large amount of money and time required for a successful H u m D e m project, new technology is not well accepted. Therefore making simple, affordable and mature technology is necessary so that projects can continue even after external funding has ceased.' H u m D e m programs determine sui tabi l i ty of technology by the cost and status of development maturi ty. T h e most readily used and implemented technology is human protection devices, metal detectors and trained dogs. Other technologies that have been tested or evaluated wi th respect to H u m D e m are: • Explosives detection using chemical sensors, biological sensors and nuclear quadrupole resonance • Opt ica l sensors such as infrared, hyper-spectral, visible and laser ranging • Electromagnetic sensors including ultra-low frequency, microwave, S A R and ground penetrating radar • Acoust ic sensors such as water-jet echo [50], impulse and ultrasound Bruschin i and Gros discuss in detail most of the above technologies in [5]. In particular, the authors suggest that G P R is a mature technology and that it is "near ready" for use in H u m D e m programs. The H u m D e m field is not receptive to new technologies, regardless of the promises of detection improvements [3]. For example, estimates for C a m b o d i a sug-gest that mine clearance is happening at a rate of 15 square kilometers per year wi th the possibil i ty of demining all land that is immediately in need for settlement, 2 (a) P M N Landmine (b) V S - 5 0 Landmine Figure 1.1: Anti-personnel landmines in situ, obtained from the Canadian Forces Landmine Action Database [20]. critical development or agricultural uses taking place in 5-10 years [2, 51]. Technol-ogy proponents suggest that these numbers could be drastically reduced, but the reality is that field tested technology is a difficult and slow process. Specifically in Cambodia, heat and moisture can dramatically reduce the lifespan or functionality of any electronics making new and untested technology dangerous to deminers. Not all applications of demining take place in such demanding environments. The declining cost and size of GPR components, such as analog/digital converters, amplifiers and high speed switches contribute to the increasing interest to imple-menting GPR as a landmine detection device. Additionally, technology improve-ments in computer processors and specialized DSP chips allow additional realtime and field-based signal processing of GPR surveys. While the processing of GPR is maturing, the limited resolution of GPR systems make the task of landmine detection difficult and ambiguous. When using GPR for HumDem tasks, the standard method of expert interpretation is generally not available. Experts are too expensive to train and support in these situations. Offline processing of GPR data is an option, which can allow computationally diffi-3 cult algorithms such as imaging and feature discrimination to be done off the field. Offline processing, however, requires a thorough understanding of noise present in G P R data. Two aspects of G P R data are examined in this thesis. The first is the ability of ID wavelet thresholding to remove Gaussian noise from raw G P R data. The second is to suppress horizontal clutter present in G P R data using thresholding with a new transform, similar to wavelets, called curvelets. Both methods attempt to separate signal from noise by means of thresholding. The received data from G P R is a one-dimensional trace that contains re-flections of a transmitted pulse. To study this, a G P R trace is synthesized from a reasonable physical model of a G P R pulse. The type of noise and amount of noise energy added is controlled allowing numerical measurements of wavelet thresholding performance. A separate case, in two-dimensions, applies to an image of received G P R traces. This image contains strong horizontal clutter events, which detract from im-age understanding. The effect of thresholding clutter events in the curvelet domain is evaluated using an image formed from synthetic G P R traces and actual G P R data from the J R C Landmine Signatures database. 1.1 Research Objectives While ID wavelet noise removal using thresholding is well established, the applica-tion to G P R is minimal. The type of wavelet decomposition, specific parameters and the choice of a wavelet basis are difficult to understand and is not well estab-lished. This research examines the most promising wavelet decomposition method for noise removal, which is the Redundant Discrete Wavelet Transform (RDWT) . 4 The hypothesis is that the method of level-constant thresholding and the number of decomposition levels used by the R D W T will have a direct effect on the ability to remove additive Gaussian white and colored noise. The second objective is to compare the ubiquitous method of subtracting an average trace from a G P R image to a proposed method that uses a threshold in the curvelet domain as a method of subtraction. This new method allows a simple noise model, like the average trace subtraction method, but without the drawbacks of global averages and subtraction. The hypothesis is that the lack of resilience of subtraction requires clutter models that are too complex, and that thresholding is more effective at suppressing horizontal clutter using a simple noise model. Additionally, both approaches to G P R processing take into account that noise removal and clutter suppression should be implemented with fast algorithms. The advantage of wavelets are that most implementations are reasonably fast. The R D W T can be implemented very easily in modern DSP chips, making it a candidate for on-line processing. 1.2 Thesis Overview This thesis first introduces the basic mechanics of G P R technology and the unique obstacles along with a selection of signal processing techniques that have been ap-plied to G P R data. A brief wavelet introduction precedes the two wavelet algorithms used in the methodology of this research. The first method, known as the redundant discrete wavelet transform has been demonstrated as a successful transform where noise can removed by the application of a simple threshold in the wavelet domain. The following section introduces a 2D "wavelet" transform, known as curvelets, which is a recent development in multi-resolution basis function decomposition and 5 similar to, but not strictly, a wavelet transform. The transform properties and con-struction is described after a breif overview. The remainder of Chapter 2 introduces the threshold operator that is used to perform noise removal. The methodology is presented in Chapters 3 and 4, for A-scan noise removal and B-scan clutter suppression, respectively. Chapter 3 introduces a method for choosing bases, the parameter that controls the number of decomposition levels in the R D W T and threshold estimation. A description of the G P R trace model explains how a physical approximation of a G P R signal is made. At this point, the noise models is introduced, first additive Gaussian white noise or A G W N , and then a filtered form of A G W N , denoted as additive Gaussian colored noise (AGCN) . The clutter suppression methodology builds on the previous material, using the synthetic G P R pulse model to form a synthetic B-scan. A brief introduction to the landmine data is made, before a novel method of clutter suppression, inspired from the success of thresholding ID signals is introduced. The remainder of Chapter 4 describes the method of evaluation and performance of the new algorithm. The results present a comparison between a proven wavelet denoising tech-nique and the more successful R D W T method. Additionally, demonstrations show the difference in the optimal mean-square-error linear filtering technique and wavelets to remove either A G W N or A G C N . New results show the relationship of input noise energy and decomposition levels along with a demonstration of compatibility of the method to the current G P R processing method of stacking. Simulations of clutter suppression are presented for the synthetic G P R image data and followed by an ex-periment on actual G P R landmine data before the thesis concludes with a discussion of results. 6 Chapter 2 Background T h e topic of this thesis is evaluating the estimation of a threshold wi th respect to ground penetrating radar. T h e concepts that are needed to perform this evaluation are a basic physical understanding of the G P R data and how wavelets can be applied. T h e motivat ion is clear: addi t ional insight into G P R processing should serve to enhance the field of landmine detection. 2.1 G r o u n d P e n e t r a t i n g R a d a r Radar operates on the principle that the range from the radar to the target can be determined by t iming the round t r ip of a t ransmit ted pulse reflected by a target. T h e assumption of radar is that the transmission medium is conducive to passing the frequency content of the pulse. T h e earth is not a homogeneous medium, therefore anticipating the frequency content to pass through the media of the earth is location dependent. G r o u n d penetrating radar overcomes this by using a wide range of frequencies, which, allows sub-surface targets to be "ranged". A pulse transmission from an antenna at a single spatial location is received 7 over the durat ion of receiver gating. T h i s resulting t ime-domain signal is called a trace, or A-scan. For G P R data to be meaningful, a G P R survey usually attempts to spatial ly sample a grid of points on the surface. T h e collection of traces from a G P R survey of several spatial locations forms a three-dimensional data set of two spatial coordinates and one of t ime called a GPR data cube. T h e method of conducting a G P R survey varies on the type of G P R . A G P R may be bi-static, mono-static or neither. In mono-static mode, the same antenna transmits and receives the G P R pulse. T h e burden falls on the transceiver electronics to properly gate the transmit and receive times and protect the receiver from the large in i t ia l radiated pulse. In bi-static mode, separate antennas are used for t ransmit and receive, but their spatial relationship is held constant. The analog in seismic jargon for this method of operation is called a common-offset gather. Each trace represents a constant angle between the transmitter and receiver. In bi-static mode, antenna polar izat ion relative to transmit and receive can improve target detectability [43]. Addi t iona l ly , common mid-point gathers can be used to determine velocity models of the subsurface. In this survey, the angle between transmit and receive antenna is increased after each pulse transmission. M a n y practical applications of G P R are being developed. H igh resolution G P R solutions are used to image the internal structure of roads and bridges for the appearance of cracks. M i l i t a r y and security applications have been suggested to image potential threats behind walls of caves and buildings [6]. Snow pack layers, internal ice structure and avalanche v i c t i m location have also been demonstrated using G P R [40, 42]. Geophysical scientists make extensive use of G P R for the study of near surface features [15, 58]. 8 2.1.1 Principles of Operation There are many variations of G P R systems, which are defined by the type of antenna and method of pulse generation. Impulse radar generates a pulse in the time domain and uses an antenna, which is non-dispersive. The other common G P R system creates a pulse in the frequency domain and may use either a dispersive antenna to transmit a "chirp" pulse or a non-dispersive antenna. This type of radar is called Frequency Modulated Continuous Wave ( F M C W ) or a popular variation in HumDem research is Stepped Frequency Continuous Wave (SFCW) [18, 53]. These radars usually fall into a broader category known as Ultra-Wideband (UWB) radar. A U W B radar has a ratio of bandwidth to center frequency that exceeds one. This is called the fractional bandwidth and is defined by [62] as A U W B impulse radar has several advantages, because it uses low average power and has a wide range of frequencies. This allows it to image the subsurface of a wide variety of materials. Research by Cherniakov [12] shows that as the fractional bandwith of a radar approached the theoretical limit of 2 it obtained maximal spec-trum efficiency. Spectrum efficiency is defined as maximizing depth of penetration to enable detection of buried objects. It was also shown that a radar with Bp > 1.1 has a high spectrum efficiency, making U W B - G P R radars an excellent choice for depth of penetration when searching for underground events. A n impulse G P R in commercial applications commonly uses a Gaussian pulse to excite an antenna, This pulse can be created by using a transistor in avalanche mode, with results in (2.1) (2-2) 9 1.5 F W H M -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Figure 2.1: P ic to r i a l representation of the F W H M parameter of a Gaussian bump. an output that is an approximation of an exponential charge and discharge cycle. A typical ly G P R specification is the F u l l W i d t h , Ha l f M a x i m u m ( F W H M ) of a pulse. T h e F W H M is a common measurement of the w id th of a gaussian pulse and is defined by FWHM = 2 C T v / 2 1 o g e 2 (2.3) The w id th of a pulse indicates bandwidth , but usually the working bandwidth of a radar is l imi ted by the antenna. The response of an antenna to the input of (2.2) is to radiate a voltage wavefield which is the derivative wi th respect to t ime, V'(t) = ~ e - t 2 / 2 a \ (2.4) A n event is defined as the filtered form of (2.4) that occurs at some point in t ime of a received G P R trace. Th i s event corresponds is the reflection of the t ransmit ted pulse from a subsurface feature and can be related to depth if the velocity of the transmission media is known. In the case of a metall ic reflector the event is the second derivative of (2.2) because a metall ic reflector acts as a second t ransmit t ing antenna. 10 G P R antenna development is crucial to successful G P R applications. T h e antenna controls the interaction of the radiated waveform wi th the surface of the ground by impedance matching. T h e antenna also determines the frequency content of a radiated pulse, which dictates spectral efficiency thereby controll ing depth of penetration and target visibi l i ty. H u m D e m researchers have developed an antenna that exhibits excellent G P R performance wi th respect to spectral efficiency and ground coupling [59]. D a t a from this laboratory G P R and antenna is used to test the algori thm presented in Chapter 4. T h e characterization of an event in a G P R trace is determined by the t iming and ampli tude of the occurrence. Due to spherical spreading of the t ransmit ted wavefront and attenuation of the subsurface, the ampli tude of the received signal for a mono-static configuration is shown by "[17] to be proport ional to f 2adsec9 / o e\ a* sec* 0 6 • [ Z - b ) T h e distance, d, is the line, in meters, normal w i t h the top of the point reflector to the surface. A s the G P R survey is performed, the angle 9 increases, which further attenuates the signal. In practice, some commercial G P R systems attempt to compensate for this attenuation using an exponential gain applied dur ing or after data collection [39]. W h i l e it is the prerogative of the user interpreting the data to apply a gain, it is important that noise removal be done beforehand to avoid non-linear changes to the data statistics. T h e determination of t iming in a homogeneous medium is based on the hor-izontal location, x, of the antenna relative to the target and the depth, d, of the target. For a point target, or a target size that is near the resolution l imi t of the radar, the round t r ip t ime is determined by twice the velocity of pulse propagation 11 times the line of sight distance to the target t= -Vx2 + d2. (2.6) v The velocity of propagation is determined by the speed of light scaled by the relative permittivity of the ground, e v = C-. (2.7) If the subsurface is a combination of media, the determination of velocity can be-come intractable. Additionally, the influence of water on e is dramatic, capable of increasing the permittivity by an order of magnitude. 2.1.2 Data formatting and visualization In order to process G P R data, there are three common formats for presentation of G P R data. Each format chooses a set of dimensions from the 3D cube of G P R data. The most common format to perform image processing and understanding is a B-scan. B-scans are are two dimensional images with time on the vertical axis and a spatial dimension on the horizontal. Most often the spatial dimension is the "down-track", which is the direction that the antenna traveled as each trace is collected. Specific spatial coordinates in a B-Scan represent an A-scan. A n A-scan is a time-domain waveform that is measured from the voltage on the antenna after an initial pulse is transmitted into the ground. Typically an A-scan is measured in milli-Volts, but may also be presented in a normalized amplitude. The amount of time measured in an A-scan varies with the G P R system in use, and depends on the center frequency of the system. A low frequency G P R can typically penetrate further into the ground and requires a longer time to allow the pulse to reflect back to the antenna. 12 Due to the nature of a G P R pulse containing mult iple events, an A-scan is described as "wiggly" and has many smooth oscillations. T h i s descriptive of A-scans also refers to a method of viewing B-scans called a wiggle plot, which is actual A -scans vert ically oriented and plotted close together. T h i s method is used most often by human operators to make interpretations of G P R data. In Figure 2.2, a pixel intensity B-scan image is presented, where each pixel color is scaled based on the ampli tude in the A-scan . Occasionally C-scans are used in image processing. A C-scan is a "depth" slice, which takes the A-scan values over a l l spatial locations to form a intensity image. T h e C-scan contains both spatial dimensions for axes, which are the down-track and cross-track. T h e cross-track dimension is the spatial direction that the antenna is moved after a down-track survey is completed. A C-scan'is a subset of a B-scan, because it is generated from the set of al l B-scans. Each B-scan is, in turn, a subset of A-scans. 2.1.3 Sources of Noise There are many noise sources present in a G P R device. D u r i n g the creation of the in i t ia l source waveform, noise may be added by interference from external sources or poor high-frequency digi tal design of electronics. After a pulse is transmitted, cross-coupling from the antennas may generate spurious, large ampli tude events prior to any reflection events. In addit ion, upon reception of a G P R signal, clock j i t ter, analog-to-digital quantization noise and amplifier noise determine a noise floor and the operating dynamic range. A t the resolution of G P R , the majori ty of signal energy comes from large impedance changes in the ground. Changes in soil water content w i l l dramatical ly 13 A-Scan at x=25cm and y=25cm B-Scan at x=25cm Amplitude y (cm) Figure 2 .2 : Selected Examples of different scan types. Each contains different per-spectives of a G P R data cube. 1 4 attenuate the signal ampli tude as it propagates into the ground. Addi t iona l ly , if the subsurface contains large amounts of near wavelength objects, like glacial t i l l , depth of penetration and target localizat ion w i l l also decrease. Collectively, events that are due to actual reflections of the source waveform, but reduce target signal reflections are referred to as clutter. Positional Noise It is important to consider that the antenna is a physical object that is never precisely positioned. In the lateral sense, the antenna posit ion is cr i t ical if any transforms are taken on the data in the X - Y dimension. B o t h the wavelet and Fourier transforms expect data to be sampled on a fixed gr id . In order to correct this noise, antenna position must be recorded and the data interpolated onto a fixed gr id . In the application of landmine detection this becomes a problem wi th handheld detectors, which can be swung in arcs or irregular lines [32]. Mos t laboratory solutions use robots or antenna jigs to stabilize and accurately determine spatial location. In addit ion to lateral movement, it is important to consider vertical and angular antenna posit ion. B o t h angular and vertical posit ioning affect the coupling of energy into the ground. A s the coupling changes, the ampli tude of the received signal varies. T h e result is that algorithms must adapt to an unknown variance in amplitude. Addi t iona l ly , vertical posit ion uncertainty creates temporal variance of the received signal. Usua l ly this noise is minimized by allowing the antenna to travel flush to the ground, but as ground changes occur, the hyperbolic assumption of a target response can be invalidated. Often this issue is lessened by the durat ion of scans, which are taken over contours that can be considered flat. 15 Sys tem Noise The most well understood system noise is thermal noise. T h e r m a l noise can be mod-eled by additive Gaussian noise. The noise is the cumulative effect of the antenna gain, amplifier noise figure and the transmission line. It can be quickly reduced by trace averaging or stacking of mult iple traces from the same spatial location. C o m -mercial G P R systems implement this in hardware and w i l l provide the user w i th the abi l i ty to average up to 32 traces. Other system noise, which is less understood and might not be well modeled by Gaussian noise comes from electronic design and susceptibility. For example, poor circuit board design of a transmitter can add noise to the generated pulse. Reception of external signals may also corrupt the signal. T h e quantization noise of the analog-to-digital converters also l imits the dynamic range of the system. T h e combination of al l sources of noise is very likely Gaussian, but not necessarily white. C l u t t e r noise Clut te r is defined by the user and s tr ic t ly speaking is not noise. Clu t te r is the cumulative effects of reflections from actual boundaries that are i l luminated by the G P R , but impede wi th analysis of the pr imary targets under study. T h e most common form of clutter in any G P R application is the ground. T h e occurrence of a reflection from the ground offers very li t t le information to a G P R image and may interfere wi th near-surface objects. Landmines are an example of a near surface target that can be impeded by the presence of an air-ground event. Other forms of clutter, w i th respect to landmines, are targets that have prop-erties very close to the desired target and can result in a false positive determination of a pr imary target. These false positives can be metal fragments, stones and other 16 debris that produce events similar to a landmine event in either ampli tude, fre-quency content and shape or a combination of the three. Generally, this form of clutter is filtered using discr iminat ion algorithms. A discr iminat ion algori thm w i l l use some unique feature of the desired target to suppress al l other possible targets that do not contain any such features. Another form of clutter, sometimes called background noise, appears as the por t ion of a signal, which does not necessarily represent a physical object, but ap-pears as a slowly varying horizontal event in a B-scan. T h e most common form is mult iple reflection that occurs between the ground and antenna, often exacer-bated by antenna r inging. Subsurface planar reflectors are sometimes grouped w i t h background noise, because of the typical ly large horizontal component. T h e terms background noise and clutter are used interchangeably in the literature and are problem specific. Clu t te r noise models are often considered additive, such that the subsurface can be represented by a t ra in of dirac functions, each of which represent either target or clutter temporal locations. The G P R signal is then convolved wi th the dirac t ra in to obtain a simulated received trace. Clu t te r is normal ly the l imi t ing noise i n G P R systems [18]. Besides target depth, clutter is the l imi t ing factor that determines the detectability of a subsurface target. 2.2 Selected Review of G P R Signal Processing A l l G P R data is processed, whether by human or machine. In a simple case of u t i l i ty finding, it is usually enough for a human to interpret u t i l i ty location based on the existence of a hyperbola in an unprocessed B-scan. In the case for humani tar ian demining, it is usually necessary to perform addit ional machine-based processing 17 and discr iminat ion to perform accurate and reliable landmine detection. G P R processing can be divided into three categories, restoration, transfor-mat ion and discr iminat ion. Restoration is concerned wi th the removal of noise, correction of system error or correction of the sensor degradation. Transformation is used as a functional tool to perform signal processing, but in the case of migrat ion is also an output. Discr iminat ion is a key aspect to landmine detection. Di sc r imi -nation uses a feature, hopefully unique, to separate noise and clutter from targets. For discr iminat ion clutter often refers to anything that is not a desired target. T h i s definition varies depending on the features used to discriminate. A typical signal processing chain for landmine detection may include sev-eral restorative techniques, possibly a transformation and always a discr iminat ion step. T h e importance of developing an excellent understanding and usefulness of restoration tools is paramount to successful landmine detection. 2.2.1 Review of A-scan Processing Methods W h e n processing a G P R A-scan, most methods focus on signal restoration. Restora-tive methods include, DC-offset removal, Gaussian noise filtering and feature en-hancement. W i t h the exception of DC-offset removal, which is sometimes performed by the equipment, most restorative techniques are chosen based on the remaining operations of the signal processing chain. Feature enhancement is usually done for the benefit of an expert interpreter. These techniques are often times carried out at the discretion of expert, where the various outputs are visual ly compared. For A-scans, these methods are frequency filtering, averaging and gain adjustment. Frequency filtering consists of applying a low, high or band-pass filter to select certain features. A s a G P R signal travels into 18 the earth it is dominated by the time-variant low pass effects, which reduces the high-frequency content of late t ime events. H i g h pass filtering is performed when it is necessary to dewow an A-scan. A "wow" signal is the low frequency effects of transmitter and receiver coupling, which produces a slow decaying signal into the A-scan [39]. Noise filtering is not a standardized process, nor is there a correct choice for all applications. T h e only ubiquitous noise removal that is performed on A-scans is stacking. Stacking, or in G P R jargon, trace averaging, is the average of repetitive scans from the same spatial location. Commerc ia l G P R systems often include this form of noise reduction as a pr imary parameter choice when performing a G P R survey. T h e purpose of stacking is to remove system noise. Other methods include Wiener filtering, which assumes the noise model is well understood and attempts to perform an inverse noise filter in the frequency domain . T h i s method is also called m i n i m u m square error filtering, because it seeks to minimize the error measure, E[/ — f]2. T h i s filter seeks to invert the degradation response, H(u) in the frequency domain using [28] 1 \H{u)\2 F(u G{to). (2.8) [H(uS) \ H { U J ) \ 2 + K_ T h e parameter K, assumes additive white Gaussian noise and can be estimated from the variance of the noise. T h e result is to suppress the noise and degradation in the corrupted signal, G(OJ) and produce an estimate of the signal. / . Wiener filtering has been applied for different special cases, which are described in [17]. Discr imina t ion can also be performed on A-scans using parametric target models, frequency features, or t ime features. Discr imina t ion using t ime features is done wi th matched filtering and is usually not well suited to G P R , because of the t ime variant filtering of the subsurface. Frequency features can be discriminated by 19 using the unique spectral content of a target to suppress other events [37]. M c C l u r e and C a r i n explored the idea of compactly representing G P R waveforms using a matching pursuits a lgori thm in [52]. T h i s method attempted to find compact atoms that represented fundamental structures of radiated and reflected electromagnetic waveforms wi th success in approximating t ime and frequency parameters of G P R signals in noise, which proved useful for discr iminat ion. 2 . 2 . 2 Review of B-scan Processing Methods T h e majori ty of B-scan processing is focused w i t h clutter suppression and discr im-inat ion. Clu t te r suppression is a difficult and subjective topic that requires careful definition of what is clutter. Horizonta l clutter is well studied, w i th current research showing moderate success at removing the air-ground event and other persistent horizontal events. A common form of horizontal clutter reduction is to use average trace sub-tract ion ( A T S ) . Given a collection of traces located spatially at points in x and y. A T S uses an average trace for a part icular B-scan located at posit ion y, 1 N fy = ~N z\2 fxy (2.9) x=l to perform a subtract ion for a given B-scan, fx,y — fx,y fy (2.10) T h i s is effective in enhancing the hyperbolas or any non-horizontal event in G P R data. T h e drawback is that a l l horizontal events are removed or attenuated, there-fore i f a target contains horizontal information A T S is not suitable. Other examples focus only on the air-ground event, such as W u [68], where the air-ground event is 20 adaptively estimated using a shifted and scaled reference signal to perform subtrac-t ion. T h e nature of subtract ion is that A T S is not robust to phase changes or ampli tude changes. If a clutter event decays to zero, A T S w i l l continue to use the same static ampli tude calculated from the average to subtract. For computat ional cr i t ical applications, it can be a useful method and is used in some discr iminat ion algorithms in [69]. Gader et. al uses a variat ion that calculates the derivative in x instead of the average, which is used in the signal processing chain previous to discr iminat ion and is highly sensitive to noise [25, 26]. T h e major problem of A T S is the instabi l i ty when subtracting coherent noise. T h i s can lead to artifacts that can decrease the usabil i ty of a data set. U l r y c h et. al studied this method of clutter suppression wi th respect to other competing methods. In part icular, the a - t r immed mean was suggested as a way to handle variable clutter, but at the addi t ional cost of careful expert decisions to choose a. T h e median was also studied, which led to [66], a s tudy on the robustness of L-moments as estimators. New research was presented on the possibil i ty of using a James-Stein estimate for noise and clutter suppression on a B-scan. T h i s method uses a shrinkage parameter to adjust an local observation based on the global mean [65]. After B-scan processing, G P R data is often times used to form a 3D image. T h i s is referred to as migration and is the process of focusing a B-scan or G P R data cube into a coherent image. T h i s technique is sometimes performed in landmine detection, but due to computat ional complexity and resolution it is not popular. It does present the most understandable output for human understanding, because the hyperbolas are collapsed and the t ime axis is recast to depth [71, 46, 27, 43]. 21 2.3 Review of Wavelet Transforms Wavelets and Wavelet transforms are an excit ing topic of research, which was in-spired by the collaboration of Mor le t and Grossmann [30]. T h e principle of wavelets is that a signal can be represented by its frequency and temporal content. T h e re-lationship between t ime and frequency, necessarily, is a tradeoff of min ima l support in one domain demanding max ima l support in the other. In the world of signals, transient events present a challenge to estimation and compactly representation [48]. If a t ime series contains a mixture of signals w i th varying amount of support in both domains, there is a fundamental l imi t to the determination of a precise t ime location and frequency content for each component. T h e elegance of wavelets is an efficient representation that attempts to balance the inherent uncertainty between time and frequency. Wavelets are often introduced as a continuous transform, which underlies the fact that each discrete representation of wavelets makes a unique sacrifice to obtain efficient algorithms and signal representations. T h e continuous wavelet transform, qualitatively, is the par t i t ioning of a signal into t ime and frequency atoms by the application of translated bandpass filters. These atoms, ideally, attempt to maximize frequency and time localizat ion. Us ing the notat ion of Ma l l a t [48], a wavelet has a zero average such that /+ o o i/j{t)dt = 0. (2.11) -oo To obtain localization i n time, the wavelet is translated by u and scaled by s, l M 0 = ^ ( ~ ) • (2-12) T h e scale parameter is usually greater than zero. T h e wavelet transform can now 22 be defined as the convolution of the wavelet with the function under study, or The sampled continuous wavelet transform should not be confused with the discrete wavelet transform. Creative solutions to use the sampled continuous wavelet transform attempt to characterize the type of discontinuity. Mallat used the wavelet modulus maxima to estimate the Lipschitz exponent to classify a singularity [47]. Hsung extended this work by approximating signals without directly identifying the singularity [38]. The author uses a modified form of wavelet modulus maxima detection to estimate an inter-scale ratio and inter-scale difference of the sampled continuous wavelet transform coefficients to estimate the signal in noise. A wavelet transform can be described as a multi-resolution decomposition of a function. Resolution may not be precisely correct in the continuous domain, but is certainly applicable in the discrete. When the continuous wavelet transform is discretized, the scale parameter defines the time/frequency resolution of the cor-responding decomposition. This discrete multi-resolution decomposition will allow efficient signal processing to approximate, compress and estimate signals. 2.3.1 Discrete Wavelet Transform The orthogonal discrete wavelet transform is a well studied implementation of wavelets in the discrete domain. It became popular amongst applied scientists when a fast filter bank approach was introduced. Alongside the fast computation and the ability of the D W T to compactly represent mixed signal functions, researchers were also able to derive many analytical results that suggested profound implications not accessible before. (2.13) 23 Amongst these discoveries was the work that influenced the field of signal processing by Donoho that suggested non-linear thresholding could achieve near opt imal signal estimation in white noise [23]. T h e substance of that c la im relies on the orthogonality of the discrete wavelet transform. T h e continuous wavelet transform is far from orthogonal and offers no imme-diate discrete solutions. To solve this problem the discrete wavelet transform makes two dramatic changes to the choice of the scale and translation parameters s and u. T h e first change is to choose s by choosing the number of scales dyadically. T h e new discrete parameter, j is s = 2j, j eZ (2.14) and the discrete translation parameter, k is dependent on the scale u = 2jk, keZ. (2.15) Th i s choice of discretization of the continuous parameters is motivated by a fast implementation algori thm, but solidly rooted in sampling and information theory. Intuitively, the dyadic part i t ion of the scale parameter in the D W T allows a higher localization in low frequency, where t ime localizat ion is necessarily min ima l . T h e dependence of k on the discrete scale parameter, j balances this uncertainty at higher frequencies when signals become highly localized in time. T h i s choice of scaling and translation forms a set of wavelet basis functions. T h e dyadic choice of j and k results in N basis functions, which is equal to the length of the signal f[n]. T h e index, m is related to each translat ion and di la t ion of the discrete wavelet basis to form g[m\. T h e discrete wavelet transform can be represented by F[m} = (f,gm), (2.16) 24 JW] *(T) +Ai+l (fe) G ] - F i + 1 (fe) Intialize the i teration wi th : A°(k) = f(n) H = Scaling F i l t e r G = Wavelet F i l te r T h e wavelet coefficients are FJ for j = 1 , . . . , J and the scaling coefficients are FJ+1(k) = AJ{n). Figure 2.3: T h e D W T is implemented by an iterative application of high pass and low pass filters that are constructed by the wavelet and scaling function. where gm is an orthogonal set of basis functions, (9m-, 9n) = C m n S m n . (2-17) Often, wavelet bases are chosen such that Cmn = 1, which makes a basis orthonor-mal . However, the D W T is better understood as an implementation of cascaded digi ta l filters, which is shown in Figure 2.3. T h e l imi ta t ion of the D W T is the lack of shift invariance. For an intuit ive understanding, consider the second coarsest scale, which only contains two coeffi-cients that represent the signal frequency and temporal content at this part icular scale. These two coefficients are sufficient for perfect reconstruction of the signal using the D W T , but if approximations or estimations in noise are made the lack of shift invariance becomes evident. Signal estimations can be performed by thresh-25 olding which chooses a part icular coefficient based on the coefficients ampli tude. Shift invariance implies that a shifted version of the signal under study w i l l produce a shift in the wavelet coefficients, while not significantly changing the ampli tude of the coefficients. W i t h o u t shift invariance, the success of thresholding the D W T coefficients to remove noise w i l l depend on the signal's location in t ime. 2.3.2 Redundant Discrete Wavelet Transform T h e Redundant Discrete Wavelet Transform ( R D W T ) is a fast, convenient method of obtaining shift invariance. T h e term R D W T is ambiguous, because the nature of obtaining coefficient redundancy is not unique. A literature review suggests that the undecimated version of the discrete dyadic wavelet transform is the R D W T . T h i s is equivalent to other names, such as the Shift-Invariant Discrete Wavelet Transform ( S I D T W ) [45], Undecimated Wavelet Transform ( U D W T ) , Stat ionary Wavelet Transform ( S W T ) [54, 56] and the original method Algorithme a Trous by Ma l l a t . These methods are not equivalent to the Cyc leSpin [13] or Shifted Orthog-onal Wavelet Transforms, which obtain redundancy by means of shifting the input signal and re-calculating the wavelet transform. A diagram of the R D W T in Figure 2.4 shows the s imilar i ty to the D W T . Each level, j , is a bandpass filter of the original signal. T h e bandwidth of each filter is chosen dyadically, which is identical to the D W T . In the orthogonal D W T , the signal is decimated at each stage of decomposition, because the exact location at lower frequencies is ambiguous. However, without decimation the R D W T coefficients fully represent the uncertainty of low frequency events w i th the same sampling resolution of the signal, thus producing l inearly dependent coefficients. T h e choice to use a redundant representation is a double-edged sword. Proven 26 ^ + 1 ( k ) Gj Hi, CP A T 2 j Gi+1 Intialize the i teration wi th : A°(k) = f (n) H° = Scal ing Fi l te r G° = Wavelet F i l t e r The wavelet coefficients are. F- 7 for j = 1 , . . . , J and the scaling coefficients are FJ+1(k) = AJ{n). Figure 2.4: Rather then decimating the signal at the output of the filters, the filters are up-sampled which preserves the frequency resolution at al l wavelet scales. 27 results for the D W T do not transfer to the R D W T , though experiments suggest that the results are similar. T h e nature of a redundant transform is that the approx-imat ion of the signal is no longer unique. Therefore, the complexity of choosing the smallest set of coefficients for a given quali ty of reconstruction becomes N P -complete. Current research by Donoho [22] and Gr ibonva l [29] suggest that if certain conditions for sparsity are met, a set of coefficients can be chosen by the min imiza t ion of a linear program. Her rmann [34] has demonstrated these results by using linear programming to mainta in high quali ty reconstruction using a sparsity constraint when choosing a set of redundant coefficients. Tropp [64] has instead used a greedy algori thm to choose a set of coefficients. Another method used by P i z u r i c a [57] uses the redundancy of the R D W T to estimate a local spatial indicator. T h i s indicator allows coefficients to be chosen across scales, exploi t ing the similarities that are present at different scales. 2.3.3 Curvelet Transform Curvelets provide an excellent representation of G P R events in B-scans for signal processing. T h e nature of curvelets is that they efficiently represent a piecewise smooth curve which contain singularities normal to the direction of smoothness. T h i s property is the precise nature of G P R data, which suggests that curvelet co-efficients can provide a localized sparse representation where non-linear operations, such as thresholding, can be performed. T h e same arguments that apply in the I D wavelet case work in this 2D representation, namely that op t imal linear denois-ing methods, in the sense of m i n i m u m risk, can be further improved by non-linear estimators using curvelets [7]. 28 The seminal curvelet papers [7] and [8] describe the nature of curvelets and their properties. A n approach to an intuit ive understanding is to compare the shape of the curvelet atoms. T h e separable 2D wavelet bases are nearly isotropic and well contained wi th in a the support of the scale, j. However, curvelets are designed such that width = length2. T h i s makes a curvelet anisotropic and requires an addi t ional orientation parameter, 9. Curvelets also require a scale parameter, j , and posit ion vector, fc, to represent the location and frequency content. T h e bui ld ing blocks for sparse representation of curves can be made using these three parameters: scale, spatial location and orientation. To insure sparsity, the curvelets are constructed using a parabolic law scaling law, which provides a natural extension to smooth curves. Therefore, unlike 2D wavelets, which tend to represent a curve as a collection of points i n the wavelet domain, curvelets minimize the number of coefficients needed to represent the curve by adapting to the scale of the curve. T h e most important result presented in [7] is the rate of energy decay, in an L-2 sense, of a incomplete reconstruction of a signal, / , using the m largest ampli tude coefficients of a a curvelet decomposition is | | / - / £ | | B C m - 2 ( l o g n ) 3 . (2.18) T h e strength of this statement is that adaptive methods, which are not necessarily tractable can approach m~2 as m —> oo, whereas wavelets approach TO-1. Despite the log term, a curvelet representation approaches the adaptive l imi t . T h e benefit is that curvelets are not adaptive and do not require the iterative computations nec-essary to adapt to a signal. W i t h the results in [9], a new technique of constructing curvelets has introduced a tight frame, which is used in this work, and the same approximation rates hold . 29 Tight frames will preserve the energy in both domains, known as Parseval's relation for the Fourier domain and echoed here in the curvelet domain as E K / ' ^ > | 2 = H / I I M M 2 ) < V / e L 2 ( ] R 2 ) , ( 2 . 1 9 ) m which makes reconstruction straight forward by using the adjoint. As expected, the reconstruction is the same for curvelets as for wavelets, which is the sum of coefficients in the curvelet domain projected onto the basis functions, f = }^(f;lrn)lm. ( 2 . 2 0 ) m The partitioning of curvelets has some similarity to wavelets in terms of scale and translation, but not in orientation. When wavelets are extended to two dimensions, the frequency partitioning does not represent the fullness of the extra dimension. In a ID signal, a singularity can be described by a point in time. Two dimensional wavelets operate under the same premise, that singularities occur at localized points in an image. However, edges are singularities that are oriented .on a line. The angular partitioning of curvelets allows these singularities to be represented by a location and orientation. Wi th a wavelet representation, this edge would be represented by a collection of points rather than a smooth line. Curvelets overcome this inefficient representation through a frequency par-tition that includes not only dyadic frequency partitioning, but also a doubling of angular partitions at every other scale. Angular localization with respect to fre-quency is subject to the same uncertainty as temporal localization. Therefore, the angular parameter, 6, provides finer resolution at higher frequencies. The result is a partitioning, which allows a piecewise-smooth edge containing a singularity to be represented by a sparse set of coefficients. 3 0 0015 . 0 01 , Figure 2.5: A n image of a single curvelet is shown. The curvelet is directional and localized wi th a definite angular component. These new atoms can be described as "needle-like" at fine scales [9], whose rectangular dimensions follow the parabolic scaling law of width = length2. W h e n these functions overlap an edge that follows a smooth curve, they produce large coefficients. T h i s produces a sparse representation for edges, which can be adapted to the hyperbolic events of G P R data. 2.4 Noise Removal using Threshold Estimations T h e foundation for wavelets is la id and the extensions to various discrete imple-mentations are described. T h e next step explores the application of wavelets and curvelets to the estimation of a signal in the presence of noise. Statisticians refer to this problem as non-parametric regression, or in engineer speak: denoising. T h e goal of noise removal w i l l be to estimate a signal, / from a noisy obser-vation d. T h e model is additive so d = f + n , where n is the noise present in the data observation. If the data is transformed into the wavelet domain using a set of orthogonal basis functions gm, by D[m] = (d,gm), (2.21) then the under lying model is also transformed such that D[m} = (f,gm) + (n,gm). (2.22) If it was possible to know the value of N(m) it would be a simple matter of sub-traction. Instead, a decision operator is introduced, fi[m], which w i l l shrink each wavelet coefficient of the data. The application is f = }~2n{m}D[m}gm, (2.23) 771 which produces an estimate, / that can be measured wi th reference to / by the the risk: Mf) = E[\\f - / | | 2 ] . (2.24) T h e m i n i m u m risk can be obtained if the nature of the noise is known in the wavelet domain using oracle estimation [48] w i th a decision operator defined as n[m] = 2 . (2.25) 1 ' \D[m]\2 + a2 T h i s holds only if the noise is Gaussian and white w i th a zero mean. T h e remark-able result shown by Donoho [23] was that the decision operator defined in (2.26) produces a risk which is nearly min ima l w i th respect to oracle estimation. T h i s non-linear threshold operator would instead make a decision to either keep a coefficient that is above a threshold, T , or k i l l i t , by mul t ip ly ing wi th a zero. More formally, hard thresholding is, r if \D[m]\>T, (2.26) if \D[m]\ < T 32 T h i s produces a lower risk in the estimate / , but higher amounts of spurious os-cillations, where / contains transients. A n alternative to hard thresholding is soft thresholding. Soft thresholding applies the same principle of keep or k i l l , but re-moves the sharp discontinuities of hard thresholding by shr inking coefficients above T. Soft thresholding is defined as 0 \{\D[m]\<T, n[m] = { D[m] -T i f D[m] > T • (2-27) D[m} + T i f D[m]<-T Donoho showed in [23] that T could be chosen by s imply knowing the stan-dard deviation of the noise and the number of points, T = a^2logeN. (2.28) T h i s result is coupled wi th the abi l i ty to estimate the standard deviation from the finest scale coefficients of D[m]. If the wavelet coefficients are arranged, such that 0 < m < N/2 then an estimate of the standard deviation is a = median|D[m]| /0.6745, for 0 < m < 7V/2[48]. (2.29) These results led to a surge of research in the application of threshold estima-tors. T h e results were promising, but a nagging problem remained: shift invariance. T h e above results hold for additive Gaussian white noise thresholded from the D W T . Due to the lack of shift invariance, the resulting estimations were excellent in the sense of risk, but lacking in smoothness. In [45, 44], L a n g et. al were inspired by R . Coifman to threshold a shift invari-ant wavelet transform to remove additive Gaussian white noise. Us ing the R D W T , they d id a performance comparison on several of the I D and 2D signals. Compar ing 33 these results to Donoho's work [21] w i th soft thresholding they were able to show that as the energy of the additive noise decreased, thresholding the R D W T began to outperform soft thresholding of the D W T . These results were demonstrated by fixing a number of the parameters available in the algori thm. Donoho and Coifman also attempted their own versions of a shift-invariant wavelet threshold algori thm. In [13], the authors introduced the Sp inCycle a lgori thm which used the results from [1] to produce a set of redundant coefficients to threshold. These results also out-performed classic D W T thresholding methods. Thresholding in the wavelet framework is summarized by theoretical and ex-perimental results present in current research. T h e questions that have not been answered in present research are the significance of choosing the number of decom-posit ion levels and i f R D W T thresholding can be applied to colored Gaussian noise using a threshold estimator. T h e next two chapters address an aspect of these ques-tions by exploring methods of evaluating and apply ing a threshold estimator to G P R data. 34 Chapter 3 Removing Noise from a G P R A-scan The objective of the following methodology is to evaluate and understand specific practical questions that have not been addressed in the application of wavelet noise removal to G P R signals. Current research suggests that wavelet denoising can be effectively performed by redundant transforms, but it is still unclear how to choose a threshold estimator, scale parameter or basis. It is beyond the scope of this thesis to explore exhaustively each choice, but instead focus is placed on an the scale parameter after a wavelet basis and threshold estimator is chosen. The hypothesis is a level-constant threshold estimator applied to the R D W T transform of a G P R signal is a suitable method for removing white and colored Gaussian noise if the scale parameter is properly chosen. There are four ingredients, which are necessary to test the suitability of the thresholded R D W T for G P R signals. Each of the following sections discusses the specific choices based on the following questions: 35 • How can a synthetic G P R pulse be generated? • How can a threshold be estimated and applied to the R D W T ? • How can a wavelet basis be chosen? • How will the scale parameter affect Gaussian noise removal of the thresholded R D W T ? In order to evaluate the suitability of R D W T thresholding a direct compar-ison to the existing theoretical results is made. The performance is juxtaposed to the choice of thresholding the R D W T and Wiener filtering on G P R signals. The metrics provide comparisons can be made to the existing literature. Specifically, the theoretical results of wavelet literature focus on the use of a risk metric, introduced in (2.24). Proven results in wavelet denoising show the decrease in a normalized risk metric, 11/-/II2 r(f) TVa 2 ~ Na2' [ 6 ' This metric assumes that zero-mean Gaussian noise, has been added to the signal / . The difference of the signal and its estimate / are normalized by the length, N of the signal and the variance cr2 of the Gaussian noise. Evaluation of the maximum scale parameter, J is presented in the familiar engineering metric of signal-to-noise ratio, SNR = 101og 1 0 ( )^ . (3.2) ° V I I / - / I I 2 / When comparing the signal, / to the estimate, / , emphasis is placed on the improve-ment of the output signal-to-noise ration, SNR0Ut relative to the signal-to-noise ratio of the noisy input data, S N R i m p = S N R o u t - S N R i n . (3.3) 36 T h i s emphasizes the changes i n performance when experimental parameters are chosen. T h e experiments are conducted in M A T L A B using R ice Univers i ty ' s i m -plementation of the R D W T from the Rice Wavelet Toolbox [31] and the wavelets are generated wi th Wavelab [24]. A description of the ingredients to the R D W T thresholding process w i th a discussion of the rationale and relevance follows. T h e parameters of R D W T thresholding are explained before the G P R and noise model is described. 3.1 Determining a Basis for A-Scan Approximation In Ma l l a t ' s Wavelet Tour of Signal Processing he labels a chapter: E S T I M A T I O N S A R E A P P R O X I M A T I O N S . A m o n g the properties of the wavelet transform, one ele-gant result is that signals are sparsely represented for a given wavelet basis function. T h i s property makes estimation of a signal in noise simple, because the correlation of the signal to the wavelet basis provides high ampli tude coefficients, which are easily chosen v i a a threshold. T h e foundation of successful estimation relies on a small number of coefficients being an effective approximation of the signal. A basis function should be chosen such that the decomposition presents a sparse set of coefficients and the signal is well approximated wi th a smal l number of coefficients. D a t a adaptive methods [14, 49] choose a "best" basis from a dict ionary of bases using a sparsity measure to make a decision. Rather than form a best-basis, this work instead uses the log-energy dis t r ibut ion measure, N E | F [ m ] | 2 l o g | F H | 2 . (3.4) m—l to choose a part icular mother wavelet. 37 When picking a basis function, the support of the bases affect sparsity, but usually at a trade-off in the smoothness of reconstruction. This trade-off is related to the to wavelet construction, which demands that the wavelet support is determined by the number of vanishing moments. The wavelet has p vanishing moments if /+oo tvip{t)dt = 0. (3.5) -oo A n intuitive understanding of vanishing moments is that a wavelet acts a differentia-tor. If a singularity is differentiable q times, then a wavelet will be able to "detect" it when the vanishing moments, p is greater. If a wavelet is chosen for the number of vanishing moments, then the minimum support is 2p — 1. For G P R signals, the underlying signal is highly differentiable and it is more important to produce large amplitude coefficients, which make signal estimation possible. Using (3.4) as a guide to determine the best choice for a basis, minimal support should increase sparsity. Table 3.1 verifies that the Daubechies wavelet with 2 vanishing moments produces the sparsest representation with a log-entropy metric. Basis Entropy Daubechies4 -302.3 Symmlet4 -294.6 Coiflet2 -290.1 Batt lel -289.5 Symmlet6 -273.5 Coiflet3 -272.6 Daubechies6 -270.4 Coifletl -268.5 Symmlet5 -266.2 Battle3 -225.8 Table 3.1: Summary of entropy measurements on selected basis functions. Table 3.1 is generated by decomposing a 2048 point synthetic G P R trace 38 with one reflection event. Each entry represents a particular basis for decomposition calculated with the Wavelab function MakeONFilter. The sparsity of a decomposition does not guarantee the ability to approx-imate a signal. However, the alternative to using a sparsity metric is to find the smallest set of coefficients for a given approximation rate. Searching for this optimal set of coefficients which satisfies a certain quality of reconstruction is an NP-complete problem. Therefore, the sparsity property becomes an efficient and reasonable as-sumption for choosing a wavelet basis. 3.2 De te rmin ing the N u m b e r of Levels of Decompos i -t ion e In addition to determining a basis function, it is important to determine a useful number of decomposition levels. The maximum number of decomposition levels, J , determines how many different wavelet scales are calculated. It is possible to decompose a signal using the Redundant Discrete Wavelet Transform to a maximum of J = log 2 N levels. However, it is not clear if a signal should be decomposed to the maximum scale. As the parameter J increases, the support of the wavelet increases to a maximum of TV when J = log 2 N. Because the efficiency of a median operator increases with N and the advantage of thresholding decreases with sparsity, it is expected that the optimal value of J is less then log 2 N. There appears to be no research that suggests an appropriate value of J . The choice of scale parameter depends on the signal, noise and wavelet basis function. The wavelet basis controls the balance between frequency and temporal localization. The amount of energy that the noise and signal under study have at particular scale, 39 j determines if threshold estimation and application is appropriate. T h e evaluation of the scale parameter compares the signal-to-noise improve-ment calculated from (3.3). T h e number of levels, additive white noise and colored noise, and the energy of the noise is varied to il lustrate the scale parameter effect. T h e measurements are performed on a synthetic G P R trace. Eva lua t ion of scale parameter experiment: 1. Generate noisy data from the addit ion of a model and Gaussian Noise. 2. Decompose to J levels using the R D W T . 3. Est imate a threshold and apply. 4. Reconstruct the signal and measure SNRimp-5. Repeat steps 1-4 for a new realization of noise w i t h larger SNRin. 6. Repeat steps 1-5 to obtain an average SNR^mp. T h e reference model is a synthetic G P R pulse that models an approximation of a physically realistic version of equation (2.4). T w o different experiments are conducted, one to examine the effects of Gaussian white noise and the second for Gaussian colored noise. T h e coloration of the noise is described in Section 3.4.2, which assumes the difficult case where the noise has nearly the same characteristics as the signal. The input noise is measured wi th (3.2) and varies from -9dB up to 30dB in 3dB increments. T h e range of SNRin is representative of the signal attenuation as it propagates into the ground, which makes late t ime reflection events decay into the noise floor of the receiver. For each value of SNRin, the number of levels is 40 varied from the m i n i m u m of one to the m a x i m u m of l o g 2 N. Each experiment is conducted twenty times and the average SN-Rimp is calculated wi th the mean. The result form a surface where an op t imal number of levels can be chosen for the part icular signal and noise type. However, because choosing the number of levels is equivalent to bandpass filtering, a more generalized understanding can be implied based on the frequency content of the signal and the noise under study. 3.3 Estimation of a Threshold from Noisy Data T h e choice of a threshold estimator should only suppress a coefficient i f it is below a predicted value of noise. For the D W T a global threshold estimator can approxi-mated for al l coefficients. T h i s global threshold is val id, because the noise remains white in the wavelet domain wi th a constant standard deviation. W h e n thresholding the R D W T , the noise does not remain white and a threshold should be estimated for each coefficient individual ly . Donoho proposes in [23] that the opt imal threshold for the D W T is, T = c V 2 1 o g e / Y , (3.6) which, for certain classes of signals, achieves nearly min imax risk wi th respect to the opt imal non-linear decision operator. T h i s result makes the removal of white noise in the D W T domain a simple matter of estimating the standard deviat ion of the noise. Us ing the equation (2.29), an estimate of a is reasonable i f the signal is sparse relative to the noise. For the D W T , the finest scale provides N/2 coefficients to estimate the standard deviation. T h e support of the wavelet at the fine scale is very narrow, which makes the signal sparse relative to the noise. For the R D W T , each scale provides JV coefficients to make an estimate, but the sparsity of the signal •41 coefficients decreases wi th j. Therefore to generate a threshold estimator requires either, advance knowledge of the noise behavior or mult iple realizations of a signal wi th iden t ica l noise distributions. Us ing a hybr id method, which neither sets a threshold for each coefficient, nor uses a global threshold for al l coefficients w i l l be called a level-constant threshold. T h i s threshold is constant for a given level and is calculated from an estimate of the standard deviation of the coefficients at level j, Tj = ajy/logeN. (3.7) Th i s is applied to the wavelet coefficients using a hard thresholding technique given in equation (2.26) to the matching level of coefficients. T h e assumption is that the decomposition is reasonably sparse and an estimate of a w i l l remove the majori ty of noise without being too aggressive. T h i s assumption is tested by the scale parameter experiments outl ined above and a calculation of risk w i t h respect to the competing methods. T h e expectation is that the scale parameter has a direct effect on the abil i ty to remove noise. A s the J increases, the decreasing sparsity introduces more bias in the coefficients making the median estimator too aggressive. T h e alternative, which is to estimate a global threshold from only the finest coefficients is sufficient and successful' for white noise. The results of L a n g [44], Ma l l a t [48] and Donoho [21] justify this method of thresholding. However, i n the case of colored noise, the average standard deviat ion w i l l vary drastically for each level, j. Therefore, a global threshold does not apply. Thus , a level-constant threshold should prove to be a reasonable signal estimator i n the presence of colored noise. T h e verification of the level-constant threshold is tested against the perfor-mance of classic wavelet denoising using the D W T and Wiener filtering. A calcula-t ion of the normalized risk for the three methods is presented as increases. T h i s 42 Level 6 Level 5 Level 4 Level 3 Level 2 Level 1 Figure 3 .1: Using the D W T , a global threshold is estimated, but is ineffective for removing colored noise, which is clearly above the green threshold. A level-constant threshold approximates a threshold that is calculated for each coefficient, which adapts to the noise coloration. 43 allows a direct comparison of the proven D W T thresholding results in [48] to the R D W T level-constant threshold method. 3.4 G P R A-Scan Model A G P R A-scan can be measured as a voltage, Vrx from a receive antenna. T h e result is a convolution of several impulse responses from the environment and the system wi th a source function, Vs, Vrx = Vs(t)*hrx(t)*htx{t)*hc{t)*hg(t)*ht{t) + n{t). [18] (3.8) T h e impulse response from cross-coupling, hc(t) is described in [18] and is one source of horizontal clutter in a B-scan. T h e impulse response of the ground, hg(t) is quite complicated as it depends on the material , moisture content and granularity. Due to the variat ion wi th time, accurate ground models should be considered uniquely for specific applicat ion. A simplification of the ground model can be constructed using a constant attenuation for a time-invariant frequency response. T h i s methodology assumes the ground is constant for al l frequencies. Refer to the work of I rv ing [41] for a frequency dependent inversion of simple ground cases. The target impulse response, ht(t) is also complex, but simple reflectors such as the ground or point reflectors can be roughly approximated as specular. T h e input signal, V s is modeled as the derivative of a Gaussian bump like in (2.4). T h i s is a simplification, because a real pulse tends to be asymmetric and long-tailed, but depends upon the transmitter design and operating characteristics. T h e impulse responses of the receiving and t ransmit t ing antenna are the same for a G P R in mono-static operation, which makes hrx(t) = htx(t). T h e development 44 of a synthetic antenna model, hn(t) is described in Section 3.4.1. Us ing the antenna model , the source pulse is filtered by hn(t) to form a synthetic t ransmit ted pulse. T h e reflection events are specular, which corresponds to a dirac function in t ime making the final G P R pulse model, Vrx = Vs{t) * hn{t) * hn{t) * 8{t - U). (3.9) T h e target locations are located at the times in the array i j , which is the round t r ip t ime from the antenna. T h e pulse must be filtered twice by the antenna, once dur ing transmission and once upon reception, thus the two antenna responses in (3.9). T h e physical parameter, F W H M , used to construct Vs is 100ns, which is approximately the F W H M of the experimental G P R demonstrated by the Be lg ium Roya l M i l i t a r y Academy H u m D e m research group [59]. T h i s compliments the an-tenna model, which is a simulation of the same experimental G P R . The simplifications and assumptions made to produce the G P R impulse model are based on method of noise removal. Wavelets can be thought of as a sin-gulari ty detector. W i t h i n l imits , the singularity can undergo linear filtering without dramatical ly affecting the abil i ty of a wavelet to efficiently represent i t . There-fore, the simplification of the ground and target to specular reflectors is reasonable. Addi t ional ly , ignoring the cross-coupling term is justified, because it should have similar frequency and temporal support as reflection events making the results on the air-ground event transferable. 45 3.4.1 Antenna Model Bart Scheers, in his doctoral thesis [60], developed a normalized antenna impulse response hN.tx — \ ~—7=htx and hN>rx = W ——T==KX (3.10) V Z « V h V a V h which describes the frequency dependent effects of the ground, fg, impedance of the system and antenna Zc, Za and the transmission coefficient T t x , r r x . The goal of the synthetic antenna model is to simulate a reasonable G P R pulse in the time domain. For the purpose of determining the suitability of the R D W T for noise removal it is not critical to understand each parameter of the model. Instead, the experiment results from an antenna gain plot from Fig. 4-11 in Scheers' ([60]) thesis are used to calculate hpj directly. Using points that estimated by sight from the gain response of Antenna 2 of the referenced figure, a piecewise model of GI(UJ) is generated. The normalized frequency response, HN(LU) is related to the gain by 4-7T Gt(u) = -^\HN(u>)\2. (3.11) To find a the coefficients for a digital filter to approximate hpj, the piecewise model is fed into a Yule-Walker regression estimator to generate the coefficients bi, cn of a digital filter, HN(z) = • (3.12) Using (3.12) it is straight forward to find the normalized impulse response and consequently filter coefficients of the theoretical pulse. In Figure 3.2 the Yule-Walker fit is compared to the piece-wise estimation. The final assumption is that h,Ntix(t) = hpf.rx(t) a n d notation is simplified to hn(t). Using this model a M A T L A B script is created to filter the synthetic input 46 Simulated Antenna Gain Frequency, ui (Ghz) Figure 3.2: Us ing the Yule-Walker method to estimate an example antenna response. pulse by the simulated antenna response. Figure 3.3 shows the result of the physical approximat ion of the G P R pulse. 3.4.2 Noise Model T h e noise model, n(t) is straight forward to model in M A T L A B . To produce Gaus-sian white noise, the M A T L A B function r andn generates an arbi t rary mat r ix values which is normal ly distr ibuted, zero mean and has a standard deviation of one. T h i s is referred to as the additive Gaussian white noise, or A G W N model . T h e A G W N model is used to verify the performance of the R D W T level-constant threshold method to the classic method using discrete wavelet thresholding. T h e colored noise model is created using the antenna model from the previous section and G W N . T h e additive Gaussian colored noise, or A G C N is filtered white 47 Simulated GPR Pulse 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Time (ps) Figure 3.3: T h i s pulse shows an typica l pulse that might be received wi th only one reflector at 1000ns. noise and is represented as an operator on the additive noise model, d(t) = f(t) + hn(t)*n(t). (3.13) The antenna model ensures that the G C N is "in-band", meaning that the noise has a similar frequency content as the signal, / . T h i s presents a challenging case for the verification of the level-constant thresholding of the noisy data, d. T h e A G C N model is not unique and therefore does not represent an ex-haustive study of A G C N noise removal. Conceptually, an A G W N model has a flat spectrum. In the context of thresholding, a flat spectrum means that if the spec-tral content of the signal is less then the noise, it is thresholded to zero. For the A G C N , the spectral content is similar to the signal. T h i s makes the estimation and suppression much more difficult. These two cases should garner a general intuit ive understanding, which w i l l allow an extension to other cases. 48 3.5 Summary T h i s chapter introduced the methodology for testing the R D W T in the presence of Gaussian colored noise. The next chapter w i l l introduce the methodology for suppressing clutter in B-scans before the final results are presented. T h e main points from this chapter are: • A simulated G P R pulse is the signal under study. • A level-constant threshold is defined. • A wavelet basis is chosen based on sparsity. o A risk and signal-to-noise metric is used to evaluate results. • A noise model for both colored and white Gaussian noise is shown. T h e experimental results use these bui ld ing blocks to evaluate the effect of choosing the m a x i m u m number of decomposition levels on the abil i ty of the R D W T to estimate a G P R signal in white or colored noise. In order to establish a baseline performance, experiments are designed to compare the D W T and R D W T in the presence of white noise, before a performance evaluation on colored noise. 49 Chapter 4 Clutter Suppression in G P R B-Scans Clut ter suppression is a common step in the signal processing chain towards the goal of landmine detection. T h e purpose is to minimize damage of the target signal event, while suppressing anything that is not related to the target. There are several reasonably successful methods of clutter reduction present in current literature. For a general overview the book wri t ten by D . Daniels, Surface-Penetrating Radar [18] and a more comprehensive review in the recent book edited by the same author, •Ground Penetrating Radar [19] are recommended. A brief overview of clutter sup-pression methods w i l l be presented here to develop a motivat ion for the proposed method. Adapt ive methods, such as the one introduced by Brooks et. al [4] uses a linear time-invariant model to perform a non-parametric regression of the system response to separate clutter. Another adaptive method of clutter suppression at-tempts to suppress mult iple reflection events, where the author uses an L - l weighted 50 estimator instead of the t radi t ional least squares filter to derive a clutter model [33]. V a n K e m p e n et. al examines the data for a set of features to form a model and suppresses clutter w i th a modified Wiener filter [67]. A non-adaptive method is presented in [63] where a pr inc ipal component decomposition is used wi th linear estimation. T h e assumption is that the target is contained wi th in the largest principle components and a linear approximat ion w i l l achieve clutter suppression. A n alternative to regression models, is discr iminat ion, where clutter is suppressed by choosing a target, which forces al l other data to be accepted as clutter. Zhao proposes a successful method of discr iminat ion using a discrete hidden Markov model, which finds both metal and non-metal landmines in G P R data [70]. T h i s work is derived from the results of Her rmann in the seismic realm [36]. Her rmann uses the parsimonious curvelet representation of seismic data to itera-tively suppress mult iple target reflections. Herrmann's work tackles the difficult problem of selectively removing events that are nearly identical to the target, but smaller in ampli tude and later in t ime [34]. Herrmann's method has also been ex-tended to removing temporal clutter using mult iple data sets from the same location, but several years apart [16] and reducing the dimensionali ty of the seismic imaging problem [35]. T h e method proposed in this work also supplements the results presented by Nuzzo and Quar t a [55], which compares r-p and wavelet transforms for effec-tive horizontal clutter suppression. Nuzzo et. al used the 2D D W T to selectively threshold only the horizontal components. T h e two main features presented in their method are: • Identify the wavelet scales that clutter appears in . 51 Or ig ina l Curvelets: 41.6 d B Wavelets: 38.4 d B Figure 4.1: T h i s example demonstrates reconstruction of a simulated hyperbola wi th only 0.5% of the curvelet and wavelet coefficients. Notice how curvelets maintain the smoothness of curve and st i l l reproduce the discontinuity across the curve. • Suppress horizontal wavelet coefficients at identified scales. T h e algori thm proposed makes several important changes to this general approach. Fi rs t , the curvelet transform is used, instead of the 2D wavelet transform. Figure 4.1 demonstrates the motivat ion of using curvelets to approximate G P R reflection events in B-scans. Second, a noise model is generated using features from the data set, which intr insical ly define the wavelet scales and horizontal clutter features. Addi t iona l ly , the noise model forms the threshold operator, which puts the algori thm into the well researched and understood domain of wavelet thresholding. L ike clutter, the noise model should be determined for each applicat ion. 52 Therefore, a synthetic clutter and target model is formed from actual landmine data. To verify the noise model , the synthetic B-scan is generated using the physical approximation of A-scans developed i n Chapter 3. T h e description of clutter suppression methodology begins wi th the syn-thetic B-scan model . Afterwards, the clutter and target models are described and the proposed algori thm is introduced. T h e chapter concludes wi th a discussion of performance metrics for evaluation of the proposed method. 4.1 B-scan Model using Synthetic G P R Pulse T h e synthetic B-scan is carefully modeled such that reasonable clutter and target approximations are made. T h i s allows a careful s tudy of the efficiency and l imi ta -tions of the proposed method of clutter suppression. In order to create a B-scan, targets are placed into a time vector and con-volved wi th the impulse response of a G P R pulse. T h e target trace is a vector of dirac functions scaled by the predicted target amplitude. For example, the ground might be a dirac function wi th unit ampli tude occurring at lOOOps. Due to the F W H M of lOOps, the sampling rate is lOps, and an index of 100 would contains the target. 4.1.1 C lu t t e r M o d e l T h e clutter model used for simulation is an approximation of clutter types that exist in a real G P R scan. Figure 4.2 shows a real G P R survey of soil w i th no targets. T h e synthetic clutter model reflects the three highlighted events of 4.2. T h e clutter is scaled to unit ampli tude for al l events. Each event is s imilar to the physical events present in Figure 4.2. The specific t ime locations of clutter 53 Figure 4.2: Th i s G P R B-scan contains examples of clutter without targets. A n -notation A shows the air-ground event, the clutter in oval B is a highly variable ampli tude event and C shows a clutter event that does not exist for the duration of the scan. 51 are not cr i t ical and arbi t rar i ly defined as: • Synthetic Event, EA' Ground , at lOOOps wi th constant ampli tude • Synthetic Event, E ^ : C lu t t e r at 1250ps wi th varying ampli tude • Synthetic Event, Ec'- C lu t te r at 1500ps which does not persist T h e ampli tude of each event is defined by a function, before it is convolved wi th the synthetic pulse, EA=1 (4.1) • ' EB =cos(2vr0.02x), 1 > x < 50 (4.2) 0, 1 > x < 24 EC = < e (x-25 ) / e 4 ) 25 > X < 29 ' (4-3) 1, 30 > x < 50 4.1.2 Target M o d e l T h e target model is generated by calculat ing a spatial hyperbola attenuated by an exponential attenuation function dependent on the distance from the antenna to the target center. T h e spatial hyperbola is calculated by creating an array of t ime shifts, t(x) = -^Z2 + (x- xt)2. for x = 0 , 0 .01 , . . . , 0 . 25 cm. (4.4) v T h e velocity is defined in (2.7), the target center, xt, is 25 cm and the target depth, Z, is the difference between target depth and ground, which is 5 cm. The target t ime is converted to a M A T L A B vector index by scaling (4.4) w i th the t ime sampling rate, lOps and rounding to the nearest integer, i d x = round ( ^ 1 0 ^) . (4.5) V l 0 x l 0 1 2 y V ; 55 (a) Point Target (b) Large Target Figure 4.3: Example of synthetic targets used in the experiments. In M A T L A B , the image is indexed by, i d x and the attenuation is calculated by (2.5) and assigned to the image array, 2 e - 2 a d B(iebc) = - - ^ - . (4.6) T h e attenuation coefficient, a was chosen to represent sandy soil, which has a low attenuation of around l O d B / m . T h e distance from the antenna to the target is found by d = \J Z2 + (x — xt)2 and the target gain function, p-2aZ G = ~ ^ n (4-7) normalizes the peak of the the hyperbola to unit amplitude. T h e normalizat ion puts the clutter, target and ground on similar ampli tude scales, which is often the case. Th i s represents a reasonable model for a point target. A n addit ional target is used, which has larger dimensions. T h i s w i l l result in a horizontal component that facilitates testing of the average trace subtraction method. It is calculated i n the same manner, except the center of the hyperbola is extended horizontally, as presented in Figure 4.3. 56 4.2 Proposed Method of Clutter Suppression A new method is proposed for suppressing horizontal clutter. T h e foundation is that the existing prevalent method, average trace subtraction ( A T S ) is often times sufficient for highlighting hyperbolas in G P R data for reasons of computat ional complexity. Therefore, new methods must be reasonably fast or especially successful in suppressing clutter. T h e new wrapped curvelet transform, which is released in CurveLab [10] is sufficiently fast to consider as an application clutter suppression. Addi t ional ly , because the curvelet domain is efficient at representing G P R traces, thresholding should prove to be competing method to A T S . The act of thresholding "mutes" an event, which meets the cri ter ia of the threshold. Unl ike subtraction, thresholding does not add energy, therefore eliminat-ing the possibil i ty of constructive interference. If an appropriate threshold can be calculated that represents the clutter contained in the B-scan then a simple scaled threshold can be applied in the curvelet domain. Th i s w i l l suppress events that are clutter like and ignore events that do not match the estimated model . T h e clutter model becomes the key to successful implementation of the proposed method. T h e clutter model for A T S is the average trace, which has two major as-sumptions. T h e first is that the resulting average ampli tude is representative of al l traces. T h i s is usually not possible for a G P R system. A m p l i t u d e varies as a result of posit ional noise and subsurface material changes. T h e second assumption is that the phase remains the same. Th i s is usually true, however, not al l horizontal clutter events are present over the window of a-B-scan." If a clutter event only exists in a port ion of a B-scan the resulting subtraction w i l l actually place an event where none exists. Th i s would be a similar case i f the clutter event changed phase 180 degrees resulting in constructive interference. 57 T h e proposed method overcomes the pit-falls of A T S by adding ampli tude and phase resilience, but retaining a simple noise model . If a quasi-accurate clutter model is generated, the phase and ampli tude resilience of the method comes from the parsimonious representation of clutter. T h i s enables clutter to be muted w i t h a high probabil i ty that the target is preserved. T h e amount of mut ing comes from a shr inking parameter, A, which must be empirical ly discovered, or controlled by an expert user to determine the tradeoff between clutter removal and target signal preservation. T h e clutter model must be quasi-accurate, which is the same assumption that A T S uses. In the average trace case, the clutter is assumed to be stationary in phase and ampli tude over the spatial dimension of the B-scan window. A n alternative clutter model is dubbed the edge model. T h e edge model uses the average of only two traces located at the edges of the B-scan. If there are M traces in a B-scan, which are column vectors, then the 1-D edge noise model is To form the threshold operator, the noise model is first transformed into an image. T h e noise model image has dimensions identical to the B-scan under study, for t ime and M for space. For each trace in the 2-D noise model image, Q, is the calculated edge model trace mul t ip l ied by a length-M vector of ones, There are several reasons that make a noise model calculated from the edge traces a reasonable attempt to model horizontal clutter. T h e first is that the clutter is highly likely to be present in the edge of a B-scan, while the target is not. T h e second is that a large target w i l l have a significant horizontal energy that would hedge = ~ ( b l + b M ) • (4.8) (4.9) 58 be detected by an noise model calculated from the average trace. However, if the clutter is not present at the edges, the edge trace model does not work. T h e proposed method of clutter suppression is summarized: 1. Take the curvelet transform of the data B-scan. 2. Calculate the edge noise model and form a B-scan noise model . 3. Take the curvelet transform of the noise model B-scan. 4. Use the noise model in the curvelet domain as a threshold. 5. Scale the threshold and apply to the data curvelet coefficients. 6. Inverse transform the thresholded data. In notation, the set of curvelet basis vectors, gm produces a set of coefficients from the data, d and noise model, Q, dm=(d,gm) (4.10) Clm={n,gm). (4.11) T h e threshold operator formed from the noise model and the scaling parameter, A, {0 if \x\ < A | Q m | , (4.12) 1 if | i | > A | f i m | . A n estimate of the the clutter free data, d is made from thresholding wi th F, and reconstructing, d = }^F{dm)gm. (4.13) m T h i s method has the advantage of being quite easy to visualize the result. T h e noise model is "muted" from the data by zeroing only the coefficients that are 59 estimated to belong to noise. T h e phase of the data is not important , because the absolute value of the coefficients is considered when thresholding. T h i s is allows incoherent estimation, because the noise model only provides a spatial location and approximate ampli tude. Coherent estimation, such as subtraction, requires careful alignment of phase and ampli tude to prevent constructive interference. T h e am-pli tude resilience of thresholding comes from control of the A parameter. A large value increases the tolerance for ampli tude variation at the cost of damaging closely spaced events. 4.3 Method of Performance Evaluation To verify that the proposed method is successful, several means of verification are presented. Us ing the synthetic data, numerical results are presented by calculating T h e t radi t ional image processing definition of P S N R uses A = 255, which is the difference between the m a x i m u m and m i n i m u m values of a 255 intensity level image. To calculate A in this applicat ion, the difference of the max imum and m i n i m u m amplitudes are taken from the input data, Numer ica l metrics are useful for comparing the amount of energy removed from the noisy data, but visual quali ty is not well represented. Therefore, the numerical results are presented to verify that clutter suppression is successful in removing the clutter energy, but visual experiments are conducted to demonstrate the difference i n artifacts. the P S N R , (4.14) A = max|ct| — min |d | . (4.15) 60 T h e visual experiments allow a direct comparison of the proposed method and A T S . Image artifacts, like target energy damage, are best viewed by comparison of the input and output images of each algori thm. A residual image, d — d is also calculated, which emphasizes the artifacts either created, or not removed by each method. 61 Chapter 5 Experimental Results and Discussion The results presented in this chapter are from two separate methods. T h e first set of results show the sui tabi l i ty of a level-constant thresholded R D W T for removing white or colored Gaussian noise. T h e second set of results are generated by the application of clutter suppression on real and synthetic data. A discussion follows each experimental section. T h e experimental results for A-scan noise address these issues: • Performance when additive noise is Gaussian and white. • Performance when additive noise is Gaussian and colored. • Evalua t ion of scale parameter and input noise energy. • Compa t ib i l i t y of proposed method to trace stacking. In addit ion to the I D results, experimental results for B-scan clutter sup-pression w i l l compare two methods, the proposed curvelet thresholding method and 02 ' average traces subtraction. T h e experiments present: • V i s u a l performance for synthetic data. • Residual images for synthetic data. • Peak S N R results for synthetic data. • V i s u a l performance for real data. 5.1 Experimental Results of Noise Removal in A-scans T h e main results of this section are presented in Table 5.1'. T h i s experiment uses the best choice for a scale parameter determined a priori. T h e input noise is additive wi th a signal-to-noise ratio of 6dB. T h e signal is a simulated G P R pulse at lOOOps containing 2048 points. Table 5.1: Summary of the average SNRirnp of the level-constant thresholded R D W T and global thresholded D W T wi th respect to white noise and colored noise. T h e results show that the global thresholded D W T can not effectively remove colored noise, while the R D W T makes significant improvements. T h i s establishes the sui tabil i ty of level-constant R D W T thresholding as a method of noise removal. The next set of experiments evaluate the specific choice of parameters and the influence on this method. W h i t e Noise Colored Noise R D W T D W T 17.8 d B 13.9 d B 13.5 d B 0.1 d B 63 Or ig ina l wi thout noise Or ig ina l w i t h A G W N : 5 d B Or ig ina l w i t h A G C N : 6 d B R D W T and A G W N : 18 d B R D W T and A G C N : 14 d B D W T and A G W N : 14 d B D W T and A G C N : 1 d B Figure 5.1: T h e R D W T visual ly removes the noise and preserves the signal using a level-constant threshold. T h e D W T wi th a global threshold can effectively remove the white noise, but not colored noise. 64 5 . 1 . 1 P e r f o r m a n c e o f G a u s s i a n W h i t e N o i s e R e m o v a l T h i s experiment uses the normalized risk measure from equation (3.1) to evaluate the performance of the R D W T to the D W T and Wiener filtering. Th i s normalized risk, | | / — f\\2/Na2, is used i n the wavelet literature to show that thresholding exceeds linear methods of noise removal as the signal length, N, increases. T h e experiment sets a constant value for a and generates the input Gaussian white noise for a simulated G P R signal w i th increasing values of N. T h i s makes SNRin of the noise lower as N increases, because the energy of white noise is deter-mined by HiVcr21|. T h e normalized risk is similar to a signal-to-noise measurement, but measures remaining noise energy on a scale of 0 to 1. A normalized risk mea-surement of one means that no noise energy is removed and a zero means a l l noise is removed, while preserving the integrity of the signal. For each value of N, the experiment is conducted 20 times and the average risk is calculated according to (3.1). T h e D W T thresholding method is performed by estimating a global threshold from the noise variance of the finest scale coefficients. Wiener filtering is performed by the function included w i t h Wavelab [24]. T h e R D W T threshold is determined by the level-constant method using an estimate of each level w i th J — 8. The results show a decay for D W T thresholding which exceeds Wiener filter-ing, confirming the results presented in wavelet literature [48]. T h e R D W T removes more noise ini t ial ly, but decays slightly slower then the D W T . T h i s shows that the R D W T is more effective then the D W T at smal l values of when using a level-constant threshold. 65 Risk decay of signal estimation in A G W N 0.08 0.07 0.06 CM 0.05 b 0.04 0.03 0.02 0.01 0 \ -\ \ \ x \ \ \ \ \ \ \ . \ \ —. • • R D W T Wiener D W T 10 11 12 13 log 2 N points 14 15 Figure 5.2: Level-constant thresholding of the noisy R D W T coefficients is effective even when the signal length is short. Level-constant thresholds are less effective when the signal length is large. 1.4 1.2 1 £ 0.8 0.6 Risk Decay for A G C N with constant a R D W T Wiener D W T 11 12 13 N points (log scale) 14 15 Figure 5.3: G l o b a l thresholding of the D W T removes no measurable noise, thus producing no risk decay. Wiener filtering is marginal ly effective for low noise, but level-constant thresholding is successful at removing colored noise. 66 5.1.2 Performance of Gaussian Colored Noise Removal T h i s experimental result demonstrates the noise removal abi l i ty of each competing method in the presence of "in-band" colored Gaussian noise. The experiment uses the normalized risk metric so that a comparison can be made to the previous results of A G W N noise removal. T h e experiment is conducted in the same manner as the previous section using a synthetic G P R pulse wi th an increasing amount of points. T h e value of a remained constant, while N increased. T h i s choice reduces the amount of noise energy as the number of points increased. T h e noise is colored by the antenna model, hn, and generated randomly for each t r ia l . T h e results represent the average normalized risk for 20 trials. There are three outstanding results in Figure 5.2. F i r s t , the D W T method of global thresholding appears to by unable to reduce the amount of noise, regardless of the number of points available. T h e R D W T method, however, s t i l l produces a decay in risk as increases in the presences of G C N . Final ly , Wiener filtering produces marginal results, while the risk for R D W T thresholding remains effective, though five times larger wi th respect to A G W N performance. 5.1.3 Dependencies of Number of Decomposition Levels T h i s experiment evaluates the dependence of the scale parameter and SNRin using the SNRimp metric. For Figure 5.4, the 2048 point simulated G P R signal is used to test the performance of the R D W T level-constant thresholding method in the presence of white noise. T h e results show that as SNRin decreases to zero, decomposing to a higher number of levels yields sl ightly better results. In the case of min ima l noise, like 39dB, 67 Figure 5.4: T h e thick black line shows the m a x i m u m S N R improvement for a given number of levels of decomposition and additive G W N . A notable feature of the figure is the sharp drop of S N R improvement when the noise energy is small . the coefficients are not sparse at high levels and show a negative improvement at high values of J for A G W N . T h e same test is performed on A G C N using the simulated G P R signal. T h e results are quite different for the larger values of SNRin. T h e opt imal choice of J remains the same, but the metric SNRimp falls off as noise energy increases. T h e choice of in-band colored noise influences this drop, because the large noise events have similar characteristics to the G P R signal and thresholding becomes ambiguous. However, it is important to note that if clutter is modeled as in-band noise, then it would be suppressed by approximately l O d B . 68 Figure 5.5: A s the input S N R of the colored noise reaches OdB, the abi l i ty to dis-tinguish noise from signal diminishes, because noise events appear nearly identical to signal events, however the opt imal choice of J remains the same. 5 . 1 . 4 T h r e s h o l d i n g M u l t i p l e A - s c a n s b e f o r e S t a c k i n g The first experiment is designed to test the efficacy of operator order. Thresholding is a non-linear operator, which means the order of operations w i l l likely produce dif-ferent results. T h e ubiquitous nature of stacking suggests that the order of thresh-olding and stacking should be examined. Figure 5.6 demonstrates the difference between the three different cases: • Perform no wavelet thresholding, just stack. • Perform wavelet thresholding before stacking. • Perform wavelet thresholding after stacking. T h e noise is A G W N applied to a 2048 synthetic G P R signal. Wavelet denoising is performed using a level-constant R D W T threshold. T h e input noise is 9dB and 69 Thresholding performs better after stacking PQ T3 a > O S-i e OT O n l y Stacking Before Stacking ' — After Stacking 12 4 8 16 Number of traces Figure 5.6: Stacking scales the input noise a by 1/VT, however, this scaling is diminished by stacking before thresholding. Regardless of the order, wavelet noise removal s t i l l provides a significant improvement. J = 8, which corresponds to the opt imal scale parameter from Figure 5.4. T h e SNRimp metric is chosen to determine the amount of noise energy removed. T h e results indicate that level-constant thresholding should be performed after stacking. Th i s is an important and fortunate feature, because most G P R data is only available post-stack. T h e experiment does not test the assumption that each trace is aligned in time and phase, which is a basic assumption of stacking. However, if only one trace is available the results show that wavelet thresholding can s t i l l provide noise removal, whereas stacking cannot be performed. T h e second experiment confirms that the number of points, - N, of a trace and the noise coloring has l i t t le effect on the results presented in Figure 5.6. T h i s set of experiments tested four cases: A G W N wi th 2 and 32 input traces and A G C N wi th 2 and 32 traces. Each case was tested as the number of points in the input 70 25 S If fe. to fe. to lOcfc O Median • • Mean — Stacked 9 10 11 12 log2 iV points (a) White noise, T = 2 40 35 5 30 '•imp { 25 fe; to 20 15 10 O Median • • • Mean " ~ Stacked 9 10 11 12 log2 Ar points (b) White noise, T = 32 O Median • Mean — •' ' ~ Stacked o p. • > p.... o . 9 10 11 12 u log 2 / v points (d) Colored noise, T = 32 10 11 log2 N points (c) Colored noise, T = 2 Figure 5.7: Stacking before thresholding is the better. If stacking is performed in the wavelet domain, the median is an improvement over the mean when T = 32. signal varied from 256 to 32,768. In addi t ion to testing the order of stacking and thresholding, a modification to the stacking operator was made. T h e median opera-tor is used to perform a different estimate of a mean trace after thresholding in the wavelet domain. T h e median operator is more suitable to determining the mean i f a probabil i ty dis t r ibut ion function is skewed. T h e nature of thresholding suggests that the efficiency of median stacking should be tested. T h e results in Figure 5.7 show that as the number of points decreases the difference between pre-stacking and post-stacking shrinks. T h e median stacking operator offers no advantage for a small number of traces, but 3dB of improvement 71 in large stacks. M i n o r differences aside, noise coloring has no effect on the order of operations. Even in the presence of colored noise it is better to pre-stack G P R data. 5.2 Discussion of Noise Removal Results T h e results show that the level-constant thresholded R D W T is useful for denois-ing a synthetic G P R signal in Gaussian noise. Several aspects are tested and the scale parameter, J , is shown to have a dramatic effect on the abil i ty to perform noise removal. Noise coloring is shown to reduce the performance of level-constant thresholding, although performance s t i l l exceeds other methods tested. T h e order of operations is shown to be best performed by finding the stacked trace before wavelet denoising. The performance difference of the level-constant thresholded R D W T between white and colored noise is due to the abi l i ty to estimate the noise in the redundant wavelet domain. Us ing the redundant transform, a different estimate of a is obtained at each level. In Figure 5.3, we can see that this allows the threshold to adapt to the coloring of the noise. W i t h the understanding that the R D W T is the application of bandpass filters, the level-constant threshold is an estimate of the standard deviation for a part icular frequency band of the additive colored noise, which allows it to adapt to the coloring of the noise. T h i s choice of level-constant thresholding has l imitations, because i f the es-t imat ion is wrong, the corresponding threshold w i l l not be effective. T h e results from testing the scale parameter imp ly that when the noise is especially low, the threshold operator removes more signal energy then noise. The extreme case is in Figure 5.4, where the result of denoising makes the corresponding signal worse then the input . 72 W h e n comparing R D W T denoising to Wiener filtering and D W T thresh-olding in Figure 5.2, the normalized risk of the R D W T shows better performance in i t ia l ly then either method. A s N increases, the difference shrinks and D W T denois-ing overtakes the R D W T method. T h i s is at t r ibuted to the choice of level-constant thresholding of the R D W T . T h e median, which is used to estimate the threshold for each level, is no longer efficient for low noise and large N. T h e level constant threshold attempts to estimate the standard deviation of the noise at the low fre-quency levels where the sparsity property no longer holds. T h i s means the threshold predicts a large value for dj, which removes more signal then noise at that level. For white noise, it would be better to use a global threshold estimated from the finest scale coefficients of the R D W T for very low noise applications. W h e n applying thresholding techniques to G P R data, results show that thresholding should be performed after trace averaging. T h e reason is that an estimate of a for pre-stacked traces w i l l be larger then post-stack. Because the threshold is dependent on a, it w i l l be larger and remove large amounts of signal energy. W h e n the traces are averaged the value of a is reduced by a factor of 1/VT. T h e corresponding estimate of a is smaller and the threshold is smaller. T h i s only holds when the dis t r ibut ion is Gaussian, because the mean is an excellent point estimator. W h e n N is increased the difference between pre- and post-stacking is re-duced, especially when T = 32. T h i s is because the estimation of the threshold is improved and quali ty of noise removal increases. However, the difference can not overcome the tremendous gain of pre-stacking, which makes thresholding more ef-fective. However, this result does not take into account that a threshold could be estimated for each coefficient in the redundant domain. A level-constant threshold 73 assumes all coefficients at a part icular level have a similar standard deviation. W h e n only one trace is available this assumption is reasonable, but as the number of traces increases better estimations of am are possible. T h e benefits of level-constant thresholding seem to outweigh the negatives. Signals w i th small number of points are better estimated by the R D W T method then either D W T or Wiener filtering. Often for real G P R systems the number of points is under 2048, which increases the advantage of level-constant thresholding in the R D W T domain. T h i s implies the level-constant R D W T threshold would be a candidate for the removal of Gaussian noise, white or colored of G P R A-scans. 5.3 Simulation of Clutter Suppression in B-scans The following experimental results compare the baseline method of average trace subtraction ( A T S ) to a the proposed method of thresholding in the curvelet domain. T h e methods are compared on synthetic data to establish an understanding of the proposed method and l imitat ions. Rea l da ta is then tested to see the visual effects of the proposed method. T h e results of the synthetic tests are summarized by Table 5.2. Th i s experi-ment compares the effectiveness of A T S to the proposed method using both a point target and a large target. W h i l e a P S N R metric is not a visual indicator of quality, it shows the amount of noise energy removed. T h e results is that the proposed method removes more clutter energy then the A T S does. T h e remaining experiments focus on visual quality. T h e output of each algori thm and the residuals for the synthetic data are compared. Fol lowing the synthetic test, a test on real landmine data is presented. 74 Orig ina l Average Trace Subtract ion Proposed Point Target 17.6 21.8 40.2 Large Target 17.6 21.3 40.1 Table 5.2: T h e p e a k - S N R of A T S shows an improvement, but much less then the proposed method for both targets. 5.3.1 Synthetic B-scan Experiments T h e visual synthetic experiments are conducted to establish an understanding of the method. T h e first experiment explores the effect of choosing a noise model . The average trace subtract ion method assumes a noise model that is an image containing the average trace at every spatial point. T h i s experiment compares visual output of A T S and then performs the proposed method also using an average trace ( A T ) model. A s a comparison, the edge noise model is also presented and the results are shown for the proposed method. T h e experiment uses a large target, which penalizes the use of an average noise model . F igure 5.8(b) shows that the target energy has been damaged, and ad-di t ional artifacts have been introduced. T h e average noise model is not a good choice for the proposed method either. Target energy is also damaged by the thresholding, but no addi t ional artifacts are introduced. T h e next experiment in Figure 5.9 uses a point target and the same clutter model to compare A T S and the the residuals. T h e residuals emphasize that A T S cannot remove clutter which varies in ampli tude or phase. However, thresholding is not subjected to that l imi ta t ion and effectively removes al l clutter events. T h e residual of the proposed method contains only artifacts at the sharp edges of the clutter. T h e ground, having constant ampli tude, is easily removed by both methods. T h e final synthetic experiment compares the large target to both methods. 75 (a) Input Noisy Model (b) Target Model (c) Avg Trace (AT) Noise Model (d) Proposed Method with AT model (e) Edge Noise Model (f) Proposed Method with edge model Figure 5.8: In 5.8(a) the original data w i th clutter is presented, but the result should be just the target as in 5.8(b). T h e proposed method is calculated wi th an average noise model , 5.8(d) and edge noise model , 5.8(f). To understand the output of the proposed method, the noise model used for the threshold estimation is shown in 5.8(c) and 5.8(e). 76 (a) Target Model (b) Input Noisy Model (c) Proposed Method Result (d) Proposed Method Residual (e) Average Trace Subtraction Result (f) Average Trace Subtraction Residual Figure 5.9: In 5.9(d) the proposed method leaves v i r tua l ly no clutter energy and the target is s t i l l present. A T S leaves tremendous artifacts, but does manage to remove the ground clutter. 77 T h i s provides results for common subsurface targets, which are larger than point targets. W i t h these common targets, using A T S shows that significant artifacts are introduced, seen in Figure 5.10(e). T h e proposed method in Figure 5.10(c) demonstrates similar results for the large target as the point target. T h e method provides near complete removal of a l l clutter and very l i t t le change to the target energy. 5.3.2 Landmine B-scan Experiments T h i s experiment uses landmine data from the J R C landmine signatures database. T h e B-scan under study represents a physical location at 25cm from a survey that done on 50cm X 50cm area wi th 1cm resolution. T h e landmine is a P M N - t y p e anti-personnel device buried at 5cms in loamy soil. T h e data can be obtained from [11] and the system is described i n [59]. T h e results of clutter suppression are shown in Figure 5.11, which use the edge noise model shown in 5.11(b) to create a threshold operator. F igure 5.11(b) is the edge noise model calculated from the data domain. A n y event that is strongly horizontal and exists near the edge w i l l be shown in the noise model . T h e result from the proposed method w i l l be to remove any events that are approximated by the calculated noise model . The result is Figure 5.11(c), which indeed shows superior clutter suppression to the average trace subtraction method. A s demonstrated before, targets that are large suffer when A T S is used to suppress clutter. There are artifacts introduced around the target as a result of the average noise model . Th i s does not occur w i th the proposed method. T h e large target remains intact and becomes the prominent feature in the result. A l so note the ground clutter is not completely removed by either method. 78 (a) Target Model (b) Input Noisy Model (c) Proposed Method Result (d) Proposed Method Residual (e) Average Trace Subtraction Result (f) Average Trace Subtraction Residual Figure 5.10: T h e large target does not effect the abil i ty of the proposed method to remove clutter. Average trace subtraction removes a port ion of target energy, adds artifacts near the target and does not remove the amplitude varying clutter. 79 Source Noise M o d e l (a) Input Data (b) Edge Noise model Adap t ive Subtract ion Average Trace Subtract ion (c) Proposed Method (d) Average Trace Subtraction Figure 5.11: T h e clutter is suppressed i n the proposed method, whereas significant artifacts are introduced in the A T S method. The noise model in 5.11(b) shows what likely events w i l l be suppressed in the result of the proposed method. 80 5.4 Discussion of Clutter Suppression Results Developing an understanding of clutter is a difficult problem, therefore simple solu-tions are sought. The computational complexity and simplicity of ATS is low, which has led to the widespread use of ATS as an initial signal processing step. Also, ATS is effective at removing the air-ground event and other strong horizontal events. In a usual G P R B-scan, the horizontal clutter is the dominant visual feature and can often impede interpretation of the data. The largest horizontal event is the ground, but multiple reflections of the ground, antenna ringing and horizontal reflectors are also significant. The advantage of using the proposed method is that events only need to be roughly approximated. This is similar to the ATS assumption that clutter can be represented by an average of all traces. The drawback of simple noise models is that they are likely to fail if the broad assumptions are invalidated. In particular, the edge noise model used in the proposed method suffers if the clutter does not exist in the edge trace. The edge model also can not handle temporal changes in the clutter. In Figure 5.11, an example where slight variation in the temporal location of the ground begins to challenge the assumptions of the edge noise model. However, the proposed method is visually more successful at removing horizontal events and minimizing artifacts then the ATS method. The point target, in Figure 5.9 is well fit to ATS clutter suppression because the horizontal component is minimal. The advantage of the noise edge model and the proposed method are evident in Figure 5.8(f) and Figure 5.10(c). The proposed method does not require the use of an average trace noise model, therefore no target damage occurs in Figure 5.10(c). When the average trace model is used on the same example, in Figure 5.8(f), no additional artifacts are introduced. The ATS method 81 struggles wi th this, because the subtraction operator is b l ind to any changes i n data. T h e proposed method uses a threshold operator, which provides a quasi-adaptive method. Th i s prevents the proposed method from adding artifacts even when the estimation of the noise model is wrong. T h e drawback of using the proposed method is that the scaling parameter A must be expert ly controlled to balance the weighting of the noise model . T h e results presented here used a weight of 2.8, which corresponds to thresholding any coefficients that exceed 3 times the noise model coefficients. It was found in these experiments that the scaling parameter only required rough adjustments and was easy to understand and control. In the case of the landmine data, a much less aggressive setting was used. Sufficient clutter suppression was obtained by only scaling the noise model 1.8 times. 82 Chapter 6 Conclusions This thesis made two primary contributions, which furthered an understanding of thresholding in the wavelet and curvelet domain. The first contribution is an evalua-tion of level-thresholding the R D W T coefficients for noise removal in a G P R A-scan. The second evaluates the use of thresholding curvelets for clutter suppression in a G P R B-scan. There are four secondary contributions that stemmed from the results of the evaluation of removing noise using level-constant thresholding in the R D W T domain: • R D W T level-constant thresholding is suitable for use on the synthetic G P R signal presented here, and likely suitable for most G P R signals. • The scale parameter, which controls the number of decomposition levels is critical to the performance of a level-constant threshold method. • The amount of noise energy present in the noisy signal affects the performance and choice of the scale parameter. 83 • T h e order of operations, if stacking is performed, is to stack first and threshold afterwards. T h i s work has buil t on the success of previous work in wavelet threshold-ing. Previous achievements presented results on global thresholding the R D W T to remove white noise. These results have extended the use of R D W T to perform non-parametric signal regression using very li t t le a priori information. T h e threshold is determined by the data, but the scale parameter must be estimated by other means. T h i s work has also presented a reasonable physical model for a G P R A -scan. Us ing this model and other pract ical aspects of G P R technology the results can interpreted and applied to real data. Specifically, the knowledge of stacking before thresholding is useful, because the R D W T denoising can be applied direct ly to current G P R data. In the realm of clutter suppression, there are three secondary contributions that stemmed from evaluating the proposed method of thresholding an edge noise model in the curvelet domain: • Clu t te r suppression can be performed by thresholding in the curvelet domain. • A n approximation of clutter using a simple model, such as the edge noise model is sufficient when using thresholding. • A m p l i t u d e and phase resilience was demonstrated using thresholding. T h i s work uses the well established theory of thresholding in the wavelet domain to re-position the problem of clutter suppression as threshold estimation problem. Because of the phase and ampli tude resilience of thresholding effective clutter suppression can be achieved by only estimating an approximate model of 84 clutter. This enables users to arbitrarily design a clutter model that suits a specific need. 6.1 Future Work This work can be extended in three different ways. The first is to perform a physical evaluation of level-constant thresholding to accurate backscatter models. Second, complex clutter models could be created to enhance the effectiveness of clutter thresholding. Finally, combining both level-constant thresholding and clutter sup-pression in the curvelet domain to perform noise' removal and clutter suppression simultaneously. The most important extension is the combination of noise removal and clutter suppression. The suitability of methods has been evaluated, therefore the combi-nation seems quite natural and exciting: Because the curvelet domain functions in generally the same way as the redundant discrete wavelet domain, this extension should be quite natural. The scale parameters will have to be carefully chosen such that signal preservation and noise removal is balanced. The evaluation of physical data with respect to an accurate radar backscatter model is important to determine if the critical features of the signal are being pre-served when noise removal takes place. In addition, a physical analysis of random clutter would allow the possibility of extending the idea of colored Gaussian noise to clutter models. The development of complex clutter models could allow temporal clutter variance to be accounted for. The possibility of an interactive noise model could allow an operator to "select" features for removal. A clutter feature could be first identified by a user and then the angle and area of influence is then chosen to create 85 a noise model and suppress the clutter feature using thresholding. 86 Bibliography [1] G . B e y l k i n . O n the Representation of Operators in Bases of Compac t ly Sup-ported Wavelets. Siam Journal on Numerical Analysis, 29(6):1716-1740, 1992. [2] I lar ia Bott igl iero. 120 Million Landmines Deployed Worldwide: Fact or Fic-tion? Leo Cooper, Pen & Sword Books, Barnsley, South Yorkshire , England , 2000. [3] J . W . Brooks . Appl ica t ions of G P R Technology to Humani ta r ian Demining Operations in Cambod ia : Some Lessons Learned. In Internation Symposium on Technology and the Mine Problem, Monterey, C A , 1998. [4] John W . Brooks , L u c van K e m p e n , and Hichem Sahl i . P r i m a r y study in adap-tive clutter reduction and buried minelike target enhancement from G P R data. In Detection and Remdiation Technologies for Mines and Minelike Targets V, volume 4038 of Proceedings SPIE, pages 1183-1192, Orlando, F A U S A , 2000. S P I E . [5] C . Brusch in i and B . Gros. A Survey of Current Technology Research for the Detection of Landmines . Sustainable Humanitarian Demining: Trends, Tech-niques and Technologies, pages 172-187, 1998. [6] G . J . B u r t o n and G . P . Ohlke. Exp lo i t a t ion of Mi l l imete r Waves for Through-W a l l Surveillance dur ing M i l i t a r y Operations i n U r b a n Terrain. Technical re-port, Roya l M i l i t a r y College of Canada, M a y 24, 2000 2000. [7] E . J . Candes and D . L . Donoho, editors. Curvelets - A suprisingly Effective Nonadaptive Representation for Objects with Edges. Curves and Surfaces. V a n -derbilt Univers i ty Press, Nashvil le , T N , 1999. [8] E . J . Candes and D . L . Donoho. Curvelets and curvilinear integrals. Journal of Approximation Theory, 113(l):59-90, 2001. 87 [9] E . J . Candes and D . L . Donoho. New tight frames of curvelets and opt imal representations of objects w i t h piecewise C-2 singularities. Communications on Pure and Applied Mathematics, 57(2):219-266, 2004. [10] Emmanue l Candes, Laurent Demanet, D a v i d Donoho, and L e x i n g Y i n g . Curve-L a b 1.0, 2005. Avai lab le from: http://www.curvelet.org/. [11] Joint Research Centre. Anti-personnel Landmine Signatures Database. A v a i l -able from: h t t p : //apl-database. j r c . i t / . [12] M . Cherniakov and L . Donskoi . Frequency band selection of radars for buried object detection. Ieee Transactions on Geoscience and Remote Sens-ing, 37(2):838-845, 1999. [13] R . R . Coifman and D . L . Donoho. Translat ion Invariant De-Noising. In Wavelets and Statistics, volume Lecture Notes in Statistics, pages 125-150. Springer-Verlag, New York , 1995. [14] R . R . Coifman and M . V . Wickerhauser. Entropy-Based Algor i thms for Best Basis Selection. Ieee Transactions on Information Theory, 38(2):713-718, 1992. [15] Lawrence B . Conyers and Dean G o o d m a n . Ground-penetrating radar : an introduction for archaeologists. A l t a M i r a Press, Walnut Creek, C A , 1997. [16] J a m i n Cr i s t a l l , M o r i t z Beyreuther, and Fel ix Her rmann . Curvelet processing and imaging: 4D adaptive subtract ion. In CSEG National Convention. 2004. [17] D . J . Daniels, D . J . Gunton , and H . F . Scott . Introduction to Subsurface Radar . lee Proceedings-F Radar and Signal Processing, 135(4):278-320, 1988. [18] D . J . Daniels and Inst i tut ion of Elec t r ica l Engineers. Surface-penetrating radar. I E E radar, sonar, navigation and avionics series ; 6. Inst i tut ion of Elect r ica l Engineers, London , 1996. [19] D . J . Daniels and Inst i tut ion of Elec t r ica l Engineers. Ground penetrating radar. I E E radar, sonar, navigation, and avionics series ; 15. Inst i tut ion of Elect r ica l Engineers, London , 2nd edit ion, 2004. [20] Canadian Department of Na t iona l Defense. Canadian Forces Landmine Database, Dec 18, 2003. Avai lable from: http://ndmic-cidnm.forces.gc.ca/index.asp?lang=e. [21] D . L . Donoho. De-Nois ing by Soft-Thresholding. Ieee Transactions on Infor-mation Theory, 41(3):613-627, 1995. 88 [22] D . L . Donoho and M . E l a d . O p t i m a l l y sparse representation in general (nonorthogonal) dictionaries v i a 1(1) min imiza t ion . Proceedings of the National Academy of Sciences of the United States of America, 100(5):2197-2202, 2003. [23] D . L . Donoho and I. M . Johnstone. Ideal Spat ia l Adap ta t ion by Wavelet Shrink-age. Biometrika, 81(3):425-455, 1994. [24] D a v i d Donoho, M a r k Reynold Duncan , X i a o m i n g Huo , and Ofer L e v i . Wavelab 802. Avai lab le from: h t t p : / / w w w - s t a t . s t a n f o r d . e d u / w a v e l a b / . [25] M . ; Y u n x i n Zhao Gader, P^D.; Mys tkowsk i . Landmine detection wi th ground penetrating radar using hidden Markov models. Geoscience and Remote Sens-ing, IEEE Transactions on, 39(6):1231-1244, 2001. [26] P. D . Gader, B . N . Nelson, H . Fr igui , G . Vail let te , and J . M . Kel le r . Fuzzy logic detection of landmines wi th ground penetrating radar. Signal Processing, 80(6):1069-1084, 2000. [27] J . Gazdag. Wave-Equat ion Mig ra t ion w i th Phase-Shift Me thod . Geophysics, 43(7):1342-1351, 1978. [28] Rafael C . Gonzalez and Richard E . Woods. Digital image processing. Prentice H a l l , Uppe r Saddle River , N . J . , 2nd edition, 2002. [29] R . Gr ibonva l and M . Nielsen. O n the strong uniqueness of highly sparse repre-sentations from redundant dictionaries. Independent Component Analysis and Blind Signal Separation, 3195:201-208, 2004. [30] A . Grossmann and J . Mor le t . Decomposi t ion of H a r d y Functions into Square Integrable Wavelets of Constant Shape. Siam Journal on Mathematical Anal-ysis, 15(4):723-736, 1984. [31] Dig i t a l Signal Processing Group . Rice Wavelet Toolbox v2.4, 2001. Avai lab le from: h t t p : / / w w w - d s p . r i c e . e d u / s o f t w a r e / r w t . s h t m l . [32] Frdr ic Guerne, Ber t rand Gros, M a r c Schreiber, and Jean-Daniel N icoud . G P R M i n e Sensors for D a t a Acquis i t ion in the F i e ld . In International Workshop on Sustainable Humanitarian Demining, Zagreb, Croa t ia , 1997. [33] A . G u i t t o n and D . J . Verschuur. Adap t ive subtract ion of multiples using the L - l - n o r m . Geophysical Prospecting, 52(l) :27-38, 2004. [34] Fel ix J . Her rmann. Curvelet imaging and processing: an overview. In CSEG National Convention. 2004. 89 [35] F . J . Herrmann and P.P. Moghaddam. Sparseness- and continuity-constrained seismic imaging with Curvelet frames. Submitted for publication, 2005. [36] F . J . Herrmann and D. J . Verschuur. 3-D Curvelet Domain multiple elimination with sparseness constraints. Submitted for publication, 2005. [37] K . C . Ho, P.D. Gader, and J .N . Wilson. Improving landmine detection using frequency domain features from ground penetrating radar. In Geoscience and Remote Sensing Symposium. 2004. IGARSS '04. Proceedings. 2004 IEEE In-ternational, volume 3, pages 1617-1620 vol.3, 2004. [38] Tai-Chiu Hsung, D.P. -K. Lun, and Wan-Chi Siu. Denoising by singularity detection. Signal Processing, IEEE Transactions on [see also Acoustics, Speech, and Signal Processing, IEEE Transactions on], 47(11):3139-3144, 1999. [39] Sensors & Software Inc. E K K O - f o r - D V L pulseEKKO 100 Manual, 2001. [40] A . Instanes, I. Lonne, and K . Sandaker. Location of avalanche victims with ground-penetrating radar. Cold Regions Science and Technology, 38(1):55-61, 2004. [41] J . D . Irving and R. J . Knight. Removal of wavelet dispersion from ground-penetrating radar data. Geophysics, 68(3):960-970, 2003. [42] C. Jaedicke. Snow mass quantification and avalanche victim search by ground penetrating radar. Surveys in Geophysics, 24(5-6) :431-445, 2003. [43] J . van der Kruk, C. P. A . Wapenaar, J . T . Fokkema, and P . ' M . van den Berg. Three-dimensional imaging of multicomponent ground-penetrating radar data. Geophysics, 68(4):1241-1254, 2003. [44] M . Lang, H . Guo, J . Odegard, C. Burrus, and R. Wells. Nonlinear Processing of a Shift Invariant D W T for Noise Reduction. In SPIE Symp. OE/Aerospace Sensing and Dual Use Photonics, Algorithm for Synthetic Aperture Radar Im-age, Orlando, F L , 1995. SPIE. [45] M . Lang, H . Guo, J . E . Odegard, C. S. Burrus, and R. O. Wells. Noise reduc-tion using an undecimated discrete wavelet transform. Ieee Signal Processing Letters, 3(1):10-12, 1996. [46] K . H . Lee and G . Q. Xie. A New Approach to Imaging with Low-Frequency Electromagnetic-Fields. Geophysics, 58(6):780-796, 1993. 90 [47] S. M a l l a t and S. Zhong. Character izat ion of. signals from multiscale edges. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 14(7):710-732, 1992. [48] S. G . M a l l a t . A wavelet tour of signal processing. Academic Press, San Diego, 2nd edition, 1999. [49] S. G . M a l l a t and Z . F . Zhang. M a t c h i n g Pursui t s w i th Time-Frequency D i c -tionaries. Ieee Transactions on Signal Processing, 41(12):3397-3415, 1993. [50] A d a m Marcus . Inventions to Lif t the Scourge of Landmines, M a y 6, 1998 1998. Avai lab le from: http://www.csmon i to r . com/du rab l e /1998 /05 /06 / f -p l4 s l . h tm . [51] Rohan M a x w e l l . Cambod ia : A Coun t ry Profile. Journal of Mine Action, 5(1), 2001. [52] M . R . M c C l u r e and L . C a r i n . Ma tch ing pursuits w i t h a wave-based dict ionary. Ieee Transactions on Signal Processing, 45(12) :2912-2927, 1997. [53] I. L . Mor row and P . van Genderen. Effective imaging of buried dielectric objects. Ieee Transactions on Geoscience and Remote Sensing, 40(4):943-949, 2002. [54] G . P . Nason and B . W . Silverman. T h e Stat ionary Wavelet Transform and some Stat is t ical Appl ica t ions . In A Antoniadis and G . Oppenheim, editors, Wavelets and Statistics, volume Lecture Notes in Statistics, pages 281-299. Springer-Verlag, New York , 1995. [55] L . Nuzzo and T . Quar ta . Improvement in G P R coherent noise attenuation using tau-p and wavelet transforms. Geophysics, 69(3):789-802, 2004. [56] J . C . Pesquet, H . K r i m , and H . Carfantan. Time-invariant or thonormal wavelet representations. Ieee Transactions on Signal Processing, 44(8):1964-1970, 1996. [57] Aleksandra P izur i ca . Image Denoising Using Wavelets and Spatial Context Modeling. P h . d , Univer is ty of Gent, 2002. [58] Jane M a r y Anastas ia Rea . Ground penetrating radar applications in hydrology. Universi ty of B r i t i s h Columbia , Vancouver, 1997. [59] B . Scheers, M . Acheroy, and A . Vander Vorst . T ime-domain simulat ion and characterisation of T E M horns using a normalised impulse response. lee Proceedings-Microwaves Antennas and Propagation, 147(6) :463-468, 2000. 91 [60] B a r t Scheers. Ultra- Wideband Ground Penetrating Radar, With Application to the Detection of Anti-Personnel Landmines. P h . d , Roya l M i l i t a r y Academy, 2001. [61] U S Department of State. Hidden Ki l l e r s : T h e Globa l Landmine Cris is . Depar-ment of State Publ ica t ion 10575, Bureau of Po l i t i c a l -Mi l i t a ry Affairs, Office of Humani ta r i an Demin ing Progams, September 1998 1998. [62] James D . Taylor . Ultra-wideband radar technology. C R C Press, B o c a Ra ton , F L , 2001. [63] S. Tjora , E . Eide , and L . L u n d h e i m . Evalua t ion of methods for ground bounce removal in G P R ut i l i ty mapping. In Ground Penetrating Radar. 2004- GPR 2004- Proceedings of the Tenth International Conference on, volume 1, pages 379-382, 2004. [64] J . A . Tropp. Greed is good: Algor i thmic results for sparse approximation. Ieee Transactions on Information Theory, 50(10):2231-2242, 2004. [65] T . J . U l rych , M . D . Sacchi, and J . M . G r a u l . Signal and noise separation: A r t and science. Geophysics, 64(5):1648-1656, 1999. [66] T . J . U l r y c h , D . R . Velis, A . D . Woodbury, and R . D . Sacchi. L-moments and C -moments. Stochastic Environmental Research and Risk Assessment, 14(1):50-68, 2000. [67] L . V a n K e m p e n , H . Sahl i , E . Nyssen, and J . Cornells . Signal processing and pattern recognition methods for radar A P mine detection and identification. In Detection of Abandoned Land Mines, 1998. Second International Conference on the (IEE Conf. Publ. No. 458), pages 81-85, 1998. [68] R . W u , A . Clement, J . L i , E . G . Larsson, M . Bradley, J . Habersat, and G . Maksymonko . Adapt ive ground bounce removal. Electronics Letters, 37(20):1250-1252, 2001. [69] X u X i a o y i n , E . L . Mi l l e r , C M . Rappapor t , and G . D . Sower. Stat is t ical method to detect subsurface objects using array ground-penetrating radar data. Geo-science and Remote Sensing, IEEE Transactions on, 40(4):963-976, 2002. [70] Y . X . Zhao, P . Gader, P . Chen , and Y . Zhang. Tra in ing D H M M s of mine and clutter to minimize landmine detection errors. Ieee Transactions on Geoscience and Remote Sensing, 41(5):1016-1024, 2003. 92 M . S. Zhdanov, P. Traynin, and J . R. Booker. Underground imaging by frequency-domain electromagnetic migration. Geophysics, 61(3):666-682, 1996. 93 Appendix A Image Plot t ing Function Due to the subjectivi ty of plot t ing intensity images, a detailed description of the plot t ing function which was used to generate B-scan plots is described. T h e im-portant features are the color scaling, image scaling and image translation. T h e function used to re-scale and translate is called r e l _ i m , shortened from, relative image. T h e purpose was to allow results to be plotted at the same intensity scale as the input . M A T L A B includes a image plot t ing function called image, which interprets an entry in a mat r ix as an index to a colormap: T h e colormap used in this thesis for plot t ing B-scans is white at the middle index of the colormap mat r ix and l inearly scales in shades of gray to black at the beginning and end of the array using the code: function [cmap] = seisgray; i f nargin <1, ncolor = 32; end tmpl = [[0:(ncolor-1)]/(ncolor-1) ; ... [0:(ncolor-1)]/(ncolor-1) ; ... [0:(ncolor-1)]/(ncolor-1)]'; 94 tmp2 = [[(ncolor-l):-1:0]/(ncolor-1) ; ... [(ncolor-1):-l:0]/31 ; ... [(ncolor-1):-l:0]/3i]'; cmap = [tmpl ; tmp2]; T h i s choice of color mapping makes no dist inct ion between negative and positive events i f the image array is translated such that the zero indexes the middle index of the colormap. T h i s is the purpose of the r e l _ i m function, which scales and translates an input image array relative to a second input: 7, Image something relative to another image '/. Usage: '/. rel_im(source, image) function im = rel_im(source, im) scale = .5/max(max(abs(source))); im = im.*scale; image(round(im.*64+32)) T h e purpose is to compare images at a common scale relative to a part icular source image. T h i s thesis uses the input image as the base scale and all other images are scaled accordingly. If the input to the function is an array, it takes the m a x i m u m value and normalizes the second argument to this value. Otherwise, a simple scalar can be used as the first argument to choose an arbi t rary scale. Note there are only 16 levels of gray to choose from, between pr int ing and the dynamic range of the color scale, some smal l ampli tude information w i l l be lost. The choice to plot the images this way was carefully made i n order to convey useful information that is not biased. It is a difficult topic to present results fairly, but every effort was made to remove this bias. 95
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Threshold estimation using wavelets and curvelets on...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Threshold estimation using wavelets and curvelets on ground penetrating radar data for noise and clutter… Harrison, Dustin 2005
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Threshold estimation using wavelets and curvelets on ground penetrating radar data for noise and clutter suppresion |
Creator |
Harrison, Dustin |
Publisher | University of British Columbia |
Date Issued | 2005 |
Description | This thesis provides an evaluation of the Redundant Discrete Wavelet Transform with application to the removal of additive white or colored Gaussian noise on a synthetic GPR signal. Special attention is given to the parameter that controls the number of decomposition levels. Evaluation is performed using a level-dependent threshold to estimate and remove noise. Results are presented using noisy synthetic Ground Penetrating Radar pulses to compare Wiener filtering and thresholding the Redundant and Non-redundant Discrete Wavelet transform. Additional results are presented on the effects of choosing a number of decomposition levels using signal-to-noise ratio measurements, which suggest the importance of choosing this parameter. Recommendations are made and supported which determine the order of thresholding before or after the practice of trace averaging. Using GPR images, an application of a novel 2D threshold model in the newly discovered curvelet domain is compared to average trace subtraction. Promising results are presented on both synthetic and actual landmine data, which shows thresholding as a viable method of clutter suppression. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2009-12-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065408 |
URI | http://hdl.handle.net/2429/16564 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2005-0467.pdf [ 6.86MB ]
- Metadata
- JSON: 831-1.0065408.json
- JSON-LD: 831-1.0065408-ld.json
- RDF/XML (Pretty): 831-1.0065408-rdf.xml
- RDF/JSON: 831-1.0065408-rdf.json
- Turtle: 831-1.0065408-turtle.txt
- N-Triples: 831-1.0065408-rdf-ntriples.txt
- Original Record: 831-1.0065408-source.json
- Full Text
- 831-1.0065408-fulltext.txt
- Citation
- 831-1.0065408.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065408/manifest