close

Вход

Забыли?

вход по аккаунту

?

Wavelet transform in denoising magnetic archaeological prospecting data.

код для вставкиСкачать
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Þ ¼ pffiffiffi 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
Документ
Категория
Без категории
Просмотров
0
Размер файла
566 Кб
Теги
data, prospective, wavelet, transform, denoising, magnetic, archaeological
1/--страниц
Пожаловаться на содержимое документа