3D Gravity Inversion on Unstructured Grids

Abstract

Compared with structured grids, unstructured grids are more flexible to model arbitrarily shaped structures. However, based on unstructured grids, gravity inversion results would be discontinuous and hollow because of cell volume and depth variations. To solve this problem, we first analyzed the gradient of objective function in gradient-based inversion methods, and a new gradient scheme of objective function is developed, which is a derivative with respect to weighted model parameters. The new gradient scheme can more effectively solve the problem with lacking depth resolution than the traditional inversions, and the improvement is not affected by the regularization parameters. Besides, an improved fuzzy c-means clustering combined with spatial constraints is developed to measure property distribution of inverted models in both spatial domain and parameter domain simultaneously. The new inversion method can yield a more internal continuous model, as it encourages cells and their adjacent cells to tend to the same property value. At last, the smooth constraint inversion, the focusing inversion, and the improved fuzzy c-means clustering inversion on unstructured grids are tested on synthetic and measured gravity data to compare and demonstrate the algorithms proposed in this paper.

1. Introduction

Conventionally, gravity forward modeling is implemented by dividing the subsurface into a large number of rectangular prisms, and the gravity anomaly is calculated by integrating the gravitational effect due to each prism. However, structured grids have limited capability when modeling structures with complex shapes (such as topographic surfaces), and a staircased configuration from structured grids results in a loss of accuracy. In comparison, unstructured grids can easily be designed to conform to arbitrary straight-edged boundaries, topography, and subsurfaces. On the other hand, the local refinement or coarsening of the mesh, in regions where higher accuracy is required or lower accuracy is acceptable, can be performed without affecting the grids outside these regions. These properties can ultimately result in higher accuracy and lower computational cost.

For gravity forward modeling with unstructured grids, several methods have been developed. Talwani and Ewing [1], and Talwani [2] used numerical integration techniques by dividing models of arbitrary shape into polygonal prisms or laminas. Barnett [3] and Okabe [4] developed techniques for gravity and magnetic modeling due to homogeneous polyhedral bodies. Zhang et al. [5] and Cai and Wang [6] used the finite-element method to calculate the gravitational response of 3D rectangular grids, while May and Knepley [7] used the fast multipole method to tackle the same problem. Jahandari and Farquharson [8] used finite-volume and finite-element methods to gravity forward modeling and their results showed that the quadratic finite-element scheme is more accurate but also more computationally demanding scheme. The flexibility of unstructured tetrahedral meshes has been proved to be advantageous when complex model geometries have to be modeled. However, the numerical schemes (such as finite-volume and finite-element methods) handling such tetrahedral meshes can lead to less-accurate results in comparison with the structured, regular grids.

Three-dimensional (3D) regularized inversion has been a popular approach in right rectangular prisms [9,10,11,12]. In regularization inversion, a regularization term enables us to incorporate more information so that a smooth or sharp model can be yielded [13]. Li and Oldenburg [9,14] used an objective function that penalize both the deviation of the recovered model from a reference model and the structural complexity in three spatial directions. A focusing method is developed to reconstruct a compact and sharper image of the subsurface [11,15,16]. A fuzzy c-means clustering algorithm [17,18,19,20] is applied to induce the recovered physical properties to gather around clusters determined by some priori information. In addition, a depth weighting function is vital to recovering the model at a proper depth on structured grids. Li and Oldenburg [9] firstly introduce a depth weighting function in the Tikhonov formulation. Portniaguine and Zhdanov [16] have introduced a nonlinear depth weighting that equalizes the sensitivity of the observed data to cells located at different depths and horizontal positions. In a mathematically less rigorous way, Commer [21] directly applied a weighting function to the gradient of data misfit functional and the results revealed a great improvement in depth resolution. However, when these methods are directly used to gravity inversion on unstructured grids, a discontinuous and hollow model is more likely to be yielded. This arises because volume distinction of unstructured grids also has a considerable effect on model distribution besides depth variation of grids. This is more pronounced when the local refinement or coarsening of the unstructured mesh is used.

