فهرست مطالب وحید ابراهیم زاده اردستانی

In recent decades, the study and recognition of liniments has become an important issue in the field of geosciences. Liniments consider could be a schedule hone for the interpretation of gravity information. Besides, it is imperative for a wide extend of topographical information. Be that as it may, elucidation results of gravity information are found to be exceedingly variable among interpreter and need consistency indeed inside a person. In this manner, modern strategies have been displayed to progress the unwavering quality of basic elucidation, and these methods help interpreters to reach similar results from the same gravity data.      Potential field methods play a fundamental role in geophysical explorations. One of the main goals in the interpretation of potential field data is to determine the location and estimate the depth of magnetic and gravity anomalies. Quantitative interpretation methods of potential field data, such as standard Euler, have always been modified in order to increase the accuracy of determining the characteristics of subsurface sources, and generalizations have been made in order to improve the reliability of the results. Standard Euler method is based on choosing the dimensions of the window and depends on the structural index. Therefore, defining a window with suitable dimensions and moving it within the data grid or profile and choosing the appropriate structural index can provide the results of this method with higher accuracy. Since the lack of accurate determination of the structural index can lead to wrong results in depth estimation, Euler's generalizations are presented to remove the structural index from the calculation process. Quality of the field data poses great influence on the Euler inversion solutions. If the data has low signal-to-noise ratio, the computational process will be masked. This issue makes it difficult to outline boundaries of the causative sources.      In this research, in order to determine the location and estimate the depth of the gravity anomaly, a generalization of the standard Euler approach using the singular value decomposition method and the ratio of the horizontal gradient to the analytical signal have been used. Also, the quantity of phase congruency was applied to the data in order to interpret the results more accurately. Compared to the standard Euler, the results show more compliance with the anomaly boundaries. Moreover, the depth estimation interval (histogram) in the improved Euler is smaller than the standard Euler. For further investigation and increasing the credibility of the interpretation, phase congruency was applied to synthetic and real data aiming at more accurately determining the location of the subsurface anomalies.

