Archaeological Prospection Archaeol. Prospect. 14, 130–141 (2007) Published online 26 February 2007 in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/arp.289 Wavelet Transform in Denoising Magnetic Archaeological Prospecting Data B.TSIVOURAKI AND G. N.TSOKAS* Laboratory of Exploration Geophysics, Department of Geophysics, School of Geology, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece ABSTRACT Magnetic measurementsin archaeological prospecting are often affected by culturalnoise having the same high-frequency content as anomalies arising from buried antiquities. Also, in many cases the microrelief of the ground surface causes a noise that is coherent, pseudorandom and periodic. The main cause ofthiskind ofnoise isthe ploughing ofthe earth.Inthispaper the efficiencyofawavelet denoising scheme is tested with respect to these types of unwanted disturbances. The proposed scheme combines the cyclospinning algorithm with a variable threshold calculated in each cycle of thealgorithm.Testsonsyntheticandrealdata showa satisfactoryperformance ofthetechniquein suppressing both the white noise and the coherent noise caused by the systematic undulations of the ground surface. Copyright # 2007 John Wiley & Sons, Ltd. Key words: magnetic; noise; filtering; wavelet; denoising; fourier Introduction In archaeological prospecting, total magnetic field or gradiometer data are often contaminated with noise and artefacts giving rise to so-called ‘cultural noise’. In general, the presence of cultural noise in data distorts the signal coming from shallow depth anomalies, resulting in poor quality of any subsequent processing. Therefore, the first step of any processing of such data is the suppression of the noise. The cultural noise is randomly distributed both spatially and in its characteristic frequencies depending on the source. Usually this kind of noise has a highfrequency content similar to that corresponding with anomalies that arise from magnetic sources * Correspondence to: G. N. Tsokas, Laboratory of Exploration Geophysics, Department of Geophysics, School of Geology, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece. E-mail: gtsokas@geo.auth.gr Copyright # 2007 John Wiley & Sons, Ltd. of archaeological interest. As a result, denoising by conventional filtering is difficult without affecting the signal’s sharp variations. Conversely, denoising in the wavelet domain by shrinking the noisy wavelet coefficients achieves significant noise suppression while preserving the signals sharp variations. Also, in many cases the microrelief of the ground surface causes a noise that is coherent. This is readily seen because a regular pattern is superimposed on the resulting images. Linear shallow ditches due to ploughing are the most common source of this kind of noise. On one hand, their effect confuses the final images and thus the anomalies do not look like a plan view of the concealed ruins. On the other hand, they make the interpretation difficult because they can be mistaken for features concealed in the subsurface. Coherent noise is pseudorandom, with mostly low-frequency content and in some cases periReceived 7 December 2005 Accepted 31 August 2006 Denoising Prospecting Data odic. Although, the separation of these features from the signal needs a special treatment, the denoising technique in the wavelet domain achieves a significant suppression of the surface microrelief noise. In this paper the discrete wavelet transformation is used for denoising total magnetic field and gradiometer data. The principles of denoising in wavelet domain are briefly stated. Synthetic and real examples are used to demonstrate the efficiency of the technique. The results show a satisfactory performance of the technique in suppressing both white noise and coherent noise caused by the systematic undulations of the ground surface. Wavelets The wavelet transformation is a practical signal analysis tool having applications in a variety of scientific and engineering areas including geophysics (Moreau et al., 1997; Fedi and Quarta, 1998; Horby et al., 1999; Ridsdill-Smith and Dentith, 1999; Billings and Herrmann, 2003; Saracco et al., 2003). It can be interpreted as a 131 decomposition of the signal in terms of wavelets that play the role of basis functions. Unlike the Fourier transformation where the basis functions are sinusoids, the wavelets are oscillatory, well localized and rapidly decaying functions and therefore suitable for analysing signals with spatially localized features. Figure 1 illustrates the oscillation of the Morlet wavelet. The wavelet basis is a family of functions constructed through dilations and translations of the analysing ‘mother’ wavelet 1 xb cab ðxÞ ¼ pﬃﬃﬃ c b 2 R; a > 0 a a where a is the scale factor and b is the translation. The wavelet coefficients of a signal f are the numbers Z wða; bÞ ¼ < f; cða; bÞ > ¼ fðxÞca;b ðxÞdx where ca;b (x) is the complex conjugate of ca,b(x). The computation of the transform is implemented for a predefined set of scales. At each scale a, the signal is compared with the scaled wavelet which can be positioned at different offsets according to the value of parameter b. By Figure 1. The Morlet wavelet. Copyright # 2007 John Wiley & Sons, Ltd. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp B. Tsivouraki and G. N. Tsokas 132 varying b the procedure may be repeated for all possible offsets of the scaled wavelet relative to the signal. Thus the wavelet transformation is a representation of the signal in the space-scale domain. Comparing the Fourier and the wavelet transforms we notice that the Fourier transform gives the maximum information in the frequency domain but completely loses positional information. The wavelet transformation decomposes a signal at a number of scales (frequency ranges) keeping the positional information. Figures. 2 and 3 illustrate the differences between the Fourier transformation of a signal with local features and its wavelet transformation. Denoising in wavelet domain Noise reduction may be performed in several ways. The most common and therefore classic denoising methods are local polynomial regression, adaptive filtering, optimal Wiener filtering and moving average filters. In the past decade the wavelet transform has also been used (Donoho and Johnstone, 1994) for cleaning noisy data by shrinking noisy wavelet coefficients via thresholding. This method is based mainly on two properties of the wavelets: they concentrate the energy of a smooth signal in a few wavelet coefficients, and the transformation of random noise is also random noise. Therefore it is reasonable to assume that small coefficients represent the noise and can be set to zero, whereas the large ones contain the signal’s energy and should be kept. The work of Donoho and Johnstone (1994) comprises the basis of a class of denoising methods in the wavelet domain. In general, the wavelet domain denoising method has the advantage that it preserves the signal’s characteristics. The performance of the method is affected by the choice of wavelet, the thresholding function and the threshold selection rule. There is no unique denoising policy covering all data types. Each data type needs an empirical selection of the proper combination of the above-mentioned parameters. The present study is based on a wavelet denoising scheme that combines the shift invariant cycle-spinning algorithm (Coifman and Donoho, 1995) and a variable threshold calculated Amplitude 1 0.5 0 -0.5 -1 0 50 100 150 200 discrete space (m ) 250 300 3 3.5 80 Amplitude 60 40 20 0 0 0.5 1 1.5 2 2.5 discrete angular frequency (radians) Figure 2. The upper part ofthe figure shows a sinusoidal signalwith two local frequencies (one follows the other in space), and the lower part depicts its Fourier transform. Note that although the Fourier transform fully characterizes the frequencies present in the data, the positions at which they occur cannot be readily inferred. The question: where does each sinusoid begin cannot be answered using this representation. Copyright # 2007 John Wiley & Sons, Ltd. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp Denoising Prospecting Data 133 Figure 3. The upper part shows the same sinusoidal signal as in Figure 2 and the lower part depicts its wavelet transformation. It is obvious that wavelets extract the same information as Fourier transform plus localization, i.e. the position at which the sinusoids occur. in each cycle of the algorithm (TsivourakiPapafotiou, 2004). For a number of circular shifts of the signal, the cycle-spinning algorithm performs the following actions: (i) denoises the shifted data; (ii) unshifts the denoised data; (iii) averages the results obtained. This method produces a reconstruction of the signal subject to much weaker Gibbs’ phenomena. The selection of the threshold is based on an estimation of the noise level in the data and the wavelet transformation of random noise with the same standard deviation. The workflow consists of the following steps: (i) generation of white Gaussian noise with same standard deviation as the data; (ii) wavelet transformation of the noise with the Coiflet1 wavelet—evaluation of the mean value m and the standard deviation s of the noise coefficients in the finest level of the transformation; Copyright # 2007 John Wiley & Sons, Ltd. (iii) select as threshold the value (m þ 2.5s). As the noise is random, repetition of the threshold selection procedure at each positional shift of the signal gives a slightly different threshold in each case. The results show a significant suppression of the noise and very good recovery of the original signal. Examples Tests on synthetic data Figure 4 shows the total intensity magnetic anomaly produced by two adjacent prismatic bodies with their upper surface buried at 1 m depth. This signal was corrupted with normal white noise with standard deviation s ¼ 2.85 giving a signal to noise ratio of 11. This kind of noise is random and Gaussian, having a power spectral density that is flat over a wide range of frequencies. The lower part of Figure 4 shows the data after application of the wavelet denoising Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp B. Tsivouraki and G. N. Tsokas 134 40 nT 20 0 -20 -40 0 10 20 30 0 10 20 30 40 50 60 70 40 50 60 70 40 nT 20 0 -20 -40 x(m) Figure 4. Denoising of synthetic magnetic data contaminated with white noise having standard deviation 2.85. In the upper part of the figure the solid line depicts the original signal and the dots depict the noisy signal. The lower part shows the coincidence of the initial signal (solid line) with the recovered one (dotted line). algorithm. Denoising was also performed using optimal Wiener filtering (although this is not shown in the figure). The quality of the signal estimation using the wavelet method is demonstrated by the high signal to noise ratio (SNR) of the recovered anomaly, as presented in Table 1. It is evident from Table 1 that the proposed scheme produces a very good recovery of the signal. The Wiener filtering also shows its optimal performance because the exact signal’s model has been used. Figure 4 shows the coincidence of the original signal and its estimation computed by the proposed scheme. As a two-dimensional example we present the denoising of the anomaly caused by a buried rectangular prismatic body with its upper surface buried at the depth of 1 grid unit (in the particular case the grid unit or sampling interval is set to 1 m). This signal was corrupted with normal white noise with standard deviation s ¼ 2. The SNR in the original synthetic data is 27 which became 29 for the wavelet domain denoised data. This is shown in Figure 5 where the original noisy data are shown in the grey scale map of the left-hand side and the denoised version of the same data in the wavelet domain are presented on the right-hand side. Total magnetic field data from the exploration of the archaeological site of Europos in northern Greece Europos was a commercial centre on the banks of the River Axios in northern Greece (Region of Macedonia). The ruins of the ancient urban centre and associated features are buried in the Table 1. Signal to noise ratio (SNR) of recovered signal Method SNR of the recovered signal MaximumValue Minimum value Wiener 12.54 Proposed algorithm 15.38 Original signal: 24.43 Recovered signal: 22.41 Original signal: 24.43 Recovered signal: 21.13 Original signal: 29 Recovered signal: 21.15 Original signal: 29 Recovered signal: 28.45 Copyright # 2007 John Wiley & Sons, Ltd. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp Denoising Prospecting Data 135 nT 80 160 Noisy 2d data SNR=27.2 0 -80 Recovered 2d signal SNR=29.04 60 50 40 30 20 10 10 20 30 40 50 60 0 10 20 30 40 50 60 5 10 m Figure 5. The noisy signal is shown on the left part and the recovered one after wavelet denoising is on the right. The initial signal was contaminated by white noise having standard deviation s ¼ 2. Evidently, the attempted suppression of the noise was performed efficiently. subsurface, near a modern village that bears the same name. The original data are shown in the left part of Figure 6. They comprise a 20 m by 20 m grid of magnetic total field data, extracted from the mesh of grids established on the site for the acquisition of geophysical data (Tsokas et al., 1994). Readings were collected at 1-m intervals along parallel profiles spaced 1 m apart. Part of the ruins of the ancient urban complex was concealed in the subsurface under this particular area of land. These ruins of foundation walls would have shown up better if a greyscale or dot-density image had been constructed for these data. Nevertheless, we chose the presentation in the form of a contour map because it demonstrates better the effect of denoising of the magnetic field. The noise was modelled as white Gaussian with standard deviation s ¼ 2. The application of the wavelet denoising scheme resulted in the distribution of the Earth’s total magnetic field shown in the right-hand side of Figure 6. For a Copyright # 2007 John Wiley & Sons, Ltd. better illustration of the method’s performance a single profile was extracted from the grid and plotted as a solid line in Figure 7. The denoised readings along the same profile are depicted by stars in the same figure. It is obvious that the wavelet method gives satisfactory results without affecting the anomalies that carry the information about the archaeological features in the subsurface. The noise artefacts have disappeared resulting in a smooth version of the field. Data from the magnetic survey of Makrygialos in northern Greece (Pieria) As a second field example we use part of the data collected from the total magnetic field survey of the archaeological site of Makrygialos in northern Greece (Region of Macedonia). The results of the survey have been reported by Tsokas et al. (1997). The site contains the remnants of a neolithic settlement which was threatened by the Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp B. Tsivouraki and G. N. Tsokas 136 Figure 6. The raw magnetic total field data over a single 20 m by 20 m grid are shown in the contour map on the left-hand side. The particular grid is part of the mesh of grids established on the ground surface for the exploration of the site, which hosts the ruins of ancient Europos in the Region of Macedonia in northern Greece (Tsokas et al., 1994). The grid is selected for demonstration of the proposed denoising scheme. Evidently the distribution of the total field suffers from high-frequency noise, but a lot of information about the subsurface situation is also present. The denoised version of the data is shown on the right, where the enhancement of the image is evident.The contours of both maps are in nT. 65 60 55 50 nT 45 40 35 30 25 20 0 2.5 5 7.5 10 x (m ) 12.5 15 17.5 20 Figure 7. The solid line depicts the total magnetic field variation along the first profile of the grid of Figure 6.The denoised data are depicted by stars. Copyright # 2007 John Wiley & Sons, Ltd. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp Denoising Prospecting Data construction of a modern motorway. Again, we used a small fraction of the whole data set consisting of a 40 m by 40 m grid. The data acquisition parameters are the same as in the case of the example of Europos. The noise was modelled as white with s ¼ 2. The map of the original data set is shown in the upper left part of Figure 8. Contour map presentation was preferred for the same reason as in the previous example. The anomalies 137 associated with two parallel ancient ditches are clearly seen. These data were denoised with the use of the proposed scheme and the results also are shown in Figure 8. Note that the weak anomalies running perpendicular to the response of the ditches have been eliminated. These anomalies are due to small undulations of the surface relief caused by ploughing. Furthermore the distribution of the field has been clearly denoised. The data along the fifth Figure 8. Raw magnetic total field data are shown in the map on the upper part. They were collected for the exploration of the archaeological site of Makrygialos in northern Greece which contains the remnants of a neolithic settlement (Tsokas et al.,1997). The magnetic signatures of two parallel ditches dominate the map. Also, high-frequency noise and microrelief-induced noise are present.The denoised version of the data is shown on the lower right-hand side. Contours of both maps are in nT. Copyright # 2007 John Wiley & Sons, Ltd. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp 138 B. Tsivouraki and G. N. Tsokas Figure 9. Plot of a single profile (fifth row) of the data shown in Figure 8. The solid line and the stars depict the raw and the denoised data respectively. profile of the grid are shown in Figure 9, both in the original and denoised form. The suppression of small noise artefacts is obvious. Gradiometer data from the prehistoric site Avgi in northwest Greece The wavelet method was also used for denoising gradiometric data from the prehistoric site Avgi in northwest Greece (Region of Macedonia). The data were subjected to a sequence of processing functions which are stated at the top of Figure 10. The outcome is shown by the image in the left part of the figure. The remnants of the foundation walls of an ancient building are clearly seen in the picture. However, ploughing of the field has created a sequence of parallel linear anomalies which distort the otherwise clear rectangular image of the target. Denoising by wavelets resulted in the image shown in the right-hand side of the figure. In comparison to conventional processing, the latter image is much enhanced. Random noise has been removed and the effect of ploughing has been significantly reduced. The effect of a ploughing track is the only visible feature of the image shown in Figure 11 (left). The particular grid was also one of those surveyed by magnetic gradiometry in Avgi. The Copyright # 2007 John Wiley & Sons, Ltd. denoised version of the distribution of the total field gradients (strictly speaking the quantity mapped is the first vertical difference of the vertical component of the field) demonstrates the filtering power of the wavelets (Figure 11 right). The image illustrates the utility of wavelet domain filtering since the effect of the ploughing track has been almost completely eliminated. However, it must be stressed that the proper selection of parameters was a matter of trial and error. Furthermore, the ploughing tracks have an approximate width of 30 cm and so the wave number content of the ploughing track anomaly was restricted to higher wave numbers. Discussion and conclusions Shift invariant denoising of magnetic data, using a variety of thresholds in the wavelet domain, tends to give an estimation of the original signal with a significant suppression of both white noise and the coherent noise caused by the systematic undulations of the ground surface. In general, the wavelet denoising method preserves the characteristics of the signal. For better performance, the method requires profiles that have a relative large dyadic (power of 2) Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp Denoising Prospecting Data 139 Figure 10. In the left part, gradiometer data are shown, which were collected at the prehistoric site of Avgi in northern Greece, denoised by conventional filtering. The processing sequence is written at the top of the figure. GEOPLOT v3.0 software (GEOSCAN RESEARCH) was used for this purpose and for the construction of the images. The data denoised in the wavelet domain are shown on the right-hand side.The noise was modelled as white Gaussian with s ¼ 2. Units are in nT/0.5 m. length. This can be achieved by symmetrical extension of the data on both sides. As a profile’s length increases, more multiresolution levels are available for thresholding. On the other hand, the fact that usually the third decomposition level is defined as the deepest one available for thresholding (Donoho and Johnstone, 1994; Kovac, 1998), requires a profile length of at least 32 measurements. Using wavelets in a denoising scheme comprises an alternative to various other filtering techniques. For instance, wave number filtering Copyright # 2007 John Wiley & Sons, Ltd. by means of the Fourier transform is a well established procedure which would yield similar results, both for suppressing random noise and high wave number coherent spurious effects. However, the advantage of wavelets lies in the fact that the spectrum calculated as an intermediate step provides both frequency and localization information. The method is relatively easy to implement and comprises a simple scheme in general. It can be used in any step of the conventional processing sequence (before any interpolation) or it can be applied directly to Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp B. Tsivouraki and G. N. Tsokas 140 nT/m -10 -8 -6 -4 -2 0 2 4 A256 15 10 10 5 5 0 0 5 0 10 5 8 10 A256 DENOISED 15 0 6 15 0 5 10 15 10 m Figure 11. Original field gradiometer data (left) from the prehistoric site Avgi in northwest Greece, and the wavelet domain denoised version of the same data (right). The particular grid (20 m by 20 m) is a small part of the explored area. The noise standard deviation was estimated as s ¼ 2. raw data. It does involve a high degree of subjectivity with respect to the selection of thresholds and the thresholding technique. Also, a suitable mother wavelet has to be chosen at the outset of the procedure. This is also a matter of experience and the level of training of the processor. Experience with archaeological data indicates that the Coiflet1 and Daubechies2 wavelets tend to give comparatively better results. Acknowledgements The authors are grateful to Professor D. Donoho for his freeware wavelet toolbox for Matlab. Special thanks to the editor Dr Chris Gaffney and the unknown reviewers for their valuable comments and suggestions. References Billings SD, Herrmann F. 2003. Automatic detection of position and depth of potential UXO using a Copyright # 2007 John Wiley & Sons, Ltd. continuous wavelet transform. In Proceedings of SPIE Vol. 5089, Detection and Remediation Technologies for Mines and Mine like Targets VII, Harmon RS, Holloway JH, Broach JT (eds). 1012– 1022. Coifman R, Donoho D. 1995. Translation Invariant De-noising. Yale University and Stanford University. Donoho D, Johnstone I. 1994. Ideal spatial adaptation via wavelet shrinkage. Biometrica 81: 425–455. Fedi M, Quarta T. 1998. Wavelet analysis for the regional-residual separation of potential field anomalies. Geophysical Prospecting 46: 507–525. Horby P, Boschetti F, Horowitz FG. 1999. Analysis of potential field data in the wavelet domain. Geophysical Journal International 137: 175–196. Kovac A. 1998. Wavelet thresholding for unequally spaced data. PhD. Thesis, University of Bristol. Moreau F, Gibert D, Holschneider M, Saracco GG. 1997. Wavelet analysis of potential fields. Inverse Problem 13: 165–178. Ridsdill-Smith TA, Dentith MC. 1999. The wavelet transform in aeromagnetic processing. Geophysics 64: 1003–1013. Saracco G, Moreau F, Mathe PE, Hermitte D. 2003. Inverse problem in archaeological magnetic Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp Denoising Prospecting Data surveys using complex wavelet transform. European Geophysical Society, Geophysical Research Abstracts 5: 13196. Tsivouraki-Papafotiou B. 2004. Geophysical explorations and geophysical signal processing. PhD thesis, University of Thessaloniki. Tsokas GN, Giannopoulos A, Tsourlos P, Vargemezis G, Tealby J, Sarris A, Papazachos CB, Savo- Copyright # 2007 John Wiley & Sons, Ltd. 141 poulou T. 1994. A large scale geophysical survey in the archaelogical site of Europos (northern Greece). Journal of Applied Geophysics 32: 85–98. Tsokas GN, Sarris A, Pappa M, Bessios M, Papazaxos CB, Tsourlos P, Giannopoulos A. 1997. A large scale magnetic survey in Makrygialos (Pieria), Greece. Archaeological Prospection 4(3): 123–128. Archaeol. Prospect. 14, 130–141 (2007) DOI: 10.1002/arp

1/--страниц