In this paper, by analyzing the gradient of objective functional, we find that the distribution of density is related to the conventional gradient formulation of objective functional. Subsequently, two mechanisms of depth weighting are compared and analyzed in detail: (1) a regularization-based strategy: applying a depth weighting to the regularization term [9,11,14,16]; and (2) gradient-based strategy: applying a depth weighting to the gradient of data misfit [21,22]. The comparative analysis indicates that the second method would improve gravity depth resolution more directly and effectively, and this improvement is less influenced by regularization parameters. Moreover, the regularization terms will have a restraining effect at all depths. Besides, we develop an improved fuzzy c-means clustering method combined with spatial constraints. Compared with previous works [17,18,19,20], it measures property distribution of inverted models in both the spatial domain and parameter domain, and it can provide a more continuous distribution of properties on unstructured grids.

To demonstrate the effectiveness of our algorithm, we tested both on a synthetic and a field dataset from the Vinton salt dome by Bell Geospace. The forward modeling is accomplished with the algorithm of Okabe [4] based on unstructured grids, which has been used in gravity researches [17,23,24]. We use TetGen [25] to generate tetrahedral grids. TetGen is an open-source software that generates the tetrahedral and Voronoï grids as required. The algorithms of TetGen are Delaunay-based, which ensures high quality for the whole grids by maximizing the minimum dihedral angle in the grids. The results show that the problem in gravity inversion with the anomaly concentration to the surface is largely solved. The inverted models show sharp boundaries and true positions.

2. Methods

The technical methods of this paper mainly include three issues: forward modeling algorithm on tetrahedral grids, Fuzzy c-means clustering constraints, and depth weighting. In this chapter, we will introduce a forward modeling algorithm on tetrahedral grids and Fuzzy c-means clustering constraints only, and the depth weighting will be described and analyzed in Section 3.

2.1. Gravity Forward Modeling Based on Tetrahedral Grids

In the right-handed Cartesian coordinate system, as shown in Figure 1, based on Barnett [3] and Okabe [4], by successively imposing thrice coordinate rotations on every edge of the polyhedral body, the gravity anomalies g of a homogeneous polyhedral body can be converted to the sum of the contributions from all edges, which can be expressed as

$$g=G\rho {\displaystyle \sum _i^n\left[\cos {\psi }_i{\displaystyle \sum _j^m{J}_j(i)}\right]}$$

(1)

In Equation (1), G is Newton’s gravitational constant, and ρ is the density of the cell. n is the facet number of the polyhedral body. m is the edge number of the ith facet. For tetrahedral grids, n = 4 and m = 3. ψ_i is a rotation angle of the ith facet, while J_j(i) is the contribution of jth edge in the ith facet, which is expressed as

$$J_j(i)={\left[\left(X\sin {\theta }_j-Y\cos {\theta }_j\right)\ln \left(X\cos {\theta }_j+Y\sin {\theta }_j+R\right)-2Z{\tan }^{-1}\frac{(1+\sin {\theta }_j)(Y+R)+X\cos {\theta }_j}{Z\cos {\theta }_j}\right]}_{(X_j,Y_j)}^{(X_{j+1},Y_{j+1})}$$

(2)

where

$$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{ccc}\cos \psi & 0& -\sin \psi \\ 0& 1& 0\\ \sin \psi & 0& \cos \psi \end{array}\right]•\left[\begin{array}{ccc}\cos \phi & \sin \phi & 0\\ -\sin \phi & \cos \phi & 0\\ 0& 0& 1\end{array}\right]•\left[\begin{array}{c}{x}^{\prime }-x_{\mathrm{obs}}\\ y^{\prime }-y_{\mathrm{obs}}\\ z^{\prime }-z_{\mathrm{obs}}\end{array}\right]$$

(3)

$$R={(X^2+Y^2+Z^2)}^{1/2}$$

(4)

In the above equations, φ, ψ, and θ are respectively three rotation angles for the jth edge, and can be calculated as follows [4],

