Introduction

Geophysical studies show ample evidences of the successful use of the resistivity method in groundwater prospecting. Many geophysical techniques have been applied to groundwater investigations with some showing more success than others. Among all the geophysical techniques, the geoelectrical resistivity method is more successful in groundwater studies. Direct current resistivity methods of geophysical exploration are in extensive use globally for aquifer mapping and estimation of aquifer parameters viz., resistivity and thickness of the subsurface earth. While applying the current to the subsurface of the earth, due to its heterogeneity, and nonlinear characteristics of various soils, rock formations and aquifer zones, it will effectively distract the applied current and gave us information about the subsurface, as the resistance varies from various contents of the subsurface earth. Due to the non linearity noises were encapsulated with the information. A better technique need to eliminate such noises. Here experimented wavelet denoising algorithm with different thresholding functions and the results proved it would be the best tool for purging the noise.

The use of Schlumberger electrode configuration for conducting vertical electrical sounding is more common. Vertical Electrical Sounding (VES) method is one of the efficient methods in identifying the aquifers in a more reliable manner. In addition to the interpretation of VES data using curve matching procedure through theoretical curves prepared by different computational methods, the inversion has been used for many years to interpret electrical resistivity data. The theoretical interpretation of geoelectrical resistivity data using Adaptive Neuro Fuzzy Inference System (ANFIS) was adopted in the research work. The processing of each component of the algorithm of the investigation essentially was qualitative followed by quantitative analysis. Thus it has been validated with different field data. The field data include, the VES data collected from Kanyakumari field areas. The results are compared with earlier results/litholog for the published data and with litholog for the Kanyakumari field areas.

Methodology

Direct current resistivity methods of geophysical exploration are in extensive use globally for aquifer mapping and estimation of aquifer parameters viz., resistivity and thickness (Kosinky and Kelly 1981; Niwas and Singhal 1981; Mazac et al. 1985; Yadav and Abolfazli 1998). The physical basis of the resistivity method is based on the relative distribution of impressed current in the earth controlled by subsurface resistivity distribution. Different computational methods were adapted by the various researchers to invert geoelectrical data. (Flathe 1955; Van Dam 1964; Mooney et al. 1966; Ghosh 1971; Srinivas et al. 2010, 2012).

In general, the characteristic sounding curves are represented in multiple layers Each of the four sets has particular properties that may be roughly classified. Types H and K have a definite minimum and maximum, indicating a bed, or beds, of anomalously low or high resistivity, respectively, at intermediate depth. Types A and Q show fairly uniform change in resistivity the first increasing, the second decreasing with depth. Obviously these curves also may be combined (Telford et al. 1990). Here in this research paper different combination of curves were attempted for interpretation.

Theoretical background

Wavelet transform plays major role in analysing the data in broad spectrum of time–frequency. Event related potential determination can be easily done through this transformation more clearly. The spatial and temporal information will provide the real time data strategy in eagle’s eye. Advantages over the conventional time–frequency analysis lies in decompose the signal of any time–frequency multiresolution functions. The admissibility condition is as follows.

ψ (t) ϵ L2 (R) satisfies certain admissibility conditions as

$$C \, \uppsi \, = \mathop \int \limits_{R} \frac{{\left| {\widehat{\varPsi }\left( \omega \right)} \right|^{2} }}{\left| \omega \right|} d\omega < \, \infty$$
(1)

\(\varPsi \left( t \right)\) is called wavelet; \(\widehat{\varPsi }\) (ω) is fourier transform of \(\varPsi \left( t \right)\). The Continuous Wavelet Transform (CWT) is defined as the sum over all time of the signal multiplied by scaled, shifted versions of the wavelet function,

$$\varPsi_{a,b} \left( t \right) = \left| a \right|^{{ - \frac{1}{2}}} \varPsi \left( {\frac{t - b}{a}} \right)$$
(2)
$$W_{f} \left( {a,b} \right) = \left\langle {f\left( t \right), \varPsi_{a,b} \left( t \right) = \varPsi_{a,b} \left( t \right) = \frac{1}{\sqrt a }\varPsi \left( {\frac{t - a}{a}} \right) a,b \in R, a > 0} \right.$$
(3)

Wavelet denoising

This application of wavelet denoising is a nonparameteric regression analysis and it works well with added gaussian noise in the true resistivity data. The noisy data have the form