In standard gravimetric correction methods, after the raw gravity data sets were corrected for drift, tide, latitude, and free-air effects to obtain free air anomalies, the effect of the mass between the reference surface and ground surface is eliminated in two steps including Bouguer and terrain corrections. But this study removes this effect in one step through the forward modeling method. To do this, two things are necessary for finding more accurate answers. First, how is the underground discretization, and to what extent a network of Digital Terrain Model (DTM) is available? Quad tree mesh accessible in Simulation and Parameter Estimation in Geophysics (SimPEG) is a very accurate and advanced meshing algorithm to discretize subsurface based on our requirements. This meshing system can choose the size of cells in the desired locations. Hence, using this flexible discretization, it is possible to define the smaller cells in borders, near the topographic region, which helps a for more precise answers. Having a dense DTM, the SRTM GeoTiff pictures are downloaded from USGS Earth explorer with 1 arc-second (90 m) resolution (https://doi.org/10.5066/F7PR7TFT), and then height information is extracted from these pictures through GeoToolkit (http://toolkit.geosci.xyz) script. Assuming a flat geoid for our study area, topography is extracted from the SRTM and the pictures are interpolated to estimate the elevation at the gravity observation points. The gravity effect of the model space (the space between the reference surface and topography) is computed via numerical forward modeling assuming a constant density (2.67 gr/cm3). This procedure is done by the Simulation module in SimPEG and is considered as the Bouguer and terrain corrections simultaneously. These corrections are subtracted from the free-air anomalies, which yields the complete Bouguer anomaly. This method is powerful in contrast to other standard methods. In standard methods, Bouguer correction considers Bouguer slab approximation. Therefore, accuracy is lost. Also, in large-scale problems, curvature correction becomes necessary. Also, terrain correction for removing the effects of the mass between the lowlands and heights of the region is inevitable. Terrain correction considers two approximations. First, it uses average height. Hence this procedure has a low precision. Secondly it divides the surrounding area into three zones (near, middle, and far) and computes the effects of middle and far zones with lower precision. Therefore, it decreases the accuracy of the results. The mentioned method is tested on 399 ground gravity data with a grid spacing of about 5 km prepared by the National Cartographical Center of Iran (NCC) in an area of about 200 km in 200 km located in parts of Central Zagros and Central Iran. The results obtained from this one-stage correction method are more accurate and less complicated in doing compared to the results of the usual procedure. Because in this method, we have no simplifying assumptions such as infinite Bouguer slab in Bouguer correction or using relative heights in terrain correction that exist in standard methods.

The flexural strength and consequently the elastic thickness (Te) of the lithosphere are often calculated utilizing the statistical relationships between gravity anomalies and topography data. The elastic thickness of the lithosphere is one of the most important parameters in determining the geological properties and characteristics of the crust as well as its rheological behavior. This prospective study was designed to investigate the use of the gravity and topography data in order to calculate the two fundamental functions, i.e. admittance and coherence, which are the basic functions for determining the elastic thickness of the lithosphere. In this research, the theory of the two objective functions is modeled. The study offers some important insights into the involved parameter and how they affected the functions. Additionally, the results show the impact of initial assumptions of different parameters on the retrieved value of the Te parameter. Taking into account the initial assumptions about parameters such as the elastic thickness () of the plate, the loading ratio () and the degree of correlation () between the surface and subsurface loading regimes applied to the plate and using spectral methods, the value of these two functions are calculated. After examining the theoretical models, the real-data analysis is conducted in order to examine the observed values and compare them with the theoretical values. Based on what was mentioned in the following, this function can be able to retrieve the elastic thickness of the lithosphere. We concluded that:One of the more significant findings to emerge from this study, by examining the theoretical models, is that the best function to determine the elastic thickness of the lithosphere is to use the Bouguer coherence function. A comparison of the results reveals that the characteristic curves of the admittance function will change drastically with the change of the involved parameters. Consequently, having complete knowledge of the structure of the region and involved parameters is essential while using free air admittance in order to determining Te. As it was illustrated in this research, the characteristic curve of the coherence function can be used to determine Te based on the roll over wavelength. For the reason that there was no significant change in the shape of characteristic curve of the functions with the initial value of the loading ratio (), this parameter had less impact on the inversion process. Having considered rare correlation between free air gravity and topography, it is reasonable to accept that the free air coherence method is less applicable in determining Te than Bouguer coherence. The findings from previous studies provide support for the key arguments.

One of the most important geophysical problems in exploration of mineral deposits is to estimate the depth of buried structure using observed gravity data. In this paper, we are trying to estimate mass anomaly depth by using one of the intelligence methods as Particle Swarm Optimization (PSO) with simple shapes as sphere, horizontal and vertical cylinder. In this modeling, two parameters of depth (z) and shape factor (q) were considered as particles and the maximum and minimum depth were used as the prior information. The method is tested for synthetic models with random noise. The method gives precise results for synthetic models contaminated with random noise which is quite acceptable and promising. This technique was also successfully applied to real data for mineral exploration. The applied real data belongs to an area with hilly topography located in the Fars province close to the Abadeh city where the barite deposit is under exploration. The method is used for a profile of real data that is provided from the residual anomalies and passed from the main detected positive anomaly in the area. The estimated depth from this method was 8.15 m which was in good agreement with the results obtained through Euler method and also outcrop‐scale observations.

The Moho discontinuity is a boundary between the crust and upper mantle that reveals the difference between them with changes in seismic velocity, density, chemical structure, and constituents. Estimating the Moho depth and studying its lateral changes is one of the important goals of geophysical studies. The current study aims to estimate the depth and topography of the Moho discontinuity in the southwestern part of the Baltic Sea, including parts of the central European system, the Trans-European Suture Zone, Caledonian Crustal Suture, and the Ringkobing-Fyn High. This area has been one of the most attractive regions for Geoscientists in the last decades due to its complicated geological structures caused by different tectonic events. For this purpose, a three-dimensional model of the crustal structures based on gravity data forward modeling in the study area has been presented. Previous seismic / non-seismic results have been used to constrain the model and reduce its degree of freedom. This model includes sedimentary sequences, crustal thickness, Moho topography, and high-velocity lower crust expansion in the region and shows the tectonic structures of the study area. This study used a combination of marine, land, and EGM2008 gravity data and modeled them with IGMAS+, Interactive Gravity and Magnetic Application System. The interactive modeling program allows the user to change the geometry as well as the density and susceptibility of the primary model and observe results quickly during the processing. In the software, the model structure could be be more user friendly by eliminating additional details and dividing the whole model into vertical sections. Our primary model consists of three main layers of sediments, crust and upper mantle. The sedimentary layer is divided into two major parts, pre-Permian and post-Carboniferous. Also, the crustal layer is divided into the upper crust and the high-density lower crust. Besides, the upper crust is composed of the upper crust of the Baltica and the upper crust of Avalonia. The last layer of the model is a part of the upper mantle. The model space consists of 16 vertical planes stretching 385 kilometers east-west with an equal distance of 15 kilometers, covering the entire study area. The initial model was developed based on seismic sections and previous models, and it has been improved using interactive forward modeling of gravity data. The result shows a good agreement between the measured and modeled Bouguer anomaly, and the Root Mean Square Error of the model is 1.12 mGal. The model correlates clearly with major tectonic units. It indicates that the Caledonian collision resulted in the amalgamation of Baltica and Avalonia is the most prominent tectonic event in the area, and the Caledonian crustal suture between them is interpreted from changes in physical parameters at crustal levels. There is a relatively thick crystalline crust in the area, and the depth of Moho discontinuity varies from 26 to 42 km. The results also indicate that the transition from the Paleozoic crust of the Central European Basin to the Precambrian crust of the Eastern European Craton occurs within the Tornquist Zone.

Zagros orogeny is one of the most active orogenic belts among the mountain ranges extending approximately 2000 kilometers from the Anatolian fault in eastern Turkey to the Minab fault in southern Iran. Concerning the importance of this region as well as the essential role of elastic thickness in controlling the rate of deformation under applied loads, determination of Te in Zagros Fold and Thrust belt has been conducted. The lithosphere's elastic thickness (Te) is a convenient measure of the flexural rigidity, which is defined as the resistance to bending under applied loads. To determine the elastic thickness of the lithosphere, the spectral admittance function is applied. We applied the load deconvolution of the admittance function between free-air gravity and topography data for estimation of Te. The Free air anomalies with a five arc-minute resolution are utilized in this study. In flexural isostatic studies, the gravity and topography data are compared with theoretical models to estimate several parameters of the lithosphere. In the simplest model, a plate has been flexed by a surface load, with the magnitude of the resulting deflection, which is governed by Te. Using the random fractal surfaces as the initial surface and subsurface loads applying at lithosphere, the lithosphere is modeled, and the post flexural gravity and topography are determined. Based on these new fields, the predicted admittance function is determined. Finally, the best-fitting Te is one that minimized the misfit between the observed and predicted functions. Additionally, the weighted misfit by the jackknife error is applied to estimate the observed admittance. The accuracy of the method is checked through synthetic modeling. Two fractal surfaces are used as the two initial surface and subsurface loads applied to the lithosphere. After calculating the corresponding gravity and topography data by the load deconvolution method, the observed and predicted admittance are estimated. The best-fitting Te will be obtained by minimizing the misfit between observed and predicted functions. After confirming the accuracy of the method in Te determination, the technique will be applied to the real data acquired from the NCC as follow. We consider a three-layered crust during the lithosphere modeling on which the internal loading is applied on the middle crust. To model the lithosphere, the global CRUST 1.0 is applied by treating the lithosphere as a three-layer crust. The 2D map of Te variations in the target area is depicted by utilizing the load deconvolution of the admittance function between free-air gravity and topography data. High-precision ground gravity data, which is more accurate than satellite data, allows us to detect more details on Te variations in the region. Based on the obtained results, the estimated range of Te in the survey region can be considered low to intermediate. This predicted range is in good accordance with the area's geology background as it is regarded as a young, active orogeny system. Te range and hence the lithosphere's predicted strength to deformation is supported by the previous studies using different geophysical and seismological studies. The mean value of Te in the area is 37±2 km. The maximum amount is detected in the Sanandaj-Sirjan zone. The overall predicted trend of Te follows the geological background of the region. Additionally, the estimated trend for Te and the strength to the applied load and deformation is in good agreement with the previous geophysical and seismological studies conducted in the region.

Gravity surveying is applied for studying geological structures, for example, basement topography underneath the sediment loads. In potential areas for hydrocarbon and groundwater resources, depth of basement can be estimated using different optimization methods, including stochastic global optimization algorithms. These methods include many functions call of the forward function, so usual forward approaches that discrete the sediment volume into a set of right rectangular prisms need too much computational time. This can be controversial issue while implementing three-dimensional stochastic inversion. In this study, 3D Cauchy-type integral as a fast forward function is applied to accelerate the gravity inversion for 3D determination of the depth to basement. Cai and Zhdanov (2015a,b) introduced this effective approach for potential fields modeling. This method in modeling the sediment-basement interface not only replaces all prisms of conventional volume approach with a gridded basement, but also uses simple mathematical terms in comparison with customary prismatic methods which include trigonometric and logarithmic expressions. Synthetic forward modeling of both of our realistic basin models assesses the validity of the forward operator. Evaluation time for one of the model basins based on the Cauchy-type integral in comparison with the prismatic method which was carried out by two different techniques of forward modeling, is 15 and 50 order lower. Implementing genetic algorithm on the gravity data, the depth of the basement was recovered. The misfit of our data achieved by the algorithm with initial population equal to 10 times of total number of parameters and carrying 700 generations, was lower than 2 mGal. Optimal values were obtained as 80% and 20% for crossover and mutation, respectively. In addition, due to the non-uniqueness of the gravity problem, the genetic algorithm uses a smoothing constraint. By fixing the optimal parameters of genetic algorithm, the optimization process is repeated to find the optimal value for the smoothing factor yielding the most accurate model based on the RMS of the reconstructed model. Results show that a smoothing factor between 0.005-0.015, reconstructs stable solutions. Besides, applying a Gaussian filter, a smoothing filter with the kernel size equal to 11×11 to the calculated depths, achieves more stable evaluations. Noisy synthetic and noise-free gravity data were inverted for one symmetric basin and the algorithm has been able to successfully reconstruct the basement. The case study area is the Aman-Abad alluvial plain (Iran) which its main parts are located in the Sanandaj-Sirjan zone in the Zagros Mountains of Iran. The suitable parameters of the genetic algorithm are found by synthetic tests to invert real gravity data to image the interface of the impermeable layer groundwater. The most common polynomial regression, i.e., degree 1 is applied to calculate residual gravity anomaly. Reconstructed depths from residual gravity anomaly match properly with gravity anomaly trend. Deep parts of the basement (as impermeable surface) have been estimated about 150m which it looks promising for groundwater resources. According to the previous gravity studies, the calculated maximum thickness of sediment is lower than 200 m and the well data specified depth of the basement is 140 m.

In this paper, three-dimensional (3D) inversion of gravity data using graph theory is used. The methodology was initially introduced by Bijani et al. (2015) and, here, we provide more details for the steps and required parameters of the algorithm. An ensemble of simple point masses are used to model a homogenous subsurface body. Then, in the presented inversion methodology, the model parameters are the Cartesian coordinates of point masses and their total mass. Consequently, the algorithm is able to reconstruct the skeleton of the subsurface body and to yield its total mass. Here, the set of point masses is associated to the vertices of a weighted full graph in which the weights are computed by the Euclidean distances separating vertices in pairs. Then, the Kruskal’s algorithm can be used to solve the Minimum Spanning Tree (MST) problem for the graph. A stabilizer, called equidistance function, is obtained using the MST, which computes the statistical variance of the distances among point masses. The function restricts the spatial distribution of points, and suggests a homogeneous distribution for the point masses in the subsurface. Here, a non-linear global objective function for the model parameters comprising data misfit term and equidistance function with balancing provided by a regularization parameter that should be minimized. A genetic algorithm (GA) is used for the minimization of the objective function. GA consists of a random search algorithm based on the mechanism of natural selection and natural genetics. Then, to solve the optimization problem in our algorithm, there is no need to calculate the derivatives of the objective function with respect to model parameters, or any matrix operation. Simulations for two synthetic examples, including a vertical and a dipping dike, demonstrate the efficiency and effectiveness of the implementation of the present algorithm. The skeleton and total mass of the bodies are estimated very accurately. We also show that although the search limits for the model parameters must be used, they are not very limitative. Even with less realistic bounds, acceptable approximations of the body are still obtained. Unlike Bijani et al. (2015) which used the L-curve method for estimating the regularization parameter, here, we present a new strategy to approximate the parameter. We demonstrate that if: 1. the equidistance function converges almost monotonically to zero with increasing numbers of generation; 2. minimum of the objective function at the final iteration becomes small; and 3. the predicted data by the reconstructed model is approximately close to observed data, then, the selected regularization parameter is nearly optimum and the results are reliable. This provides a suitable and inexpensive methodology for estimating the regularization parameter. The method is tested on gravity data from the Mobrun ore body, north east of Noranda, Quebec, Canada. The anomaly is associated with a massive body of base metal sulfide, mainly pyrite, which has displaced volcanic rocks of middle Precambrian age (Grant and West, 1965). With application of the algorithm, a skeleton of the body is obtained which extends about 350 m in the east direction, and shows a maximum extension of 200 m in depth.

In this study, Simulated Annealing algorithm is applied on Rayleigh wave group velocities to image the density variations and shear and compressional wave velocities structure of the crust and upper-mantle of Makran subduction zone. Based on previous studies, surface wave dispersion measurements are primarily sensitive to seismic shear wave velocities. However, it has been proved that the sensitivity to compressional wave velocity is significantly smaller than the sensitivity to shear wave velocity. Also the sensitivity function for the density is smaller than the one for the shear wave velocity. Therefore, shear wave velocity variations are mainly the model parameters in surface wave dispersion analysis. Simulated Annealing is a probabilistic technique for finding the global optimum of a given function. It is especially useful to approximate global optimization in a large search space. The Simulated Annealing method like the Monte-Carlo method, samples the whole model space and can avoid getting stuck in local minima.To evaluate calculation efficiency and effectiveness of Simulated Annealing algorithm, two noise-free and two noisy (10% of white Gaussian noise) synthetic data sets are firstly inverted. Then, a real data from Makran region is inverted to examine the applicability and robustness of the proposed approach on real surface wave data.
In next step, gravity data inversion was applied with a priori information based on surface wave analysis results to obtain Moho depth variations and crustal density structure. The reason for using gravity data set is that surface waves group velocity is sensitive to average velocity variations and has a good lateral sensitivity, whereas gravity anomaly is sensitive to depth variations of discontinuities and has a good vertical resolution.
Our results show that the Moho depth across the Makran subduction zone increases from the Oman seafloor and Makran forearc setting to the volcanic arc. Generally, the crust in the western Makran is thicker than the eastern part and the maximum crustal thickness in the Makran region reaches 46 to 48 km below the Taftan-Bazman volcanos. The Moho map clearly depicts the western edge of the Makran subduction zone, where the Minab fault (representing the eastern edge of the Hormuz Straits) marks the boundary between the thick continental crust of the Arabian plate and the thin oceanic crust of the Oman Sea. Our results show clearly that the high-velocity slab of the Arabian plate subducts northwards beneath the low-velocity overriding lithosphere of Lut block in the western Makran and Helmand block in the eastern Makran.

Summary: First and second vertical gradients are widely used in the interpretation of gravity data. Vertical gradients are sensitive to noise. Accuracy of vertical gradient calculation directly effects the accuracy of interoperations. Therefore, accurate and without noise calculation of vertical gradient is vital. The most common method to calculate vertical gradients is to use Fourier transform. Low noise in the gravity data causes that the vertical gradients, calculated by Fourier transform, have severe noise. In this research, we have used discrete cosine transform (DCT) to calculate vertical gradients. Results of DCT and Fourier transform are completely equal when the gravity data are noise free, but in the case of noisy data, DCT has better performance than FFT. This improvement is investigated by using signal to noise ratio (SNR). The SNR of the results of DCT compared to Fourier transform is larger, therefore less noise enters in the calculation of vertical gradients by using DCT. We have tested these two transforms on the synthetic data containing Gaussian noise. First and second vertical gradients are calculated by DCT and Fourier transform. The results have shown that DCT in comparison with Fourier transform is less sensitive to noise. Moreover, these two transforms are used for calculating first and second vertical gradients of gravity data obtained from Safo manganese mine. The results have shown that less noise enters in vertical gradient map obtained using DCT. Edge detection of anomalies is one of the usage of gradients in the interpretation of gravity data. Analytic Signal has been used for edge detection of anomalies in the cases of real and synthetic gravity data. Vertical gradient of analytic signal calculated by DCT, compared to Fourier transform, has less noise and better quality. Introduction: The first and second vertical gradients are used to distinguish the difference between two adjacent anomalous bodies, reduce the effects of interference of the amplitudes of anomalies, separate the local field superimposed on the background determine, and to determine the location and dimensions of the anomalies. The gradient data are used in direct interpretation and inversion and inputs of many interpretations. Vertical gradients are more sensitive to noise than bouguer map, and second vertical gradients is more sensitive to noise than first vertical gradients. Methodology and Approaches: Generally, the relation between cosine and sine transforms with Fourier transform is equal to ‘F=C + iS’. Where F, denotes Fourier transform and C and S denote cosine and sine transforms, respectively. If the function is positive or constant, imaginary part of the Fourier transform is equal zero. In other words, Fourier transform and DCT will be equal. It can be assumed that noise is a function that is added to gravity signal. Because of the nature of noise, Fourier transform and DCT of noise will not be equal. Then, the signal to noise ratios for DCT and Fourier transform of noisy gravity signals are not the same. Hence, less noise enters the gradient map when calculating vertical gradients using DCT. Results and Conclusions: In this research, we have used DCT to calculate vertical gradients. We analytically have proven that DCT in noise-free data has equal results with Fourier transform results. However, the DCT of noisy data have shown less noise enters the data compared to Fourier transform. We have examined this issue by using synthetic data containing Gaussian noise. A comparison between vertical gradients obtained by DCT and Fourier transform indicates that DCT is less sensitive to noise. We have obtained similar results by using real data. The results of analytic signal calculated by DCT and Fourier transform in both cases of real and synthetic data have been compared. Edge detection calculated by using DCT has shown less noise due to less noise in the DCT and also has better quality compared to using Fourier transform.

The aim of this research is to develop methods to derive sharp images of subsurface, which show the geometry of the structures in the region of interest. Structural inversion of gravity data -deriving robust images of the subsurface by delineating litho-type boundaries using density anomalies- is an important goal in a range of exploration settings (e.g., ore bodies, salt flanks). In this paper, an L1-norm based inversion approach is investigated. The density distribution of the subsurface is modeled on a uniform grid of cells. The model parameters are the density of each cell that is inverted by minimizing the L1-norm objective function using linear programming (LP) while satisfying a priori density constraints. LP method for structural inversion has two advantages: 1. It offers a natural way to incorporate a priori information regarding the model parameters. 2. It gives a subsurface image with sharp boundaries (structure). The inversion quality depends on a good priori estimation of the minimum depth of the anomalous body.

Depth estimation of geological structures is one of the most important objectives in geophysical studies. Euler deconvolution (standard Euler) is a well-known method in the depth estimation. Based on standard Euler, various methods are introduced to reduce error of the depth estimation. In this study we have used a new method called Euler RDAS. This method is based on the standard Euler. In this method derivatives of analytic signal and first vertical are used in Euler equation. Applying derivatives of analytic signal for depth estimation is better than the analytic signal. There is no problem of choosing structural index in this method which increases accuracy in the depth estimation. To examine the performance of this method, depths of several synthetic models are estimated and their results are compared to that of the standard Euler. In all models, results of the RDAS Euler shows fewer errors in high depth model in comparison to that of the standard Euler. This method was tested on synthetic data with noise. RDAS Euler sensitive to noise due to usage of high-degree derivatives is less than the standard Euler. Study shows that if the noise in the data is reduced by methods such as filter upward, it can provide an appropriate estimate of the depth. In this study, the depth of gravity anomalies caused by the masses of hematite located in Kerman, has estimated using standard Euler and RDAS Euler. Upward continued by 3 m has been used for reducing noise in this data. There are 3 possible hematite masses in residual map of the study area. The minimum of high depth anomalies indicates masses of hematite, which is calculated using RDAS Euler and was about 5 meters and a maximum of high depth is about 20 meters. The minimum depth obtained using standard Euler for this anomaly is about 5 meters and the maximum depth is about 40 meters. Responses of RDAS Euler is more compatible that of standard Euler with the boundary anomalies and has smaller vertical interval which can be a criterion for more precise solutions of the RDAS Euler. For further examination, the gravity data of hematite mine is used with inverse modeling of Camacho method. Minimum and maximum upper depths obtained for these anomalies are 5 to 35 meters, respectively. In addition, the modeling results is compared with the results of depth estimation of Euler. To this, 10 points in the anomalies area are pointed and the calculated depth of these points using standard Euler, RDAS Euler and modeling are shown. Root mean square error (RMS) between Euler’s and modeling results is calculated. In comparison of the results in different methods, which are not standard, results of two methods that have the lowest RMS error is considered as the selection criteria. RMS for standard Euler with MATLAB code, Geosoft, and RDAS Euler are equal to 11.43, 8.5, and 2.66 respectively. The results of two methods among these three methods which are used to estimate the depth of hematite masses are close hence they can be more reliable results. Therefore, due to fewer errors of RMS of RDAS Euler and modelling results are more accurate than that of the standard Euler.

Enhancement on the edges of the causative source is a tool in the interpretation of potential field data. There are many methods for recognizing the edges, most of which involve high pass filters based on derivatives of potential field data. In this paper, a new edge detection method is introduce, called the enhanced mathematical morphology (EMM) filter for interpretation of field data. The EMM filter uses the ratio of the erosion of the total horizontal derivative to the dilation of the total horizontal derivative to recognize the edges of the sources, and can display the edges of the bodies simultaneously. Edge detection of the potential field data has been widely used as a significant tool for geophysical exploration technologies, which can delineate the horizontal locations of causative sources. Normally, various high-pass filters are used to recognize the edges of the potential field data (Evjen 1936; Fedi and Florio 2001; Verduzco et al. 2004; Cooper and Cowan 2006; Cooper and Cowan 2011; Ma and Li 2012; Ma 2013).
The EMM filter uses the ratio of the erosion of THD to the dilation of THD to recognize the edges of the source. Mathematical morphology was developed by Matheron and Serra in 1964 (SERRA, 1983); which is an image analysis and recognition tool. The structuring element (SE) is a basic operator in mathematical morphology, used to interact with an image and to draw conclusions about how a shape fits or misses the shapes in the image. SE consists of a matrix of 0s and 1s that can have any arbitrary shape and size. The basic operations of mathematical morphology are dilation and erosion. Dilation is defined as the maximum value in the window ascertained by the SE. Erosion is defined as the minimum value in the window ascertained by the SE. The EMM filter is expressed as (Lili et al., 2013):
where imerode (F,SE) and imdilate (F,SE) represent the erosion and dilation of the THD, respectively.
In this paper, a new relationship is presented for EMM filter that is tested on synthetic data with and without noise as well as the real potential field data in Qom salt dome. The EMM method successfully delineates the edges of the causative sources, which gives better resolution of the deeper source than other filters, and can display the edges of the bodies in a more centralized way. In this article, a new relationship is defined for the EMM filter as:
The EMM filter was used to recognize the edges of the sources. It can display the edges of the shallow and deep bodies simultaneously. The EMM filter does not require the computation of vertical derivatives, which makes this method computationally stable. The EMM filter is tested on synthetic, and real potential field data in Qom salt dome and the edge detection was done with reasonable results.

Summary To determine the location of subsurface objects from the potential field maps, edge detection methods are used. This methods are mainly derivatives of potential field. For edge detection, several methods or filters such as analytic signal, tilt angle, total horizontal derivative of tilt angle and normalized total horizontal derivative are commenly used in geophysical interpretation. In this regard, the method of normalized anisotropy variance (NAV) for edge detection is investigated in this paper. This filter directly do not use high-order derivatives, and also, vertical derivative in calculations, thus, stability of noisy data is preserved in this method. One of the effective parameters in the calculation of this method is the optimal window number that is obtained by maximum of variation of normalized anisotropy variance and window number. In order to assess performance of the NAV-method in the edge detection of gravity anomalies, cubic models used with positive and negative density contrast and different depth, and the NAV-method and also the above-mentioned filters were applied to synthetic models and real gravity data of Safo Manganese mine. Consequently, acceptable results from applying the NAV method for edge detection of gravity anomalies were obtained.
Introduction A variety of techniques based on horizontal and vertical derivatives of potential field data have been developed as effective tools for the edge detection of potential field anomalies. Location of maximum horizontal derivative can be used an indicator of the location of the anomaly’s edges. The first vertical derivative (FVD) is positive over the source, zero-crossing value of FVD over the edge, and outside of a vertical sided source is negative. The peaks of the horizontal derivatives indicate the edges and zero over the body. One of the most important problems in using derivative-based edge detection techniques is that the noise in the potential field data increases during the process. Analytic signal is one of the filters that employs both horizontal and vertical derivatives to determine the boundary of the anomalies. Local phase filters are set of filters based on horizontal and vertical derivatives of potential field data that are used in geophysical interpretation. Tilt angle, total horizontal derivative of tilt angle and normalized total horizontal derivative (TDX) are conventional local phase filters that are sensitive to noise. To achieve the best result for edge detection of gravity anomalies, we introduce normalized anisotropy variance or NAV method that reduces sensitivity to noise.
Methodology and Approaches To determine the edge of subsurface bodies, NAV edge detection filter has been applied to gravity data. The NAV method has been investigated as a tool to detect the edge of anomalies.In this regard, its ability to determine the edge of synthetic models from gravity data, in a noisy and free-noise state has been demonstrated in this paper. Anisotropy scale and window number are two important parameters on the value of the NAV. The anisotropy scale in anomalies separation in the vicinity of each other, is an effective parameter. An empirical approach to select these parameters has been applied. Window number (M) is associated with data quality, for noise-free data M=3 is acceptable. For data that are contaminated with noise, the value of M is determined by variation of the maximum value of NAV with window number. Since the NAV does not use any of higher-order derivatives directly, it is less sensitive to noise than the other methods. In MATLAB platform, a computer program has been prepared and its outputs are extracted.
Results and Conclusions In this research, capability of NAV method in estimating boundaries of synthetic and real potential field anomalies has been investigated. Vertical derivative that creates complexity for the interpretation, has not been used in this method; hence this method or filter for deep resources is more stable. Similar to tilt angle filter, the zero value of NAV corresponds to the edges of anomalies. The value of the NAV filter is positive over the body and negative out of it; therefore, area of the body can be detected from other parts. The results also show that this method is less sensitive to noise and the boundaries that are determined for the model, are in good accordance with reality. For the model that a small prism superposed on a big prism, the NAV could detect the edges of the superposed prism with high accuracy.

The gravity method is one of the geophysical tools used for geological, engineering and environmental investigations where the detection of geological boundaries, cavities, subsurface karstic features, subsoil irregularities, or landfills are essential. In higher accuracy measurements, the microgravity method has been widely and successfully used for locating and monitoring subsurface materials.
Since microgravity methods measure gravity variations at the surface, they are directly influenced by the density distribution in the subsurface and particularly by the presence of formation material, which may create a mass deficit relative to the density of the surrounding terrain. In many cases, deep or small-scale heterogeneities generating low-amplitude anomalies can be detected and the reliability of further interpretation requires highly accurate measurements which are carefully corrected for any quantifiable disturbing effects. The main purpose of the research, that was conducted in small part of a dam site, is to determine the quality and type of subsurface structures in location of tunnel construction. Study area for collecting microgravity data was located at a small part, considered for construction of Siah Bisheh dam, road Tehran to Chalous. Position of microgravity stations were over a tunnel path which in some parts encountered with collapsing structures. The study area was part of Alborz Mountains. Geology formation(Shemshak formation), consisting of lime beds together with igneous rocks which are severely affected by fractures. Data were collected along 13 profiles with separating distance of 15 m. The stations distance and number of data were 15 m and 148 respectively. Bouguer gravity anomaly was calculated after making corrections such as earth tide, free air, Bouguer, topography and terrain effects. The regional effect obtained using a program that is written in FORTRAN to fit orthogonal and orthonormal polynomials on the observed data and then residuals were estimated. Three negative anomalies were distinguishable in residual gravity map. Data of these anomalies are modeled with a 3-D inversion approach using GROWTH 2.0 software. The GROWTH 2.0 is an inversion tool which enables the user to obtain, in a nearly automatic and non-subjective mode, a 3D model of the subsurface density anomalies based on the observed gravity anomaly data. The current version of the tool has been developed from an earlier code (Camacho et al., 2002). In a nearly automatic approach, the software provides a 3-D model informing on the location and shape of the main structural building blocks of the subsurface structures. Then densities contrast of these anomalies was estimated. Result of the inversion was a 3-D distribution of densities contrast. To show this distribution of the densities contrast, the horizontal and vertical sections at different depth and different horizontal positions were selected and interpreted. From these sections it is indicated that the effective depth of the data, for identifying martial of subsurface structures from the inversion, is about 50 m. In the sections, areas with low densitig contrasts are related to the fractured limestone and those with high contrast ones are related to the compact limestone or igneous rocks. Existence of igneous and lime rocks that have more density and compactness, increase the quality of the structures in the path of the tunnel construction. Areas including fractured limestone, with lower density, decrease the quality of the structure and increase the risk of water permeability and collapsing in the path of the tunnel construction. Thus by interpreting of the results of the microgravity data inversion, areas with high and low compactness and good and bad quality rocks for tunnel construction are recognized, those are related to the fractured or karstic limestone and limestone and igneous rocks. Also boundaries of these formations where densitig contrasts vary suddenly, are related to the existence of faults.

Summary: In this paper, 3D inversion of gravity data for determination of subsurface density distribution is made using geostatistical co-kriging method. Co-kriging is a mathematical interpolation and extrapolation tool. It uses the spatial correlation between the secondary variables and a primary variable to improve the estimation of the primary variable at un-sampled locations. The Co-kriging method gives weights to the data so as to minimize the estimation variance (the co-kriging variance). In this paper, the primary variable is density, (estimated by ρ*) and the secondary variable is gravity g. For determination of kernel matrix, the subsurface area is divided into large number of rectangular blocks of known sizes and positions. The unknown density contrast of each prism is the parameter that should be estimated. In addition, the weighting matrix has also been used in order to improve the depth resolution. Preconditioned conjugate gradient method has been used for inversion. The computer program has been written in MATLAB and tested on synthetic gravity of a rectangular prism model. The results indicate that the geometry and density of the reconstructed model are close to those of the original model. The gravity data acquired in an area, which includes concealed manganese ore bodies (Safoo mine site), in northwest of Iran. The results show a density distribution in the subsurface from the depth of about 5 to 35-40 m. These results are in good agreement with the results of the borehole drilled in the site.
We may encounter two problems in gravity data inversion: non-uniqueness and non stability of solutions. The first one occurs for two reasons: The first reason is known as the theoretical ambiguity of the unknown nature of potential theory. The second reason is known as algebraic uncertainty, which is considered when the number of parameters is greater than the number of observations. The second problem may occur because of bad condition (ill-condition) of the kernel matrix and the presence of noise in the data. For finding a unique and stable solution, constraints should be considered in the objective function, and then, the new objective function, which is replaced the initial objective function, should be minimized.
Methodology and Approaches: For determination of kernel matrix, the subsurface of the survey area is divided into a large number of rectangular prisms of known sizes and positions. The unknown density contrast of each prism is the parameter to be estimated. This kind of parameterization is flexible for the reconstruction of the subsurface model, but generates more unknown model parameters than observations (here N Results and By applying co-kriging stochastic algorithm on synthetic data in states of without and with random noise, good results for the density and depth of the model have been achieved. By applying the method on actual data from Safoo manganese mine site, the results obtained for the depth and density of the subsurface ore body are in good agreement with the results of drilling implemented in the mine site. Furthermore, as a result of applying the method on the data, general shape of the subsurface ore body was well determined.

We apply two forward methodologies in order to study density and susceptibility structure of the crust and upper mantle. The study area is a profile crossing the Zagros collision zone located as margin of Eurasia-Arabia converging plates. Gravity modeling focusing on lithospheric structure is performed in thermodynamic framework in which chemical composition is important and provides an understanding of deep layers in lithosphere like Moho and Lithosphere-Asthenosphere Boundary. Results on the crustal thickness show minimum values beneath the Arabia Platform and Central Iran (42–43 km), and maximum values beneath the Sanandaj Sirjan zone (SSZ; 55–63 km). Results on the lithosphere thickness a long profile also indicate that the Arabian lithosphere is approximately 220 km thick, toward North West of Iran especially below the Central Iran rises up to 90 km. In the profile (central Zagros), lithosphere thinning occurs in wider region, from the Zagros fold thrust belt to the Sanandaj Sirjan zone. Our results are based on application of average Proterozoic mantle compositions in modeling beneath the Arabian Platform, Mesopotamian Foreland Basin and Iranian Plateau. After rough estimation of upper crust via integrated modeling by elevation, gravity and geoid data, the distribution of density and magnetic susceptibility values allows us to perform a study in crustal scale. Afterwards, determination of the homogenous blocks with the same density and susceptibility, the geometry to different crustal layers including sediments, upper, middle and lower crust deep to Moho boundary were refined in crust-scale study based on regional model in lithospheric scale. Presence of Main Zagros Fault is a bold point in our modeling which leads to better fit of gravity data.

In this paper the 3D inversion of gravity data using two different regularization methods, namely Tikhonov regularization and truncated singular value decomposition (TSVD), is considered. The earth under the survey area is modeled using a large number of rectangular prisms, in which the size of the prisms are kept fixed during the inversion and the values of densities of the prisms are the model parameters to be determined. A depth weighting matrix is used to counteract the natural decay of the kernel, so the inversion obtains reliable information about the source distribution with respect to depth. To generate a sharp and focused model, the minimum support (MS) constraint is used, which minimizes the total area with non zero departure of the model parameters from a given a priori model. Then, the application of iteratively reweighted least square algorithm is required to deal with non-linearity introduced by MS constraint. At each iteration of the inversion, a priori variable weighting matrix is updated using model parameters obtained at the previous iteration. We use the singular value decomposition (SVD) for computing Tikhonov solution, which also helps us to compare the results with the solution obtained by TSVD. Thus, the algorithms presented here are suitable for small to moderate size problems, where it is feasible to compute the SVD. In Tikhonov regularization method, the optimal regularization parameter at each iteration is obtained by application of the parameter-choice method. The method is based on the statistical distribution of the minimum of the Tikhonov function. For weighting of the data fidelity by a known Gaussian noise distribution on the measured data and, when the regularization term is considered to be weighted by unknown inverse covariance information on the model parameters, the minimum of the Tikhonov functional becomes a random variable that follows a distribution. Then, a Newton root-finding algorithm can be used to find the regularization parameter. For truncated SVD regularization, the Picard plot is used to find a suitable value of truncation index. In math literature, a plot of singular values together with SVD and solution coefficients is often referred to as Picard plot. To test the algorithms, a density model which consists of a dipping dike embedded in a uniform half-space is used. The surface gravity anomaly produced by this model is contaminated with three different noise levels, and are used as input for introduced inversion algorithms. The results indicate that the algorithms are able to recover the geometry and density distribution of the original model. In general, the reconstructed model is more sparse using TSVD method as compare with Tikhonov solution. This especially happens for high noise level, where there is an important difference between two solutions. In this case, while TSVD produces a sparse model, the solution of Tikhonov regularization is not sparse. Furthermore, the number of iterations, which is required to terminate the algorithms, is more for TSVD as compare with Tikhonov method. This feature, along with automatic determination of regularization parameter, makes the implementation of the Tikhonov regularization method faster than TSVD. The inversion methods are used on real gravity data acquired over the Gotvand dam site in the south-west of Iran. Tertiary deposits of the Gachsaran formation are the dominant geological structure in this area, and it is mainly comprised of marl, gypsum, anhydrite and halite. There are several solution cavities in the area so that relative negative anomalies are distinguishable in the residual map. A window of residual map consists of 640 gridded data, which includes three negative anomalies, that is selected for modeling. The reconstructed models are shown and compare with results obtained by bore holes.

Optimization methods such as Simulated Annealing (SA), Genetic Algorithm (GA), Ant Colony Optimization (ACO) and Partial Swarm Optimization (PSO) are very popular for solving inverse problems. Gravity inversion is one of the fields of geophysical exploration that mainly used for mining and oil exploration. In this research optimization method is considered for gravity inversion of a sedimentary basin to find out its geometry. A new method for 3D inversion of gravity data is developed. In such optimization methods before anything else having a forward model is necessary. This model is the relationship between Bouguer gravity anomalies and a combination of prisms. The gravity anomalies of the density interface are generated by equating the material below the interface to a series of juxtaposing rectangular blocks. The stochastic optimization method that is used for solving inverse problem is Simulated Annealing. With a try and error method at first the cost function of a lot of choices is measured. Then the minimum value between them is used as the primary set of variables for introducing to model that is the start configuration of model variables. After that the repeating algorithm starts until the cost function reaches to as small as possible value by considering the geophysical constraints of the place that is being studied. Finally the geometry of prisms that is the depth of each prism is achieved. This geometry shows the situation of anomaly source. For better inversion regional anomaly is used with introducing some unknown constants to the model. The values of unknown parameters of model are extracted in a continuous range that is obtained from priori information of the region. At first the algorithm was used in solving a synthetic problem that is defined by developing random data in an arbitrary region, so with using forward model again and again finally the cost function reaches to its minimum value. For evaluating the success or failure of the algorithm, the contour map of depths that was used in forward model is compared with map that is produced after inversion and discrepancies was considered. The good results of synthetic problem lead to implementing the algorithm for real gravity data from Aman Abad region in Arak city. For this purpose at first the gravity data must be changed to Bouguer gravity anomalies by some reduction. After that the repeated algorithm is implemented and finally the results are compared with some priori information which is obtained by madding boreholes around the region. This priori information claims that the minimum and maximum depth of prisms is between 70 meters to 120 meters. Results of the algorithm are compatible with this information and show the power of algorithm in solving this kind of inverse problem. It must be mentioned that without using this priori information the results of gravity inversion are not unique and many different interpretations can be concluded even interpretations that cannot be accepted in geophysical view. As a norm to find the ability of algorithm in solving this problem the cost function can be helpful which shows the amount of error. The maximum value of it was 0.5 mgal in some regions in this research.

In this paper, a new formula is applied for the calculation of gravity anomalies from a cylinder model representing a geological body. Compared to conventional methods, this new development allows the cylinder to be freely oriented in space. 38 gravity forward models are produced each of which anomaly attributes are calculated. Then a linear relationship is established between attributes and source parameters as a new formula. These attributes are relations which calculated from coordinates of special parts of residual gravity anomaly curve. Using this linear relationship, source parameters can be then estimated by gravity anomaly attributes. Linear relationship is obtained by least-square method and minimizing differences. Consequently, Compared to previous methods, this new development considers multiple factors that have impact on geophysical observations (some neglected in previous studies) and more variables are considered such as dip angle, strike direction, size, depth to the top, and density of the cylinder. These parameters are important for determination of the geometric of subsurface geological body. Based on a series of forward modeling using the new formula, a multiple linear regression system has been developed. The multiple linear regression method relates the variations of residual gravity anomaly which is changed by variations source parameters. Based on previous studies, we suppose that the shape as well the amplitude of the gravity anomaly will change with the changes of cylinder occurrence. We use the multiple linear regression to examine if there is a linear relationship between each parameter of the cylinder and a series of attributes from the gravity anomaly. Actually, in this algorithm, we assessed the variability of a residual gravity anomaly as a function of the source occurrence. In previous study, seven most significant parameters of a cylinder-like body were identified, which influence the shape of the gravity anomaly as well its intensity. We build up a linear regression system that contains six linear relationships between the six most significant parameters of a cylinder-like body and seven attributes from the gravity anomaly. Those equations allow the user to estimate the multiple parameters of an elongated geological body simultaneously under the constraint of gravity observations. Integrating those calculations into gravity surveys helps in making a drilling decision. In this article, a shaft situated in Siah bisheh dam has been considered as a case study to verify proposed formula. This case study was a part of The Siah Bishe Pumped Storage project. The project was located in the Alborz mountain range, 125 km north of Tehran. The site can be reached on the main Chalus road, connecting Tehran with the Caspian Sea. The project area lies in the southern part of the Paleozoic- Mesozoic Central Range of the alpine Alborz mountain chain. The rock sequences in the project area consist of massive limestones, detrital series (sandstones, shales) and volcanic rocks of Permian formations, Triassic dolomites and Jurassic formations with black shales and sandstones. Several tectonic faults are crossing the project alignment. In this area, main purpose of gravity surveying was exploration of collapse zones. The gravity data used in this study come from gravity department of institute of geophysics in IRAN. The spatial resolution of the original gravity observations was 15 m between stations. The Bouguer anomaly grid was then interpolated to a spacing of 3 m using a minimum curvature gridding algorithm. The residual Bouguer anomaly is obtained after removing a first order polynomial trend to the Bouguer anomaly. Moreover, its engineering and geological parameters have been calculated. Collapsed zone in abovementioned shaft has been determined and illustrated using microgravity data. We applied the new method to residual gravity anomaly curve of collapsed zone. Finally, this zone has also been mathematically modeled using inclined cylinder.

In this paper the inversion of gravity data using L1–norm stabilizer is considered. The inversion is an important step in the interpretation of data. In gravity data inversion, the goal is to estimate density and geometry of the unknown subsurface model from a set of known observation measured on the surface. Commonly, rectangular prisms are used to model the subsurface under the survey area. The unknown density contrasts within each prism are the parameters which should be estimated. The inversion of gravity data is an example of underdetermined and ill-posed problem, i.e. the solution can be non-unique and unstable. Thus, in order to find an acceptable solution regularization should be imposed. Solution is usually obtained by minimizing a global objective function consisting of two terms, data misfit and the regularization term. Data misfit measures how well an obtained model can reproduce the observed data. Usually, it is assumed noise in gravity data is Gaussian, therefore a L2–norm measure of the error between observed and predicted data is well suited for data misfit. There are several choices for a stabilizer, depends on type of features one wants to see from inverted model. A typical choice is a L2 –norm of a low-order differential operator applied to the model, which also a priori information and depth weighting can be incorporated (Li and Oldenburg, 1996). In this case the objective function is quadratic, then minimization of the function results a linear system to be solved. However, the models recovered in this way are characterized by smooth feature which are not always consistent with the real geological structures. There are situations in which the sources are localized and separated by sharp, distinct interfaces. To deal with this problem, during last decades, researchers have proposed a few types of stabilizer. Last and Kubik (1983) presented a compactness criterion for gravity inversion that seeks to minimize the area (or volume in 3D) of the causative body. Portniaguine and Zhdanov (1999) based on this stabilizer, who named the minimum support (MS), developed the minimum gradient support (MGS) stabilizer. For both constraint, the regularization term can be written as the weighted L2–type norm of the model. Therefore, the problem of the minimization of the objective function can be treated same as conventional Tikhonov functional. The only difference is that a priori variable weighting matrix for model parameters incorporated in the regularization term. Thus the Iteratively Reweighted Least Square (IRLS) algorithm is required to solve the problem. Other possibility for stabilizer is the minimization of the L1-norm of model or gradient of model, the latter indicates total variation regularization. The L1–norm stabilizer allows occurrence of large elements in the inverted model among mostly small values. Therefore, it can be used to obtain sharp boundaries and blocky features. Although the L1–norm stabilizer has favorable properties, in reconstruction of sparse models, its numerical implementation in a minimization problem can be difficult because its derivatives with respect to an element is not defined at zero. To overcome this difficulty, in this paper, the L1–norm stabilizer is approximated by a reweighted L2 –norm term. The algorithm is extended to gravity inverse problem, which needs depth weighting and other priori information to be included in the objective function. For estimating the regularization parameter, which balances between two terms of objective function, the Unbiased Predictive Risk Estimator (UPRE) method is used. The solution of the resulting objective functional is found using Generalized Singular Value Decomposition (GSVD), also provides for efficient determination of the regularization parameter at each iteration. Simulation using synthetic data of a dipping dike demonstrates that the method is capable to reconstruct focused image, boundaries and slop of the reconstructed model are close to those of the original model. The method is applied on gravity data acquired over the Gotvand dam site, in the south-west of Iran. The results show rather good agreement with those obtained from the boreholes.

The main objective of interpretation of acquired gravity data on the Earth''s surface is to determine the contrasts in density or shape/dimension of mass anomalies. Interpretation of gravity data can be done through an inversion process. In this research، a block model has been considered for the subsurface anomalous mass. By considering a constant initial density (about 2. 6 gr/cm3) for all blocks and by using inversion method، distribution of density of the anomalous mass was estimated and interpreted. In this research، Occam method is used to invert 246 gravity data collected in 2007. Results of the gravity data inversion show sufficient fit between observed and calculated gravity data. Using this inversion method، distribution of density in the subsurface layers related to sediments and basement are estimated in this area. Since there is a density contrast between sedimentary layers and basement، the estimated density distribution can help to explore the lithology of formations as well as the discontinuities in them. Densities less than 2 gr/cm3 in horizontal and vertical sections obtained from the inversion are attributed to the alluviums. The depth of these sediments، which include sand، silt and clay of different percentages، is estimated to be less than about 200 m. Unequal density distribution along the layers is taken to indicate fractures. In fact، these fractures are associated with part of the Tabarteh fault in this area، which caused numerous earthquakes (but less than 5 Richters in magnitude) around the Arak and Dawood Abad cities in past years.

In this paper the 3D inversion of gravity data is considered. The goal is to reconstruct models of subsurface density distribution using a set of known gravity observations measured on the earth surface. The subsurface under the survey area is divided into large number of rectangular blocks of known sizes and positions. The unknown density contrasts within each prism define the parameters to be estimated. This kind of parameterization is flexible for the reconstruction of the subsurface model, but requires more unknown model parameters than observations (here N << M, where N is the number of data and M is the number of model parameters). The final density distribution will be obtained by minimizing a global objective function consists of data misfit and a regularization term. The inverse problem is solved in data space, which needs inverse of matrix with N×N dimension, as compared with M×M dimension system in model space inversion. This methodology was used by Pilkington (2009) in 3D inversion of magnetic data. To solve the resulting set of linear equation, the conjugate gradient method is used. Combination of data-space method with conjugate gradient leads to keep the storage and computational time to a minimum. The iteratively-defined regularization matrix, which is used in objective function, is a combination of three diagonal matrix; namely depth weighting, compactness and hard constraint matrices. The compactness constraint was introduced in Last and Kubik (1983) and developed in Portniaguine and Zhdanov (1999), who used term "minimum support stabilizer", is considered here to produce models with non-smooth features. It is a suitable and well-known constraint for identifying geologic structures which have material properties that vary over relatively short distances. The depth weighting matrix, introduced in Li and Oldenburg (1998), is used in regularization term to counteract the natural decay of the kernel with depth. The hard constraint allows us to incorporated priori geological and geophysical information into inversion process. While depth weighting and hard constraint matrices both are independent of the iteration index, the compactness depends on iterations. In order to recover a feasible image of the subsurface, realistic lower and upper density bounds are imposed during the inversion process. The computer program is written in MATLAB and tested on synthetic data produced by a model consists of two cubes. The cubes have same dimension and density, but located at different depths. The results indicate that the algorithm is efficient to handle large-scale gravity inverse problems. For the shallow cube the geometry and density of the reconstructed model are close to those of the original model, but for the deeper body the resolution decrease and a smooth image of subsurface obtained. The gravity data acquired over the Safo mining camp in the north-west of Iran, which is well-known for manganese ores, are used as a real modeling case. The results show a density distribution in the subsurface from about 5 to 35-40 m in depth and about 35 m extent in the x direction, which are close to those obtained by bore-hole drilling on the site.

One of the most important exploration problems in geophysics is to estimate the geophysical parameters from the observed or residual gravity anomaly related to a buried structure, such as depth, amplitude coefficient and geometrical shape factor. The gravity anomaly expression produced by a simple geometrically shaped model (sphere or cylinder) can be represented by an appropriate analytical formula. Several interpretative methods have been developed to interpret gravity field data assuming a fixed simple geometricalmodel such as a sphere, a horizontal cylinder or a vertical cylinder. In most cases, these methods consider the geometrical shape factor of the buried body to be a priori assumed, and the depth variable may thereafter be obtained by graphical methods applied to the residual anomaly. However, only a few methods have been developed to determine the shape of the buried structure from the residual gravity anomaly. Consequently, the accuracy of the results obtained by these methods depends on the accuracy within which the residual anomaly can be separated from the observed gravity anomaly. In this study, a new and simple method has been developed to estimate the depth, amplitude coefficient and geometrical shape factor of a buried structure from the observed (composite) or residual gravity anomaly related to a cylinder or sphere-like structure. The method is based on nonlinearly constrained mathematical modeling and also stochastic optimization approaches. This method consists of three main steps: The first step is oriented to formulate a nonlinearly constrained optimization model (NCOM) which mathematically describes the geophysical gravity problem related to the studied structure. The (NCOM) model is to optimize a mathematical objective function on an unbounded subset (defined by mathematical inequalities constrains in which the geophysical parameters are generally surmised to satisfy) contained in the free geophysical parameters. Ignoring these mathematical constrains probably leads to general error estimations of the parameters. In this research, the objective function was taken as the statistical likelihood function which depends on the deviations between the observed and synthetic points and also on the number of observations. The second step is directed to suggest an interior penalty function to transform the (NCOM) model into a nonlinearly unconstrained optimization one (NUOM). The goal of using the penalty function is to eliminate the constraints of the (NCOM) model and make them reactive in a new target function of the (NUOM) model. The target function of the (NUOM) model considers both the objective function of the (NCOM) model and the suggested interior penalty function. The third step is to solve the (NUOM) model by the adaptive simulated annealing algorithm, a stochastic approach, well-known for optimizing numerical functions of several real decision variables. The obtained solution of the (NUOM) model includes the geophysical gravity parameters of the studied structure such as: depth, amplitude coefficient and shape factor. A statistical analysis has been carried out to demonstrate the accuracy and the precision of the suggested interpretative method. We applied this method to some theoretical synthetic examples in order to evaluate the precision of the suggested method. We also used the method to estimate the mentioned parameters for the gravity anomaly of the Abadeh site. The obtained results had an appropriate agreement with other methods.