$$\cos \phi =-\frac{S_{yz}}{{\left(S_{yz}^2+S_{zx}^2\right)}^{1/2}}$$

(5)

$$sin\phi =-\frac{S_{zx}}{{\left(S_{yz}^2+S_{zx}^2\right)}^{1/2}}$$

(6)

$$cos\psi =-\frac{S_{xy}}{{\left(S_{yz}^2+S_{zx}^2+S_{xy}^2\right)}^{1/2}}$$

(7)

$$\sin \psi =\frac{S_{yz}^2+S_{zx}^2}{{\left(S_{yz}^2+S_{zx}^2+S_{xy}^2\right)}^{1/2}}$$

(8)

$$\cos \theta =\frac{X_{j+1}-X_j}{{\left[{(X_{j+1}-X_j)}^2+{(Y_{j+1}-Y_j)}^2\right]}^{1/2}}$$

(9)

$$\sin \theta =\frac{Y_{j+1}-Y_j}{{\left[{(X_{j+1}-X_j)}^2+{(Y_{j+1}-Y_j)}^2\right]}^{1/2}}$$

(10)

In the above Equations (5)–(8), Syz, Sxz, and Sxy are the projected areas of the ith facet respectively onto the y-z-, z-x-, and x-y- planes. It should be pointed out that the vertices of the facet must be arranged clockwise in the X-Y-plane before Equation (3) is calculated. However, the generation of tetrahedral meshes is automatic and vertices are arranged and stored randomly. Therefore, we introduce a judgment and rearrangement of vertices before the gravity contribution of each facet is calculated. For more details about the rotations and judgment, please see Appendix A.

2.2. Fuzzy C-Means Clustering Inversion with Spatial Constraints

The inverse problem is solved by minimizing the Tikhonov parametric functional [13]. Fuzzy c-means clustering method (FCM) is used as the regularization term [18,19]. The objective function of gravity inversion is formulated as follows:

$$\mathsf{\Phi }={\mathsf{\Phi }}_{\mathrm{d}}+\lambda {\mathsf{\Phi }}_{\mathrm{m}}={\Vert Fm-d_{\mathrm{obs}}\Vert }^2+\lambda {\displaystyle \sum _k^p{\Vert W_{C_k}^{f/2}(m-C_k)\Vert }^2}$$

(11)

where F is a forward modeling operator. m are density parameters of all cells. d_obs are observed data. ${\mathsf{\Phi }}_{\mathrm{d}}$ is a data misfitting term, and ${\mathsf{\Phi }}_{\mathrm{m}}$ is a regularization term to constrain model structure. λ is a trade-off parameter. ${\mathsf{\Phi }}_{\mathrm{m}}={\displaystyle \sum _k^p{\Vert W_{C_k}^{f/2}(m-C_k)\Vert }^2}$ is the clustering term to classify the physical property values into distinct clusters $C=\left\{C_k|k=1,2,\cdots ,p\right\}$. In this paper, we assume that p and C are determined from a priori information. The parameter f, known as fuzzification parameter, is generally set to 2.0 [26]. Here, W_Ck is a diagonal matrix that measures the degree of the physical property values belonging to the kth cluster,

$$W_{Ck}=diag\left[\frac{{(m-C_k)}^{\raisebox{1ex}{$-2$}\!\left/ \!\raisebox{-1ex}{$(f-1)$}\right.}}{{\displaystyle \sum _k^p{(m-C_k)}^{\raisebox{1ex}{$-2$}\!\left/ \!\raisebox{-1ex}{$(f-1)$}\right.}}}\right]$$

(12)

Φ_m can effectively classify the physical property values into distinct clusters in the parameter domain, however, it has no effect on the spatial distribution of density. As a result, the property spatial distributions may be discontinuous and discrete. So, on structured grids, Sun and Li [18,19] added the smoothness regularization term to promote continuous in the spatial domain. However, it is difficult to balance the smoothness and clustering term and reconstruct a continuous model on unstructured grids.

