Enhancing Geophysical Signals Through the Use of Savitzky-Golay filtering method

Resumen

En la cuenca de Sidi Chennane en Marruecos, existen inclusiones estériles de caliche que dificultan la extracción de las rocas fosfatadas y que son difíciles de detectar. La resistividad de los caliches excede un valor de 200 Ω-m contra 150 Ω-m para la roca fosfatada. Se realizó un trabajo de prospección eléctrica de tipo Schlumberger en una zona de 50 ha en la cuenca de fosfatos de Oulad Abdoun. Se obtuvieron modelos del perfil geológico mediante el filtro de Savitzky-Golay, y se logró localizar las inclusiones de caliche y las estimaciones de cálculo de las reservas de los fosfatos se encuentran muy circunscritas.

Abstract

In the Sidi Chennane basin exploitation inclusions of sterile hardpan—so-called “derangements” — are hard to detect and interfere with phosphate extraction. A Schlumberger resistivity survey over an area of 50 hectares was carried out. The Savitzky-Golay filtering method was used as a tool for denoising the data.

Savitzky-Golay (S-G) filters are one of the filters which can smoothen out the signal without much destroying its original properties. Despite their exceptional features, they were rarely used so far in the geophysical signal processing. The aim of this paper is the investigation of their properties in detail from the geoelectrical signal processing aspect. The experimental results indicate that S-G filter is better for denoising geoelectrical signal. Models of the geology were successfully obtained from Savitzky-Golay filtering method, which help mapping the phosphate deposit inclusions and the estimations of phosphate reserves were improved and better constrained.

Introduction

Morocco is the world's third largest phosphate producer, after the USA and China. Total mine production recorded by the Ministry of Energy and Mines in 2003 was 29.39 Mt so more than 75% of world reserves. Four major phosphate basins are now known and are being exploited, three of which are located in central-northern Morocco. The four main deposits of phosphate are: the Oued Eddahab basin situated in Sahara, the central ganntour basin near Youssoufia, the Meskala basin at east of Essaouira and the Oulad Abdoun basin situated near Khouribga. The existence of morocco sedimentary phosphate rock has been known since 1908 in the Meskala basin, but it had not generated significant interest until the discovery, in 1917, of the Oulad Abdoun basin.

The geological investigations carried out in Sidi Chennane phosphatic deposit in the Oulad Abdoun basin revealed a phenomenon of sterile hardpan, these sterile bodies are formed by accumulations of silicified limestones or by limestone blocks within an argillaceous matrix (Michard, 2008) which interfere with phosphate extraction and their resistivity is higher than the phosphate-rich mineral resistivity. The application of the electric prospection methods constitutes a suitable means to map these sterile bodies in order to establish a model of their distribution and would permit the definition of these structures before the mining front reaches them.

Geological framework and methodology of work

The Ouled Abdoun basin is the largest phosphate basin in Morocco. It's located about 100 km in south-east of Casablanca. The phosphate deposits of Ouled Abdoun area belongs to the western Moroccan Meseta, commonly considered being stable. The local sedimentary deposits resulting from a large transgression occurred in mid-Cretaceous. It consists of (Kchikach; 2002) : marly limestone and gypsum of Cenomanian, Turonian white limestones, Senonian marl and yellow marly limestones, phosphatic series dated from Maastrichtian to Ypresian and Lutetian calcareous Thersitean slab. The Neogene continental deposits cover locally the marine series.

The loose phosphatic Levels exploited are cited according to their succession (Kchikach; 2002) Layer III: Maastrichtian, Layer II: Paleocene, Layers I and 0, “forrow” A and B: Ypresian. They are typically separated by phosphatic indurated limestone benches, more or less important and more or less regular called “infill”, II / I layers infill, III / II layers infill, etc.

Data acquisition

