نشریه پژوهش های ژئوفیزیک کاربردی
سال پنجم شماره 2 (پیاپی 10، پاییز و زمستان 1398)

Seismic data can present a remarkably good image of the subsurface, and therefore, seismic methods have found considerable attention in oil and gas exploration industry. Wavelet or source deconvolution is one of the most important procedures in seismic processing used to increase the time resolution of the seismic sections while requires a reliable wavelet. The accuracy of the wavelet depends severely on the complexity of the wavelet phase. In this paper, through a smooth estimation of the wavelet amplitude spectrum, we go for obtaining the impulse response of the earth via a phase retrieval algorithm. Despite the conventional deconvolution methods, here just the Fourier amplitude spectrum information of the data is inverted as a phase retrieval problem. In the next step, deconvolution of the recorded impulse response from the data leads to a better estimation of the wavelet with any desired phase spectrum. Therefore, the presented algorithm is considered as a ''blind deconvolution'' method.

In statistical seismic deconvolution, the wavelet is estimated from the data; however, it is easier to estimate its amplitude spectrum, and the phase is usually missed or is very inaccurate. This makes deconvolution of mixed-phase wavelets more problematic. It has still remained a challenge among geophysicists to estimate a reasonable seismic wavelet from the data to perform deconvolution efficiently. In this paper, as it is going to be illustrated, the proposed phase retrieval algorithm can be used for deconvolution of mixed-phase wavelets.

A reflection seismogram, after some specified processing steps, can be regarded as a convolution of the source wavelet with the reflectivity series and some additive noise. For obtaining the reflectivity series describing the earth, an appropriate wavelet is needed for deconvolution. The reflectivity model is obtained only by the amplitude spectrum of the observed data, looking for a solution whose predicted amplitude spectrum matches the observations of amplitude spectrum, up to a constant considered for errors in the data. Actually, it refers to an amplitude-only inversion problem, which reconstructs the Fourier phase of a signal from the Fourier amplitude. Obtaining the reflectivity model through the phase retrieval algorithm is an ill-posed inverse problem and has to be solved through regularization method. The problem is solved for a reflectivity series based on the fast iterative shrinkage/thresholding algorithm (FISTA) allowing extra constraints while preserving the computational simplicity. In multichannel deconvolution for improving temporal sparsity while preserving the lateral continuity of the estimation, we define a combined regularization function based on sparsity and second-order total variation. Split Bregman algorithm is used to solve the corresponding proximity function. Then, a wavelet with any desired phase spectrum is estimated using the obtained reflectivity.

In this paper, we proposed a deconvolution algorithm which just needed a smooth approximation of the source wavelet amplitude spectrum. The desired performance of the proposed phase retrieval method on the numerical and field seismic examples confirmed its efficiency by enhancing the resolution of the seismic section and obtaining the accurate reflectivity model. The approach for solving deconvolution discussed here had no limitations for the phase of the extracted wavelet and could obtain wavelets having complex structures with aan acceptable accuracy.

Definition of a suitable cost function is a critical step in solving a problem based on inverse theory, and minimization of the designed cost function can also be very challenging. Non-linear cost functions are usually solved by using iterative algorithms where selection of the parameters in each iteration has a profound effect on the speed of convergence and the quality of the final solution. Therefore, a scheme for automatic determination of the parameters in a general framework can be highly effective in solving inverse problems arising in applied geophysics. Seismic nonstationary deconvolution is a highly ill-conditioned problem, which is required to be solved when improving the vertical resolution of the data is intended. The solution of this problem can be very challenging specifically when a high-resolution solution is desired and when the contaminant noise in the data is non-Gaussian or spike-like. In this paper, we consider this problem assuming that the seismic wavelet and the medium quality factor (Q) are known. Specifically, we consider the minimization of a general cost function for solution of seismic non-stationary deconvolution. The iteratively reweighted least squares (IRLS) algorithm is a common technique for solving this kind of problems in geophysics. However, automatic determination of the regularization parameter in each iteration of IRLS is not an easy task and making the algorithm inefficient. Here, we extend the recently developed method, called iterative reweighted and rrefined least squares (IRRLS) method, for treating seismic deconvolution and propose a new scheme based on the secant method for automatic update of the regularization parameter in each iteration.

