Inversion of Amplitude from the 2-D Analytic Signal of Self-Potential Anomalies

2.1. Forward modeling

Following Nabighian [71], the 2-D analytical signal amplitude (ASA) or the total gradient (TG) of SP anomaly is given as.

E1

where V_xand V_yare the horizontal and vertical derivatives of the SP anomalies, respectively. The 2-D ASA of few of the idealized geobodies can be estimated following the general expression of SP anomaly ASA (x)for idealized causative bodies (horizontal/vertical cylinder, sphere, inclined sheet, etc.) at any point on the surface of the earth (Figure 1) [70, 71, 74, 75], which is given by the equation:

E2

where kis the amplitude coefficient/factor related to the physical properties of the source, zis the depth from the surface to the top of the body (sphere, cylinder, sheet, line of poles, point pole), x₀ (i = 1, …,N) is the horizontal position coordinate on the surface, and qis the shape factor. The shape factors (q) for horizontal/vertical cylinder, sphere, inclined sheet, line of poles, and point poles are 1, 1.5, 1, 0.5, and 1, respectively. The total derivation of the ASA or TG anomaly can be found in various literatures [71, 72, 73, 74, 75, 81].

The above equation can be used to interpret single structure. However, in order to interpret the numerous structures, the expression (2) can be rewritten as [36]:

E3

where ASA_j(x_i) is the SP anomaly at x_ilocation for jth body and Mis the number of bodies.

2.2. Global optimization

Several global optimization approaches have been effectively used in interpretation of different geophysical data (e.g., [35, 36, 37, 75, 77, 82, 83, 84, 85, 86]). In the present study, an alternative method of simulated annealing (SA) which is called as very fast simulated annealing (VFSA) was applied in the present study to elucidate the total gradient SP anomaly. VFSA is different from SA in terms of faster cooling schedule due to its finer Cauchy probability distribution for the arbitrary selection of every model parameter. VFSA takes any value in a model space, whereas SA does it in a predefined model space and hence the resolution increases for VFSA. Moreover, in VFSA optimization process, it does not recall all models and hence needs very small memory [35, 87, 88]. In every geophysical inversion, the main objective is to minimize the error function or the misfit. In the present work, the misfit (φ) between the observed and calculated/model response was used for SP data interpretation because of the fact that the objective function gets affected (increases) near zero crossing of SP anomaly (after [78]).

where N is the number of data point,

For the present inversion of SP data using VFSA optimization process, different parameters such as initial temperature, cooling schedule, number of iterations, and number of moves per temperature were taken as 1.0, 0.4, 2000, and 50. To find out the global/optimum solution, probability density function (PDF) was taken within 60.65% limit and ambiguity study has also been carried out based on the techniques developed by Mosegaard and Tarantola [89], Sen and Stoffa [90]. The details of the inversion process can be found in different literatures such as Sen and Stoffa [77], Sharma [86], Sharma and Biswas [78], Biswas [82], and Biswas [5]. Because in the earlier studies, the VFSA algorithm has established a very good performance as a powerful optimization method for estimating multidimensional and multimodal error functions. Hence, this optimization strategy is applied for SP data enhancement and regularization. The present VFSA algorithm for interpretation of ASA or TG of SP anomaly was carried out in Windows 8 environment using MS FORTRAN Developer studio on a simple desktop PC with Intel Core i7 processor. For each step of optimization, an overall of 10⁶ forward computations (2000 iteration × 50 number of moves × 10 VFSA runs) were accomplished and accepted models were stored in memory. The total time required (not CPU time) to compute a sole inversion is 35 s.

3.1. Parameter search range

An appropriate selection of initial guess values or the search range for every model parameters is a significant objective of any inversion approach. For the quantitative elucidation of SP anomaly from ASA, the horizontal position of the source is found from the highest anomaly, the depth from half width [71], and the size of the amplitude from highest amplitude and depth. Srivastava and Agarwal [70] stated that the most stable parameter for interpretation of SP anomaly is the horizontal location, which changes very little. Moreover, it was discussed earlier that limiting the shape factor also gives the consistent results in terms of ambiguity and error [35]. Hence, in order to find the accurate or the true value of each model parameter, initially the search ranges for each model parameters were kept wide so as to find out the most probable solution. Next, the search spaces were reduced to find the more appropriate results. When the model parameter gives the utmost results and very least error, the location (x₀) and the shape factor (q) were fixed to their real/true values to further reduce the error in the amplitude and the depth so as to find the actual depth of the subsurface bodies and reduce ambiguity in the final interpretation.

