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

In the present study, analytic signal amplitude (ASA) or total gradient (TG) inversion of self-potential anomalies has been carried out using very fast simulated annealing (VFSA) global optimization technique. The results of VFSA optimization demonstrate the applica- tion and efficacy of the proposed method for idealized synthetic hypothetical models and real single and multiple geological structures. The model parameters deciphered here are the amplitude coefficient ( k ), horizontal location ( x 0 ), depth of the body ( z ), and shape ( q ). Inversion of the model parameter suggests that constraining the horizontal location and the shape factor offers the most reliable results. Investigation of convergence rate, histogram, and cross-plot examination suggest that the interpretation method developed for the self- potential anomalies is stable and the model parameters are within the estimated ambiguity. Inversion of synthetic noise-free and noise-corrupted data for single structures and mul- tiple structures in addition to real field information exhibits the viability of the method. The model parameters estimated by the present technique were in good agreement with the real parameters. The method has been used to invert two field examples (Sulleymonkoy anomaly, Ergani, Turkey, Senneterre area of Quebec, Canada) with application of subsur - face mineralized bodies. This technique can be very much helpful for mineral or ore bodies investigation of idealized geobodies buried within the shallow and deeper subsurface.


Introduction
The self-potential (SP) technique has an important significance in mineral and ore explorations [1][2][3][4][5][6][7][8]. The method has an extensive range of application, viz., mining industries [9][10][11][12], sulfide, graphite exploration, and groundwater exploration [13,14], study of groundwater The main objective of the present work is to interpret the analytic signal amplitude (ASA) or the total gradient (TG) of SP anomaly over different idealized or causative body, which satisfies the Laplace's condition. It is well known that the analytical signals derived from SP anomaly are correctly known as the TG [71][72][73][74]. The interpretation of ASA or TG derived from SP has been sparsely done in the present literature. It was only interpreted using ant colony optimization (ACO) [75]. In the current work, very fast simulated annealing (VFSA) was applied to decide the different model parameters associated with idealized subsurface bodies for ASA or TG from SP anomalies. The present inversion method has an advantage over other approaches for its pliability and its proficiency to converge toward global optima. The method has an ability to avoid from getting stuck in local minima and it has very high resolution, faster calculation as well as less memory without negotiating the resolution of the model parameters [2,3,35,[76][77][78][79][80]. Moreover, it does not need a priori information for the elucidation of SP anomalies.
In this study, inversion using VFSA algorithm was performed with the help of synthetic noise-free and noise-corrupted synthetic data using single and multiple structures and two field data from Sulleymonkoy anomaly, Ergani, Turkey, Senneterre area of Quebec, Canada. The results from the present method were compared with other well-established SP anomaly interpretations, such as ACO techniques. The present method used as a comprehensive method for quantitative elucidation of SP anomalies derived after various subsurface idealized geobodies.

Forward modeling
Following Nabighian [71], the 2-D analytical signal amplitude (ASA) or the total gradient (TG) of SP anomaly is given as. (1) where V x and V y are 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: (2) where k is the amplitude coefficient/factor related to the physical properties of the source, z is the depth from the surface to the top of the body (sphere, cylinder, sheet, line of poles, point pole), x 0 (i = 1, …,N) is the horizontal position coordinate on the surface, and q is 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]: where ASA j (x i ) is the SP anomaly at x i location for jth body and M is the number of bodies.

Global optimization
Several global optimization approaches have been effectively used in interpretation of different geophysical data (e.g., [35-37, 75, 77, 82-86]). In the present study, an alternative method

Minerals
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, and are the ith observed and model responses and and are the maximum and minimum values of the observed response, respectively.
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 6 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.

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 0 ) 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.

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 0 , 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 −4 . Then, a statistical mean model was calculated utilizing models that have misfit lesser than 10 −4 , 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 −2 . 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.

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 0 , 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 3b shows 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 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 0 , 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.  demonstrates 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 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 0 , 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 Figure 4. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (10% random) data for sphere.  Table 1. Actual model parameters, search range, and interpreted mean model for noise-free, 10% random noise with uncertainty for sphere.  (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (20% random) data for sphere.   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.

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  Table 3. Actual model parameters, search range, and interpreted mean model for noise-free, 20% random noise with uncertainty for horizontal/vertical cylinder. Figure 8. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (10% random) data for cylinder.
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 cylinderlike 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 Figure 9. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (20% random) data for cylinder.   corrupted in the models mentioned above. Figure 11c and d shows the histogram for noisecorrupted 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.

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 Figure 12. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (10% random) data for sheet.  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.

Model parameters
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) Figure 14. Fittings between the observed and model responses for sheet-(a) noise-free synthetic data; (b) noise-free synthetic data; (c) 10% random noisy synthetic data; and (d) 20% random noisy synthetic data. Figure 13. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −4 for noise-free data) (green), and models with PDF > 60.65% (red) for noise-free data, (b) cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisy data) (green), and models with PDF > 60.65% (red) for noisy (20% random) data for sheet.
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 noisecorrupted data. Histogram and cross-plots show alike as revealed in other examples for single structures. However, it is not presented here for brevity.

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  Table 6. Actual model parameters, search range, and interpreted mean model for noise-free, 20% random noise with uncertainty for sheet. Figure 15. Fittings between the observed and model responses for multiple structure-(a) noise-free synthetic data; (b) 20% random noisy synthetic data.

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:

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  Table 7. Actual model parameters, search range, and interpreted mean model for noise-free and 20% random noise with uncertainty for multiple structure.
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 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. Figure 17. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noise-free data) (green), and models with PDF > 60.65% (red) for field data-single body, (b-d) crossplots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noise-free data) (green), and models with PDF > 60.65% (red) for field data-multiple bodies-Sulleymonkoy anomaly, Ergani, Turkey.  Inversion of Amplitude from the 2-D Analytic Signal of Self-Potential Anomalies 21

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.

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. Figure 20. (a) Cross-plots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noise-free data) (green), and models with PDF > 60.65% (red) for field data-single body, (b-d) crossplots between amplitude coefficient (k), depth (z), shape location (x 0 ) for all models having misfit < threshold (10 −2 for noisefree data) (green), and models with PDF > 60.65% (red) for field data-multiple bodies-Senneterre area of Quebec, Canada.  2.3 ± 0.3 Table 9. Search range and interpreted mean model for Senneterre area of Quebec, Canada.

Model parameters
Inversion of Amplitude from the 2-D Analytic Signal of Self-Potential Anomalies 25 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.