According to the convolutional model of the Earth, a seismic signal can be modeled as convolution of the source generated wavelet (w) with the Earth impulse response. The Earth impulse response contains the reflectivity information of the layer boundaries and the elasticity effects of the medium such as attenuation, absorption, etc. The aim of non-stationary seismic deconvolution is to recover information about subsurface from non-stationary seismic signals by solving a multi objective non-linear optimization. Solving this optimization without any prior information about w and attenuation (Q) will be impossible. There are some methods for estimating w and Q from surface seismic data. In the case of known Q and w, the problem changes to a linear optimization with non-linear cost function. The IRLS is a popular method in geophysics for solving this kind of non-linear cost functions but it can be time consuming when there is no information about the noise bound. An alternative IRRLS has been proposed that solves non-linear cost functions iteratively. Furthermore, an automatic algorithm has been developed for updating the regularization parameter at each iteration.

We have proposed a method, based on the IRRLS and secant method of root finding, for high resolution constrained non-stationary deconvolution. Numerical examples from simulated results confirmed that the proposed method is not sensitive to the initial parameters and provide high-resolution estimates of seismic reflectivity series.

The host rock in Chumalu area includes the Eocene volcanic rocks and the Oligocene monzonite intrusive mass that have undergone propylitic-carbonate, silica, and silica-argillic hydrothermal alterations. The hydrothermal activities resulting from the injection of intrusive mass have produced two types of mineralization in the region: the first mineralization type is a silica-shear zone having N70E trend, and the second type is in the form of a bunch of veins with an approximate trend of northwest to southeast, located in the north of the first mineralization type. The mineralization is of the form of vein, veinlet, massive, scattered and replaced, and follows the trend of the faults. The ore contains lead, zinc and fluorine, along with copper, gold and silver, which are formed of chalcopyrite, pyrite, galena, sphalerite, cerussite, calcite and fluorite minerals. Based on the results of the geological studies conducted in this study, the mineralization characteristics of Chumalu area are of the type of lead and zinc hydrothermal deposits. Induced polarization (IP) and resistivity data acquisition were first designed and carried out using rectangular array to determine the anomalous zones, and then, using dipole- dipole array to determine the lateral and depth limits of these anomalies. Then, modeling of IP and resistivity data was done using smooth inversion method. For better representation of the results, all two-dimensional (2D) sections were combined and as a result, a three-dimensional (3D) model was presented. Consequently, an anomaly with an approximate east-west orientation in the south of the area has been extended on the silica-shear zone, and some anomalies with an approximate north-south orientation in the north of the area, in relation to the fractures of the basaltic deposits, which can be linked with metallic mineralization in the area. Anomalous areas with high chargeability and medium to low resistivity values have been determined. High chargeability values for some anomalies in the north of the area are probably due to the presence of scattered pyrites in depth. Finally, 6 points were proposed for drilling.

Application of IP and resistivity geophysical methods and their integration with geological studies is of great importance in the exploration of sulfide metal deposits. The Chumalu area is located 70 km northwest of Zanjan in the Tarom-Hashtjin metallogenic zone, which has possible potential for metal reserves and mineral resources. The purpose of this research is to collect and process raw IP and resistivity data using the smooth inverse modeling, interpretation of the results based on geological studies and finally determining the anomalous limits, depth and thickness of probable mineral masses in the area and at the end, suitable locations for drilling are suggested.