3.2. Synthetic examples

The VFSA inversion technique was utilized considering synthetic noise-free and noise-corrupted data (10 and 20% random noise) for self-potential anomaly for various subsurface structures derived from the ASA anomaly. At first, every model parameter was inverted for every data. Synthetic hypothetical data were produced utilizing Eq. (2) for different idealized geobodies, and 10 and 20% random noise is corrupted to the synthetic data. VFSA inversion was applied utilizing noise-free and noisy data to retrieve the genuine model parameters and analyze the impact of noise on the deciphered model parameters. In general, a reasonable search range/space for every model parameter was chosen and one VFSA run was performed. Thereafter, the best possible convergence of every model parameter was studied (k, x₀, z, and q) and misfit.

Next, to acquire the mean model, 10 VFSA runs were executed. At that point, histograms were constructed from the retrieved models whose misfit error was lesser than 10⁻⁴. Then, a statistical mean model was calculated utilizing models that have misfit lesser than 10⁻⁴, which also exist inside one standard deviation. Besides, cross-plots were also developed and analyzed to check whether the model parameters were inside the high PDF region (60.65%). Further, noises were also added in the data and the process was repeated again where the misfit error was lower than 10⁻². Finally, the comparison between the observed and model responses was shown for every model. This process was applied for each hypothetical synthetic noise-free, noisy, and field cases.

3.2.1. Model 1 (sphere)

Inversion of the SP data was performed as said above utilizing synthetic noise-free data for sphere-like structure. Figure 2 demonstrates the convergence example for every model parameter. Figure 3a demonstrates the histogram for every model parameter (k, x₀, and z). The histogram uncovers that all the parameters of the body can be very much interpretable after inversion. Moreover, very less ambiguity in the interpretation of the main three model parameters was found. The cross-plots analysis (Figure 4a) demonstrated that the model parameters were near their true value (green). The final estimated model parameters were within the ambiguity limits and inside the high PDF region (red). The fittings between the observed and model response are shown in Figure 6a. The deduced model parameters and mean model are shown in Table 1. Next, another model was also selected to see the variation in the amplitude, location, and depth. Inversion was repeated the same way as discussed above. Figure 3bshows the histogram for the three model parameters. Figure 5a demonstrates the cross-plots for this model as well. The fittings between the observed and model response are shown in Figure 6b, and the elucidated parameters and mean model are shown in Table 2.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	−5000	−10,000 to 0	−5002.2 ± 27.4	−4907.7 ± 152.2
x0 (m)	200	0–300	200.0 ± 0.0	200.0 ± 0.0
z(m)	20	0–50	20.0 ± 0.0	19.9 ± 0.3
q	1.5	0–2	1.5 ± 0.0	1.5 ± 0.0
Misfit			6.6 × 10−9	9.4 × 10−5

Table 1.

Actual model parameters, search range, and interpreted mean model for noise-free, 10% random noise with uncertainty for sphere.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	−7500	−10,000 to 0	−7497.4 ± 34.6	−7037.4 ± 151.6
x0 (m)	250	0–300	250.0 ± 0.0	250.0 ± 0.0
z(m)	30	0–50	30.0 ± 0.0	29.9 ± 0.3
q	1.5	0–2	1.5 ± 0.0	1.5 ± 0.0
Misfit			4.2 × 10−8	7.1 × 10−4

Table 2.

Actual model parameters, search range, and interpreted mean model for noise-free, 20% random noise with uncertainty for sphere.

3.2.2. Model 2 (horizontal and vertical cylinder)

Another synthetic example for horizontal and vertical cylinder-like body was taken for inversion of the SP data. Figure 7a shows the histogram for every model parameter (k, x₀, and z). The histogram shows that all the parameters of the cylindrical structure can be very much elucidated. The cross-plots analysis (Figure 8a) demonstrates that the model parameters were near their true value (green) and the final expected model parameters were within the ambiguity limits and within high PDF region (red). Observed and model response fits are shown in Figure 10a. The elucidated parameters and mean model are shown in Table 3. Similarly, for the sphere-like structure, another model was selected to see the variation in the model parameters. Figure 7b demonstrates the histogram for the three model parameters. Figure 9ademonstrates the cross-plots of this model. Observed and model responses are shown in Figure 10b, and the elucidated parameters and mean model are shown in Table 4.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	−7500	−10,000 to 0	−7497.9 ± 23.7	−7357.3 ± 122.8
x0 (m)	200	0–300	200.0 ± 0.0	200.0 ± 0.0
z(m)	30	0–50	30.0 ± 0.0	31.1 ± 0.4
q	1.0	0–2	1.0 ± 0.0	1.0 ± 0.0
Misfit			5.2 × 10−11	9.8 × 10−4