$$y_{i} = g\left( {x_{i} } \right) + e_{i}$$
(4)

Here \(g\left( {x_{i} } \right)\) is the true function to be experimented and \(e_{i}\) which represents the noise added to the signal. This true data which additive noise is transformed to a denoised data which is of the form

$$d^{*} = d + \in$$
(5)

where \(d\) is the transform off the true data and \(\in\) is the transform of the noise.

Hard thresholding

$$\begin{aligned} D\left( {U, \lambda } \right) = U \quad for\, all\, \left| {U > \lambda } \right| \hfill \\ = \, 0 \quad{\text{ Otherwise}} \hfill \\ \end{aligned}$$
(6)

Soft thresholding

$$D\left( {U, \lambda } \right) = sgn\left( U \right){ \hbox{max} }\left( {0, \left| U \right| - \lambda } \right)$$
(7)

Universal threshold

$$\lambda^{u} = \sigma \sqrt {2logn}$$
(8)
$$\sigma = MAD \left( {d_{J - 1} } \right) = median_{i} \left( {\left| {d_{J - 1, i} - median_{k} \left( {d_{J - 1, k } } \right)} \right|} \right)$$
(9)

Minimax threshold

It is better to know some prior information about the data, but it is difficult to understand the distribution of complex data. Explicitly, the obtained geoelectrical apparent resistivity data from different geologic conditions will suffer different kind of noises which could not be easily discriminated. Since there is no homogeneous model to fix the distribution, the data is never in uniform behaviour. For the data it is good to find a best estimator whose maximal risk is minimal among all estimators. The prior information of the signal set is \(\varTheta\). But this set doesnot specify the probability distribution of data in \(\varTheta\). The more prior information smaller the set \(\varTheta\) (Mallat 2009). For the data \(f \in \varTheta\) with noise as given in Eq. (4). The risk estimation \(\tilde{F} = DX\) is

$$r \left( {D,\, f} \right) = E\left\{ {DX - f^{2} } \right\}$$
(10)

The probability distribution of data in set \(\varTheta\) is unknown; the precise risk cannot be estimated. Only the range can be determined. The maximum risk of this range is given by (Donoho and Johnstone 1998)

$$r \left( {D, \varTheta } \right) = \sup E\left\{ {DX - f^{2} } \right\} \quad {\text{with}}\quad f \in \varTheta$$
(11)

In this estimation of minimax function, the lower bound risk of above equation with all possible with operators D is represented as

$$r_{n } \left( {D, \varTheta } \right) = Inf r \left( {D \in \varTheta } \right) \quad {\text{with}}\quad D \in O_{n}$$
(12)

where \(O_{n}\) represents the sets of all operators

Rigrsure thresholding

It is a soft threshold evaluator of unbiased risk. W = [\(\omega_{1} , \omega_{2,} \ldots \ldots \ldots ..\omega_{N}\)] is a vector consists of the square of wavelet coefficients from small to large. Select the minimum value \(r_{b} \left( {b^{th} r} \right)\) from risk vector, which is given as,

$$R \, = \left\{ {r_{i} } \right\}_{i = 1,2, \ldots \ldots .N} = \frac{{\left[ {{\text{N}} - 2{\text{i }} + \left( {{\text{N }} - {\text{i}}} \right) \omega_{i} + \mathop \sum \nolimits_{k = 1}^{l} \omega_{k} } \right]}}{N} .$$
(13)

as the risk value. The selected threshold is

$$\lambda = \sigma \sqrt {\omega_{b} }$$
(14)

where, \(\omega_{b}\) is the b th squared wavelet coefficient (coefficient at minimum risk) chosen from the vector and σ is the standard deviation of the noisy signal.

Heursure thresholding

If the signal to noise ratio is very small, the sure method estimation is not good to choose. The Heuristic threshold is given by,