This research is based on field and laboratory studies. In order to conduct geological studies, alteration, mineralization and determination of the genesis of the mineralization, 8 thin sections and 9 polished sections of surface and subsurface samples from the mineralization ore and host rock were prepared and studied. Moreover, field operations of IP and resistivity data acquisition were first made in the form of 4 networks using rectangular array in order to determine the anomalous limits with a current line of 800 meters long, electrode spacing of 20 meters, and the distance between the lines 50 meters. Then, 6 lines using dipole-dipole array were surveyed in order to explore the lateral and depth limits of these anomalies with electrode spacing of 20 meters. Then, modeling of IP and resistivity data was made by 2D smooth inversion method using RES2DINV and ZondRes2D software packages. Furthermore, all 2D sections were combined using RockWorks software and a 3D model was presented. Ultimately, possible mineralization limits were identified and suitable sites for drilling were proposed.

Based on the geological studies, lead and zinc mineralization in olivine basalt volcanic host rock and part of the monzonitic intrusive mass in Chumalu area occurs during hydrothermal processes and in the form of veins along fault structures. Chumalu mineralization is of hydrothermal lead and zinc type. Based on geophysical and geological studies, anomalies have been introduced that have high chargeability and medium to low resistivity, as probable mineralization. Finally, 6 locations were proposed for drilling.

One of the most important aspects of mineral deposit exploration is depth estimation values of the mineral masses. Gravity method is used widespread for detection of mineral deposits. A new approach is presented in order to interpret residual gravity anomalies due to simple geometrically shaped bodies such as horizontal cylinder, vertical cylinder, and sphere. This approach is mainly based on using feed forward modular neural network (MNN) inversion for estimating the shape factor, depth, and the amplitude coefficient. The sigmoid function has been used as the activation function in the MNN inversion. The new approach has been tested first on synthetic data from different models using only one well-trained network. The results of this approach show that the parameter values estimated by the modular inversion are almost identical to the true parameters. Furthermore, noise analysis has been made. The inversion of noisy data produces satisfactory results for the data up to 5% of random noise. The reliability of this approach is demonstrated for real gravity field anomalies taken over a chromite deposit near Sabzevar City, Khorasan Province, Iran.

Forward modeling plays an important role in gravity data interpretation. Gravity data interpretation aims mainly to estimate the depth and location of the causative target. It is known that the gravity data interpretation is non-unique where different subsurface causative targets may yield the same gravity response (anomaly); however, a priori information about the geometry of the causative target may lead to a unique solution (Roy et al., 2000; Aboud et al., 2004). Neural networks (NNs) provide means to build mathematical models that relate input data to desired output data. The neural networks do not know the physics of the forward problem; they have only catalogs of the input/output pairs of the forward mapping that have been fed to it. In this paper, MNN inversion is used mainly to compute the depth and the shape factor of the causative target from a gravity anomaly. NNs can offer a unique solution, especially for noisy data, when acknowledge of a task is not available or unknown nonlinearity between input and output may exist (Bhatt and Helle, 2002; Al-Garni, 2010).

NNs can be considered as universal approximation which can approximate any function in terms of its variables. Generally, a NN is fed by a training set of a group of examples from which it learns to estimate the mapping function described by the example patterns. NNs algorithms may be divided into two main groups, which are supervised (associative) learning and unsupervised (self-organization) learning. The supervised learning is based on desired outputs. During the training, the NN tries to match the outputs with the desired values. In unsupervised learning, the method is not given any target value where the desired output of the network is unknown. During the training, the network performs some kind of data compression such as dimensionality reduction or clustering. The NN inversion that has been used for training is based on the MNN architecture. A MNN is characterized by a series of independent NNs moderated by some intermediary. Each independent NN serves as a module (local expert) and operates on separate inputs to accomplish some subtask of the task that the network wishes to implement (Azam, 2000). The outputs of the modules are mediated by an integrated unit called gating network, which does not permit to feed information back to modules.

NN inversion of gravity data over simple geometric shaped bodies such as sphere, horizontal cylinder, and vertical cylinder has been investigated in this paper. MNN inversion has been used in order to obtain three parameters: shape factor, depth, and amplitude coefficient. This approach has been tested first on synthetic data using only one welltrained network, and then, on a field example taken from Sabzevar area, Iran. The results show the upper and bottom depths of the ore body are about 8 m and 32 m, respectively.