Electrical sounding is a method to investigate the change in earth resistivity with depth at a particular location. Horizontal electrical profiling is a method to determine lateral variations in earth resistivity within a limited depth range. Traditionally, arrangements using four electrodes (two current-transmitting electrodes and two voltage-sensing electrodes) are used for either vertical soundings or horizontal profiling. For vertical soundings, the electrodes are arranged symmetrically according to a center, with increasing distances between electrodes used to explore deeper depths. In the profiling mode, the distance between the potential and current dipoles (a dipole consists of a pair of like electrodes) is maintained while the array is moved along the profile line for mapping lateral changes (Chouteau, 2001; Ouadif, 2011).

When using profiling techniques, an estimate of the earth resistivity is calculated by using the well-known relation between resistivity, an electric field, current density (called Ohm's Law), the geometry and spacing of the current and potential electrodes. When the earth is not homogeneous neither isotropic, this estimated value is called the apparent resistivity, ρ_upp, which is an average of the true resistivity in the measured section of the earth:

Where K is the geometrical factor that depends on the electrodes arrangement.

To cover all the zones being able to be disturbed, we carried out, during the geophysical prospection in a parcel of 50 ha, 41 vertical electrical soundings then 5151 resistivity measurements as horizontal profiling (SYSCAL-R2 resistivity instrument IRIS Instruments) using the well-known Schlumberger array, in order to map the spatial distribution of the sterile hardpan inclusions. The 5151 stations of the resistivity survey is a compilation of 51 profiles spaced at 20 m. For every profile, there were 101 stations separated by 5 m.

The VES (some of them have been done in old trenches of exploitation where these bodies are visible) helped us: (1) to determine the intrinsic resistivity for the different terms of the phosphatic series (Baba, 2012), indeed, the apparent resistivity values above the sterile bodies are between 200 and 250 Ωm and the observed values on normal phosphatic series are around 150 Ωm, (2) to select the appropriate length of the measurement device for the profiling survey. Anomalies are in the upper part of the phosphate layer 25 meters to 60 meters of depth, leading us to choose 120 m for the maximum device length of the resistivity profiling survey.

Data Analysis

To map resistivity contrasts, field data, apparent resistivity measurements obtained from horizontal profiling with 40 m, 80 m and 120 m device length, are plotted on a map of the surveyed area, using Surfer desktop software. Kriging interpolation was used in order to create maps of the corresponding geophysical data.

Iso-resistivity map was prepared and interpreted in terms of resistivity and thickness of sub-surface layer and resistivity results were correlated with the existing lithology:

• •

  The kriged maps of electrical resistivity distribution AB = 40 denote a higher mean resistivity distribution representing the calcareous Thersitean slab effect.

• •

  The kriged maps of electrical resistivity distribution AB = 80 and AB = 120 show anomalous zones of directions NE-SW located at an average depth of 60 m. Thus, the majority of anomalies affect only the subsurface of the phosphate series.

Denoising geoelectrical data

Noise is often a significant issue in electrical resistivity data, field data rarely comes clean, thus robust processing methods is required, in order to remove non-geological noise and extract useful information from raw data.

As a result of the above processes, the maps resulting by the interpolation of the geophysical measurements often has poor quality, containing high percentages of random or systematic noise which hinder the valuable information related to the subsurface targets which are here essentially the inclusions. Therefore, the suppression of the noise levels and the enhancement of the signals carrying the useful information is an important process in any processing approach.

Savitzky-Golay

Savitzky and Golay (1964), interested in smoothing of noisy data obtained from chemical spectrum analyzers, demonstrated that fitting a polynomial to a set of input samples and then evaluating the resulting polynomial at a single point within the approximation interval is equivalent to discrete convolution with a fixed impulse response (Savitzky, 1964).

The least squares digital polynomial smoothing filters, popularized by Savitzky-Golay are widely used for smoothing and differentiation signal processing in many fields namely in Biomedical signals usually known as non-stationary such as electroencephalogram (EEG) and electrocardiogram (ECG) (Azami 2012), in elastography and magnetocardiogram and mainly in absorption spectroscopy (Hargittai 2005). ECG signal can potentially corrupted by various types of noise which lead to incorrect the diagnosis. Many types of filters are available to smoothing the noisy ECG signal Golay is one of the filter which can smoothen out the signal without destroying its original properties[(Hassanpour 2007; Orfanidis 2010; Schafer 2011; Hassanpour et al. 2012).

The filters are constructed to fit a particular polynomial to a windowed portion of the signal, with the least squares method, and then replace the central point of the window with the value of the polynomial at that point to produce a smoothed output. The Saviztky-Golay filtering can be considered as a generalized moving average.

The polynomial can be represented as:

Where ρ is the corresponding apparent resistivity data vector, x_i represent the north coordinate of a gidded point of the resistivity map.

To design these filters, we have first to decide on the length of the filter k, order of polynomial p, order of derivative n and the size of the smoothing window N where N is odd, and N ≥ p + 1. The coefficients of a Savitzky-Golay filter, when applied to a signal, perform a polynomial is fitted to points N = N_r + N_l + 1 of the signal. Nr and Nl are signal points in the right and signal points in the left of a current signal point, respectively.

To find these polynomial coefficients, we solve a least-squares inversion:

Where

We can find the vector of polynomial coefficients via least-squares solution of the matrix equation

The matrix inversion expresses each of the polynomial coefficients in c as a linear combination of the rows of (M^t M)⁻¹M^tρ.

Importantly, the value of the polynomial at point ρ₀ is simply given by c₀ (since all other values in the polynomial are zero). The Savitzky-Golay filter of derivative order 0 is given by the (time-reversed) coefficients of the central row of (M^t M)⁻¹M^tρ. In general, the coefficients for the Savitzky-Golay filter of derivative order n will be given by the (time-reversed) row n + 1 of the matrix inverse.

FiltersResults From Denoised Data

To achieve a high level of smoothing without attenuating the extrema in the data, a Savitzky-Golay filter with a quartic (order 4) polynomial and the mean centering was applied.

To map resistivity contrasts, field data, apparent resistivity measurements obtained from smoothed data are plotted as a map of the surveyed area.

The resistivity maps generated by the Savitzky-Golay smoothed data allow enhancing the qualitative interpretation of iso-apparent resistivity, in figure 8, 9 and 10, the postulated anomalies appear clearly delineated by the Savitzky-Golay filtered resistivity maps, thus a high order polynomial (n = 4) allows a high level of smoothing without attenuation of data features.

Figure 8.

Maps of the apparent resistivity, raw and smoothed data; for Schlumberger electrodes spacing AB = 40 m

Figure 9.

Maps of the apparent resistivity, raw and smoothed data; for Schlumberger electrodes spacing AB = 80 m

Figure 10.

Maps of the apparent resistivity, raw and smoothed data; for Schlumberger electrodes spacing AB = 120 m.

Conclusion

A detailed geoelectric resistivity survey in a parcel of 50 ha was executed using the horizontal electrical profiling approach with Schlumberger array. Iso-apparent resistivity maps were constructed at electrode spacings equal 40, 80 and 120 m in order to study the lateral variations in the geoelectric behavior in phosphatic series of Sidi Chennane.

In this paper, the Savitzky-Golay smoothing method is used to denoise geophysical data to detect anomalous zones of Sidi Chennane phosphatic series, their localisation would permit the mining engineers to get around them during the exploitation.

The Savitzky-Golay filter, method of least-squares-fit smoothing and differentiation of digital data, is versatile and extremely simple to use.

The qualitative interpretation of the smoothed resistivity maps allows defining resistivity contrast, consequently we have delimited the crossing dominate area from a “normal” into a “disturbed” area. We have also observed a best similarity to original signal and slightly denoised it; indeed, this filter is often preferred because, when it is appropriately designed to match the waveform of an oversampled signal corrupted by noise, it tends to preserve the width and height of peaks in the signal waveform. The efficiency of Savitzky-Golay noise reduction filtering estimation technique is shown to be significantly better.