Table 3.

Actual model parameters, search range, and interpreted mean model for noise-free, 20% random noise with uncertainty for horizontal/vertical cylinder.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	−1000	−2000 to 0	−999.6 ± 3.3	−963.1 ± 22.6
x0 (m)	250	0–300	250.0 ± 0.0	250.0 ± 0.0
z(m)	20	0–50	20.0 ± 0.0	19.6 ± 0.4
q	1.0	0–2	1.0 ± 0.0	1.0 ± 0.0
Misfit			5.9 × 10−9	1.6 × 10−4

Table 4.

Actual model parameters, search range, and interpreted mean model for noise-free, 10% random noise with uncertainty for horizontal/vertical cylinder.

3.2.3. Model 3 (inclined sheet)

A synthetic example for 2-D inclined sheet-type body was also taken for inversion of the SP data. Figure 11a shows the histogram for every model parameter (k, x₀, and z). The cross-plots (Figure 12a) also demonstrate that the model parameters are near their actual value (green) and the final plausible model parameters lie inside the estimated ambiguity limits and inside high PDF region (red). Observed and model response fits are shown in Figure 14a. The elucidated parameters and mean model are given in Table 5. Another model was selected to see the variation in the model parameters. Figure 11b displays the histogram for the three model parameters. Figure 13a shows the cross-plots of this model. Observed and model response fits are shown in Figure 14b, and the elucidated parameters and mean model are shown in Table 6.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	100	0–200	100.1 ± 0.4	96.3 ± 1.2
x0 (m)	250	0–300	250.0 ± 0.0	250.0 ± 0.0
z(m)	40	0–50	40.0 ± 0.0	40.6 ± 0.4
q	1.0	0–2	1.0 ± 0.0	1.0 ± 0.0
Misfit			4.7 × 10−7	8.1 × 10−5

Table 5.

Actual model parameters, search range, and interpreted mean model for noise-free, 10% random noise with uncertainty for sheet.

Model parameters	Actual value	Search range	Mean model (noise-free)	Mean model (noisy data)
k(mV)	500	0–1000	500.0 ± 1.6	493.5 ± 10.5
x0 (m)	200	0–300	200.0 ± 0.0	200.0 ± 0.0
z(m)	20	0–50	20.0 ± 0.0	19.5 ± 0.3
q	1.0	0–2	1.0 ± 0.0	1.0 ± 0.0
Misfit			2.0 × 10−9	6.2 × 10−4

Table 6.

Actual model parameters, search range, and interpreted mean model for noise-free, 20% random noise with uncertainty for sheet.

3.3. Noise analysis

To see the effectiveness of the present inversion results, the synthetic noise-free data must be corrupted with different degrees of noises and reinterpreted using the inversion method to check its performance and robustness of the method. Hence, 10 and 20% random noise was added to the data for sphere, cylinder, and sheet-type structure, and the procedure was repeated again to examine the effect of noise. Figure 3c and d displays the histogram for noise-corrupted sphere-like structure. Figures 4b and 5b show the cross-plots for noisy data. Tables 1 and 2 show the interpreted mean model for noisy data and the fits between the observed and model response for noisy model are shown in Figure 6c and d. For horizontal and vertical cylinder-like structure, again noise was corrupted with the two models mentioned above. Figure 7c and d shows the histogram for noise-corrupted horizontal/vertical-like structure. Figures 8b and 9b show the cross-plots for noisy data. Tables 3 and 4 show the interpreted mean model for noisy data and the fits between the observed and model response for noisy model are shown in Figure 10c and d. Also, for 2-D inclined sheet-type structure, the same amount of noise was corrupted in the models mentioned above. Figure 11c and d shows the histogram for noise-corrupted inclined sheet-type structure. Figures 12b and 13b illustrate the cross-plots for noisy data. Tables 5 and 6 show the interpreted mean model for noisy data and the fits between the observed and model response for noisy model are shown in Figure 14c and d. It can be seen from the study of histogram and cross-plots that the model parameters are interpreted very precisely, and it also advocates that the appraised model parameters for all structures are inside the ambiguity limits and within high PDF region. This also suggests that the inversion methodology developed for the elucidation of SP anomalies can precisely decide every model parameter and even if the data are highly corrupted with different degrees of noises.

3.4. Effect of complicated structure