Thus, we combine spatial constraints into the FCM term to help recover spatial continuity of physical property. We reset component wise of m as follows:

$${m^{\prime }}_j=\left(m_j+{\displaystyle \sum _l^{q_j}{m}_l}\right)/\left(q_j+1\right)$$

(13)

where q_j is set to be the adjacent cells number of jth cell, and it is usually equal to 4 except for the cells at boundaries. ${m^{\prime }}_j$ denotes the average of properties in the jth cell and its adjacent cells. The new Φ_m can measure the inverted model in both spatial domain and parameter domain. When the new Φ_m is minimized, it induces the properties of the jth cell and its adjacent cells to gather a certain cluster C_k. So, a continuous property distribution in spatial domain and a focused property distribution in parameter domain can be obtained.

In Equation (11), depth weighting function is not involved as it does in gravity inversion based on structured grids. Compared with structured grids, cell volumes are no longer the same as a constant. The dual effects of volume and depth difference of grids would result in a discontinuous and hollow model. While, the local refinement or coarsening of the mesh is an important advantage of unstructured grids, which will inevitably lead to different mesh volumes. So, it’s crucial to solve this problem before gravity inversion based on unstructured grids used widely in field data. We will discuss the depth weighting in the next section.

3. Analysis of Depth Weighting on Unstructured Grids

In gravity inversion, the depth resolution is a serious issue. Different kinds of depth weighting have been proposed to counteract the decay of kernel function for improving the depth resolution. Typically, Li and Oldenburg [9,27] introduced a weighting function of cell depth,

$$W_{\mathrm{L}\&\mathrm{O}}=\frac{1}{{(z+z_0)}^{\beta /2}}$$

(14)

where z is the cell depth, β is usually 2 for gravity inversion, z₀ is a constant.

Besides, Portniaguine and Zhdanov [16] proposed a weighting function as the integrated sensitivity matrix, i.e.,

$$W_{\mathrm{Zh}}=\mathrm{diag}{\left({\displaystyle \sum _i{F}_{ij}^2}\right)}^{1/2}$$

(15)

where F_ij denotes the element of the kernel function F.

On structured grids, the above two depth weighting functions are widely used in gravity inversions. While on unstructured grids, Figure 2 shows the variations of above 2 depth weighting functions with the depth.

Results in Figure 2 are calculated over an area of 500 m × 500 m × 500 m composed of tetrahedral grids. z₀ is chosen to be 20 in Equation (14). As volume of grids are different for tetrahedral meshes, the depth weighting $W_{\mathrm{Zh}}$ is normalized by cell volumes. In Figure 2, similar tendencies can be seen from two weightings. $W_{\mathrm{L}\&\mathrm{O}}$ is only related to depth, and all cells at the same depth have the same weightings. However, $W_{\mathrm{Zh}}$ are different for cells at the same depth, especially near the surface, which is caused by difference of cell volume. On unstructured grids, if cells at the same depth have the same weightings, it will result in that the cells with larger volumes have the priority to change their densities first to fit the observed data. So based on unstructured grids, $W_{\mathrm{Zh}}$ is better than $W_{\mathrm{L}\&\mathrm{O}}$ to improve resolution of gravity inversion.

Stabilized gradient-based inversion methods prove to be efficient for 3D imaging problems with highly parameterized models. Based on the principle of gradient-based inversion methods, setting the gradient to near zero can minimize the objective function efficiently [21]. The gradient of Equation (11) can be expressed as

$$\nabla \mathsf{\Phi }=\nabla {\mathsf{\Phi }}_{\mathrm{d}}+\lambda \nabla {\mathsf{\Phi }}_{\mathrm{m}}=F^T(Fm-d_{\mathrm{obs}})+\lambda {\displaystyle \sum _i^p{W}_z{W}_{Ci}^f(m-C_i)}$$

(16)

Firstly, let’s focus on the gradient of data misfit

$$\nabla {\mathsf{\Phi }}_{\mathrm{d}}=F^{\mathrm{T}}(Fm-d_{\mathrm{obs}})$$