Radio magnetotelluric (RMT) and ground penetrating radar (GPR) are known as the near-surface geophysical methods in groundwater investigations. The RMT method provides information about the variation of the electrical resistivity of 50 m of the uppermost part of the ground. High-resolution structural information can be extracted from the GPR processed sections of the very shallow ground. Combining the obtained data using these two methods lead to valuable results on the identification of near-surface layers and structures. In this study, we propose a new constraint for the two dimensional (2D) inversion of the RMT data. We have investigated a known aquifer located in Heby, Sweden, to assess the constrained inversion results using a joint interpretation approach. RMT and GPR surveys have been carried out along two survey lines having the lengths of 870 m and 550 m, respectively. The results show that thick saturated zones are distinguished quite well either in the joint interpretation results or when using the constrained inversion approach. In such cases, the main problem is to locate the water table in the inverted RMT sections. Imposing smooth regularization in the inversion results turns rather sharp boundaries into the gradual transition zone in the final resistivity models. Thus, using the GPR common-offset (CO) reflections as constraints in the inversion of the RMT can recover the water table as a sharp interface in the RMT inverted model. Thin saturated zone has not been recognized in the RMT sections, due to low resolution of the RMT method. For verification of the results, we have evaluated a synthetic model with similar physical properties to the study area. In such circumstances, the results need to be improved either in the joint interpretation or the constrained inversion approach using CO sections. Hence, harder constraints through our proposed scheme have been incorporated into the inversion routine to detect a thin aquifer and achieve a more realistic model.

The RMT and GPR methods are among the most useful non-invasive methods, which can provide continuous data for groundwater exploration. The RMT method due to its limited range of frequencies (10-250 KHz) has low resolution, especially at very shallow depth, and the GPR method itself suffers from its limited penetration depth. Hence, it seems that combining the modeling results of these two methods leads to a more accurate anomaly definition. Reflection (seismic or GPR) data are usually used as constraints in electromagnetic data inversion. Although all reflectors in seismic and GPR sections are not attributed to the distinct resistivity contrasts, in GPR they are mainly related to the dielectric contrast or may occur due to the thin layers embedded in homogenous geological formations. Thus, we propose an alternative scheme to incorporate interfaces with distinct resistivity contrast in the RMT data inversion.

Using all GPR reflections as constraints in the RMT data inversion may cause some artifacts in the final inverted model. In low clay content formations, such as clean sand and gravel formations, dielectric constant and resistivity are mainly related to the volumetric water content. Therefore, we propose a new structural constraint based on the assumption that the resistivity and water content contrasts occur at the same boundaries. To establish this constraint, we have used common mid-point (CMP) velocity analysis as well as the combination of Topp’s and Archie’s relationships. As a result, an initial resistivity model has been deduced from the CMP velocity analysis that can be used as a priori information in the RMT data inversion.

Thick saturated zones (having thicknesses of more than 10 m) have been distinguished quite well by applying smooth constraint inversion of the RMT data as the joint interpretation of The RMT and GPR data leads to a reasonable outcome in this regard. Although sharp boundaries are mapped as gradual interfaces in the inverted resistivity section of the RMT data, such interfaces are recovered well by incorporating the GPR result as a priori information in the constrained inversion of the RMT data. The water table at a depth of 10 to 20 m, and consequently, the saturated zone is resolved well in this constrained inversion method. It correlates to the borehole log information. On the other hand, thin saturated layers could not be distinguished in the RMT sections due to its low resolution. It means that the water table at a depth of 10 to 15 m is not mainly detected when only the determinant mode data are used. In such areas, the constrained inversion of the RMT data using the water table location deduced from the CO GPR data also fails. However, we have incorporated harder constraints through the model covariance matrix and prior information in our proposed constrained inversion routine. Using this approach, a local thin aquifer has been recognized well. Furthermore, our proposed technique can be used in the inversion of other electric and electromagnetic data.