$$\lambda = \left\{ {\begin{array}{ll} \lambda_{1} &\quad A > B \\ \hbox{min} \left( {\lambda_{1} ,\lambda_{2} } \right) & \quad A \ge B \end{array} } \right.$$
(15)

where \(A \, = \frac{S - N}{N}\) and \(B \, = \left( {log_{2} N} \right)^{3/2} \sqrt N\). The N is length of wavelet coefficient vector and is the sum of squared wavelet coefficients given as

$$S \, = \mathop \sum \limits_{i = 1}^{N} \omega_{i}^{2} .$$
(16)

A small threshold which gives noisy output and large threshold can cut significant part of signal thus the information contained in the signal may be lost.

Adaptive Neuro Fuzzy inversion

The role model for Soft Computing is the Human Mind. Unlike the conventional or hard computing, it is tolerant of imprecision, uncertainty and partial truth. It is also tractable, robust, efficient and inexpensive. Apart from neural networks and fuzzy logic, ANFIS is hybrid approach network architecture integrates the concepts of both ANN and FL. Thus the fundamental knowledge was supplied by the ANNs and the fuzzified membership function rules thus evaluate the real structure of the earth more clearly.

In order to make full use of soft computing for intelligent subsurface characterization, it is important to note that the knowledge based system has been framed carefully and the implementation of the hybrid systems should aim to improve prediction and its reliability. At the same time, the improved systems should contain small number of sensitive user-definable model parameters and use less CPU time. The future development of hybrid systems should incorporate various inter disciplinary knowledge and maximize the amount of useful information extracted between data types so that reliable extrapolation could be obtained. In recent years, considerable attention has been devoted to the use of hybrid neural-network/fuzzy-logic approaches as an alternative for pattern recognition, statistical and mathematical modeling. It has been shown that neural network models can be used to construct internal models that recognize fuzzy rules.

Neuro-fuzzy modeling is a technique for describing the behavior of a system using fuzzy inference rules within a neural network structure. The model has a unique feature in that it can express linguistically the characteristics of a complex nonlinear system. Geoelectrical resistivity method interpretation can be done with various tools. ANFIS network was developed by Jang 1993. An adaptive network is a multilayer feed forward network in which each node performs a particular function on giving the input to the network. Electrical resistivity data collected using the vertical electrical sounding (VES) method is used for training the data set. The data collected from the field will be the apparent resistivity data, while on interpretation the subsurface parameters viz., resistivity and thickness of the individual layers are obtained. Trained data set will be the reference data for interpreting the subsurface layer parameters of the earth. This dataset will be used to train ANFIS by adjusting the membership function parameters that best models this data.

ANFIS architecture

The basic structure of the fuzzy inference system that maps input characteristics to input membership functions, input membership function to rules, rules to a set of output characteristics, output characteristics to output membership functions, and the output membership function to a single-valued output or a decision associated with the output.

ANFIS system consists of 5 layers, layer symbolized by the box is a layer that is adaptive. Meanwhile, symbolized by the circle is fixed. Each output of each layer is symbolized by O 1,i with i is a sequence of nodes and 1 is the sequence showing the lining. Here is an explanation for each layer, namely:

Layer 1

Serves to raise the degree of membership

$$O_{1,i} = \, \mu_{A} \left( x \right),i = \, 1,2.$$
(17)

and

$$O_{1,i} = \, \mu_{B} \left( y \right),i = \, 1,2.$$
(18)

with x and y are the input for the i-th node

$$\mu A(x) = 1/[1 + (det(x - ci)/ai)^{ \wedge } 2bi],$$

by {ai, bi and ci} are the parameters of membership function or called as a parameter premise.

Layer 2

Serves to evoke firing-strength by multiplying each input signal.

$$O_{2,i} = \, w_{i} = \, \mu_{A} \left( x \right)x\mu_{B} \left( y \right), \, i \, = 1, \, 2.$$
(19)

Layer 3

Normalize the firing strength

$$O_{3,i} = \overline{{w_{i} }} = \frac{{w_{i} }}{{w_{1} + w_{2} }},i = 1,2.$$
(20)

Layer 4

Calculating the output based on the parameters of the rule consequent {p i , q i and r i }

$$O_{4,i} = \overline{{w_{i} }} f_{i} = \left( {p_{i} x + q_{i} y + r_{i} } \right).$$
(21)

Layer 5

Counting the ANFIS output signal by summing all incoming signals will produce

$$\mathop \sum \limits_{i} \overline{{w_{i} }} f_{i} = \frac{{\mathop \sum \nolimits_{i} w_{i} f_{i} }}{{\mathop \sum \nolimits_{i} w_{i} }}.$$
(22)

Functionally, there are almost no constants on the node functions of an adaptive network except piecewise differentiability.

Results and discussion

Resistivity sounding data were obtained from certain regions of Kanyakumari district (8°6′31.79″ N; 77°30′50.29″ E). Synthetic noise was added to the data with 2, 5, 10 and 15 % for which Rigrsure, Heursure, Sqtwolog, Minimax wavelet thresholding functions were applied and the results were experimented. The experimental results of VES 1 with certain noise percentage using Daubechies decomposition of various levels was tabulated (Table 1). The experimental analysis were done by using MATLAB software (MATLAB R 2013b). Besides noting down the performance analysis of various methods, the effect of decomposition levels also has been checked. For comparing and analysing the performance the true resistivity and depth data were compared with the litholog data available near to the data location.

$${\text{Average error in true Resistivity }} = mean\left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} Actual \,resistivity - Obtained \,resistivity} \right)$$
$${\text{Average error in depth }} = mean\left( {\frac{1}{N}\mathop \sum \limits_{i = 1}^{N} Actual \,depth - Obtained \,depth} \right)$$

