Threshold Estimation Using Wavelets and Curvelets on Ground Penetrating Radar Data for Noise and Clutter Suppression by Dustin Harrison B . S . , S o u t h D a k o t a School of M i n e s and Technology, 2000 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR T H E DEGREE OF Master of Applied Science in T H E F A C U L T Y OF G R A D U A T E STUDIES (Electrical and C o m p u t e r Engineering) The University of British Columbia M a y 2005 Â© D u s t i n H a r r i s o n , 2005 Abstract T h i s thesis provides an evaluation of the R e d u n d a n t Discrete Wavelet Transform w i t h a p p l i c a t i o n to the removal of additive white or colored G a u s s i a n noise on a synthetic G P R signal. Special attention is given to the parameter that controls the number of decomposition levels. E v a l u a t i o n is performed using a level-dependent threshold to estimate and remove noise. Results are presented using noisy synthetic G r o u n d P e n e t r a t i n g R a d a r pulses to compare W i e n e r filtering and thresholding the R e d u n d a n t 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, w h i c h suggest the importance of choosing this parameter. Recommendations are made and supported w h i c h determine the order of thresholding before or after the practice of trace averaging. U s i n g G P R images, an application of a novel 2 D threshold m o d e l i n the n e w l y discovered curvelet d o m a i n is compared to average trace subtraction. Promising results are presented on b o t h synthetic and actual landmine data, w h i c h shows thresholding as a viable m e t h o d 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 2.1 7 Ground Penetrating Radar 7 2.1.1 Principles of Operation 9 2.1.2 Data formatting and visualization iii 12 2.1.3 2.2 2.3 2.4 Sources of Noise 13 Selected R e v i e w of G P R S i g n a l Processing 17 2.2.1 R e v i e w of A - s c a n Processing M e t h o d s 18 2.2.2 R e v i e w of B - s c a n Processing M e t h o d s 20 R e v i e w of Wavelet Transforms 22 2.3.1 Discrete Wavelet Transform 23 2.3.2 R e d u n d a n t Discrete Wavelet T r a n s f o r m 26 2.3.3 Curvelet Transform 28 Noise R e m o v a l using T h r e s h o l d E s t i m a t i o n s 31 3 Removing Noise from a G P R A-scan 3.1 D e t e r m i n i n g a Basis for A - S c a n A p p r o x i m a t i o n 37 3.2 D e t e r m i n i n g the N u m b e r of Levels of D e c o m p o s i t i o n 39 3.3 E s t i m a t i o n of a T h r e s h o l d from N o i s y D a t a 41 3.4 G P R A-Scan Model 44 3.4.1 Antenna Model 46 3.4.2 Noise M o d e l .' 3.5 4 35 ' Summary 49 Clutter Suppression in G P R B-Scans 4.1 47 50 B - s c a n M o d e l using Synthetic G P R P u l s e 53 4.1.1 Clutter Model 53 4.1.2 Target M o d e l 55 4.2 P r o p o s e d M e t h o d of C l u t t e r Suppression 57 4.3 M e t h o d of Performance E v a l u a t i o n 60 iv 5 E x p e r i m e n t a l R e s u l t s a n d Discussion 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 5.4 6 62 Discussion of Clutter Suppression Results Conclusions 6.1 78 81 83 Future Work 85 Bibliography Appendix A : 87 Image P l o t t i n g F u n c t i o n v 94 List of Tables 3.1 E n t r o p y measurement of different basis functions 38 5.1 Effect of noise coloring on SNR 63 5.2 P S N R of proposed m e t h o d and average trace s u b t r a c t i o n u s i n g the R D W T and D W T . . . imp vi 75 List of Figures 1.1 Example Landmines in situ 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 4.2 Example clutter 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 SNR 68 imp 3 52 '. of R D W T level-constant thresholding for A G W N vii 54 5.5 SNR 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 imp of R D W T level-constant thresholding for A G C N target 69 77 5.10 Comparison of residuals between proposed method and A T S for large target 79 5.11 Proposed method demonstrated on P M N landmine data. . . . . . . . vm 80 List of Abbreviations AGCN Additive Gaussian Colored Noise AGWN Additive Gaussian White Noise ATS Average Trace Subtraction DSP - D i g i t a l Signal Processing- DWT Discrete Wavelet Transform FWHM Full W i d t h Half Maximum GCN Gaussian Colored Noise GPR Ground Penetrating Radar GWN 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 â€¢ Dr. M . R . Ito provided me w i t h funding and academic advisement w h i c h facilitated m y completion. â€¢ D r . Felix H e r r m a n n teaches an i n s p i r i n g course on wavelets, w h i c h motivated this research and his great knowledge of the topic gave me the techniques and insight to finish. â€¢ T h e U n i v e r s i t y of B r i t i s h C o l u m b i a also saw enough potential i n m y determination and skills to provide a d d i t i o n a l funding for this project. â€¢ E m m a n u e l Candes, L a u r e n t Demanet, D a v i d D o n o h o , L e x i n g Y i n g deserve gracious thanks for m a k i n g C u r v e L a b available for research. â€¢ Joint Research Centre for h a v i n g the brilliant foresight to facilitate h u m a n i t a r i a n l a n d m i n e detection research by the creation of the J R C L a n d m i n e s 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. DUSTIN The University of British Columbia May 2005 x HARRISON T o m y wife, D o r a L i for her patience a n d love a n d m y sister, A l i s s a for showing me I can do it xi Chapter 1 Introduction T h e r e are millions of landmines w o r l d w i d e , w h i c h are difficult to find and remove. T h e concept of humanitarian demining, w h i c h is sometimes referred to as H u m D e m , is to provide simple, sustainable and affordable m e t h o d of detection and removal of landmines. T h e groups that implement these programs are often non-governmental organizations w h o d i d not place the mines and generally have l i t t l e information on mine type or placement locations. T h e most c o m m o n and h i g h l y visible aspect of d e m i n i n g are the activities that i n poor countries devastated by war. T h e objective of H u m D e m groups i n these countries are to develop programs w i t h i n 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. T h e sustainability objective of m a n y H u m D e m programs is critical to imple- m e n t i n g a m e t h o d of l a n d m i n e clearance that can be performed by locally trained d e m i n i n g groups. P r e v i o u s estimations state that 100 m i l l i o n or more landmines are currently laid i n the w o r l d , however that number has been reduced to 60 m i l l i o n or less [61]. T h e number of landmines, danger of detection a n d the removal procedures increase the cost and t i m e requirements to nearly impossible levels for poor 1 countries. T h e cost of detection and removal is also related to the amount of t i m e a l a n d m i n e is i n the g r o u n d . Delays i n l a n d m i n e clearance causes vegetation coverage, w h i c h increases the challenge to detect minefields and mine locations. D u e 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 m a k i n g simple, affordable a n d mature technology is necessary so that projects can continue even after external funding has ceased.' H u m D e m programs determine s u i t a b i l i t y of technology by the cost and status of development maturity. T h e most readily used a n d implemented technology is h u m a n protection devices, m e t a l detectors and trained dogs. O t h e r technologies that have been tested or evaluated w i t h respect to H u m D e m are: â€¢ Explosives detection using chemical sensors, biological sensors and nuclear quadrupole resonance â€¢ O p t i c a l sensors such as infrared, hyper-spectral, visible and laser ranging â€¢ Electromagnetic sensors i n c l u d i n g ultra-low frequency, microwave, S A R and ground penetrating radar â€¢ A c o u s t i c sensors such as water-jet echo [50], impulse and ultrasound B r u s c h i n i and G r o s discuss i n detail most of the above technologies i n [5]. In particular, the authors suggest that G P R is a m a t u r e technology a n d that it is "near ready" for use i n H u m D e m programs. T h e H u m D e m field is not receptive to new technologies, regardless of the promises of detection improvements [3]. F o r example, estimates for C a m b o d i a suggest that mine clearance is h a p p e n i n g at a rate of 15 square kilometers per year w i t h the possibility of d e m i n i n g all l a n d that is i m m e d i a t e l y i n 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]. Technology 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 G P R components, such as analog/digital converters, amplifiers and high speed switches contribute to the increasing interest to implementing GPR as a landmine detection device. Additionally, technology improvements in computer processors and specialized DSP chips allow additional realtime and field-based signal processing of GPR surveys. While the processing of G P R is maturing, the limited resolution of G P R 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 diffi3 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 I D 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 reflections 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 image 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 I D wavelet noise removal using thresholding is well established, the application 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 established. This research examines the most promising wavelet decomposition method for noise removal, which is the Redundant Discrete Wavelet Transform ( R D W T ) . 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 D S P 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 applied 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 2 D "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 construction 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 ( A G C N ) . 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 I D 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 technique 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 experiment 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 w i t h respect to ground penetrating radar. T h e concepts that are needed to perform this evaluation are a basic physical u n d e r s t a n d i n g of the G P R d a t a and how wavelets can be applied. T h e m o t i v a t i o n is clear: a d d i t i o n a l insight into G P R processing should serve to enhance the field of l a n d m i n e detection. 2.1 Ground Penetrating Radar R a d a r operates on the p r i n c i p l e that the range from the radar to the target can be determined by t i m i n g the r o u n d t r i p of a t r a n s m i t t e d pulse reflected by a target. T h e assumption of radar is that the transmission m e d i u m is conducive to passing the frequency content of the pulse. T h e earth is not a homogeneous m e d i u m , therefore anticipating the frequency content to pass t h r o u g h the m e d i a 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 d u r a t i o n of receiver gating. T h i s resulting t i m e - d o m a i n signal is called a trace, or A - s c a n . F o r G P R d a t a to be meaningful, a G P R survey usually attempts to s p a t i a l l y sample a grid of points o n the surface. T h e collection of traces from a G P R survey of several spatial locations forms a three-dimensional d a t a set of two spatial coordinates a n d one of t i m e called a GPR data cube. T h e m e t h o d of c o n d u c t i n g a G P R survey varies o n the type of G P R . A G P R m a y be bi-static, mono-static or neither. I n mono-static mode, the same antenna transmits a n d receives the G P R pulse. T h e b u r d e n falls o n the transceiver electronics to p r o p e r l y gate the t r a n s m i t a n d receive times a n d protect the receiver from the large i n i t i a l radiated pulse. I n bi-static mode, separate antennas are used for t r a n s m i t a n d receive, b u t their s p a t i a l relationship is held constant. T h e analog in seismic jargon for this m e t h o d of operation is called a common-offset gather. E a c h trace represents a constant angle between the transmitter a n d receiver. I n bi-static mode, antenna p o l a r i z a t i o n relative to t r a n s m i t a n d receive c a n improve target detectability [43]. A d d i t i o n a l l y , c o m m o n m i d - p o i n t gathers c a n be used to determine velocity models of the subsurface. I n this survey, the angle between transmit a n d receive antenna is increased after each pulse transmission. M a n y p r a c t i c a l applications of G P R are being developed. H i g h resolution G P R solutions are used to image the internal structure of roads a n d bridges for the appearance of cracks. M i l i t a r y a n d security applications have been suggested to image potential threats b e h i n d walls of caves a n d buildings [6]. Snow pack layers, internal ice structure a n d avalanche v i c t i m l o c a t i o n have also been demonstrated using G P R [40, 42]. G e o p h y s i c a l scientists make extensive use of G P R for the s t u d y 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 (2.1) 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 spectrum 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, (2-2) This pulse can be created by using a transistor in avalanche mode, with results in 9 1.5 FWHM -2 -1.5 -1 -0.5 0 1 0.5 1.5 2 F i g u r e 2.1: P i c t o r i a l representation of the F W H M parameter of a G a u s s i a n b u m p . an output that is an a p p r o x i m a t i o n of an exponential charge and discharge cycle. A t y p i c a l l y G P R specification is the F u l l W i d t h , H a l f M a x i m u m ( F W H M ) of a pulse. The F W H M is a c o m m o n measurement of the w i d t h of a gaussian pulse and is defined by FWHM = 2 C T v /21og 2 e (2.3) T h e w i d t h of a pulse indicates b a n d w i d t h , but usually the w o r k i n g b a n d w i d t h of a radar is l i m i t e d by the antenna. The response of an antenna to the i n p u t of (2.2) is to radiate a voltage wavefield w h i c h is the derivative w i t h respect to time, V'(t) = ~e- t 2 / 2 a \ (2.4) A n event is defined as the filtered form of (2.4) t h a t occurs at some point i n t i m e of a received G P R trace. T h i s event corresponds is the reflection of the t r a n s m i t t e d pulse from a subsurface feature a n d can be related to d e p t h if the velocity of the transmission m e d i a is k n o w n . I n the case of a metallic reflector the event is the second derivative of (2.2) because a metallic reflector acts as a second t r a n s m i t t i n g 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 w i t h the surface of the ground by impedance matching. T h e antenna also determines the frequency content of a radiated pulse, w h i c h dictates spectral efficiency thereby c o n t r o l l i n g depth of penetration and target visibility. H u m D e m researchers have developed an antenna that exhibits excellent G P R performance w i t h respect to spectral efficiency and ground c o u p l i n g [59]. D a t a from this l a b o r a t o r y G P R and antenna is used to test the a l g o r i t h m presented i n C h a p t e r 4. T h e characterization of an event i n a G P R trace is determined by the t i m i n g and a m p l i t u d e of the occurrence. D u e to spherical spreading of the t r a n s m i t t e d wavefront and attenuation of the subsurface, the a m p l i t u d e of the received signal for a mono-static configuration is shown by "[17] to be p r o p o r t i o n a l to f a* sec* 0 / o e\ 2adsec9 6 â€¢ [ Z - b ) T h e distance, d, is the line, i n meters, n o r m a l 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, w h i c h further attenuates the signal. In practice, some c o m m e r c i a l G P R systems a t t e m p t to compensate for this attenuation using an exponential gain applied d u r i n g or after d a t a collection [39]. W h i l e it is the prerogative of the user interpreting the d a t a to a p p l y a gain, it is i m p o r t a n t that noise removal be done beforehand to avoid non-linear changes to the d a t a statistics. T h e determination of t i m i n g i n a homogeneous m e d i u m is based on the horizontal location, x, of the antenna relative to the target a n d the depth, d, of the target. For a point target, or a target size that is near the resolution l i m i t of the radar, the r o u n d t r i p time is determined by twice the velocity of pulse propagation 11 times the line of sight distance to the target t= -Vx v 2 + d. (2.6) 2 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 become 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 D u e to the nature of a G P R pulse containing m u l t i p l e events, an A - s c a n is described as "wiggly" and has m a n y s m o o t h oscillations. T h i s descriptive of A-scans also refers to a m e t h o d of v i e w i n g B-scans called a wiggle plot, w h i c h is actual A scans vertically oriented and plotted close together. T h i s m e t h o d is used most often by h u m a n operators to make interpretations of G P R data. I n F i g u r e 2.2, a pixel intensity B - s c a n image is presented, where each pixel color is scaled based on the a m p l i t u d e i n the A - s c a n . Occasionally C-scans are used i n image processing. A C-scan is a "depth" slice, w h i c h takes the A - s c a n values over a l l spatial locations to form a intensity image. T h e C-scan contains b o t h spatial dimensions for axes, w h i c h are the downtrack 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 all B-scans. E a c h B - s c a n is, i n t u r n , a subset of A-scans. 2.1.3 Sources of Noise T h e r e are m a n y noise sources present i n a G P R device. D u r i n g the creation of the i n i t i a l source waveform, noise may be added by interference from external sources or poor high-frequency digital design of electronics. After a pulse is transmitted, crosscoupling from the antennas may generate spurious, large a m p l i t u d e events prior to any reflection events. In a d d i t i o n , u p o n reception of a G P R signal, clock jitter, analog-to-digital quantization noise and amplifier noise determine a noise floor and the operating d y n a m i c range. A t the resolution of G P R , the m a j o r i t y of signal energy comes from large impedance changes i n the ground. Changes i n soil water content w i l l d r a m a t i c a l l y 13 A-Scan at x=25cm and y=25cm B-Scan at x=25cm Amplitude y (cm) F i g u r e 2.2: Selected E x a m p l e s of different scan types. E a c h contains different perspectives of a G P R d a t a cube. 1 4 attenuate the signal a m p l i t u d e as it propagates into the g r o u n d . A d d i t i o n a l l y , if the subsurface contains large amounts of near wavelength objects, like glacial t i l l , depth of penetration and target localization 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 i m p o r t a n t to consider that the antenna is a physical object that is never precisely positioned. I n the lateral sense, the antenna position is critical if any transforms are taken on the d a t a i n the X - Y dimension. B o t h the wavelet and Fourier transforms expect d a t a to be sampled on a fixed g r i d . I n order to correct this noise, antenna position must be recorded and the d a t a interpolated onto a fixed g r i d . I n the application of l a n d m i n e detection this becomes a p r o b l e m w i t h h a n d h e l d detectors, w h i c h can be s w u n g in arcs or irregular lines [32]. M o s t laboratory solutions use robots or antenna jigs to stabilize and accurately determine spatial l o c a t i o n . In a d d i t i o n to lateral movement, it is i m p o r t a n t to consider vertical and angular antenna p o s i t i o n . B o t h angular a n d vertical p o s i t i o n i n g affect the c o u p l i n g of energy into the g r o u n d . A s the c o u p l i n g changes, the a m p l i t u d e of the received signal varies. T h e result is that algorithms must adapt to an u n k n o w n variance in a m p l i t u d e . A d d i t i o n a l l y , vertical position uncertainty creates t e m p o r a l variance of the received signal. U s u a l l y this noise is m i n i m i z e d by allowing the antenna to travel flush to the ground, but as g r o u n d changes occur, the h y p e r b o l i c assumption of a target response can be invalidated. Often this issue is lessened by the d u r a t i o n of scans, w h i c h are taken over contours that can be considered flat. 15 System Noise T h e most well understood system noise is t h e r m a l noise. T h e r m a l noise can be m o d eled by additive G a u s s i a n noise. T h e 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 m u l t i p l e traces from the same spatial l o c a t i o n . C o m mercial G P R systems implement this i n hardware and w i l l provide the user w i t h the ability to average up to 32 traces. O t h e r system noise, w h i c h is less understood and might not be well modeled by G a u s s i a n noise comes from electronic design and susceptibility. F o r example, poor circuit b o a r d design of a transmitter can a d d 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 limits the d y n a m i c range of the system. T h e combination of all sources of noise is very likely G a u s s i a n , but not necessarily white. C l u t t e r noise C l u t t e r is defined by the user and s t r i c t l y speaking is not noise. C l u t t e r is the cumulative effects of reflections from actual boundaries that are i l l u m i n a t e d by the G P R , but impede w i t h analysis of the p r i m a r y targets under study. T h e most c o m m o n form of clutter i n any G P R a p p l i c a t i o n is the ground. T h e occurrence of a reflection from the ground offers very little information to a G P R image and may interfere w i t h near-surface objects. L a n d m i n e s are an example of a near surface target that can be impeded by the presence of an air-ground event. O t h e r forms of clutter, w i t h respect to landmines, are targets that have properties very close to the desired target and can result i n a false positive determination of a p r i m a r y target. These false positives can be metal fragments, stones and other 16 debris that produce events similar to a l a n d m i n e event i n either a m p l i t u d e , frequency content and shape or a c o m b i n a t i o n of the three. Generally, this form of clutter is filtered using d i s c r i m i n a t i o n algorithms. A d i s c r i m i n a t i o n a l g o r i t h m w i l l use some unique feature of the desired target to suppress all other possible targets that do not contain any such features. A n o t h e r form of clutter, sometimes called background noise, appears as the p o r t i o n of a signal, w h i c h does not necessarily represent a physical object, b u t appears as a slowly v a r y i n g horizontal event i n a B - s c a n . T h e most c o m m o n form is m u l t i p l e reflection that occurs between the ground and antenna, often exacerbated by antenna r i n g i n g . Subsurface planar reflectors are sometimes grouped w i t h background noise, because of the t y p i c a l l y large h o r i z o n t a l component. T h e terms background noise and clutter are used interchangeably i n the literature and are problem specific. C l u t t e r noise models are often considered additive, such that the subsurface can be represented by a t r a i n of dirac functions, each of w h i c h represent either target or clutter t e m p o r a l locations. T h e G P R signal is then convolved w i t h the dirac t r a i n to o b t a i n a simulated received trace. C l u t t e r is n o r m a l l y the l i m i t i n g noise i n G P R systems [18]. Besides target depth, clutter is the l i m i t i n g 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 d a t a is processed, whether b y h u m a n or machine. In a simple case of u t i l i t y finding, it is usually enough for a h u m a n to interpret u t i l i t y location based on the existence of a h y p e r b o l a i n an unprocessed B - s c a n . In the case for h u m a n i t a r i a n demining, it is usually necessary to perform a d d i t i o n a l machine-based processing 17 and d i s c r i m i n a t i o n to perform accurate and reliable l a n d m i n e detection. G P R processing can be d i v i d e d into three categories, restoration, transfor- m a t i o n and d i s c r i m i n a t i o n . R e s t o r a t i o n is concerned w i t h the removal of noise, correction of system error or correction of the sensor degradation. T r a n s f o r m a t i o n is used as a functional t o o l to perform signal processing, b u t i n the case of m i g r a t i o n is also an output. D i s c r i m i n a t i o n is a key aspect to l a n d m i n e detection. D i s c r i m i nation uses a feature, hopefully unique, to separate noise and clutter from targets. For d i s c r i m i n a t i o n clutter often refers to a n y t h i n g that is not a desired target. T h i s definition varies depending on the features used to discriminate. A t y p i c a l signal processing chain for l a n d m i n e detection may include several restorative techniques, possibly a transformation and always a d i s c r i m i n a t i o n step. T h e i m p o r t a n c e of developing an excellent understanding and usefulness of restoration tools is paramount to successful l a n d m i n e detection. 2.2.1 Review of A-scan Processing Methods W h e n processing a G P R A - s c a n , most methods focus on signal restoration. Restorative methods include, DC-offset removal, G a u s s i a n noise filtering and feature enhancement. W i t h the exception of DC-offset removal, w h i c h is sometimes performed by the equipment, most restorative techniques are chosen based on the r e m a i n i n g 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 visually compared. For A-scans, these methods are frequency filtering, averaging and gain adjustment. Frequency filtering consists of a p p l y i n g 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 d o m i n a t e d by the time-variant low pass effects, w h i c h reduces the high-frequency content of late time events. H i g h pass filtering is performed w h e n it is necessary to dewow an A - s c a n . A "wow" signal is the low frequency effects of t r a n s m i t t e r a n d receiver c o u p l i n g , w h i c h produces a slow decaying signal into the A - s c a n [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 o n A-scans is stacking. S t a c k i n g , or i n G P R jargon, trace averaging, is the average of repetitive scans from the same s p a t i a l l o c a t i o n . C o m m e r c i a l G P R systems often include this form of noise reduction as a p r i m a r y parameter choice when performing a G P R survey. T h e purpose of stacking is to remove system noise. O t h e r methods include W i e n e r filtering, w h i c h assumes the noise m o d e l is well understood a n d attempts to perform an inverse noise filter i n the frequency d o m a i n . T h i s m e t h o d is also called m i n i m u m square error filtering, because it seeks to m i n i m i z e the error measure, E[/ â€” f] . 2 response, H(u) T h i s filter seeks to invert the degradation i n the frequency d o m a i n using [28] F(u 1 [H(uS) \H{u)\ 2 \H{UJ)\ 2 + K_ G{to). (2.8) T h e parameter K, assumes additive w h i t e G a u s s i a n noise a n d can be estimated from the variance of the noise. T h e result is to suppress the noise and degradation i n the corrupted signal, G(OJ) a n d produce an estimate of the signal. / . W i e n e r filtering has been applied for different special cases, w h i c h are described i n [17]. D i s c r i m i n a t i o n can also be performed on A-scans using parametric target models, frequency features, or time features. D i s c r i m i n a t i o n using time features is done w i t h matched filtering a n d is usually not well suited to G P R , because of the t i m e variant filtering of the subsurface. Frequency features can be d i s c r i m i n a t e d 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 c o m p a c t l y representing G P R waveforms using a m a t c h i n g pursuits a l g o r i t h m i n [52]. T h i s m e t h o d a t t e m p t e d to find compact atoms that represented fundamental structures of radiated and reflected electromagnetic waveforms w i t h success i n a p p r o x i m a t i n g time and frequency parameters of G P R signals i n noise, w h i c h proved useful for d i s c r i m i n a t i o n . 2.2.2 Review of B-scan Processing Methods T h e m a j o r i t y of B - s c a n processing is focused w i t h clutter suppression and d i s c r i m i n a t i o n . C l u t t e r suppression is a difficult and subjective topic that requires careful definition of w h a t is clutter. H o r i z o n t a l clutter is well studied, w i t h current research showing moderate success at removing the air-ground event and other persistent horizontal events. A c o m m o n form of horizontal clutter reduction is to use average trace subt r a c t i o n ( A T S ) . G i v e n a collection of traces located spatially at points i n x and y. A T S uses an average trace for a particular B - s c a n located at position y, 1 fy = ~N N z\2 f y x (2.9) x=l to perform a s u b t r a c t i o n for a given B - s c a n , fx,y â€” fx,y fy (2.10) T h i s is effective i n enhancing the hyperbolas or any non-horizontal event i n G P R data. T h e drawback is t h a t a l l horizontal events are removed or attenuated, therefore i f a target contains horizontal information A T S is not suitable. O t h e r examples focus o n l y 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 subtraction. T h e nature of s u b t r a c t i o n is t h a t A T S is not robust to phase changes or a m p l i t u d e changes. If a clutter event decays to zero, A T S w i l l continue to use the same static a m p l i t u d e calculated from the average to subtract. For c o m p u t a t i o n a l critical applications, it can be a useful m e t h o d and is used i n some d i s c r i m i n a t i o n algorithms i n [69]. G a d e r et. al uses a v a r i a t i o n that calculates the derivative i n x instead of the average, w h i c h is used i n the signal processing chain previous to d i s c r i m i n a t i o n and is h i g h l y sensitive to noise [25, 26]. T h e major p r o b l e m of A T S is the i n s t a b i l i t y w h e n s u b t r a c t i n g coherent noise. T h i s can lead to artifacts t h a t can decrease the u s a b i l i t y of a d a t a set. U l r y c h et. al studied this m e t h o d of clutter suppression w i t h respect to other c o m p e t i n g methods. In particular, the a - t r i m m e d mean was suggested as a way to handle variable clutter, but at the a d d i t i o n a l cost of careful expert decisions to choose a. T h e median was also studied, w h i c h led to [66], a s t u d y on the robustness of L-moments as estimators. N e w research was presented on the possibility of using a James-Stein estimate for noise and clutter suppression on a B - s c a n . T h i s m e t h o d uses a shrinkage parameter to adjust an local observation based on the global mean [65]. After B - s c a n processing, G P R d a t a is often times used to form a 3 D image. T h i s is referred to as migration and is the process of focusing a B - s c a n or G P R d a t a cube into a coherent image. T h i s technique is sometimes performed i n l a n d m i n e detection, b u t due to c o m p u t a t i o n a l c o m p l e x i t y and resolution it is not p o p u l a r . It does present the most understandable o u t p u t for h u m a n understanding, because the hyperbolas are collapsed and the t i m e axis is recast to depth [71, 46, 27, 43]. 21 2.3 Review of Wavelet Transforms Wavelets and Wavelet transforms are an e x c i t i n g topic of research, which was i n spired by the collaboration of M o r l e t and G r o s s m a n n [30]. T h e p r i n c i p l e of wavelets is that a signal can be represented by its frequency and t e m p o r a l content. T h e relationship between t i m e and frequency, necessarily, is a tradeoff of m i n i m a l support i n one d o m a i n d e m a n d i n g m a x i m a l support i n the other. I n the w o r l d of signals, transient events present a challenge to estimation and c o m p a c t l y representation [48]. If a time series contains a m i x t u r e of signals w i t h v a r y i n g amount of support i n b o t h domains, there is a fundamental l i m i t to the determination of a precise time 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 i n t r o d u c e d 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 p a r t i t i o n i n g of a signal into time a n d frequency atoms by the application of translated bandpass filters. These atoms, ideally, attempt to m a x i m i z e frequency and time localization. U s i n g the notation of M a l l a t [48], a wavelet has a zero average such that +oo / i/j{t)dt = 0. (2.11) -oo T o o b t a i n 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 t h a n zero. T h e wavelet transform can now 22 be defined as the convolution of the wavelet with the function under study, or (2.13) 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 corresponding 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. 23 A m o n g s t these discoveries was the work t h a t influenced the field of signal processing by D o n o h o that suggested non-linear thresholding could achieve near o p t i m a l signal estimation i n white noise [23]. T h e substance of that c l a i m relies on the orthogonality of the discrete wavelet transform. T h e continuous wavelet transform is far from orthogonal and offers no i m m e diate discrete solutions. T o solve this p r o b l e m the discrete wavelet transform makes two d r a m a t i c changes to the choice of the scale and translation parameters s a n d 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 j eZ s = 2, j (2.14) and the discrete translation parameter, k is dependent on the scale u = 2 k, keZ. j (2.15) T h i s choice of discretization of the continuous parameters is motivated by a fast implementation algorithm, but solidly rooted i n s a m p l i n g and information theory. Intuitively, the dyadic p a r t i t i o n of the scale parameter i n the D W T allows a higher localization i n low frequency, where time l o c a l i z a t i o n is necessarily m i n i m a l . The dependence of k on the discrete scale parameter, j balances this uncertainty at higher frequencies when signals become h i g h l y localized i n 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 i n N basis functions, w h i c h is equal to the length of the signal f[n]. T h e index, m is related to each translation and d i l a t i o n of the discrete wavelet basis to form g[m\. T h e discrete wavelet transform can be represented by F[m} = (f,g ), m 24 (2.16) JW] *(T) +Ai+ (fe) l G] -F i + 1 (fe) Intialize the iteration w i t h : AÂ°(k) = f(n) H = Scaling F i l t e r G = Wavelet F i l t e r T h e wavelet coefficients are F J for j = 1 , . . . , J and the scaling coefficients are F (k) = A {n). J+1 J F i g u r e 2.3: T h e D W T is implemented by an iterative application of h i g h pass and low pass filters that are constructed by the wavelet and scaling function. where g m is an orthogonal set of basis functions, (9m-, 9n) = C Often, wavelet bases are chosen such that C m n mn mal. S m n . (2-17) = 1, which makes a basis orthonor- However, the D W T is better understood as an implementation of cascaded d i g i t a l filters, w h i c h is shown i n F i g u r e 2.3. T h e l i m i t a t i o n of the D W T is the lack of shift invariance. For an intuitive understanding, consider the second coarsest scale, w h i c h o n l y contains two coefficients that represent the signal frequency and t e m p o r a l content at this p a r t i c u l a r scale. These two coefficients are sufficient for perfect reconstruction of the signal using the D W T , but if approximations or estimations i n noise are made the lack of shift invariance becomes evident. Signal estimations can be performed by thresh- 25 o l d i n g w h i c h chooses a p a r t i c u l a r coefficient based on the coefficients amplitude. Shift invariance implies that a shifted version of the signal under study w i l l produce a shift i n the wavelet coefficients, while not significantly changing the amplitude 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 i n time. 2.3.2 Redundant Discrete Wavelet Transform T h e R e d u n d a n t Discrete Wavelet Transform ( R D W T ) is a fast, convenient m e t h o d of o b t a i n i n g shift invariance. T h e t e r m R D W T is ambiguous, because the nature of o b t a i n i n g coefficient redundancy is not unique. A literature review suggests that the undecimated version of the discrete dyadic wavelet transform is the RDWT. T h i s is equivalent to other names, such as the Shift-Invariant Discrete Wavelet Transform ( S I D T W ) [45], U n d e c i m a t e d Wavelet T r a n s f o r m ( U D W T ) , Wavelet T r a n s f o r m ( S W T ) [54, 56] and the original m e t h o d Algorithme Stationary a Trous by M a l l a t . These methods are not equivalent to the C y c l e S p i n [13] or Shifted O r t h o g onal Wavelet Transforms, w h i c h o b t a i n redundancy by means of shifting the i n p u t signal and re-calculating the wavelet transform. A d i a g r a m of the R D W T i n F i g u r e 2.4 shows the s i m i l a r i t y to the D W T . E a c h level, j , is a bandpass filter of the original signal. T h e b a n d w i d t h of each filter is chosen dyadically, w h i c h 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, w i t h o u t decimation the R D W T coefficients fully represent the uncertainty of low frequency events w i t h the same s a m p l i n g resolution of the signal, thus p r o d u c i n g linearly dependent coefficients. T h e choice to use a redundant representation is a double-edged sword. P r o v e n 26 ^ + 1 (k) G j A T 2j Hi, CP Gi +1 Intialize the iteration w i t h : AÂ°(k) = f (n) HÂ° = S c a l i n g F i l t e r GÂ° = Wavelet F i l t e r T h e wavelet coefficients are. F- for j = 7 F (k) J+1 = 1 , . . . , J and the scaling coefficients are A {n). J F i g u r e 2.4: R a t h e r then decimating the signal at the o u t p u t of the filters, the filters are up-sampled w h i c h preserves the frequency resolution at all 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 t h a t the approxi m a t i o n of the signal is no longer unique. Therefore, the c o m p l e x i t y of choosing the smallest set of coefficients for a given q u a l i t y of reconstruction becomes N P complete. C u r r e n t research by D o n o h o [22] a n d G r i b o n v a l [29] suggest that if certain conditions for sparsity are met, a set of coefficients can be chosen by the m i n i m i z a t i o n of a linear program. H e r r m a n n [34] has demonstrated these results b y using linear p r o g r a m m i n g to m a i n t a i n high quality reconstruction using a sparsity constraint when choosing a set of redundant coefficients. T r o p p [64] has instead used a greedy a l g o r i t h m to choose a set of coefficients. A n o t h e r m e t h o d 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, e x p l o i t i n g the similarities that are present at different scales. 2.3.3 Curvelet Transform Curvelets provide an excellent representation of G P R events i n B-scans for signal processing. T h e nature of curvelets is that they efficiently represent a piecewise smooth curve w h i c h contain singularities n o r m a l to the direction of smoothness. T h i s property is the precise nature of G P R data, w h i c h suggests that curvelet coefficients can provide a localized sparse representation where non-linear operations, such as thresholding, can be performed. T h e same arguments that a p p l y i n the I D wavelet case work i n this 2 D representation, namely that o p t i m a l linear denoising methods, i n the sense of m i n i m u m risk, can be further improved by non-linear estimators using curvelets [7]. 28 T h e seminal curvelet papers [7] and [8] describe the nature of curvelets and their properties. A n approach to an intuitive understanding is to compare the shape of the curvelet atoms. T h e separable 2 D wavelet bases are nearly isotropic and well contained w i t h i n a the support of the scale, j. However, curvelets are designed such that width = length . T h i s makes a curvelet anisotropic and requires an a d d i t i o n a l 2 orientation parameter, 9. Curvelets also require a scale parameter, j , and position vector, fc, to represent the location and frequency content. T h e b u i l d i n g blocks for sparse representation of curves can be made using these three parameters: scale, spatial location and orientation. T o insure sparsity, the curvelets are constructed using a parabolic law scaling law, w h i c h provides a n a t u r a l extension to smooth curves. Therefore, unlike 2 D wavelets, w h i c h tend to represent a curve as a collection of points i n the wavelet d o m a i n , curvelets m i n i m i z e the number of coefficients needed to represent the curve by a d a p t i n g to the scale of the curve. T h e most i m p o r t a n t result presented i n [7] is the rate of energy decay, i n an L - 2 sense, of a incomplete reconstruction of a signal, / , using the m largest a m p l i t u d e coefficients of a a curvelet decomposition is ||/-/Â£||BCm- (logn) . 2 (2.18) 3 T h e strength of this statement is that adaptive methods, w h i c h are not necessarily tractable can approach m~ 2 as m â€”> oo, whereas wavelets approach TO . Despite -1 the log term, a curvelet representation approaches the adaptive l i m i t . T h e benefit is that curvelets are not adaptive and do not require the iterative computations necessary to adapt to a signal. W i t h the results i n [9], a new technique of constructing curvelets has introduced a tight frame, w h i c h is used i n this work, and the same a p p r o x i m a t i o n rates h o l d . 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/IIMM )< 2 V / e L 2 ( ] R 2 ) , (2.19) 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.20) 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 I D 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. .on a line. However, edges are singularities that are oriented The angular partitioning of curvelets allows these singularities to be represented by a location and orientation. W i t h 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 partition that includes not only dyadic frequency partitioning, but also a doubling of angular partitions at every other scale. Angular localization with respect to frequency 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. 30 0015 . 0 01 , F i g u r e 2.5: A n image of a single curvelet is shown. T h e curvelet is directional and localized w i t h 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 = length . W h e n 2 these functions overlap an edge that follows a s m o o t h curve, they produce large coefficients. T h i s produces a sparse representation for edges, w h i c h can be adapted to the h y p e r b o l i c events of G P R data. 2.4 Noise Removal using Threshold Estimations T h e foundation for wavelets is laid a n d the extensions to various discrete implementations are described. T h e next step explores the a p p l i c a t i o n of wavelets a n d curvelets to the estimation of a signal i n the presence of noise. Statisticians refer to this p r o b l e m as non-parametric regression, or i n engineer speak: denoising. T h e goal of noise removal w i l l be to estimate a signal, / from a noisy observation d. T h e m o d e l is additive so d = f + n , where n is the noise present i n the d a t a observation. If the d a t a is transformed into the wavelet d o m a i n using a set of orthogonal basis functions g , m by D[m] = (d,g ), (2.21) m then the u n d e r l y i n g m o d e l is also transformed such that D[m} = (f,g ) + (n,g ). m m If it was possible to know the value of N(m) (2.22) it w o u l d be a simple m a t t e r of sub- traction. Instead, a decision operator is introduced, fi[m], w h i c h w i l l shrink each wavelet coefficient of the data. T h e application is f = }~2n{m}D[m}g , m (2.23) 771 w h i c h produces an estimate, / that can be measured w i t h reference to / by the the risk: Mf) = E[\\f - / | | ] . (2.24) 2 T h e m i n i m u m risk can be obtained if the nature of the noise is k n o w n i n the wavelet d o m a i n using oracle estimation [48] w i t h a decision operator defined as n[m] = ' \D[m]\ 1 2 2 + a . 2 (2.25) T h i s holds o n l y if the noise is G a u s s i a n and white w i t h a zero mean. T h e remarkable result shown by D o n o h o [23] was that the decision operator defined i n (2.26) produces a risk w h i c h is nearly m i n i m a l w i t h respect to oracle estimation. T h i s nonlinear threshold operator w o u l d instead make a decision to either keep a coefficient that is above a threshold, T , or k i l l it, by m u l t i p l y i n g w i t h a zero. M o r e formally, h a r d thresholding is, r if \D[m]\>T, (2.26) if \D[m]\ < T 32 T h i s produces a lower risk i n the estimate / , but higher amounts of spurious oscillations, where / contains transients. thresholding. A n alternative to h a r d thresholding is soft Soft thresholding applies the same p r i n c i p l e of keep or k i l l , but re- moves the sharp discontinuities of h a r d thresholding by s h r i n k i n g coefficients above T. Soft thresholding is defined as 0 \{\D[m]\<T, n[m] = { D[m] -T i f D[m] > T D[m} + T if â€¢ (- ) 2 27 D[m]<-T D o n o h o showed i n [23] that T could be chosen by s i m p l y k n o w i n g the stand a r d deviation of the noise and the number of points, T = a^2log N. e (2.28) T h i s result is coupled w i t h the a b i l i t y to estimate the standard deviation from the finest scale coefficients of D[m]. 0 < m < N/2 If the wavelet coefficients are arranged, such that 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 i n the application of threshold estimators. T h e results were p r o m i s i n g , but a nagging p r o b l e m remained: shift invariance. T h e above results h o l d for additive G a u s s i a n white noise thresholded from the D W T . D u e to the lack of shift invariance, the resulting estimations were excellent i n the sense of risk, but lacking i n smoothness. In [45, 44], L a n g et. al were inspired by R . C o i f m a n to threshold a shift invariant wavelet transform to remove additive G a u s s i a n white noise. U s i n g the R D W T , they d i d a performance comparison on several of the I D and 2 D signals. C o m p a r i n g 33 these results to D o n o h o ' s work [21] w i t h 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 i n the a l g o r i t h m . D o n o h o and C o i f m a n also attempted their o w n versions of a shift-invariant wavelet threshold a l g o r i t h m . In [13], the authors i n t r o d u c e d the S p i n C y c l e a l g o r i t h m w h i c h used the results from [1] to produce a set of redundant coefficients to threshold. These results also outperformed classic D W T thresholding methods. T h r e s h o l d i n g i n the wavelet framework is s u m m a r i z e d by theoretical and experimental results present i n current research. T h e questions that have not been answered i n present research are the significance of choosing the number of decomposition levels and i f R D W T thresholding can be applied to colored G a u s s i a n noise using a threshold estimator. T h e next two chapters address an aspect of these questions by e x p l o r i n g methods of evaluating and a p p l y i n g 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 RDWT? In order to evaluate the suitability of R D W T thresholding a direct comparison 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, r(f) 11/-/II 2 TVa ~ Na ' 2 2 [ 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 cr of the Gaussian noise. 2 Evaluation of the maximum scale parameter, J is presented in the familiar engineering metric of signal-to-noise ratio, SNR = 101og 10 ( ^) . (3.2) Â° VII/-/II / 2 When comparing the signal, / to the estimate, / , emphasis is placed on the improvement of the output signal-to-noise ration, relative to the signal-to-noise ratio SNR 0Ut of the noisy input data, S N R i m p = SNR 36 o u t - SNR i n . (3.3) T h i s emphasizes the changes i n performance w h e n experimental parameters are chosen. T h e experiments are conducted i n M A T L A B using R i c e U n i v e r s i t y ' s i m plementation of the R D W T from the R i c e Wavelet T o o l b o x [31] and the wavelets are generated w i t h Wavelab [24]. A description of the ingredients to the RDWT thresholding process w i t h a discussion of the rationale and relevance follows. The parameters of R D W T thresholding are explained before the G P R and noise m o d e l is described. 3.1 Determining a Basis for A-Scan In M a l l a t ' s Wavelet Tour of Signal Processing Approximation he labels a chapter: ESTIMATIONS 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 elegant result is that signals are sparsely represented for a given wavelet basis function. T h i s property makes estimation of a signal i n noise simple, because the correlation of the signal to the wavelet basis provides h i g h a m p l i t u d e coefficients, w h i c h 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 a p p r o x i m a t e d w i t h a s m a l l number of coefficients. D a t a adaptive methods [14, 49] choose a "best" basis from a d i c t i o n a r y of bases using a sparsity measure to make a decision. R a t h e r t h a n form a best-basis, this work instead uses the log-energy d i s t r i b u t i o n measure, N E|F[m]| log|FH| . 2 2 mâ€”l to choose a particular mother wavelet. 37 (3.4) 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 t ip{t)dt = 0. v / -oo (3.5) A n intuitive understanding of vanishing moments is that a wavelet acts a differentiator. 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 Daubechies4 Symmlet4 Coiflet2 Battlel Symmlet6 Coiflet3 Daubechies6 Coifletl Symmlet5 Battle3 Entropy -302.3 -294.6 -290.1 -289.5 -273.5 -272.6 -270.4 -268.5 -266.2 -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 approximate 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 assumption for choosing a wavelet basis. 3.2 D e t e r m i n i n g the N u m b e r of Levels of D e c o m p o s i tion 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 N levels. However, it is not clear if a signal should be decomposed to the 2 maximum scale. As the parameter J increases, the support of the wavelet increases to a maximum of TV when J = log N. Because the efficiency of a median operator 2 increases with N and the advantage of thresholding decreases with sparsity, it is expected that the optimal value of J is less then log N. 2 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 a p p l i c a t i o n is appropriate. T h e evaluation of the scale parameter compares the signal-to-noise improvement calculated from (3.3). T h e number of levels, additive white noise a n d colored noise, and the energy of the noise is varied to illustrate the scale parameter effect. T h e measurements are performed on a synthetic G P R trace. E v a l u a t i o n of scale parameter experiment: 1. Generate noisy d a t a from the a d d i t i o n of a m o d e l and G a u s s i a n Noise. 2. Decompose to J levels using the R D W T . 3. E s t i m a t e 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 6. Repeat steps 1-5 to o b t a i n an average SNR . in SNR^mp. T h e reference m o d e l is a synthetic G P R pulse that models an a p p r o x i m a t i o n of a physically realistic version of equation (2.4). T w o different experiments are conducted, one to examine the effects of G a u s s i a n white noise and the second for G a u s s i a n colored noise. T h e coloration of the noise is described i n Section 3.4.2, w h i c h assumes the difficult case where the noise has nearly the same characteristics as the signal. T h e i n p u t noise is measured w i t h (3.2) and varies from - 9 d B up to 3 0 d B i n 3dB increments. T h e range of SNRi n is representative of the signal attenuation as it propagates into the ground, w h i c h makes late time reflection events decay into the noise floor of the receiver. For each value of SNRi , n 40 the number of levels is 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 N. 2 conducted twenty times and the average SN-Ri mp E a c h experiment is is calculated w i t h the mean. T h e result form a surface where an o p t i m a l number of levels can be chosen for the particular signal and noise type. However, because choosing the number of levels is equivalent to bandpass filtering, a more generalized understanding can be i m p l i e d 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 o n l y suppress a coefficient i f it is below a predicted value of noise. For the D W T a global threshold estimator can approximated for all coefficients. T h i s global threshold is valid, because the noise remains white i n the wavelet d o m a i n w i t h 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 i n d i v i d u a l l y . D o n o h o proposes i n [23] that the o p t i m a l threshold for the D W T is, T = cV21og /Y, e (3.6) which, for certain classes of signals, achieves nearly m i n i m a x risk w i t h respect to the o p t i m a l non-linear decision operator. T h i s result makes the removal of white noise i n the D W T d o m a i n a simple matter of estimating the standard d e v i a t i o n of the noise. U s i n g the equation (2.29), an estimate of a is reasonable i f the signal is sparse relative to the noise. F o r 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, w h i c h makes the signal sparse relative to the noise. F o r the R D W T , each scale provides JV coefficients to make an estimate, but the sparsity of the signal â€¢41 coefficients decreases w i t h j. Therefore to generate a threshold estimator requires either, advance knowledge of the noise behavior or m u l t i p l e realizations of a signal w i t h i d e n t i c a l noise distributions. U s i n g a h y b r i d m e t h o d , w h i c h neither sets a threshold for each coefficient, nor uses a global threshold for all 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 = a /log N. jy (3.7) e T h i s is applied to the wavelet coefficients using a h a r d thresholding technique given in equation (2.26) to the m a t c h i n g level of coefficients. T h e assumption is t h a t the decomposition is reasonably sparse and an estimate of a w i l l remove the m a j o r i t y of noise w i t h o u t being too aggressive. T h i s assumption is tested b y the scale parameter experiments outlined above and a calculation of risk w i t h respect to the c o m p e t i n g methods. T h e expectation is that the scale parameter has a direct effect on the ability to remove noise. A s the J increases, the decreasing sparsity introduces more bias i n the coefficients m a k i n g the m e d i a n estimator too aggressive. T h e alternative, which is to estimate a global threshold from o n l y the finest coefficients is sufficient and successful' for white noise. T h e results of L a n g [44], M a l l a t [48] and D o n o h o [21] justify this m e t h o d of thresholding. However, i n the case of colored noise, the average standard d e v i a t i o n w i l l vary drastically for each level, j. Therefore, a global threshold does not apply. T h u s , 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 performance of classic wavelet denoising using the D W T and W i e n e r filtering. A calculation of the n o r m a l i z e d risk for the three methods is presented as 42 increases. T h i s 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 i n [48] to the R D W T level-constant threshold m e t h o d . 3.4 G P R A-Scan Model A G P R A - s c a n can be measured as a voltage, V rx from a receive antenna. T h e result is a convolution of several impulse responses from the environment a n d the system w i t h a source function, V rx V, s = V (t)*h (t)*h {t)*h {t)*h (t)*h {t) s rx tx c g + n{t). t T h e impulse response from cross-coupling, h (t) c [18] (3.8) is described i n [18] and is one source of horizontal clutter i n a B - s c a n . T h e impulse response of the ground, h (t) g is quite complicated as it depends on the material, moisture content and granularity. D u e to the variation w i t h time, accurate ground models should be considered uniquely for specific a p p l i c a t i o n . A simplification of the ground m o d e l can be constructed using a constant for a time-invariant frequency response. is constant for all frequencies. attenuation T h i s methodology assumes the g r o u n d Refer to the work of I r v i n g [41] for a frequency dependent inversion of simple ground cases. T h e target impulse response, h (t) t is also complex, but simple reflectors such as the ground or point reflectors can be roughly a p p r o x i m a t e d as specular. The i n p u t signal, V is modeled as the derivative of a G a u s s i a n b u m p like i n (2.4). T h i s s is a simplification, because a real pulse tends to be a s y m m e t r i c and long-tailed, but depends u p o n the transmitter design and operating characteristics. T h e impulse responses of the receiving and t r a n s m i t t i n g antenna are the same for a G P R i n mono-static operation, w h i c h makes h (t) rx 44 = h (t). tx T h e development of a synthetic antenna m o d e l , h (t) n is described i n Section 3.4.1. U s i n g the antenna m o d e l , the source pulse is filtered by h (t) to form a synthetic t r a n s m i t t e d pulse. n T h e reflection events are specular, w h i c h corresponds to a dirac function i n t i m e m a k i n g the final G P R pulse m o d e l , V rx = V {t) * h {t) s n * h {t) n * 8{t - U). (3.9) T h e target locations are located at the times i n the array i j , w h i c h is the r o u n d t r i p time from the antenna. T h e pulse must be filtered twice by the antenna, once d u r i n g transmission a n d once u p o n reception, thus the two antenna responses i n (3.9). The physical parameter, F W H M , used to construct V s is 100ns, w h i c h is a p p r o x i m a t e l y the F W H M of the experimental G P R demonstrated by the B e l g i u m R o y a l M i l i t a r y A c a d e m y H u m D e m research group [59]. T h i s compliments the antenna m o d e l , w h i c h is a s i m u l a t i o n of the same experimental G P R . The simplifications a n d assumptions made to produce the G P R impulse m o d e l are based on m e t h o d of noise removal. Wavelets can be thought of as a singularity detector. W i t h i n limits, the singularity can undergo linear filtering w i t h o u t d r a m a t i c a l l y affecting the ability of a wavelet to efficiently represent it. There- fore, the simplification of the g r o u n d and target to specular reflectors is reasonable. A d d i t i o n a l l y , ignoring the cross-coupling t e r m is justified, because it should have similar frequency a n d t e m p o r a l support as reflection events m a k i n g 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=h tx and hN = W â€”â€”T==KX >rx V Â« Vh Z V a (3.10) V h which describes the frequency dependent effects of the ground, f , impedance of the g system and antenna Z , Z c a 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 G (u) t = -^\H (u>)\ . (3.11) 2 N 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, H (z) N = â€¢ (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 YuleWalker fit is compared to the piece-wise estimation. The final assumption is that h (t). n h,N (t) = hpf. x(t) tix r a n d notation is simplified to 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) F i g u r e 3.2: U s i n g the Y u l e - W a l k e r m e t h o d to estimate an example antenna response. pulse by the simulated antenna response. F i g u r e 3.3 shows the result of the physical a p p r o x i m a t i o n of the G P R pulse. 3.4.2 Noise Model T h e noise model, n(t) is straight forward to m o d e l i n M A T L A B . T o produce G a u s sian w h i t e noise, the M A T L A B function r a n d n generates an a r b i t r a r y m a t r i x values w h i c h is n o r m a l l y d i s t r i b u t e d , zero mean and has a standard deviation of one. T h i s is referred to as the additive G a u s s i a n w h i t e noise, or A G W N m o d e l . T h e A G W N model is used to verify the performance of the R D W T level-constant threshold m e t h o d to the classic m e t h o d using discrete wavelet thresholding. T h e colored noise m o d e l is created using the antenna m o d e l from the previous section and G W N . T h e additive G a u s s i a n colored noise, or A G C N is filtered w h i t e 47 Simulated GPR Pulse 500 1000 1500 2000 2500 3000 Time (ps) 3500 4000 4500 5000 F i g u r e 3.3: T h i s pulse shows a n t y p i c a l pulse that might be received w i t h only one reflector at 1000ns. noise a n d is represented as a n operator o n the additive noise m o d e l , d(t) = f(t) + hn(t)*n(t). (3.13) T h e antenna m o d e l ensures that the G C N is " i n - b a n d " , m e a n i n g 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. The A G C N m o d e l is n o t unique a n d therefore does not represent a n ex- haustive s t u d y of A G C N noise removal. Conceptually, a n A G W N m o d e l has a flat s p e c t r u m . I n the context of thresholding, a flat s p e c t r u m means that i f the spect r a l content of the signal is less then the noise, it is thresholded t o zero. F o r the A G C N , the spectral content is similar t o the signal. T h i s makes the estimation a n d suppression m u c h more difficult. These two cases should garner a general i n t u i t i v e understanding, w h i c h w i l l allow a n extension t o other cases. 48 3.5 Summary T h i s chapter i n t r o d u c e d the methodology for testing the R D W T i n the presence of G a u s s i a n colored noise. T h e next chapter w i l l introduce the methodology for suppressing clutter i n B-scans before the final results are presented. The 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 m e t r i c is used to evaluate results. â€¢ A noise m o d e l for b o t h colored a n d white G a u s s i a n noise is shown. T h e experimental results use these b u i l d i n g blocks to evaluate the effect of choosing the m a x i m u m number of decomposition levels on the ability of the R D W T to estimate a G P R signal i n white or colored noise. I n order to establish a baseline performance, experiments are designed to compare the D W T and R D W T i n the presence of w h i t e noise, before a performance evaluation on colored noise. 49 Chapter 4 Clutter Suppression i n G P R B-Scans C l u t t e r suppression is a c o m m o n step i n the signal processing chain towards the goal of l a n d m i n e detection. T h e purpose is to m i n i m i z e damage of the target signal event, while suppressing a n y t h i n g that is not related to the target. T h e r e are several reasonably successful methods of clutter reduction present i n current literature. F o r a general overview the book w r i t t e n by D . Daniels, Surface-Penetrating Radar [18] and a more comprehensive review i n the recent b o o k edited by the same author, â€¢Ground Penetrating Radar [19] are recommended. A brief overview of clutter suppression methods w i l l be presented here to develop a m o t i v a t i o n for the proposed method. A d a p t i v e methods, such as the one i n t r o d u c e d by B r o o k s et. al [4] uses a linear time-invariant m o d e l to perform a non-parametric regression of the system response to separate clutter. A n o t h e r adaptive m e t h o d of clutter suppression attempts to suppress m u l t i p l e reflection events, where the author uses an L - l weighted 50 estimator instead of the t r a d i t i o n a l least squares filter to derive a clutter m o d e l [33]. V a n K e m p e n et. al examines the d a t a for a set of features to form a m o d e l and suppresses clutter w i t h a modified W i e n e r filter [67]. A non-adaptive m e t h o d is presented i n [63] where a p r i n c i p a l component decomposition is used w i t h linear estimation. T h e assumption is that the target is contained w i t h i n the largest principle components and a linear a p p r o x i m a t i o n w i l l achieve clutter suppression. A n alternative to regression models, is d i s c r i m i n a t i o n , where clutter is suppressed b y choosing a target, w h i c h forces all other d a t a to be accepted as clutter. Zhao proposes a successful m e t h o d of d i s c r i m i n a t i o n using a discrete hidden M a r k o v model, w h i c h finds b o t h metal and non-metal landmines i n G P R d a t a [70]. T h i s work is derived from the results of H e r r m a n n i n the seismic realm [36]. H e r r m a n n uses the parsimonious curvelet representation of seismic d a t a to iteratively suppress m u l t i p l e target reflections. H e r r m a n n ' s work tackles the difficult p r o b l e m of selectively removing events that are nearly identical to the target, but smaller i n a m p l i t u d e and later i n time [34]. H e r r m a n n ' s m e t h o d has also been extended to removing t e m p o r a l clutter using m u l t i p l e d a t a sets from the same l o c a t i o n , b u t several years apart [16] and reducing the dimensionality of the seismic i m a g i n g p r o b l e m [35]. T h e m e t h o d proposed i n this work also supplements the results presented by N u z z o and Q u a r t a [55], w h i c h compares r-p and wavelet transforms for effective horizontal clutter suppression. N u z z o et. al used the 2 D D W T to selectively threshold only the horizontal components. T h e two m a i n features presented i n their m e t h o d are: â€¢ Identify the wavelet scales that clutter appears i n . 51 Original C u r v e l e t s : 41.6 d B Wavelets: 38.4 d B F i g u r e 4.1: T h i s example demonstrates reconstruction of a simulated h y p e r b o l a w i t h o n l y 0.5% of the curvelet and wavelet coefficients. N o t i c e how curvelets m a i n t a i n the smoothness of curve and still reproduce the discontinuity across the curve. â€¢ Suppress horizontal wavelet coefficients at identified scales. T h e a l g o r i t h m proposed makes several i m p o r t a n t changes to this general approach. F i r s t , the curvelet transform is used, instead of the 2 D wavelet transform. F i g u r e 4.1 demonstrates the m o t i v a t i o n of using curvelets to approximate G P R reflection events i n B-scans. Second, a noise m o d e l is generated using features from the d a t a set, w h i c h intrinsically define the wavelet scales and horizontal clutter features. A d d i t i o n a l l y , the noise m o d e l forms the threshold operator, w h i c h puts the a l g o r i t h m into the well researched and understood d o m a i n of wavelet thresholding. L i k e clutter, the noise m o d e l should be determined for each a p p l i c a t i o n . 52 Therefore, a synthetic clutter and target m o d e l is formed from actual l a n d m i n e data. T o verify the noise m o d e l , the synthetic B - s c a n is generated using the physical a p p r o x i m a t i o n of A-scans developed i n C h a p t e r 3. T h e description of clutter suppression methodology begins w i t h the synthetic B - s c a n m o d e l . Afterwards, the clutter and target models are described and the proposed a l g o r i t h m is introduced. T h e chapter concludes w i t h a discussion of performance metrics for evaluation of the proposed m e t h o d . 4.1 B-scan Model using Synthetic G P R Pulse T h e synthetic B - s c a n is carefully modeled such that reasonable clutter and target approximations are made. T h i s allows a careful s t u d y of the efficiency and l i m i t a tions of the proposed m e t h o d of clutter suppression. In order to create a B - s c a n , targets are placed into a time vector and convolved w i t h 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 a m p l i t u d e . For example, the ground might be a dirac function w i t h u n i t a m p l i t u d e o c c u r r i n g at lOOOps. D u e to the F W H M of lOOps, the s a m p l i n g rate is lOps, and an index of 100 w o u l d contains the target. 4.1.1 Clutter Model T h e clutter m o d e l used for s i m u l a t i o n is an a p p r o x i m a t i o n of clutter types that exist i n a real G P R scan. F i g u r e 4.2 shows a real G P R survey of soil w i t h no targets. T h e synthetic clutter m o d e l reflects the three highlighted events of 4.2. T h e clutter is scaled to u n i t a m p l i t u d e for all events. E a c h event is similar to the physical events present i n F i g u r e 4.2. T h e specific time locations of clutter 53 F i g u r e 4.2: T h i s G P R B - s c a n contains examples of clutter w i t h o u t targets. A n notation A shows the air-ground event, the clutter i n oval B is a highly variable a m p l i t u d e event and C shows a clutter event that does not exist for the d u r a t i o n of the scan. 51 are not critical and a r b i t r a r i l y defined as: â€¢ Synthetic Event, EA' G r o u n d , at lOOOps w i t h constant a m p l i t u d e â€¢ Synthetic Event, E ^ : C l u t t e r at 1250ps w i t h v a r y i n g a m p l i t u d e â€¢ Synthetic Event, Ec'- C l u t t e r at 1500ps w h i c h does not persist T h e a m p l i t u d e of each event is defined by a function, before it is convolved w i t h the synthetic pulse, E =1 (4.1) A â€¢ ' E B =cos(2vr0.02x), 1 > x < 50 1 > x < 24 0, E C = < (x-25)/ 4 e e 25 > X < 29 ) ' (4-3) 30 > x < 50 1, 4.1.2 (4.2) Target M o d e l T h e target model is generated by calculating a spatial h y p e r b o l a attenuated by an exponential attenuation function dependent on the distance from the antenna to the target center. T h e spatial h y p e r b o l a is calculated by creating an array of time shifts, t(x) = -^Z 2 + (x- x) . for x = 0 , 0 . 0 1 , . . . , 0 . 2 5 c m . 2 t (4.4) v T h e velocity is defined i n (2.7), the target center, x , is 25 c m and the target depth, t Z, is the difference between target depth and ground, w h i c h is 5 c m . T h e target time is converted to a M A T L A B vector index by scaling (4.4) w i t h the time s a m p l i n g rate, lOps and r o u n d i n g to the nearest integer, i d x = round ( ^ ^2 ) . Vl0xl0 y 55 110 (4.5); V (a) Point Target (b) Large Target F i g u r e 4.3: E x a m p l e of synthetic targets used i n 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 -2ad B(iebc) = - - ^ - . e (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 . found by d = \J Z 2 + (x â€” x ) 2 t T h e distance from the antenna to the target is a n d the target gain function, -2aZ p G = ~ ^ n (-) 4 7 normalizes the peak of the the h y p e r b o l a to unit a m p l i t u d e . T h e n o r m a l i z a t i o n puts the clutter, target and ground on similar a m p l i t u d e scales, w h i c h is often the case. T h i s represents a reasonable model for a point target. A n a d d i t i o n a l target is used, w h i c h has larger dimensions. T h i s w i l l result i n a horizontal component that facilitates testing of the average trace subtraction m e t h o d . It is calculated i n the same manner, except the center of the h y p e r b o l a is extended horizontally, as presented i n F i g u r e 4.3. 56 4.2 Proposed Method of Clutter Suppression A new m e t h o d is proposed for suppressing horizontal clutter. T h e foundation is that the existing prevalent m e t h o d , average trace subtraction ( A T S ) is often times sufficient for h i g h l i g h t i n g hyperbolas i n G P R d a t a for reasons of c o m p u t a t i o n a l complexity. Therefore, new methods must be reasonably fast or especially successful in suppressing clutter. T h e new w r a p p e d curvelet transform, w h i c h is released i n C u r v e L a b [10] is sufficiently fast to consider as an a p p l i c a t i o n clutter suppression. A d d i t i o n a l l y , because the curvelet d o m a i n is efficient at representing G P R traces, thresholding should prove to be c o m p e t i n g m e t h o d to A T S . T h e act of thresholding "mutes" an event, w h i c h meets the c r i t e r i a of the threshold. U n l i k e subtraction, thresholding does not a d d energy, therefore eliminating the possibility of constructive interference. If an appropriate threshold can be calculated that represents the clutter contained i n the B - s c a n then a simple scaled threshold can be applied i n the curvelet d o m a i n . T h i s w i l l suppress events that are clutter like and ignore events that do not m a t c h the estimated m o d e l . T h e clutter model becomes the key to successful i m p l e m e n t a t i o n of the proposed m e t h o d . T h e clutter m o d e l for A T S is the average trace, w h i c h has two major assumptions. T h e first is that the resulting average a m p l i t u d e is representative of all 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 positional noise and subsurface material changes. T h e second assumption is that the phase remains the same. T h i s is usually true, however, not all h o r i z o n t a l clutter events are present over the w i n d o w of a-B-scan." If a clutter event o n l y exists i n a p o r t i o n of a B - s c a n the resulting subtraction w i l l actually place an event where none exists. T h i s w o u l d be a similar case i f the clutter event changed phase 180 degrees resulting i n constructive interference. 57 T h e proposed m e t h o d overcomes the pit-falls of A T S by a d d i n g a m p l i t u d e and phase resilience, but retaining a simple noise m o d e l . If a quasi-accurate clutter m o d e l is generated, the phase and a m p l i t u d e resilience of the m e t h o d comes from the parsimonious representation of clutter. T h i s enables clutter to be muted w i t h a high p r o b a b i l i t y that the target is preserved. T h e amount of m u t i n g comes from a s h r i n k i n g parameter, A, w h i c h must be e m p i r i c a l l y 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 i n phase and a m p l i t u d e over the spatial dimension of the B-scan w i n d o w . A n alternative clutter model is d u b b e d the edge model. T h e edge m o d e l uses the average of only two traces located at the edges of the B - s c a n . If there are M traces i n a B - s c a n , w h i c h are c o l u m n vectors, then the 1-D edge noise model is hedge = ~ ( l b + b ) M â€¢ (4.8) T o 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 time and M for space. For each trace i n the 2-D noise model image, Q, is the calculated edge model trace m u l t i p l i e d by a l e n g t h - M vector of ones, (4.9) T h e r e are several reasons that make a noise m o d e l calculated from the edge traces a reasonable attempt to m o d e l horizontal clutter. T h e first is t h a t the clutter is h i g h l y likely to be present i n 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 w o u l d 58 be detected by an noise m o d e l calculated from the average trace. However, if the clutter is not present at the edges, the edge trace m o d e l does not work. T h e proposed m e t h o d of clutter suppression is s u m m a r i z e d : 1. Take the curvelet transform of the d a t a B - s c a n . 2. C a l c u l a t e the edge noise m o d e l and form a B-scan noise m o d e l . 3. Take the curvelet transform of the noise m o d e l B - s c a n . 4. Use the noise m o d e l i n the curvelet d o m a i n as a threshold. 5. Scale the threshold and a p p l y to the d a t a curvelet coefficients. 6. Inverse transform the thresholded data. In notation, the set of curvelet basis vectors, g produces a set of coefficients m from the data, d and noise model, Q, d =(d,g ) m (4.10) m Cl ={n,g ). m (4.11) m T h e threshold operator formed from the noise m o d e l 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 w i t h F, and reconstructing, d = }^F{d )g . m m (4.13) m T h i s m e t h o d has the advantage of being quite easy to visualize the result. T h e noise m o d e l is "muted" from the d a t a b y zeroing o n l y the coefficients that are 59 estimated to belong to noise. T h e phase of the d a t a is not i m p o r t a n t , because the absolute value of the coefficients is considered when thresholding. T h i s is allows incoherent estimation, because the noise m o d e l o n l y provides a spatial location and approximate a m p l i t u d e . Coherent estimation, such as subtraction, requires careful alignment of phase and a m p l i t u d e to prevent constructive interference. T h e am- plitude resilience of thresholding comes from control of the A parameter. A large value increases the tolerance for a m p l i t u d e variation at the cost of d a m a g i n g closely spaced events. 4.3 Method of Performance Evaluation T o verify that the proposed m e t h o d is successful, several means of verification are presented. U s i n g the synthetic data, numerical results are presented by calculating the P S N R , (4.14) T h e t r a d i t i o n a l image processing definition of P S N R uses A = 255, w h i c h 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. T o calculate A i n this a p p l i c a t i o n , the difference of the m a x i m u m and m i n i m u m amplitudes are taken from the i n p u t data, A = max|ct| â€” m i n | d | . (4.15) N u m e r i c a l metrics are useful for c o m p a r i n g the amount of energy removed from the noisy data, b u t visual quality is not well represented. Therefore, the n u m e r i c a l results are presented to verify that clutter suppression is successful i n removing the clutter energy, b u t visual experiments are conducted to demonstrate the difference i n artifacts. 60 T h e v i s u a l experiments allow a direct comparison of the proposed m e t h o d and A T S . Image artifacts, like target energy damage, are best viewed b y comparison of the i n p u t a n d output images of each a l g o r i t h m . A residual image, d â€” d is also calculated, w h i c h emphasizes the artifacts either created, or not removed by each method. 61 Chapter 5 Experimental Results and Discussion T h e results presented i n this chapter are from two separate methods. T h e first set of results show the s u i t a b i l i t y of a level-constant thresholded R D W T for removing white or colored G a u s s i a n 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 - s c a n noise address these issues: â€¢ Performance w h e n additive noise is G a u s s i a n and white. â€¢ Performance w h e n additive noise is G a u s s i a n and colored. â€¢ E v a l u a t i o n of scale parameter and i n p u t noise energy. â€¢ C o m p a t i b i l i t y of proposed m e t h o d to trace stacking. In a d d i t i o n to the I D results, experimental results for B - s c a n clutter suppression w i l l compare two methods, the proposed curvelet thresholding m e t h o d and 02 ' average traces s u b t r a c t i o n . T h e experiments present: â€¢ V i s u a l performance for synthetic data. â€¢ R e s i d u a l images for synthetic data. â€¢ P e a k 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 m a i n results of this section are presented i n T a b l e 5.1'. T h i s experiment uses the best choice for a scale parameter determined a priori. T h e i n p u t noise is additive w i t h a signal-to-noise ratio of 6 d B . T h e signal is a simulated G P R pulse at lOOOps containing 2048 points. RDWT DWT W h i t e Noise C o l o r e d Noise 17.8 d B 13.5 d B 13.9 d B 0.1 d B Table 5.1: S u m m a r y of the average SNRi of the level-constant thresholded R D W T and global thresholded D W T w i t h respect to white noise and colored noise. rnp 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 s u i t a b i l i t y of level-constant R D W T thresholding as a m e t h o d of noise removal. The next set of experiments evaluate the specific choice of parameters and the influence on this m e t h o d . 63 O r i g i n a l w i t h o u t noise Original with A G W N : 5 d B Original with A G C N : 6 d B R D W T a n d A G W N : 18 d B R D W T a n d A G C N : 14 d B D W T a n d A G W N : 14 d B D W T and A G C N : 1 d B F i g u r e 5.1: T h e R D W T visually removes the noise and preserves the signal using a level-constant threshold. T h e D W T w i t h a global threshold can effectively remove the w h i t e noise, b u t not colored noise. 64 5.1.1 Performance of Gaussian W h i t e Noise Removal T h i s experiment uses the n o r m a l i z e d risk measure from equation (3.1) to evaluate the performance of the R D W T to the D W T a n d W i e n e r filtering. T h i s n o r m a l i z e d risk, | | / â€” f\\ /Na , 2 2 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 i n p u t G a u s s i a n white noise for a simulated G P R signal w i t h increasing values of N. SNRi n T h i s makes of the noise lower as N increases, because the energy of white noise is deter- m i n e d by HiVcr 1|. T h e normalized risk is similar to a signal-to-noise measurement, 2 but measures r e m a i n i n g noise energy on a scale of 0 to 1. A n o r m a l i z e d risk measurement of one means that no noise energy is removed and a zero means a l l noise is removed, w h i l e 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 m e t h o d 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]. The R D W T threshold is determined by the level-constant m e t h o d using an estimate of each level w i t h J â€” 8. T h e results show a decay for D W T thresholding w h i c h exceeds W i e n e r filtering, confirming the results presented i n wavelet literature [48]. T h e R D W T removes more noise initially, 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 s m a l l values of constant threshold. 65 when using a level- Risk decay of signal estimation in A G W N 0.08 0.07 - RDWT Wiener DWT 0.06 CM b 0.05 0.04 0.03 0.02 \ \ \ \ x \ \ \ \ \ \ \ . \ 0.01 0 \ 10 â€”. 11 12 13 log N points â€¢â€¢ 14 15 2 F i g u r e 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. Risk Decay for A G C N with constant a 1.4 1.2 RDWT Wiener DWT 1 Â£ 0.8 0.6 11 12 13 N points (log scale) 14 15 F i g u r e 5.3: G l o b a l thresholding of the D W T removes no measurable noise, thus p r o d u c i n g no risk decay. W i e n e r filtering is m a r g i n a l l y 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 a b i l i t y of each c o m p e t i n g m e t h o d i n the presence of "in-band" colored G a u s s i a n noise. T h e experiment uses the normalized risk metric so that a comparison can be made t o the previous results of A G W N noise removal. T h e experiment is conducted i n the same manner as the previous section using a synthetic G P R pulse w i t h a n increasing amount of points. T h e value of a remained constant, w h i l e N increased. T h i s choice reduces t h e amount of noise energy as the number of points increased. T h e noise is colored b y the antenna model, h , and generated r a n d o m l y for each t r i a l . T h e results represent the average n normalized risk for 20 trials. T h e r e are three o u t s t a n d i n g results i n F i g u r e 5.2. F i r s t , the D W T m e t h o d of global thresholding appears to by unable t o reduce the amount of noise, regardless of the number of points available. a decay i n risk as T h e R D W T m e t h o d , however, s t i l l produces increases i n the presences of G C N . F i n a l l y , W i e n e r filtering produces m a r g i n a l results, while the risk for R D W T thresholding remains effective, though five times larger w i t h respect t o A G W N 5.1.3 performance. Dependencies of Number of Decomposition Levels T h i s experiment evaluates the dependence of the scale parameter and SNRi n the SNRi mp using metric. For F i g u r e 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 m e t h o d i n the presence of white noise. T h e results show that as SNRi n decreases t o zero, decomposing t o a higher number of levels yields slightly better results. I n the case of m i n i m a l noise, like 3 9 d B , 67 F i g u r e 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 remains the same, but the metric SNRi mp SNRi . T h e o p t i m a l choice of J n falls off as noise energy increases. The 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 i m p o r t a n t to note that if clutter is modeled as in-band noise, then it would be suppressed by a p p r o x i m a t e l y l O d B . 68 F i g u r e 5.5: A s the i n p u t S N R of the colored noise reaches OdB, the ability to distinguish noise from signal diminishes, because noise events appear nearly identical to signal events, however the o p t i m a l choice of J remains the same. 5.1.4 Thresholding Multiple A-scans before Stacking T h e first experiment is designed to test the efficacy of operator order. T h r e s h o l d i n g is a non-linear operator, w h i c h means the order of operations w i l l likely produce different results. T h e ubiquitous nature of stacking suggests that the order of thresho l d i n g and stacking should be examined. F i g u r e 5.6 demonstrates the difference between the three different cases: â€¢ P e r f o r m no wavelet thresholding, just stack. â€¢ Perform wavelet thresholding before stacking. â€¢ P e r f o r m 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 i n p u t noise is 9 d B a n d 69 Thresholding performs better after stacking ' â€” O n l y Stacking Before Stacking After Stacking PQ T3 a > O S-i e OT 12 4 8 16 N u m b e r of traces F i g u r e 5.6: S t a c k i n g scales the i n p u t noise a by 1/VT, however, this scaling is d i m i n i s h e d by stacking before thresholding. Regardless of the order, wavelet noise removal still provides a significant improvement. J = 8, w h i c h corresponds to the o p t i m a l scale parameter from F i g u r e 5.4. SNRi p m The 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. T h i s is an i m p o r t a n t and fortunate feature, because most G P R d a t a is o n l y available post-stack. T h e experiment does not test the assumption that each trace is aligned i n time and phase, w h i c h is a basic assumption of stacking. However, if o n l y one trace is available the results show that wavelet thresholding can still 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 l e effect on the results presented i n F i g u r e 5.6. T h i s set of experiments tested four cases: A G W N w i t h 2 and 32 i n p u t traces a n d A G C N w i t h 2 and 32 traces. E a c h case was tested as the number of points i n the i n p u t 70 25 40 O Median â€¢ â€¢ Mean Stacked 35 â€” S 'â€¢imp { 5 If fe; fe. to to lOcfc O Median â€¢ â€¢ â€¢ Mean " ~ Stacked 30 25 20 15 9 10 11 log iV points 10 12 9 10 11 12 log A points r 2 2 (b) White noise, T = 32 (a) White noise, T = 2 O Median â€¢ Mean ' ~ Stacked â€” â€¢' o fe. to > 10 11 p.... 9 o. 10 p. â€¢ 11 12 u l o g / v points log N points 2 2 (d) Colored noise, T = 32 (c) Colored noise, T = 2 F i g u r e 5.7: Stacking before thresholding is the better. If stacking is performed i n the wavelet d o m a i n , the median is an improvement over the mean when T = 32. signal varied from 256 to 32,768. In a d d i t i o n to testing the order of stacking and thresholding, a modification to the stacking operator was made. T h e m e d i a n operator is used to perform a different estimate of a mean trace after thresholding i n the wavelet d o m a i n . T h e median operator is more suitable to determining the mean i f a p r o b a b i l i t y d i s t r i b u t i o n function is skewed. T h e nature of thresholding suggests that the efficiency of m e d i a n stacking should be tested. T h e results i n F i g u r e 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 3 d B of improvement 71 i n large stacks. M i n o r differences aside, noise coloring has no effect on the order of operations. E v e n i n 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 denoising a synthetic G P R signal i n G a u s s i a n noise. Several aspects are tested and the scale parameter, J , is shown to have a d r a m a t i c effect on the ability to perform noise removal. Noise coloring is shown to reduce the performance of level-constant thresholding, although performance still exceeds other methods tested. T h e order of operations is shown to be best performed by finding the stacked trace before wavelet denoising. T h e performance difference of the level-constant thresholded R D W T between w h i t e and colored noise is due to the a b i l i t y to estimate the noise i n the redundant wavelet d o m a i n . U s i n g the redundant transform, a different estimate of a is obtained at each level. In F i g u r e 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 a p p l i c a t i o n of bandpass filters, the level-constant threshold is an estimate of the standard deviation for a p a r t i c u l a r frequency b a n d 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 i m i t a t i o n s , because i f the est i m a t i o n is wrong, the corresponding threshold w i l l not be effective. T h e results from testing the scale parameter i m p l y that when the noise is especially low, the threshold operator removes more signal energy then noise. T h e extreme case is i n F i g u r e 5.4, where the result of denoising makes the corresponding signal worse then the i n p u t . 72 W h e n c o m p a r i n g R D W T denoising to W i e n e r filtering and D W T thresh- o l d i n g i n F i g u r e 5.2, the n o r m a l i z e d risk of the R D W T shows better performance i n i t i a l l y then either m e t h o d . A s N increases, the difference shrinks and D W T denoising overtakes the R D W T m e t h o d . T h i s is a t t r i b u t e d to the choice of level-constant thresholding of the R D W T . T h e median, w h i c h 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 frequency levels where the sparsity property no longer holds. T h i s means the threshold predicts a large value for dj, w h i c h removes more signal then noise at that level. For white noise, it w o u l d 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 a p p l y i n g thresholding techniques to G P R data, results show that thresholding s h o u l d 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 w h e n the d i s t r i b u t i o n 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 reduced, especially when T = 32. T h i s is because the estimation of the threshold is improved and quality of noise removal increases. However, the difference can not overcome the tremendous gain of pre-stacking, w h i c h makes thresholding more effective. However, this result does not take into account that a threshold could be estimated for each coefficient i n the redundant d o m a i n . A level-constant threshold 73 assumes all coefficients at a p a r t i c u l a r 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 a m are possible. T h e benefits of level-constant thresholding seem to outweigh the negatives. Signals w i t h s m a l l number of points are better estimated by the R D W T m e t h o d then either D W T or W i e n e r filtering. Often for real G P R systems the number of points is under 2048, w h i c h increases the advantage of level-constant thresholding i n the R D W T d o m a i n . T h i s implies the level-constant R D W T threshold w o u l d be a candidate for the removal of G a u s s i a n noise, white or colored of G P R A-scans. 5.3 Simulation of Clutter Suppression in B-scans T h e following experimental results compare the baseline m e t h o d of average trace subtraction ( A T S ) to a the proposed m e t h o d of thresholding i n the curvelet d o m a i n . T h e methods are compared on synthetic d a t a to establish an understanding of the proposed m e t h o d and l i m i t a t i o n s . R e a l d a t a is then tested to see the visual effects of the proposed m e t h o d . T h e results of the synthetic tests are s u m m a r i z e d by Table 5.2. T h i s experiment compares the effectiveness of A T S to the proposed m e t h o d using b o t h 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 m e t h o d removes more clutter energy then the A T S does. T h e r e m a i n i n g experiments focus on v i s u a l quality. T h e o u t p u t of each a l g o r i t h m and the residuals for the synthetic d a t a are compared. synthetic test, a test on real l a n d m i n e d a t a is presented. 74 F o l l o w i n g the Original Average Trace S u b t r a c t i o n Proposed P o i n t Target 17.6 40.2 Large Target 17.6 21.8 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 m e t h o d for b o t h targets. 5.3.1 Synthetic B-scan Experiments T h e v i s u a l synthetic experiments are conducted to establish an understanding of the m e t h o d . T h e first experiment explores the effect of choosing a noise m o d e l . T h e average trace s u b t r a c t i o n m e t h o d assumes a noise m o d e l 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 m e t h o d also using an average trace ( A T ) model. A s a comparison, the edge noise m o d e l is also presented and the results are shown for the proposed m e t h o d . T h e experiment uses a large target, w h i c h penalizes the use of an average noise m o d e l . F i g u r e 5.8(b) shows t h a t the target energy has been damaged, and add i t i o n a l artifacts have been i n t r o d u c e d . T h e average noise m o d e l is not a good choice for the proposed m e t h o d either. Target energy is also damaged by the thresholding, but no a d d i t i o n a l artifacts are i n t r o d u c e d . T h e next experiment i n F i g u r e 5.9 uses a point target a n d 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 w h i c h varies i n a m p l i t u d e or phase. However, thresholding is not subjected to t h a t l i m i t a t i o n a n d effectively removes all clutter events. The residual of the proposed m e t h o d contains only artifacts at the sharp edges of the clutter. T h e ground, h a v i n g constant a m p l i t u d e , is easily removed by b o t h methods. T h e final synthetic experiment compares the large target to b o t h 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 F i g u r e 5.8: I n 5.8(a) the original d a t a w i t h clutter is presented, but the result s h o u l d be just the target as i n 5.8(b). T h e proposed m e t h o d is calculated w i t h an average noise m o d e l , 5.8(d) and edge noise m o d e l , 5.8(f). T o understand the output of the proposed m e t h o d , the noise m o d e l used for the threshold estimation is shown i n 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 F i g u r e 5.9: I n 5.9(d) the proposed m e t h o d leaves v i r t u a l l y no clutter energy and the target is still present. A T S leaves tremendous artifacts, but does manage to remove the ground clutter. 77 T h i s provides results for c o m m o n subsurface targets, w h i c h are larger t h a n point targets. W i t h these c o m m o n targets, using A T S shows that significant artifacts are introduced, seen i n F i g u r e 5.10(e). T h e proposed m e t h o d i n F i g u r e 5.10(c) demonstrates similar results for the large target as the point target. T h e m e t h o d provides near complete removal of a l l clutter and very little change to the target energy. 5.3.2 Landmine B-scan Experiments T h i s experiment uses l a n d m i n e d a t a from the J R C l a n d m i n e signatures database. T h e B - s c a n under s t u d y represents a physical location at 25cm from a survey that done on 50cm X 50cm area w i t h 1cm resolution. T h e l a n d m i n e is a P M N - t y p e anti-personnel device b u r i e d at 5cms i n loamy soil. T h e d a t a can be obtained from [11] and the system is described i n [59]. The results of clutter suppression are shown i n F i g u r e 5.11, w h i c h use the edge noise m o d e l shown i n 5.11(b) to create a threshold operator. F i g u r e 5.11(b) is the edge noise m o d e l calculated from the d a t a d o m a i n . A n y event that is strongly horizontal and exists near the edge w i l l be shown i n the noise m o d e l . T h e result from the proposed m e t h o d w i l l be to remove any events that are a p p r o x i m a t e d by the calculated noise m o d e l . T h e result is F i g u r e 5.11(c), w h i c h indeed shows superior clutter suppression to the average trace subtraction m e t h o d . A s demonstrated before, targets that are large suffer w h e n A T S is used to suppress clutter. T h e r e are artifacts i n t r o d u c e d around the target as a result of the average noise m o d e l . T h i s does not occur w i t h the proposed m e t h o d . T h e large target remains intact and becomes the prominent feature i n the result. A l s o note the g r o u n d clutter is not completely removed by either m e t h o d . 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 F i g u r e 5.10: T h e large target does not effect the ability of the proposed m e t h o d to remove clutter. Average trace s u b t r a c t i o n removes a p o r t i o n of target energy, adds artifacts near the target a n d does not remove the a m p l i t u d e v a r y i n g clutter. 79 Source Noise M o d e l (a) Input Data (b) Edge Noise model Adaptive Subtraction Average Trace S u b t r a c t i o n (c) Proposed Method (d) Average Trace Subtraction F i g u r e 5.11: T h e clutter is suppressed i n the proposed method, whereas significant artifacts are i n t r o d u c e d i n the A T S m e t h o d . T h e noise model i n 5.11(b) shows what likely events w i l l be suppressed i n 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 solutions are sought. The computational complexity and simplicity of A T S is low, which has led to the widespread use of A T S as an initial signal processing step. Also, A T S 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 A T S 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 A T S method. The point target, in Figure 5.9 is well fit to A T S 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 A T S method 81 struggles w i t h this, because the subtraction operator is b l i n d to any changes i n data. T h e proposed m e t h o d uses a threshold operator, w h i c h provides a quasi-adaptive m e t h o d . T h i s prevents the proposed m e t h o d from a d d i n g artifacts even w h e n the estimation of the noise m o d e l is w r o n g . T h e drawback of using the proposed m e t h o d is t h a t the scaling parameter A must be e x p e r t l y controlled to balance the weighting of the noise m o d e l . The results presented here used a weight of 2.8, w h i c h corresponds to thresholding any coefficients that exceed 3 times the noise m o d e l coefficients. It was found i n these experiments that the scaling parameter o n l y required rough adjustments and was easy to understand and control. aggressive setting was used. In the case of the l a n d m i n e data, a much less Sufficient clutter suppression was obtained by o n l y scaling the noise m o d e l 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 evaluation 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 b u i l t on the success of previous work i n wavelet thresholding. P r e v i o u s achievements presented results on global thresholding the R D W T to remove w h i t e noise. These results have extended the use of R D W T to perform nonparametric signal regression using very little a priori information. T h e threshold is determined by the data, b u t the scale parameter must be estimated by other means. T h i s work has also presented a reasonable physical m o d e l for a G P R A scan. U s i n g this m o d e l and other p r a c t i c a l 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 d i r e c t l y to current G P R data. In the realm of clutter suppression, there are three secondary contributions that stemmed from evaluating the proposed m e t h o d of thresholding an edge noise m o d e l i n the curvelet d o m a i n : â€¢ C l u t t e r suppression can be performed by thresholding i n the curvelet d o m a i n . â€¢ A n a p p r o x i m a t i o n of clutter using a simple model, such as the edge noise m o d e l 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 i n the wavelet d o m a i n to re-position the p r o b l e m of clutter suppression as threshold estimation problem. Because of the phase and a m p l i t u d e resilience of thresholding effective clutter suppression can be achieved by o n l y estimating an approximate m o d e l 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 suppression 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 combination 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 preserved 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 m o d e l and suppress the clutter feature using thresholding. 86 Bibliography [1] G . B e y l k i n . O n the Representation of Operators i n Bases of C o m p a c t l y Supp o r t e d Wavelets. Siam Journal [2] I l a r i a B o t t i g l i e r o . tion? 120 Million on Numerical Landmines Analysis, Deployed 29(6):1716-1740, 1992. Worldwide: Fact or Fic- L e o C o o p e r , P e n & S w o r d B o o k s , Barnsley, S o u t h Y o r k s h i r e , E n g l a n d , 2000. [3] J . W . B r o o k s . A p p l i c a t i o n s of G P R Technology to H u m a n i t a r i a n D e m i n i n g Operations i n C a m b o d i a : Some Lessons L e a r n e d . In Internation on Technology and the Mine Problem, Symposium Monterey, C A , 1998. [4] J o h n W . B r o o k s , L u c van K e m p e n , and H i c h e m S a h l i . P r i m a r y s t u d y i n adaptive clutter reduction and b u r i e d minelike target enhancement from G P R data. In Detection and Remdiation volume 4038 of Proceedings Technologies SPIE, for Mines and Minelike Targets V, pages 1183-1192, O r l a n d o , F A U S A , 2000. SPIE. [5] C . B r u s c h i n i a n d B . G r o s . A Survey of C u r r e n t Technology Research for the Detection of L a n d m i n e s . Sustainable niques and Technologies, Humanitarian Demining: Trends, Tech- pages 172-187, 1998. [6] G . J . B u r t o n and G . P . O h l k e . E x p l o i t a t i o n of M i l l i m e t e r Waves for T h r o u g h W a l l Surveillance d u r i n g M i l i t a r y O p e r a t i o n s i n U r b a n T e r r a i n . Technical report, R o y a l M i l i t a r y College of C a n a d a , M a y 24, 2000 2000. [7] E . J . Candes a n d D . L . D o n o h o , editors. Nonadaptive Representation Curvelets - A suprisingly Effective for Objects with Edges. Curves and Surfaces. V a n - derbilt U n i v e r s i t y Press, Nashville, T N , 1999. [8] E . J . Candes a n d D . L . D o n o h o . Curvelets a n d curvilinear integrals. of Approximation Theory, 113(l):59-90, 2001. 87 Journal [9] E . J . C a n d e s a n d D . L . D o n o h o . N e w tight frames of curvelets a n d o p t i m a l representations of objects w i t h piecewise C - 2 singularities. Communications Pure and Applied Mathematics, on 57(2):219-266, 2004. [10] E m m a n u e l Candes, L a u r e n t Demanet, D a v i d D o n o h o , a n d L e x i n g Y i n g . C u r v e L a b 1.0, 2005. A v a i l a b l e from: h t t p : / / w w w . c u r v e l e t . o r g / . [11] Joint Research Centre. A n t i - p e r s o n n e l L a n d m i n e Signatures Database. Avail- able from: h t t p : / / a p l - d a t a b a s e . j r c . i t / . [12] M . C h e r n i a k o v a n d L . D o n s k o i . buried object detection. Frequency b a n d selection of radars for Ieee Transactions on Geoscience and Remote Sens- ing, 37(2):838-845, 1999. [13] R . R . C o i f m a n a n d D . L . D o n o h o . Wavelets and Statistics, T r a n s l a t i o n Invariant D e - N o i s i n g . In volume L e c t u r e Notes i n Statistics, pages 125-150. Springer-Verlag, N e w Y o r k , 1995. [14] R . R . C o i f m a n a n d M . V . Wickerhauser. E n t r o p y - B a s e d A l g o r i t h m s for Best Basis Selection. Ieee Transactions on Information [15] Lawrence B . Conyers a n d D e a n G o o d m a n . introduction for archaeologists. Theory, 38(2):713-718, 1992. Ground-penetrating radar : an A l t a M i r a Press, W a l n u t Creek, C A , 1997. [16] J a m i n C r i s t a l l , M o r i t z B e y r e u t h e r , a n d Felix H e r r m a n n . C u r v e l e t processing and i m a g i n g : 4 D adaptive s u b t r a c t i o n . I n CSEG National Convention. 2004. [17] D . J . Daniels, D . J . G u n t o n , a n d H . F . Scott. I n t r o d u c t i o n to Subsurface R a d a r . lee Proceedings-F Radar and Signal Processing, 135(4):278-320, 1988. [18] D . J . Daniels a n d I n s t i t u t i o n of E l e c t r i c a l Engineers. Surface-penetrating radar. I E E radar, sonar, navigation a n d avionics series ; 6. I n s t i t u t i o n of E l e c t r i c a l Engineers, L o n d o n , 1996. [19] D . J . Daniels a n d I n s t i t u t i o n of E l e c t r i c a l Engineers. Ground penetrating radar. I E E radar, sonar, navigation, a n d avionics series ; 15. I n s t i t u t i o n of E l e c t r i c a l Engineers, L o n d o n , 2 n d edition, 2004. [20] C a n a d i a n Forces Department Landmine of Database, National Dec 18, Defense. 2003. Canadian Available from: http://ndmic-cidnm.forces.gc.ca/index.asp?lang=e. [21] D . L . D o n o h o . D e - N o i s i n g b y Soft-Thresholding. Ieee Transactions mation Theory, 41(3):613-627, 1995. 88 on Infor- [22] D . L . D o n o h o and M . E l a d . O p t i m a l l y sparse representation (nonorthogonal) dictionaries v i a 1(1) m i n i m i z a t i o n . Proceedings Academy of Sciences of the United States of America, i n general of the National 100(5):2197-2202, 2003. [23] D . L . D o n o h o and I. M . Johnstone. Ideal S p a t i a l A d a p t a t i o n by Wavelet S h r i n k age. Biometrika, 81(3):425-455, 1994. [24] D a v i d D o n o h o , M a r k R e y n o l d D u n c a n , X i a o m i n g H u o , and Ofer L e v i . W a v e l a b 802. A v a i l a b l e from: h t t p : / / w w w - s t a t . s t a n f o r d . e d u / wavelab/. [25] M . ; Y u n x i n Zhao Gader, P^D.; M y s t k o w s k i . L a n d m i n e detection w i t h g r o u n d penetrating radar using h i d d e n M a r k o v models. Geoscience ing, IEEE Transactions and Remote Sens- on, 39(6):1231-1244, 2001. [26] P. D . Gader, B . N . Nelson, H . F r i g u i , G . V a i l l e t t e , and J . M . K e l l e r . logic detection of landmines w i t h ground penetrating radar. Signal Fuzzy Processing, 80(6):1069-1084, 2000. [27] J . G a z d a g . W a v e - E q u a t i o n M i g r a t i o n w i t h Phase-Shift M e t h o d . Geophysics, 43(7):1342-1351, 1978. [28] Rafael C . Gonzalez and R i c h a r d E . W o o d s . Digital image processing. Prentice H a l l , U p p e r Saddle R i v e r , N . J . , 2nd edition, 2002. [29] R . G r i b o n v a l and M . Nielsen. O n the strong uniqueness of h i g h l y sparse representations from redundant dictionaries. Independent Blind Signal Separation, Component Analysis and 3195:201-208, 2004. [30] A . G r o s s m a n n and J . M o r l e t . D e c o m p o s i t i o n of H a r d y F u n c t i o n s into Square Integrable Wavelets of C o n s t a n t Shape. Siam Journal on Mathematical Anal- ysis, 15(4):723-736, 1984. [31] D i g i t a l Signal Processing G r o u p . R i c e Wavelet T o o l b o x v2.4, 2001. A v a i l a b l e 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] F r d r i c Guerne, B e r t r a n d G r o s , M a r c Schreiber, and J e a n - D a n i e l N i c o u d . G P R M i n e Sensors for D a t a A c q u i s i t i o n i n the F i e l d . I n International Sustainable Humanitarian Demining, Workshop on Zagreb, C r o a t i a , 1997. [33] A . G u i t t o n and D . J . Verschuur. A d a p t i v e s u b t r a c t i o n of multiples using the L - l - n o r m . Geophysical Prospecting, 52(l):27-38, 2004. [34] Felix J . H e r r m a n n . C u r v e l e t i m a g i n g and processing: an overview. In National Convention. 2004. 89 CSEG [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. L u n , 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 p u l s e E K K O 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 groundpenetrating 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 Image, Orlando, F L , 1995. SPIE. [45] M . Lang, H . Guo, J . E . Odegard, C . S. Burrus, and R . O. Wells. Noise reduction 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 a n d S. Z h o n g . Pattern Analysis C h a r a c t e r i z a t i o n of. signals from multiscale edges. 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. A c a d e m i c Press, S a n Diego, 2nd edition, 1999. [49] S. G . M a l l a t a n d Z . F . Z h a n g . M a t c h i n g P u r s u i t s w i t h T i m e - F r e q u e n c y D i c tionaries. Ieee Transactions [50] A d a m on Signal Processing, Marcus. Landmines, Inventions May 6, 41(12):3397-3415, 1993. to 1998 Lift the 1998. Scourge Available of from: http://www.csmonitor.com/durable/1998/05/06/f-pl4sl.htm. [51] R o h a n M a x w e l l . C a m b o d i a : A C o u n t r y Profile. Journal of Mine Action, 5(1), 2001. [52] M . R . M c C l u r e a n d L . C a r i n . M a t c h i n g pursuits w i t h a wave-based d i c t i o n a r y . Ieee Transactions on Signal Processing, 45(12) :2912-2927, 1997. [53] I. L . M o r r o w a n d P . v a n Genderen. Effective i m a g i n g of b u r i e d dielectric objects. Ieee Transactions on Geoscience and Remote Sensing, 40(4):943-949, 2002. [54] G . P . N a s o n a n d B . W . S i l v e r m a n . T h e S t a t i o n a r y Wavelet Transform a n d some S t a t i s t i c a l A p p l i c a t i o n s . I n A A n t o n i a d i s a n d G . O p p e n h e i m , editors, Wavelets and Statistics, volume L e c t u r e Notes i n Statistics, pages 281-299. Springer-Verlag, N e w Y o r k , 1995. [55] L . N u z z o a n d T . Q u a r t a . Improvement i n G P R coherent noise attenuation using tau-p a n d wavelet transforms. Geophysics, 69(3):789-802, 2004. [56] J . C . Pesquet, H . K r i m , a n d H . C a r f a n t a n . T i m e - i n v a r i a n t o r t h o n o r m a l wavelet representations. Ieee Transactions [57] A l e k s a n d r a P i z u r i c a . Modeling. on Signal Processing, Image Denoising Using 44(8):1964-1970, 1996. Wavelets and Spatial Context P h . d , U n i v e r i s t y of Gent, 2002. [58] Jane M a r y A n a s t a s i a R e a . Ground penetrating radar applications in hydrology. U n i v e r s i t y of B r i t i s h C o l u m b i a , Vancouver, 1997. [59] B . Scheers, M . Acheroy, a n d A . V a n d e r V o r s t . T i m e - d o m a i n s i m u l a t i o n a n d characterisation of T E M horns using a normalised impulse response. Proceedings-Microwaves Antennas and Propagation, 91 147(6) :463-468, 2000. lee [60] B a r t Scheers. Ultra- Wideband the Detection Ground Penetrating of Anti-Personnel Landmines. Radar, With Application to P h . d , Royal M i l i t a r y Academy, 2001. [61] U S D e p a r t m e n t of State. H i d d e n K i l l e r s : T h e G l o b a l L a n d m i n e C r i s i s . Deparment of State P u b l i c a t i o n 10575, B u r e a u of P o l i t i c a l - M i l i t a r y Affairs, Office of H u m a n i t a r i a n D e m i n i n g P r o g a m s , September 1998 1998. [62] James D . T a y l o r . Ultra-wideband FL, radar technology. C R C Press, B o c a R a t o n , 2001. [63] S. T j o r a , E . E i d e , a n d L . L u n d h e i m . E v a l u a t i o n of methods for ground bounce removal i n G P R u t i l i t y m a p p i n g . I n Ground Penetrating 2004- Conference Proceedings of the Tenth International Radar. 2004- GPR on, volume 1, pages 379-382, 2004. [64] J . A . T r o p p . G r e e d is good: A l g o r i t h m i c results for sparse a p p r o x i m a t i o n . Ieee Transactions on Information Theory, 50(10):2231-2242, 2004. [65] T . J . U l r y c h , M . D . Sacchi, a n d J . M . G r a u l . Signal a n d 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 . W o o d b u r y , a n d R . D . Sacchi. L - m o m e n t s a n d C moments. Stochastic 68, Environmental Research and Risk Assessment, 14(1):50- 2000. [67] L . V a n K e m p e n , H . Sahli, E . Nyssen, a n d J . Cornells. Signal processing a n d pattern recognition methods for radar A P mine detection a n d identification. I n 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, a n d G. Maksymonko. A d a p t i v e ground bounce removal. Electronics Letters, 37(20):1250-1252, 2001. [69] X u X i a o y i n , E . L . Miller, C M . R a p p a p o r t , a n d G . D . Sower. S t a t i s t i c a l m e t h o d to detect subsurface objects using array ground-penetrating radar data. Geoscience and Remote Sensing, IEEE Transactions on, 40(4):963-976, 2002. [70] Y . X . Zhao, P . G a d e r , P . C h e n , a n d Y . Z h a n g . T r a i n i n g D H M M s of m i n e a n d clutter to m i n i m i z e l a n d m i n e detection errors. Ieee Transactions and Remote Sensing, 41(5):1016-1024, 2003. 92 on Geoscience 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 P l o t t i n g Function Due to the s u b j e c t i v i t y of p l o t t i n g intensity images, a detailed description of the p l o t t i n g function w h i c h was used to generate B - s c a n plots is described. The im- portant features are the color scaling, image scaling and image translation. function used to re-scale and translate is called r e l _ i m , shortened from, The relative image. T h e purpose was to allow results to be plotted at the same intensity scale as the i n p u t . M A T L A B includes a image p l o t t i n g function called image, w h i c h interprets an entry i n a m a t r i x as an index to a colormap: T h e colormap used i n this thesis for p l o t t i n g B-scans is white at the m i d d l e index of the colormap m a t r i x and linearly scales i n 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 m a p p i n g makes no d i s t i n c t i o n between negative and positive events i f the image array is translated such that the zero indexes the m i d d l e index of the colormap. T h i s is the purpose of the r e l _ i m function, w h i c h scales and translates an i n p u t image array relative to a second i n p u t : 7, Image something r e l a t i v e 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 c o m m o n scale relative to a p a r t i c u l a r source image. T h i s thesis uses the i n p u t image as the base scale and all other images are scaled accordingly. If the i n p u t 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 a r b i t r a r y scale. N o t e there are o n l y 16 levels of gray to choose from, between p r i n t i n g and the d y n a m i c range of the color scale, some s m a l l a m p l i t u d e information w i l l be lost. T h e 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
Page Metadata
Item Metadata
Title | Threshold estimation using wavelets and curvelets on ground penetrating radar data for noise and clutter suppresion |
Creator |
Harrison, Dustin |
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 |
Graduation Date | 2005-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | 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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065408/manifest