Edge detection of causative bodies is important in the interpretation of potential field data. There are many methods that can be employed to detect and enhance the edges. In recent years, many filters have been presented based on potential field gradient tensor data. Because of using nine signal components, these methods have a very high accuracy compared to previous proposed methods in this regard. In this paper, normalized horizontal modulus (NHM) method, based on the potential field gradient tensor data, has been proposed. This method makes use of the modulus of potential field gradient tensor to normalize the horizontal modulus of potential field gradient tensor. The minimum value of NHM specifies the edge of magnetized structures. The proposed method with very high degree of precision delineates the edges of anomalies and does not show distraction. This new filter is tested on synthetic data and finally, it has been applied to the aeromagnetic data of the Varzaghan area, and as a result, the location of faults in the area has been determined with high accuracy.

The potential-field gradient tensors are the second derivatives of potential-field data. Since the potential-field gradient tensor data contains nine signal components, their interpretation allows a high resolution and detailed investigation of geological structures. The methods based on the potential field gradient tensor (PGT) matrix use curvature or eigenvalue of the PGT matrix or directional methods of the PGT matrix. In the directional methods of the PGT matrix, in order to provide a more detailed map of the subsurface, Oruç and Keskinsezer (2008) have defined the directional tilt angles filters. Mikhailov et al. (2007) and Beiki (2010) have proposed the directional analytic signal to delineate the edges. However, this technique cannot display the edge of amplitude size of different anomalies simultaneously. Yuan and Yu (2015) have introduced second order directional analytic signal method and then, proposed a normalization method, which can display the large and small amplitude edges simultaneously. Yuan et al. (2016) have proposed horizontal directional theta (ED) method as an edge detection method based on the gravity gradient tensor that can weaken the above defects. It has higher resolution compared to the previous filters for delineating edges. However, for complex geological situations, the mentioned methods have some restrictions for edge detection. In this paper, to overcome these restrictions, NHM method has been proposed, and compared to other methods produces more detailed results.

The NHM is defined as follows: where HM is the horizontal modulus of the potential field gradient tensor and M is the modulus of potential field gradient tensor. The minimum value of NHM specifies the edges of magnetized structures. The NHM method beside tilt angle, total horizontal derivative of the tilt angle (THDR) and ED methods has been applied on synthetic magnetic data, and their results have compared. Moreover, the NHM is applied on aeromagnetic data from Varzeghan area.

The results of applying the tilt angle, THDR, and NHM methods on synthetic magnetic data have shown that the THDR method has been succeeded in determining the position of the bodies, however, the edges of the deeper bodies have not been recognized clearly. The tilt angle method can recognize the edges of the shallow and deep bodies simultaneously but with low accuracy. The ED method possesses high precision and high accuracy in identifying the edges, but in its results, some distortions can be observed. The NHM method of the total intensity data can display the edges of the bodies more accurately while no distortion is seen in the results of the NHM method. Therefore, the NHM filter, compared to other filters, produces more detailed results. The proposed method has been applied on the aeromagnetic data from Varzeghan area, and as a result, the recognized edges of the geological structures are found to be precise and clear.