(17)

It is well known that the gravity forward modeling operator F decays quickly as depth increases [9,28]. Therefore, when multiplying F^T with the data misfit Fm–d_obs, $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ will be less sensitive to depth, so densities of shallow cells will first change to fit d_obs. As a result, the inverse model will more easily concentrate to the surface in the traditional gradient-based inversions. So, from this perspective, the key to improve the depth resolution is how to counteract the decay of F^T in $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$.

In a conventional way, imposing depth weighting into Φ_m can improve depth resolution [9,14,16]. Then the gradient of the objective function can be expressed as

$$\nabla \mathsf{\Phi }=\nabla {\mathsf{\Phi }}_{\mathrm{d}}+\lambda \nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$$

(18)

where ${{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$ denotes the regularization term after depth weighting imposed.

Then the weighted $\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$ decays quickly and becomes less sensitive to depth as $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ does. It encourages that properties of shallow cells will first satisfy the constraint of ${{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$, such that the property differences of shallow cells in smooth inversion or the relative properties of shallow cells in focus inversion tends to zero. As a result, this feature encourages the properties of shallow cells to remain unchanged and counteracts the influence of F^T on shallow cells. So the depth resolution can be improved to some extent, but the improvement is affected by the balance of $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ and $\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$, i.e., trade-off parameter λ. With a small λ, $\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$ is not powerful enough to offset the decay of F^T in $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$, so the recovered model may be shallower than the true model. While with an overlarge λ, the recovered model may be deeper than the true model. What’s worse, as $\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$ decays quickly with increasing depth, its effect on deep cells is weakened.

Another solving method is to directly apply a depth weighting function to $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ instead of Φ_m in a mathematically less rigorous way [21]. In this way, the gradient of objective function can be expressed as

$$\nabla \mathsf{\Phi }=\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{d}}+\lambda \nabla {\mathsf{\Phi }}_{\mathrm{m}}=W_z{F}^{\mathrm{T}}(Fm-d_{\mathrm{obs}})+\lambda \nabla {\mathsf{\Phi }}_{\mathrm{m}}$$

(19)

where $W_z$ is a depth weighting function that increases quickly as depth increases, which is just adverse to the decay of F^T in $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$. Compared with Equation (18), $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ is multiplied by $W_z$ in Equation (19) instead of a depth weighted ${{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$. Therefore, $W_z$ can counteract the effect of F^T in $\nabla {\mathsf{\Phi }}_{\mathrm{d}}$ directly and effectively. Besides, different with $\nabla {{\mathsf{\Phi }}^{\prime }}_{\mathrm{m}}$, $\nabla {\mathsf{\Phi }}_{\mathrm{m}}$ is not related to depth, which means that λ would have less effect on the improvement of depth resolution as it does in Equation (18).

From the perspective of objective function gradient, above analyses indicate that the gradient-based strategy in Equation (19) is more effective to improve depth resolution than the method shown in Equation (18).

4. Results

To demonstrate the effectiveness of our algorithm, we tested both on a synthetic and a field dataset from the Vinton salt dome by Bell Geospace [29] and compared the improved FCM clustering inversion with focusing inversion and conventional FCM clustering inversion in the examples. In the following examples, the open-source software TetGen [25] is used to generate tetrahedral grids. The conjugate gradient method is implemented to solve the gravity inverse problems.

4.1. Inversions of Simulation Data

The synthetic gravity data were generated for a 3D homogeneous structure with a relative density of 1 g/cm³, as shown in Figure 3. The top bury depth is 100 m. The gravity data were calculated at the Earth’s surface on a grid with 50 m spacing and a total of 320 locations. We discretized the survey area (1000 m × 800 m × 600 m) into 237,298 tetrahedral grids with maximum cell volume of 5000 m³ when using the program TetGen [25]. Two depth weighting strategies with different inversion methods were respectively applied to recover the synclinal structure. A uniform relative density model of 0 g/cm³ was used as the initial model. A constraint on the upper-lower limit of density [0 g/cm³, 1.0 g/cm³] was used in these experiments. We carried out the FCM clustering gravity inversion with C = {0 g/cm³, 1 g/cm³} and p = 2. The inversion was taken as converged when the data misfit or relative model norm update was below their respective predetermined thresholds.

The inversion results of two weighting strategies without regularization terms (λ = 0) are shown in Figure 4, where the section passing the center in y-direction is shown only. The dashed boxes in the figure outline the true boundaries of the anomalous body. With no model structure constraints, the result with a conventional strategy in Equation (16) (as shown in Figure 4a) demonstrates a concentration near the surface, while it is improved by using gradient-based strategy in Equation (17) (as shown in Figure 4b). Then using different λ, we implemented smooth inversions, and the improved FCM clustering inversions with two strategies. Respectively, the results of two inversions are respectively shown in Figure 5 and Figure 6. It should be pointed out that, in order to better perform the influence of λ, we set λ as a constant in different inversions. In FCM clustering inversion, λ is a constant. From Figure 4, Figure 5 and Figure 6, one can find that the depth resolution improvement of the conventional strategy is closely related to the regularization parameter. Increasing λ can improve the depth resolution to some extent, while an improper λ will result in an incorrect depth, such that a shallow model was yielded with a small λ (in Figure 5a and Figure 6a). What’s more, since the regularization terms decay as depth increases, model constraints are weakened and the inverted models become blurred and fuzzy at depth in FCM clustering inversions (as shown in Figure 5a,b and Figure 6a,b). When λ further increases, the data fit becomes deteriorated. On the contrary, in Figure 5c,d and Figure 6c,d, the model is well recovered and the top depth is well determined with gradient-based strategy in two kinds of inversion methods. It indicates that the depth resolution is less related to the regularization parameter. Besides, the constraints of regularization terms are not weakening when depth increases. These results correspond with the above analyses of depth resolution.

In Figure 7, focusing inversion and conventional FCM clustering inversion are implemented with gradient-based strategy to compare with the improved FCM clustering inversion. We set a large initial λ, and then used the c to gradually reduce it. One can find that all these three methods can reconstruct a compact model. However, from Figure 7a,b, it is shown that focusing inversion and conventional FCM clustering inversion yield the inverted models into segments or with holes on unstructured grids. Focusing inversion is to reconstruct the model with minimum volume, and it has no constraint in spatial domain. In conventional FCM clustering inversion, one may increase smooth constraint to yield a continuous model, while the constraint of FCM clustering is weakened and the inverted model will not be focusing and compact. But in Figure 6, the models are more compact and no holes inside with the improved FCM clustering method. That is because the new method induces cells and their adjacent cells to tend to the same value except for the constraint of smooth regularization term.

Table 1. Inversion result list with different methods, strategies, and parameters.

Result	Inversion Method	Strategy	Trade-Off Parameter
Figure 4a	NaN	Conventional strategy	0
Figure 4b	NaN	Gradient-based strategy	0
Figure 5a	Smooth inversion	Conventional strategy	1 × 10−2
Figure 5b	Smooth inversion	Conventional strategy	5 × 10−2
Figure 5c	Smooth inversion	Gradient-based strategy	1 × 101
Figure 5d	Smooth inversion	Gradient-based strategy	5 × 101
Figure 6a	New FCM inversion	Conventional strategy	1 × 10−3
Figure 6b	New FCM inversion	Conventional strategy	3 × 10−3
Figure 6c	New FCM inversion	Gradient-based strategy	3 × 100
Figure 6d	New FCM inversion	Gradient-based strategy	5 × 100
Figure 7a	Focusing inversion	Gradient-based strategy	method of Lelièvre et al. [17]
Figure 7b	Conventional FCM inversion	Gradient-based strategy	method of Lelièvre et al. [17]

shows a clear perspective of inversion results with different methods, strategies, and trade-off parameters. The comparative analysis indicates that the second method would improve gravity depth resolution more directly and effectively, and this improvement is less influenced by regularization parameters. Moreover, the regularization terms will have a restraining effect at all depths.

4.2. Field Data Inversion

The survey area Vinton Dome is located in southwestern Louisiana, USA. We used a subset of the survey area for this study that covers approximately 16 km² in the middle of the survey area, which was acquired and processed by Bell Geospace in 2008 [29]. Before we inverted the data, a second-order trend surface was removed to get rid of background, and an upward continuation of 100 m was done to eliminate the shallow local anomalies. The processed gravity data is displayed in Figure 8.

Inversion area is 4 km × 4 km × 1 km with 169,949 cells and 1681 survey stations. The station interval was 100 m. The upper-lower limit of density was assumed to be [0 g/cm³, 0.55 g/cm³] and clusters C = {0 g/cm³, 0.55 g/cm³} was used in our inversions. Smooth inversion, conventional and improved FCM clustering inversion are implemented with a gradient-based strategy. We specifically point the x-, y-and z-axis to the East, North, and vertically downwards. Figure 9, Figure 10 and Figure 11 show the slices of the recovered model on x-y-, x-z-, and y-z-planes of different inversion methods.

Many surveys and researches have been done at the Vinton Dome. According to Coker et al. [30], the Vinton salt dome is composed of a well-defined massive cap rock extending above the salt rock. This cap rock is formed by gypsum and anhydrite that is embedded in sediments characterized by intercalated layers of sandstone and shale.

These researches indicate that the cap rock has a top depth of 160 m and a density of 2.75 g/cm³, while the underlying media have a density of 2.20 g/cm³ [31]. Besides, the sediment around slat dome has a density of about 2.0 g/cm³–2.37 g/cm³. Thus, we assumed the relative density of cap rock to be 0.55 g/cm³, and the observed gravity data were predominantly caused by the cap rock. From previous work in this area [30,31,32,33], it is known that the extension of the anomalous body is about 1200 m in the west-east direction, 900 m in the north-south direction, 200 m in z-direction, and the top depth is about 175 m.

The result of smooth inversion is blur and fuzzy in Figure 9, and the 3D extension size of the cap rock is much larger than the real geological model. In both Figure 10 and Figure 11, the maximum relative density is 0.55 g/cm³, and a more compact model is yielded. The extensions of the anomalous body are well consistent with the known information. But it is more internal continuous by using the improved FCM clustering inversion in Figure 11.

In brief, the results show that the problem in gravity inversion with the anomaly concentration to the surface is largely solved. The inverted models show sharp boundaries and true positions.

5. Conclusions

In gravity inversion based on unstructured grids, the dual effects of volume and depth difference of grids would result in a discontinuous and hollow model. However, the local refinement or coarsening of the mesh, an important advantage of unstructured grids, will inevitably lead to different mesh volumes. So, it’s crucial to solve this problem before gravity inversion based on unstructured grids used widely in field data.

From the view of the objective function gradient, we analyzed the cause of poor depth resolution and two strategies of depth weighting in gradient-based inversion methods. It indicated that the key to improve the depth resolution is how to counteract the decay of kernels in the data misfit gradient. The data-based strategy can improve the depth resolution more directly and effectively than the conventional strategy.

Besides an improved FCM clustering combined with spatial constraints is developed to obtain a continuous density distribution in the spatial domain and a focused property distribution in the parameter domain. It will reconstruct a more continuous inverted model on unstructured grids. With different regularization parameters, two weighting strategies are compared by applying smooth and the new FCM clustering inversions on synthetic data sets. The results show that the depth resolution of the gradient-based strategy is independent of regularization parameters and regularization terms does not weaken with depth increases. Compared with focusing inversion and conventional FCM clustering inversion, the improved FCM clustering inversion can provide a more continuous model.

At last, we implemented the gradient-based strategy and the improved FCM clustering inversion on field data sets. The results show that a compact model can be better consistent with the known information than smooth inversion. But it is more internal continuous by using the improved FCM clustering inversion.