It is significant to mention that that, in nature (field examples), it is very difficult to get an idealized geobody or structure. Moreover, the structures are mostly corrupted with different degree of noises. Also, there might be multiple structures within the subsurface of different types. In those cases, the anomaly from multiple structures will be very difficult to interpret. To test whether the inversion method can accurately identify the multiple structures, three different structures were taken and forward responses have been computed using Eqs. (2)and (3). Single structure and multiple structures from the forward modeling are shown in Figure 15a. Inversion of the synthetic data was carried out the same way as for the single model examples. Fits between the observed and model response are also shown in Figure 15a. Moreover, to investigate the effect of noisy data, only 20% random noise was corrupted in the data and inversion was repeated. Figure 15b displays the fits between the observed and model response for noisy data. Table 7 shows the interpreted model parameters for each structure derived from the multiple anomalies. Both noise-free and noisy data are shown in Table 7. It can be seen from Table 7 that the misfit is quite less for noise-free and noise-corrupted data. Histogram and cross-plots show alike as revealed in other examples for single structures. However, it is not presented here for brevity.

3.5. Field example

To demonstrate the efficiency of the present method, two field examples from Turkey and Canada were taken for interpretation. The field anomaly of ASA/TG from SP anomaly was taken from the earlier published literatures as described below:

3.5.1. Sulleymonkoy anomaly, Ergani, Turkey

This field example was taken from the Sulleymonkoy SP anomaly, Ergani, Turkey [91]. The ASA or TG anomaly was derived following Srivastava and Agarwal [70] (Figure 18). The anomaly was interpreted by many workers using different interpretation methods [38, 92, 93, 94, 95, 96]. The SP anomaly data were interpreted using the VFSA global optimization technique in this study. The TG anomaly shows a large magnitude peak anomaly with two small anomalies. In the present case, initially the main peak anomaly was interpreted considering a single structure. The depth obtained by the present study for single body was found to be 34.5 m on the horizontal location of 66.3 m, and the structure was interpreted as a cylindrical body. Histogram plot for single body (Figure 16a) also shows that the model parameters were precisely determined. Analysis of cross-plots (Figure 17a) also shows that the assessed model parameters were within the ambiguity limits. Next, considering multiple structures, the inversion process was repeated again for three different structures considering three peak values. Analysis of histogram plot for multiple bodies (Figure 16b–d) also shows that the model parameters were precisely determined. Investigation of cross-plots (Figure 17b–d) also shows that the estimated model parameters were within the uncertainty limits and more precise than considering single body using the whole TG anomaly. However, the depth obtained for the peak anomaly was found to be 28.5 m on the location of 63.9 m considering multiple structures. The depth obtained by other workers such as Yungul [91] as 38.8 m, Bhattacharya and Roy [94] as 40.0 m, Agarwal [92] as 30.1 m, Sundararajan and Srinivas [96] as 36 m, and Srivastava and Agarwal [95] as 28.9 m. Moreover, the estimated uncertainty was less considering multiple structures rather than a single body. Table 8 shows a comparison between the different model parameters and misfit. Figure 18 shows the comparison between the observed TG anomaly and the model response.

Model parameters	Single body		Multiple bodies
	Search range	Present method B1	Search range	Present method B1	Search range	Present method B2	Search range	Present method B3
k(mV)	1000–200,000	81281.0 ± 957.5	1000–100,000	60888.2 ± 574.0	100–50,000	24244.9 ± 442.6	1000–1,000,000	935524.0 ± 12801.3
x0 (m)	50–100	66.3 ± 0.2	50–100	63.9 ± 0.2	50–200	115.8 ± 0.2	50–200	152.3 ± 0.2
z(m)	0–50	34.5 ± 0.3	0–50	28.5 ± 0.2	0–50	30.3 ± 0.3	0–100	40.0 ± 0.2
q	0–2	1.0	0–2	1.0	0–2	1.0	0–2	1.0
Misfit		5.6 × 10−3		9.8 × 10−4		8.5 × 10−5		1.4 × 10−3

Table 8.

Search range and interpreted mean model for Sulleymonkoy anomaly, Ergani, Turkey.

3.5.2. Senneterre area of Quebec, Canada