Assessing a time-frequency representation of signal with an acceptable timefrequency resolution, and for specific purposes in different applied studies, has always been a challenge for signal processing researchers. In case of seismic data, using a time-frequency representation with high resolution will yield a higher precision in processing and interpretational applications of timefrequency analysis of data. In the most of time-frequency analysis methods, a form of smoothing is used for generating time-frequency map, which it causes energy dissipation in time-frequency plane and decreasing the resolution. Reassignment is an efficient technique for compensating this issue and increasing the resolution. It can provide a high time-frequency resolution through moving and concentrating the energy distribution in the time-frequency plane to true location. Reassignment has been applied to various time-frequency analysis methods and its performance has been presented in different researches. In this paper, the reassigned Stransform as a new development on S-transform to provide higher time and frequency resolution is utilized to extract some seismic attributes. The performance of the method in providing an acceptable time-frequency resolution is shown by testing on synthetic non-stationary chirp and seismic signals. As a seismic application, the reassigned S-transform is utilized in studying low frequency shadows through time-frequency analysis of seismic data set acquired on a hydrocarbon reservoir. For this purpose, some time-frequency attributes including single-frequency, instantaneous amplitude, instantaneous dominant frequency and sweetness factor are extracted by this method. The results show that the reassigned S-transform can provide much higher energy concentration rather than standard S-transform, and the events and anomalies can be interpreted with more precision due to their better time and space resolution in attribute sections.

The time-frequency analysis methods are among the most common signal and image processing techniques in different applied fields of electric engineering, mechanical engineering, geoscience, etc. Time-frequency methods are employed in seismic data processing and interpretation applications for denoising, attenuation estimation, deconvolution, hydrocarbon detection, channels and faults visualization and so on.
There are several time-frequency analysis methods. One of the main reasons of developing new time-frequency methods is to reach higher time-frequency resolution. The reassignment is one of the successful approaches in this field. The mission of reassignment method (RM) is to move the energy distribution of the time-frequency plane to true location. Through this way, a precise distribution of instantaneous frequency has been provided for any time sample. Reassigning is also carried out in time direction. The RM has been applied in several time-frequency methods such as wavelet transform, Wigner-Ville distribution, Gabor transform and S-transform. In this paper, the reassigned Stransform has been studied in seismic data time-frequency analysis. The method has been utilized for detection of low frequency shadows in a seismic dataset to locate probable gas reservoir.

In this paper, the performance of reassigned S-transform has been studied by its application on synthetic chirp signal and seismic trace. The results show that the method is capable of providing a well-concentrated time-frequency maps. As an application in real seismic data, the method has been utilized for studying the low frequency shadows related to probable gas bearing zones. This approach extracts some attributes including single frequency, instantaneous amplitude, instantaneous dominant frequency and sweetness factor, through time-frequency analysis of the data. The results show that the reassigned S-transform can provide higher time and space resolution, and thus, the events and anomalies can be interpreted more precise compared to standard S-transform results.

In this paper, we examine the applicability of a pseudo two-dimensional (2D) inversion technique, called 2D-DOInv, on a layered earth model, to a frequency-domain helicopter-borne electromagnetic (FDHEM) dataset. In this scheme, a one-dimensional (1D) inversion is modified with 2D Occam’s smoothness constraints between 1D models of adjacent sites in addition to the vertical smoothing. By incorporating the vertical and horizontal weighting factors in the regularization matrix, we obtain more stable solution and geologically more realistic results. Unlike 1D Occam’s inversion, in the 2DDOInv algorithm, the data of all stations along a flight line are simultaneously inverted by minimizing a common objective function. In this inversion algorithm, we are able to incorporate the inequality constraints. The inversion scheme can be parallelized using multiple processors in a single computer. To validate the algorithm, we consider synthetic responses generated over known 2D targets, a buried valley structure and a two-layer earth containing heterogeneous overburden. In comparison to the 1D Occam’s inversion, the 2D-DOInv algorithm estimates pseudo 2D cross section of subsurface resistivity structure, and efficiently reduces the effects of the multidimensional modeling cost and data noise. Finally, this inversion is applied on the real data in Kalat-e-Reshm area of Semnan Province, Iran. The resulting inverted parameters using the proposed algorithm correspond reasonably close to the known geology and to the results from electrical resistivity tomography (ERT) data inversion.

