Geophysical Journal International Geophys. J. Int. (2017) 208, 546–560 Advance Access publication 2016 November 2 GJI Gravity, geodesy and tides doi: 10.1093/gji/ggw413 Improving compact gravity inversion using new weighting functions Mohammad Hossein Ghalehnoee, Abdolhamid Ansari and Ahmad Ghorbani Department of Mining and Metallurgical Engineering, Yazd University, P.O. Box 89195–741, Yazd, Iran. E-mail: mhghalehnoee@gmail.com Accepted 2016 November 2. Received 2016 October 31; in original form 2016 January 16 Key words: Numerical solutions; Inverse theory; Gravity anomaly and Earth structure. 1 I N T RO D U C T I O N The geophysical inverse problem is an attempt to obtain information about the Earth interior from measurements of physical entities in the near surface region. When the inverse problem demands too much information from the geophysical data, its solution becomes either non-unique or unstable. Problems of this type are called illposed problems (Barbosa & Silva 1994). Many authors introduced solutions to overcome this problem in the last three decades. Green (1975) found the closest model to the initial model, Li & Oldenburg (1996, 1998) counteracted the decreasing sensitivity of the cells with depth by weighting it with an inverse function of depth. Smoothness or roughness of physical parameters distribution that control gradients in spatial direction used by Pilkington (1997) and Li & Oldenburg (1996, 1998). Boulanger & Chouteau (2001) employed a 3-D inversion algorithm that implements several constraints, including minimum distance, flatness, smoothness and compactness. Last & Kubik (1983), in particular, developed a method called ‘compact inversion’ explaining the observed anomaly by structures of minimum volume. Their strategy is to minimize the area (or volume in three dimensions) of the model. This is equivalent to maximizing its compactness. This method strongly depends on the number of model parameters (prisms), and by increasing the model 546  C parameters, the anomalies tend to concentrate near the surface. Guillen & Menichetti (1984) applied an approach which includes the search for solutions minimizing the moment of inertia with respect either to the centre of gravity or to an axis of a given dip line passing through it. The method works properly for a single gravity source, but the problem is dealing with multi-source and complicated anomalies which do not lie in one point or one axe. Barbosa & Silva (1994) generalized the compact inversion method to compact along several axes using Tikhonov’s regularization. They improved the method offered by Guillen & Menichetti (1984) for multi-source and complicated anomalies. The newest gravity compact method is the so called interactive inversion that was developed by Silva & Barbosa (2006), which estimates the location and geometry of several density anomalies. They simplified their old method (Barbosa & Silva 1994) for computational performing. The method is suitable for multi-source and even complicated anomalies depending on the quantity and quality of ‘a priori’ information. Silva et al. (2009, 2011) also developed interactive inversion successfully for 3-D gravity data closest to pre-specified geometric elements such as axes and points. One of the advantages of compact inversion over other methods is its simplicity which renders a lower computational calculation. Its second advantage is an easier constraining of upper and lower bounds. But this method has its own limitations since it highly The Authors 2016. Published by Oxford University Press on behalf of The Royal Astronomical Society. Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 SUMMARY We have developed a method to estimate the geometry, location and densities of anomalies coming from 2-D gravity data based on compact gravity inversion technique. Compact gravity inversion is simple, fast and user friendly but severely depends on the number of model parameters, that is, by increasing the model parameters, the anomalies tend to concentrate near the surface. To overcome this ambiguity new weighting functions based on density contrast, depth, and compactness models have been introduced. Variable compactness factors have been defined here to get either a sharp or a smooth model based on the depth of the source or existence of prior information. Depth weighting derived from one station of gravity data whereas the effect of gravity data is 2-D and 3-D. To compensate this limitation an innovating weighting function namely kernel function has been introduced which multiplies with weight and compactness matrixes to yield a general model weighting function. The method is tested using three different sets of synthetic examples: a body at various depths (20, 40, 80 and 140 m), two bodies at the same depth but various distances to estimate lateral resolution and three bodies with negative and positive density contrast in different depths. The method is also applied to three real gravity data of Woodlawn massive sulphide body, sulphides mineralization of British Colombia and iron ore body of Missouri. The method produces solutions consistent with the known geologic attributes of the gravity sources, illustrating its potential practicality. Improving compact gravity inversion 547 depends on the number of data and parameters. If the model parameters choose a very large number compared with the number of the data, the anomaly may unhopefully concentrate at the surface or shallower than the real depth. To solve this problem, one may decide to opt a less number of data parameter, then choosing a large dimension of prisms. In this situation the recognition of the true geometry of the ore body may be hard and, in some cases the correct anomaly position with large dimensions of prisms is also not possible. The generalized compact and interactive inversion strongly need ‘a priori’ information to yield an accurate estimation. A priori information refers to geological information, well logs and previous inversions. When there is little information about gravity sources, modelling should be done with other inversion methods. The method introduced by Li & Oldenburg (1998), which is largely used, can be applied with or without geological constraints in two and three dimensional measurements (Williams 2008). This method has many coefficients and functions that should be estimated referring to the kernel of potential field data and existence or absence of prior information. This procedure is time consuming and makes the method rather unfriendly. We developed a new approach for interpreting 2-D gravity anomalies produced by multiple and complex gravity sources at any depths. New weighting functions have been used based on the depth of the prisms (Fig. 1), kernel matrix and compactness weighting. The inversion can be made acceptable with and without prior information for either sharp or smooth model. 2 METHODOLOGY M  Ai j m j + ei i = 1, 2, . . . , N ; j=1 j = 1, 2, . . . , M (1) where m j is the density of the jth prism, ei is the noise associated with ith data point, and Ai j is kernel matrix, the elements of which represent the influence of the jth prism on the ith gravity value (Last & Kubik 1983). The data equation can be written in matrix notation g = Am + e The inversion method here is linear, like other linear geophysical inversion: given, the observed gravity data (g), finding a density distribution ‘m’ which explains ‘g’, with a certain noise level. The solution of the system in kth iteration, can be in the leastsquare problem in matrix notation (Menke 1989):  −1 −1 −1 −1 A T AWm(k) A T + µWe(k) δg k (3) δm k+1 = Wm(k) m k+1 = m k + δm k+1 , (2) (4) where Wm is model weighting matrix, We is a noise weighting matrix in kth iteration, both of which are diagonal, and δg (k) = (g obs − Am (k) ). Mu (µ) is damping factor to get rid of matrix singularity, having the positive small value and depending on the noise level of the data points. The less value of damping factor refers to the less noise level of the data points. We introduce a general weighting function including the compactness weighting, depth weighting and kernel weighting matrices in kth iteration: Wm(k) −1 = Wz Wc(k) Wa . The method used here is based on 2-D model that is illustrated in Fig. 1. For this model the gravity effect at the ith data point is given by gi = Figure 2. The kernel matrix (A) for: (a) one station (A5j ), and (b) summation N the whole station ( i=1 Ai j ). The filled triangles indicate the data points. (5) Noise weighting matrix We , can be simply written as Last & Kubik (1983):   −1 −1 (6) We(k) = diag AWm(k) AT . 2.1 Depth weighting Gravity data, like other potential field data, have no inherent resolution. When minimizing m2 = ∫ m 2 dv, final model tend to concentrate near the surface regardless of the true depth of the causative bodies. This arises because the constructed model is a linear combination of kernels that decay rapidly with depth (Li & Oldenburg 1998). Many attempts have been done by authors as Mendonca & Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 1. The 2-D model, showing data point i and prism j. 548 M.H. Ghalehnoee, A. Ansari and A. Ghorbani Silva (1994, 1995) used N × N diagonal normalization matrix D which is given by ⎛ ⎞−1/2 M  Ai2j ⎠ . Dii = ⎝ (7) the kernels. A depth weighting has been presented here appropriate for 2-D gravity data: wz j = −1 Dg. (8) Eq. (8) concentrate the anomalies near the earth’s surface due to the ambiguity or non-uniqueness problem. This mathematical solution provides little information about the true structure. Stocco et al. (2009) presented a depth weighting function that, multiplied by the values of the kernel under a measured point, gives a constant value equal to 1: −1/2 , z j + zi Wz = diag (wz ) , And eq. (3) becomes w j = A1, j β (9) where A is the kernel and j the index of the prisms under the first measured point. This function is absolutely objective and only depends on the kernel. They applied this function for compact inversion of magnetic data. By using this function, the anomalies at shallow depth can be clearly demonstrated. But with increasing the depth or the geometric complexity of the anomalies, this approach yields unreliable results. Both of the above approaches are unsuccessful; therefore, the problem of concentrating density at the surface can be solved by introducing a better depth weighting to counteract the natural decay of (10) i=1 j=1 m = A T D DG A T D + µI N  (11) where z j is the depth of the jth cell, z i is the measurement height for each data from surface, and β is a constant value. The value of β is usually chosen to reproduce the exponential decay of the gravity or magnetic response of a sphere with distance (Williams 2008). Since gravitational effect decays with the inverse distance square, Li & Oldenburg (1998) used β = 2 for gravity data. Accordingly, the same value of β has been chosen in this research. 2.2 Compactness weighting The compactness functional (minimum area or minimum support) was first introduced by Last & Kubik (1983), where the authors suggest seeking the source distribution with the minimum area (or minimum volume in 3-D case) to model the anomaly. This concept is illustrated as: if d and h are the prism dimensions (Fig. 1), a definition of area for 2-D model is (Last & Kubik 1983)  α M m j    α . (12) area = dh m j  + ε j =1 This leads to choose the compactness weighting function  α = m j  + ε −1 wcj (13) Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 3. First synthetic example, a causative 2-D body lay at the top depths of 20, 40, 80 and 140 m with the density contrast of 1 gr cm−3 . Improving compact gravity inversion 549 Wc(k) = diag wc−1 , (14) where α is compactness factor and the parameter ε is a small number of 10−7 which is introduced to provide stability as m j → 0. This weight is further developed by Portniaguine & Zhdanov (1999) who used the term ‘minimum structure’. They added a reference model (mo ) to eq. (13); therefore, the new compactness function is yielded  α −1 wcj = m j − m oj  + ε. (15) The reference model mo may be a general background model that is estimated from previous investigations, or it could be the zero model. The value of α has been chosen as a constant value by many authors. Last & Kubik (1983) and Vatankhah et al. (2014) presented α = 2; Guillen & Menichetti (1984), Barbosa & Silva (1994), Silva & Barbosa (2006), Silva et al. (2009, 2011) and Grandis & Dahrin (2014) have chosen α = 1. In this study, the value of α is variable and depends on the situations coming from the depth of causative bodies and prior information. When the source lies at large depths, the large value of α can be opted and vice-versa. If there is a subsurface geological information such as drilling or known minimum and maximum bounds, the value of compactness factor can also be large to have a sharp model. Generally, if there are not any subsurface information or little prior information, at the first investigation, the inversion begins with a few value of compactness as α ≤ 0.5. After drilling operations or any geological and geophysical studies based on first stage, the second stage (i.e. constrained inversion) is done, using a more large value of compactness factor. Finally, it should be pointed out that the appropriate range of value is 0 ≤ α ≤ 2 for 2-D gravity inversion. 2.3 Kernel weighting Since depth weighting function is completely one dimensional (e.g. for one point measurement) based on Li & Oldenburg (1996), a more appropriate weighting function should be considered. Fig. 2(a) shows kernel behaviour for fifth station g5 with respect to depth Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 4. The model at top depths 20 and 40 m (shallow); (a,b) inversion of first synthetic data with compactness factor of 0.5, (c,d) inversion with compactness factor of 1, (e,f) inversion with compactness factor of 2; the rectangle indicates the location of true anomaly. 550 M.H. Ghalehnoee, A. Ansari and A. Ghorbani (A5j ). As shown the value of kernel decays rapidly and spherically with depth. To estimate the whole data decay, we decided to add up  Nkernels in each station. Fig. 2(b) shows the result of summation i=1 Ai j , which depicts the whole decay in 2-D gravity measurement (in this example the data points are 12 with a distance of 10 m). In this situation the decay is deferent from the case for one station; hence, a new weighting function shall be introduced. Both potential data measurements and the effect of decay of the potential field response are always two and three dimensional; therefore, using depth weighting is not sufficient. This problem arises especially when some gravity sources lie at deferent depths near each other. In this situation the deep anomaly may model deeper than the real location. Consequently, a new weighting function, named kernel weighting function, is introduced here wa j = N  i=1 Ai j (16) Wa = diag (wa ) , (17) The exact expression of kernel matrix (A) for 2-D gravity data has been simplified by Last & Kubik (1983). 2.4 Inversion procedure For the beginning, if there is no priori information, the reference model is chosen m o = 0, then Wm (eq. 5) and We (eq. 6) all of which are calculated, based on value of compactness factor. After that, we will have iteration procedure with least-squares solution as eqs (3) and (4). The stopping criteria in inversion algorithms is usually based on fit between the observed and calculated data by the produced model. Typical fit or misfit estimator is the following N 2 obs − gical i=1 gi , (18) misfit = sqrt N obs 2 i=1 gi Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 5. The model at top depths 80 and 140 m (deep): (a,b) inversion of first synthetic data with compactness factor of 0.5, (c,d) inversion with compactness factor of 1, (e,f) inversion with compactness factor of 2; the rectangle indicate the location of true anomaly. Improving compact gravity inversion 551 Figure 7. Inversion of two body with the distance of 50 m and the top depth of 40 m data with standard deviation of 0.01 mGal random Gaussian noise added using: (a) compactness factor of 1, and (b) compactness factor of 2; the rectangles indicate the location of true anomaly. Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 6. Inversion of two body with the distance of 70 m and the top depth of 40 m data with standard deviation of 0.01 mGal random Gaussian noise added using: (a) compactness factor of 1, and (b) compactness factor of 2; the rectangles indicate the location of true anomaly. 552 M.H. Ghalehnoee, A. Ansari and A. Ghorbani Figure 9. Inversion of two body with the distance of 10 m and the top depth of 40 m data with standard deviation of 0.01 mGal random Gaussian noise added using: (a) compactness factor of 1, and (b) compactness factor of 2; the rectangles indicate the location of true anomaly. Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 8. Inversion of two body with the distance of 30 m and the top depth of 40 m data with standard deviation of 0.01 mGal random Gaussian noise added using: (a) compactness factor of 1, and (b) compactness factor of 2; the rectangles indicate the location of true anomaly. Improving compact gravity inversion 553 Figure 11. (a) Unconstrained inversion of third synthetic data (noise-free) with compactness factor of 0.5, (b) Constrained inversion with compactness factor of 1. where the superscripts ‘obs’ and ‘cal’ denote observed (or measured) and calculated data respectively. Besides, the approach explained here, does not depend strongly on the number of parameters (M), unlike compact inversion; it is advised to have fewer number of parameters by increasing the prism dimensions as far as possible thus having the most appropriate model and fast computational running. 2.5 Bound constraints To produce a physically meaningful model, refer to prior information obtained from geological maps, well-loggings, and laboratory tests, the density contrast of each prism must satisfy m min( j) ≤ m j ≤ m max( j) (19) If at an iteration mj exceeds one of its bounds, then it will be fixed at the violated bound. Instead of being calculated by eqs (9) and (10), the corresponding weight wj will be set at a very large value. The large value of wj will force the prism density contrast to be ‘frozen’ at the violated boundary, at least temporarily during the first few iterations. This way, the response of the modified parameter estimates will not fit the observations, which will in turn trigger the necessity for further parameter perturbation as a function of the misfit (Guillen & Menichetti 1984; Barbosa & Silva 1994; Grandis & Dahrin 2014). Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 10. Third synthetic example, three body with the density contrast of 1 and −1 gr cm−3 at various depths. 554 M.H. Ghalehnoee, A. Ansari and A. Ghorbani Figure 13. (a) Gravity data of Woodlawn ore body, New South Wales (Templeton 1981); Observed data (star marks) with calculated data (solid line), (b) unconstrained inversion result with α = 0, (c) constrained inversion result using α = 1, and (d) α = 2; ore body is located with black lines and dashed circles is the deep part of the anomaly. 3 I N V E R S I O N O F S Y N T H E T I C D ATA To illustrate the efficiency of the proposed approach, the algorithm was tested using three synthetic examples including examples for testing the depth resolution, examples for testing the lateral resolution and an example to test the negative and positive anomaly together. The first example is shown in Fig. 3, an anomaly in four top various depths as 20, 40, 80 and 140 m depths with the density contrast of 1.0 gr cm−3 . The inversion was started using a zero background model with the 5 m prism size. The anomaly is contaminated with pseudorandom Gaussian noise with zero mean and a standard deviation of 0.01 mGal. The model with damping factor of 0.01 and compactness factor of 0.5, 1 and 2 after 30 iterations are depicted in Figs 4 and 5. These models were imposed minimum and maximum bounds or constrained inversion (e.g. 0 ≤ m j ≤ 1). The models of depths 20, 40 and 80 m are coincident with the real anomaly in all value of α, but the deepest anomaly (the top depth of 140 m) is not shown properly by low value of compactness factor and is discernable only using α = 2. All of the models has sharp boundary in largest value of α which Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 12. (a) Unconstrained inversion of third synthetic data with standard deviation of 0.05 mGal random Gaussian noise added using compactness factor of 0.5, (b) Constrained inversion with compactness factor of 1. Improving compact gravity inversion 555 coincident with the real anomaly, excepting the two deep sources, because they lie deep, referring to the length of the profile and size (area or volume) of the gravity source, but the results are acceptable. The models showed in Figs 4 and 5 proved that compactness factor can be chosen depending on the depth of anomaly. In the deeper anomalies, the compactness factor shall be larger, and in shallower anomalies, the sharp model can be yielded using smaller compactness factor. The second example is two anomalies in various distances together (70, 50, 30 and 10 m distance) with the density contrast of 1.0 gr cm−3 to test lateral resolution. The data were added by zero mean and a standard deviation of 0.01 mGal Gaussian random noise. The inversion was started using a model with the prism size of 5 m and reference model of zero in various compactness factor. The model has been performed using damping factor of 0.1 and compactness factor of 1 and 2, the results of which are shown in Figs 6–9. These models are also constrained inversion (e.g. 0 ≤ m j ≤ 1). In Figs 6 and 7 anomalies are so far apart that they are easily recog- nizable even by a small compactness factor. But in Figs 8 and 9 that anomalies are closer, separation of anomalies, especially with small compactness factor is difficult or indistinguishable. However, with high compactness factor (α = 2), all anomalies can be detected except the nearest anomalies that cannot be resoluble in this level of noise (Fig. 9). The sources in Fig. 9 are so close that this distance is not recognizable even in the surface measurements with 10 m data point intervals. The third model is rather more complicated with three separated bodies with the negative and positive density at various depths (Fig. 10). The data are firstly supposed noise-free and inversion started using zero background with a model include 60×30 prisms, and dimensions of 5 m after 30 iterations (Fig. 11). Damping factor is chosen as small as 10−2 and compactness factor is 0.5 to have a smooth model (Fig. 11a). Inverted model is analogous to its original, excepting the deeper source, which has less density contrast here. By putting α = 1, to yield a sharp anomaly and imposing bounds (−1 ≤ m j ≤ 1), another model (Fig. 11b) was obtained. This model Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 14. (a) Gravity data of southeastern British Colombia ore body, Canada (after Brock 1973); observed data (star marks) with calculated data (solid line), (b) geological cross-section; the solid lines indicate drilling positions, (c) unconstrained inversion result using α = 0.0, and (d) inversion result with α = 0.0 and positivity constrained. 556 M.H. Ghalehnoee, A. Ansari and A. Ghorbani is largely similar to its true model especially with a high compact factor. The data were also modelled again with severe noise added as zero mean and a standard deviation of 0.05 mGal random noise using α = 0.5 for smooth and α = 1 for sharp and constrained inversion (Figs 12a and b). To avoid unwanted and false anomalies due to noise, the value of damping factor was increased hundreds times in comparison with noise-free data; therefore, the reasonable value of µ = 1 was chosen. The final model shown in Fig. 12, can be regarded as a satisfying model considering the high level of noise added. According to synthetic examples modelled by present algorithm, the method can be applicable even in complex and multi-sources anomalies at any depth and with any geometry, and has sufficient sensitivity and accuracy to model safe and noisy data. The inversion method is applied to three real gravity data. The first field example is taken from Templeton (1981), the gravity measurement was made in a line over the Woodlawn massive sulphide ore body, New South Wales, Australia (Fig. 13). The Woodlawn ore body is located in a well-populated, pastoral area of southeastern New South Wales, Australia. It consists of a shallow, massive sulphide deposit of complex shape. The Woodlawn ore body is best geologically classified as a stratiform volcanogenic copperlead-zinc deposit occurring within acid volcanic host rocks. The Woodlawn ore body is lenticular bodies in a wide zone of acid volcanics, pyroclastics and sediments which dip west. These rocks are considered to be transgressively invaded by thick lenses of dolerite which are conformable in places but break across the stratigraphy in other places (Whiteley 1981). An average density of the massive sulphides ore is of about 3.9 gr cm−3 and the density of the footwall ore is 2.9 gr cm−3 , which is the average density of the host rocks, and the target density contrast is 1 gr cm−3 . At the first stage, it is supposed that no information exists then inversion is started using a background of zero with 30×15 prisms of 5 m size, α = 0 and µ = 0.7 (Fig. 13b). Although the model is fully smooth, but its dip is rather similar with true ore body. Modeling was performed again with α = 1 and µ = 0.3 and imposed bounds of −0.01 ≤ m j ≤ 1 to have a more sharp result (Fig. 13c). As shown in Figs 13(b) and (c), the result of inversion is a vein with the dip a little less than 45◦ to the west. At the bottom of the layer from the depth of 100 m, the dip will be reduced. The anomaly is located at the depth of 15 m which this depth consists with the Gossan zone. Hematite, goethite, limonite and other iron oxides which produce Gossan decrease the density of original ore body near surface; therefore the density contrast is close to the background density contrast (i.e. 0 gr cm−3 ) at the top of the layer obtained by inversion. To model the deep part of the ore body, the compactness factor is increased to 2 (Fig. 13d). As shown, the deep anomaly indicates the deep part of the ore body rather admissible. Although gravity measurements are not spaced closely enough to delineate every small detail of the gravity field, and hence, the observations appear to have a random geological noise, the result obtained by inversion is acceptable enough. The second example is the data measured over sulphides mineralization, southeastern British Colombia, Canada (Brock 1973). Fig. 14, shows the geological section and related surface gravity data. The inversion was made with no bound constraints using α = 0.0 and µ = 1 with 45×15 prisms of 20 m size and null back- Figure 15. Inversion results of British Colombia ore body gravity data using bound constraints of 0 ≤ m j ≤ 1.5 and various compactness factor: (a) α = 0.5, (b) α = 1.0, (c) α = 1.5, and (d) α = 2.0. ground (e.g. mo = 0) regarding topography. The result shown in Fig. 14(c) is completely unconstrained inversion and there are both negative and positive density contrast. At the centre of the profile an anomaly of 0.65 gr cm−3 is seen smoothly which is located deeper in comparison with the real anomaly and less than real density contrast. However, produced anomaly shows a fairly poor estimate of location of the ore body and geometry of the model has almost no resemblance to the ore body. The result is rather inconsistent with our geologic information but needs more works such as drilling. Fortunately, six wells were being drilled along the profile being investigated. Our modelling was renewed according to the information obtained from drilling and detailed geological studies. Fig. 14(d), is a new model with the same compactness factor and imposing positivity bound constraints. This model is better than before but the obtained anomaly is bigger than real ore body and the value of gravity contrast is remained at 0.65 gr cm−3 . To obtain a sharp model, and even true geometry of the ore body, increasing compactness factor is necessary. Fig. 15 shows the inversion results using bound constraints of 0 ≤ m j ≤ 1.5 and various compactness factor from 0.5 to 2.0. By increasing the value of compactness factor, a more compact model is yielded as the best model with α = 2, which is obviously improved and very similar to the location and boundaries of real ore body. Fig. 16 shows the Bouguer anomaly and drilling tracks along a profile perpendicular to the strike of the Northwest Ore Body at Iron Mountain Mine, Missouri (Barbosa & Silva 1994). The host rocks are Precambrian igneous rocks consisting of pyroclastics, dacites, Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 4 I N V E R S I O N O F R E A L D ATA Improving compact gravity inversion 557 Figure 17. Residual Bouguer anomaly of the Northwest Ore Body (top); observed data are shown in stars and calculated data in solid line (after Barbosa & Silva 1994), unconstrained inversion result using α = 0.0 (bottom). rhyolites, and andesite porphyry. Paleozoic sediments consisting of a basal conglomerate, sandstones and limestones. The ore veins include hematite and magnetite respectively 76 per cent and 4 per cent, and 20 per cent of gangue minerals such as quartz, garnet and calcite (Murphy & Ohle 1968). Barbosa & Silva (1994) have removed the regional anomaly from the data graphically to enhance residual Bouguer anomaly. Fig. 17 shows the result of inversion using null background and the prisms with 10 m dimensions and full smoothness of α = 0 with no bounding imposed. The result shows a vertical anomaly Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 16. Bouguer gravity profile (above) and geological cros-section (below) of the Northwest Ore Body at Iron Mountain Mine, Missouri. After Leney (1966). 558 M.H. Ghalehnoee, A. Ansari and A. Ghorbani Figure 19. Residual Bouguer anomaly of the Northwest Ore Body (top); observed data are shown in stars and calculated data in solid line (after Barbosa & Silva 1994), bound constrained inversion result (0 ≤ m j ≤ 1.0) using α = 1.0 (bottom). with gravity contrast of 0.55 gr cm−3 which completely differs with true ore body. The ore body consists of two parts: a shallow horizontal and a deep vertical part which can be regarded as a complicated case study to invert. Therefore, choosing a suitable inversion parameter (e.g. compactness factor) is challengeable. By the above explanation, the data inverted again imposing bound constraints of 0 ≤ m j ≤ 1.0 and compactness factor of α = 0.5, α = 1.0, α = 1.5, and α = 2.0 in Figs 18–21, respectively. The result in Fig. 18 shows an anomaly which differs from true body, but the result is better than model in Fig. 17. By increasing compactness factor to α = 1.0, the second deep anomaly slightly is being shown (Fig. 19). Because the bodies are not too deep, the best result is shown when α = 1.5 (Fig. 20). In this model, two anomalies have density contrast of 1.2 gr cm−3 , which is near the true density and lies in the position rather agreeing with actual ore veins. When using α = 2.0 (Fig. 21), the result tends to be overcompact. In this situation, the yielded density contrast is more than true density contrast (which is 1.5 times more than real density contrast), but the geometry is rather consistent with the true geology. Consequently, depth weighting, kernel weighting and compactness weighting are effective to model deep anomalies acceptably either for sharp or smooth ones. Compactness factor plays a significant role in inverting sharp or smooth model which can be determined based on the existence or absence of prior information at least maximum and minimum bounds. 5 C O N C LU S I O N S We have modified the regular compact gravity inversion method based on developing new weighting functions. New weighting functions have been used based on the depth of the prisms, kernel matrix Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 18. Residual Bouguer anomaly of the Northwest Ore Body (top); observed data are shown in stars and calculated data in solid line (after Barbosa & Silva 1994), bound constrained inversion result (0 ≤ m j ≤ 1.0) using α = 0.5 (bottom). Improving compact gravity inversion 559 Figure 21. Residual Bouguer anomaly of the Northwest Ore Body (top); observed data are shown in stars and calculated data in solid line (after Barbosa & Silva 1994), bound constrained inversion result (0 ≤ m j ≤ 1.0) using α = 2.0 (bottom). and compactness weighting. The method is simple and can be easily run on a personal computer just a few seconds. The algorithm was tested by constrained and unconstrained models in synthetic and real data. It should be pointed out from synthetic and real data that the method explained here is quite general and reliable, accurate and precise, and can be applied in the inversion of any gravity data. The method also has an acceptable resolution in the depth and lateral direction. Besides, this algorithm could be developed to 3-D model and magnetic sources easily. REFERENCES Barbosa, V.C.F. & Silva, J.B.C., 1994. Generalized compact gravity inversion, Geophysics, 59(1), 57–68. Boulanger, O. & Chouteau, M., 2001. Constraints in 3D gravity inversion, Geophys. Prospect., 49, 265–280. Brock, J.S., 1973. Geophysical exploration leading to the discovery of Faro deposit, Can. Inst. Min. Metall. Bull., 66(738), 97–116. Grandis, H. & Dahrin, D., 2014. Constrained two-dimensional inversion of gravity data, J. Math. Fundam. Sci., 46(1), 1–13. Green, W.R., 1975. Inversion of gravity profiles by use of a Backus-Gilbert approach, Geophysics, 40, 763–772. Guillen, A. & Menichetti, V., 1984. Gravity and magnetic inversion with minimization of a specific functional, Geophysics, 49, 1354– 1360. Last, B.J. & Kubik, K., 1983. Compact gravity inversion, Geophysics, 48, 713–721. Leney, G.W., 1966. Field studies in iron ore geophysics, in Mining Geophysics, Vol. 1, pp. 391–417, eds Hansen, D.A., Heinrichs, W.E., Jr., Holmer, R.C., MacDougall, R.E., Rogers, G.R., Sumner, J.S. & Ward, S.H., Soc. Expl. Geophys. Li, Y.G. & Oldenburg, D.W., 1996. 3-D inversion of magnetic data, Geophysics, 61, 394–408. Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Figure 20. Residual Bouguer anomaly of the Northwest Ore Body (top); observed data are shown in stars and calculated data in solid line (after Barbosa & Silva 1994), bound constrained inversion result (0 ≤ m j ≤ 1.0) using α = 1.5 (bottom). 560 M.H. Ghalehnoee, A. Ansari and A. Ghorbani Silva, F.J.S., Barbosa, C.F. & Silva, J.B.C., 2009. 3D gravity inversion through an adaptive-learning procedure, Geophysics, 74(3), I9–I21. Silva, F.J.S., Barbosa, C.F. & Silva, J.B.C., 2011. Adaptive learning 3D gravity inversion for salt-body imaging, Geophysics, 76(3), I49–I57. Stocco, S., Godio, A. & Sambuelli, L., 2009. Modelling and compact inversion of magnetic data, Comput. Geosci., 35, 2111–2118. Templeton, R.J., 1981. Gravity surveys at Woodlawn, in Geophysical Case Study of the Woodlawn Ore Body, New South Wales, Australia, pp. 485– 494, ed. Whiteley, R.J., Pergamon Press. Vatankhah, S., Ardestani, V.E. & Renaut, R.A., 2014. Automatic estimation of the regularization parameter in 2D focusing gravity inversion: application of the method to the Safo manganese mine in the northwest of Iran, J. geophys. Engin., 11(4), doi:10.1088/1742-2132/11/4/045001. Whiteley, R.J., 1981. Geophysical Case Study of the Woodlawn Orebody New South Wales, Australia, Pergamon Press, 559 pp. Williams, N.C., 2008. Geologically-constrained UBC-GIF Gravity and Magnetic inversion with examples from the Agnew-Wiluna Greenstone Belt, Western Australia, PhD thesis, The University of British Colombia, Faculty of Geophysics. Downloaded from https://academic.oup.com/gji/article/208/1/546/2417079 by guest on 12 June 2022 Li, Y.G. & Oldenburg, D.W., 1998. 3-D inversion of gravity data, Geophysics, 63(1), 109–119. Mendonca, C.A. & Silva, J.B.C., 1994. The equivalent data concept applied to the interpolation of potential field data, Geophysics, 59(5), 722–732. Mendonca, C.A. & Silva, J.B.C., 1995. Interpolation of potential field data by equivalent layer and minimum curvature: a comparative analysis, Geophysics, 60(2), 399–407. Menke, W., 1989. Geophysical Data Analysis: Discrete Inverse Theory, Revised edition, Academic Press. Murphy, J.E. & Ohle, E.L., 1968. The iron mountain mine, iron mountain, Missouri, in Ore Deposits of the United States, 1933–1967, Vol. 1, pp. 287–302, ed. Ridge, J.D., Am. Inst. Min., Metallurg., Petr. Eng., Inc. Pilkington, M., 1997. 3D magnetic imaging using conjugate gradients, Geophysics, 62, 1132–1142. Portniaguine, O. & Zhdanov, M.S., 1999. Focusing geophysical inversion images, Geophysics, 64, 874–887. Silva, J.B.C. & Barbosa, V.C.F., 2006. Interactive gravity inversion, Geophysics, 71(1), J1–J9.