Where N represents the number of layers.

Table 1 Experimental result of Daubechies wavelet functions with different levels of decomposition for VES 1 data
Full size table

Based on the result universal threshold (Sqtwolog) shows the best performing threshold function of wavelet denoising for geoelectrical data inversion. Thus the VES data 1, 2, 3, 4, 5, 6, 7 were interpreted and the results were shown in figures. Figures 1 and 2 represents the VES 1 inverted model and error contribution of resistivity and thickness respectively. Figure 3 shows the interpreted model for VES 3 data and the corresponding error contribution is shown in Fig. 4. Table 2 shows the litholog for corresponding sounding data. Figures 5 and 6 shows the model with error contribution of VES 3. Figures 7, 9, 11 and 13 shows the interpreted layer model for sounding data 4, 5, 6 and 7 respectively. Figures 8, 10, 12 and 14 represents the corresponding error contribution of each VES data respectively. Error contribution graph shows the average error of each layer for true resistivity and depth taken into account. Approximate sharing of errors for each layer is shown in graphs.

Fig. 1
figure 1

VES 1 inverted layer model with average error

Full size image
Fig. 2
figure 2

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 1

Full size image
Fig. 3
figure 3

VES 2 inverted layer model with average error

Full size image
Fig. 4
figure 4

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 2

Full size image
Table 2 Vertical electrical sounding locations with lithology
Full size table
Fig. 5
figure 5

VES 3 inverted layer model with average error

Full size image
Fig. 6
figure 6

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 3

Full size image
Fig. 7
figure 7

VES 4 inverted layer model with average error

Full size image
Fig. 8
figure 8

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 4

Full size image
Fig. 9
figure 9

VES 5 inverted layer model with average error

Full size image
Fig. 10
figure 10

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 5

Full size image
Fig. 11
figure 11

VES 6 inverted layer model with average error

Full size image
Fig. 12
figure 12

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 6

Full size image
Fig. 13
figure 13

VES 7 inverted layer model with average error

Full size image
Fig. 14
figure 14

Error contribution map for true resistivity and thickness in comparison with actual litholog of VES 7

Full size image

Figure 15 shows the sample denoising application of wavelet denoising response to the noisy data. The correlation coefficient of wavelet denoising with respect to the information losing phenomena for VES 1 is shown in Table 3. The results obtained from Table 3 clearly portray the denoising capability of wavelet algorithm. Above 10–15 % the information gets distracted and gets lost finally when it reaches the noise level above 15 %.

Fig. 15
figure 15

Wavelet denoising application of noisy data using Daubechies wavelet functions with level 4 decomposition

Full size image
Table 3 Noise level estimation with respect to the correlation coefficient
Full size table

Conclusions

In this paper, performance of various thresholding methods of wavelet viz., Sqtwolog, Rigrsure, Heusure, Minimax is applied for denoising geoelectrical sounding dataset. Vertical Electrical Sounding data obtained from Kanyakumari district was tested and verified with actual litholog section. The effect of various decomposition levels were investigated on the performance of various thresholding functions. The results show that as level of decompostion increases from 2 to 15 % in correspondence to the information losing criteria during the wavelet denoising algorithm. As from the above studies it reveals that sqtwolog universal thresholding proves to be giving better result than the other algorithm which has been taken under consideration. This research work accompanies the wavelet denoising procedure which has been one of the best pre-processing algorithms for attempting to any signal decomposition or inversion techniques.