The FDHEM applications are industrially feasible, as long as there is a fairly fast algorithm, yet accurate enough for inversion of a tremendous amount of survey data to model near surface resistivity variations. Currently, the only way to invert such large amounts of data in the field is by utilizing a 1D approximation. The 1D inversion often leads to a series of stitched 1D models that can be misleading, since real data inherently consist information about more than one dimension. Unfortunately, multidimensional modelling has been limited by the computational expense. In this paper, we present a 2D smooth regularized inversion based on a 1D forward modelling with a particular pair of regularization coefficients that simultaneously guarantees the stability and best-fit criteria. All data and 1D models along a flight line are inverted together, giving pseudo 2D sections of the subsurface.

Here, the forward modeling routine is based on the 1D assumption in which the thicknesses of layers are fixed and increase logarithmically with depth. To implement the inversion procedure of 2D-DOInv, the 2D model is discretized with M by N blocks of constant resistivity, where the parameter vector is of length M, the data vector is of length N. In this inversion, the model regularization function is improved to smooth the resistivity model by incorporating separate vertical and horizontal smoothing factors in the regularization matrix R. To further stabilize the inverse problem, we introduce the inequality constraint using the transformation of the model parameter vector. The algorithm can parallelize computations using multiple processors on the same computer to make the inversion scheme very large scale and quick. 1D and 2D damped Occam’s inversions, which are developed in MATLAB environment, are first applied on the synthetic data sets obtained from two 2D models. Compared to the 1D damped Occam’s inversion, 2D-DOInv recovers fairly accurate 2D synthetic models from noisy data. Finally, the 2D-DOInv approach is applied on the field data.

In this paper, we have developed a pseudo 2D inversion in MATLAB environment. This inversion technique, named 2D-DOInv, is used to invert FDHEM data. 2D constrained inversion of both synthetic and field data certainly improves the 1D Occam’s inversion results of quasi-layered earth structures, although the misfit is higher. Even in case of noisy data, we can filter out the influence of the noise using our smooth regularized inversion technique, and enhance the resolution of the subsurface resistivity images. Thus, this technique can suppress both the measurements and the processes errors. For non-layered earth structures, e.g. the two-layer model comprising of a heterogeneous overburden that is common in weathered crystalline terrains; however, the reconstructed 2D pseudo section will be blurred as a result of excessive horizontal smoothing. Through this method in which the horizontal constraint is chosen to be relatively weak in the top layer, 2D-DOInv still improves the inversion results and mitigates 2D/3D effects. The example of measured data sets from a DIGHEM survey carried out over Kalat-e-Reshm area demonstrates the capabilities of the algorithm reasonably. We can see reasonable correlations between the pseudo 2D resistivity models from FDHEM data inversion and ERT 2D resistivity models.

Edge detection is a fast and qualitative interpretation method to achieve information from potential field (e.g. gravity) anomalies. In the edge detection methods, the separation of overlapping amplitudes of anomalies and accuracy of edge detection are very important. There are various methods for edge detection. Most of these methods are based on the gradients of the potential field data. The gradients are sensitive to noise. Statistical methods have been used to increase the accuracy of edge detection. Normalized standard deviations (NSTD) and correlation coefficient of multidirectional standard deviations (CCMS) are among these methods.

There are several edge detection methods based on gradients of data. Each of these methods has some strengths and some weaknesses. In the selection of these methods for a particular case, simplicity and better performance are considered. These methods include: total horizontal derivative (THD), Theta angle, Tilt angel, hyperbolic tilt angle (HTA) and a new method based on the gradients, called normalized total horizontal derivative (NTHD). In addition, the semi-statistical method of NSTD and statistical method of CCMS are among these methods that have been explained in this paper. The NSTD method is obtained from the standard deviation of the gradients, however, the CCMS method does not use the gradients. This method is completely a statistical method, which is based on correlation coefficient and standard deviation.

In this paper; after examining the above-stated edge detection methods, they have been applied on both synthetic and real data. The performances of these methods are compared in the presence of noisy data, overlapping amplitudes of anomalies and their accuracies in edge detection.