A total gradient of SP anomaly (Figure 6a) over a massive sulfide ore deposits in the Senneterre area of Quebec, Canada [97], was considered for this study. The anomaly was interpreted using enhanced local wave number technique by Srivastava and Agarwal [95] and regularized inversion by Mehanee [7] without considering TG. However, the anomaly was reinterpreted by Srivastava et al. [75] using ACO technique. Srivastava et al. [75] considered six anomalies from the TG. This anomaly was also elucidated using the present inversion method to retrieve the model parameters considering four peak anomalies, which were quite distinct. Investigation of histogram plot for multiple bodies (Figure 19a–d) also shows that the model parameters were precisely determined. Examination of cross-plots (Figure 20a–d) also shows that the appraised mean model parameters were within the ambiguity limits. The depth obtained by using the VFSA method was found to be 10.5, 4.6, 3.7, and 4.3 m, respectively. The depths obtained by Srivastava et al. [75] for the four bodies are 5, 4.2, 3.8, and 4.3 m, respectively. The assessed parameters in this study are in respectable agreement with the other work. Moreover, in the present work, all four anomalies were interpreted as multiple structures and not independently. The estimated model parameters are shown in Table 9 along with misfit. A comparison between the observed anomaly and model responses are shown in Figure 21.

Model parameters	Search range	Present method B1	Search range	Present method B2	Search range	Present method B3	Search range	Present method B4
k(mV)	10–500	122.1 ± 0.7	100–2000	1568.4 ± 13.9	0–100	85.4 ± 0.6	0–100	69.5 ± 0.3
x0 (m)	100–200	147.2 ± 0.1	100–200	166.1 ± 0.1	200–250	232.0 ± 0.1	200–250	240.9 ± 0.1
z(m)	0–20	10.5 ± 0.1	0–10	4.6 ± 0.1	0–10	3.7 ± 0.1	0–10	4.3 ± 0.1
q	0–2	0.5	0–2	1.0	0–2	0.5	0–2	0.5
Misfit		3.6 × 10−3		2.2 × 10−3		3.4 × 10−3		5.5 × 10−3
Srivastava et al. [75] (ACO) x0 and z0 are the coordinate of the source location in unit of data spacing (0.61 m).
	Search space	TGM peak 1	Search space	TGM peak 3	Search space	TGM peak 5	Search space	TGM peak 6
k(mV)	1–6250	157 ± 17	1–42,000	6998 ± 72	1–15,000	2201 ± 36	1–13,000	944 ± 11
x0 (m)	220–240	143.2 ± 2.0	265–285	167.4 ± 2.1	370–390	231.2 ± 1.2	380–410	242.5 ± 3.1
z(m)	1–50	5.0 ± 0.2	1–50	4.2 ± 0.2	1–50	3.8 ± 0.3	1–50	4.3 ± 0.1
β(=2q)	0.2–4	1.3 ± 0.1	0.2–4	2.2 ± 0.3	0.2–4	2.2 ± 0.3	0.2–4	2.3 ± 0.3

Table 9.

Search range and interpreted mean model for Senneterre area of Quebec, Canada.

3.6. Conclusions

It is well known that the local search inversion or the optimization has faster convergence rate. However, while trying to search the global optima, it can be trapped in the local minima. Hence, selection of initial guess is very important for global optimization studies. In case of complicated structure, absence of a priori information also hampers the final solution. However, global optimization methods search the best possible solution and try to find out the exact solution within a multidimensional model space. Moreover, apart from global optimization method, a statistical method can also be applied to find out the ambiguity associated with the final results.

In this study, an effort was made to examine the relevance and adequacy of VFSA on the parameter appraisals from ASA or TG of SP anomalies. In this method, the experimental studies were executed utilizing hypothetically derived data and field data. The elucidation of the amplitude coefficient, location, depth, and shape of a subsurface structure from TG anomaly can be very much established utilizing the present technique. Synthetic data tests were performed utilizing both noise-free and noisy information sets derived from idealized geobodies. The present work reveals, while interpreting every single model parameter (amplitude coefficient, location, depth, shape) together, the VFSA method produces excellent results. Moreover, multiple model bodies were interpreted very efficiently.

The subsequent histogram and cross-plots investigation proposes that the acquired parameters were within the high probability and less uncertainty. This is also supported by the rapid and stable convergence rate of the present inversion method. The viability of this approach has been effectively demonstrated employing noise-free and noisy data. The suitability of this technique for useful application in mineral investigation has been additionally shown on two field cases (Sulleymonkoy anomaly, Ergani, Turkey, Senneterre area of Quebec, Canada). The field data were presumed to be formed due to different idealized geological bodies as mentioned in the published literature. The technique can be utilized to understand numerous structures from the anomaly. The assessed inverse parameters for the field examples were observed to be in fair agreement with other alternate techniques such as PSO.

It is noteworthy to mention that the present work does not get affected by wide search range. Other studies suggested that the search range must be carried out by the workers to understand the effect of search range and the solution.