The results of applying the above-stated edge detection methods on the synthetic data show that the gradient-based edge detection methods are sensitive to noise, depths of anomalies and overlapping amplitudes of anomalies. The NTHD, NSTD and CCMS methods are less sensitive to noise than the other edge detection methods. These methods detect anomalies with different depths and separate anomalies with overlapping amplitudes. In all of these methods, as the depths of anomalies increase, the accuracy of edge detection decrease. This study show that the CCMS method has the best result when applied on the synthetic data. Furthermore, applying the CCMS method on the real data yields better results in comparison with the other edge detection methods. The results of edge detection by this method have been shown on the bouguer map. Thus, this method reduces complexities of edge detection that can be useful for the interpreter.

During wet mineral processes in mines, a large amount of wastewater is produced and stored in tailing dams. In Miduk copper mine, Shahr-e-Babak, Kerman, Iran, there are tailing dams and wastewater dams as a result of mineral processing. The wastewater dam normally is constructed in order to separate and to reserve the water content of tailings, and finally, to reuse it for mineral processing. In this study, we have used the geophysical methods (electrical resistivity tomography (ERT) and self-potential (SP)) to investigate the seepage from the wastewater dam of the Miduk copper mine. The ERT surveys show a low electrical resistivity zone (less than 20 Ωm). It can be as a result of a clayey alteration zone or fractured zone. In addition, SP measurements show a positive anomaly in the same area due to the seepage from the dam walls. It can be due to an electrokinetic source origin. The observed natural source in the same zone, confirmed the geological, hydrogeological and geophysical interpretations. These results show the ability of ERT and SP methods to investigate the seepage zones related to wastewater dam walls.

Mining like other industrial activities can affect the environment negatively. Tailing dams have long been associated with mining activities and have formed a major negative impact. Excessive seepage in the foundation of the dam threatens the integrity of the structure. Geophysics has been used intensively in environmental problems, in particular in the study of mine tailings problems. In fact, geophysical surveys especially electrical resistivity tomography (ERT) and self-potential (SP) methods constitute a comprehensive methodology for assessment of anomalous seepage conditions by detecting, mapping and monitoring the zones.

The field survey to detect the seepage in Miduk wastewater dam was composed of a geophysical survey and a hydrogeological study. The hydrological study was performed in order to identify water table, surface water and groundwater quality, and to understand the relation between water collected in dam and water resources. Such a study is very common in environmental geochemical investigations. The ERT survey was carried out along 2 survey lines, called P1 and P2, using the dipole–dipole array with an electrode spacing of 20 m in the survey line P1 that has a total length of 400 m, however, the survey along the survey line P2 has been carried out using an electrode spacing of 5 m. The survey line P2 is 100 m long. A maximum n value of 11 was used for both P1 and P2 ERT profiles. All data were acquired by a WDDS-2 resistivity meter using multi-electrodes cable. The SP survey along SP1 and SP2 survey lines was carried out on the northeast abutment of the dam. The SP measurements were taken using two non-polarizing electrodes (Cu-CuSO4), one electrode was located at the base station as a stationary electrode and the other moving along the desired line at pre-fixed stations. The base-station was chosen at a point convenient for operation but away from the expected anomaly. Electrodes polarization was controlled between measurements using two constant points in the site. In order to identify the water table, in a part of the hydrological study, water depth in 6 observation wells, located downstream of the water retention dam, was measured. In the second part of the hydrological study, hydrochemical investigation was carried out by taking 17 water samples collected from 15 stations. The hydrogeochemical characteristics of water were obtained through physicochemical analysis of the water samples.

The ERT survey shows a low electrical resistivity zone (less than 20 Ωm). It can be the result of a clayey alteration zone or fractured zone. The SP measurements show a positive anomaly in the same area due to the seepage from the dam wall. According to the hydrochemical study, the water type of this area is bicarbonate-chlorate while the water type of one borehole and the dam water sample are both sulfate-chlorate. The obtained information from the survey area confirms the geological, the hydrogeological and the geophysical interpretations