Abstract
Modeling plays an important role in engineering exploration. However, the traditional rectangular mesh cannot meet the requirement because the actual shape of the model may be very complex. Therefore, this paper chooses the Delaunay grid for modeling and applies this method to magnetic numerical simulation. At the same time, we make use of the advantages of the Delaunay grid and choose to use grid encryption technology in the boundary and complex area of the model. Then, the formula for calculating the magnetic anomaly of tetrahedron is deduced in detail, and the magnetic anomaly with surface fluctuation is calculated. For a synthetic data model, Delaunay grids with different densities are used to model and calculate surface anomalies. The two results are compared with the analytical solution of the model. The results show that the method has high accuracy. Then, the method is applied to the actual geological body modeling in the Jinchuan mining area. According to local needs, a three-dimensional model with inhomogeneous grid density is established and the surface magnetic field of the model is calculated. Finally, the simulation data are compared with the measured data. The results show that the Delaunay grid modeling method has strong applicability.
1. Introduction
The information of the Earth’s magnetic field plays an important role in the research of the Earth’s structure. The magnetic forward problem and the inversion problem are the core research contents of the magnetic data. The former mainly studies the magnetic anomaly characteristics and distribution caused by field sources or geological bodies such as different shapes, occurrences, and field source magnetic susceptibility, and the calculation process is forward calculation; the latter mainly calculates the geology based on the distribution of magnetic anomalies, the magnetic susceptibility distribution of the body, or the field source, and the calculation process is inversion calculation [1]. In the process of solving geophysical problems, the mathematical modeling of the target geological body is a very important link [2].
In the study of gravity and magnetic calculation and seismic imaging, the widely used structural elements are generally quadrilateral in two-dimensional space and generally hexahedron in three-dimensional space [3–6].
When constructing a simple geological body, this structural unit is indeed very effective, but once the shape of the geological body is relatively complicated, the effect of simply using a quadrilateral or a hexahedron is not ideal. If the structural unit is encrypted, the calculation amount is increased.
Compared to the commonly used quadrilateral or hexahedral structural elements, triangular (two dimensional) or tetrahedral (three dimensional) is more flexible and convenient when defining complex geologic bodies. We use the Delaunay method to fill the two with triangular (tetrahedron) as the basic structural unit. Dimensional (three dimensional) space, eliminates the structural constraints of nodes in the structural mesh, and the controllability of nodes and elements is good, so the boundary can be processed well, which is suitable for simulating real complex shapes.
It is very easy to form a complex shape geological body using polygonal elements. The degree of approximation depends on the number of polygon faces and the choice of vertices. It also has a very good effect on the correction of terrain. Figure 1 shows four methods for simulating surfaces. We can see the effect of different grid discretization [7].
A nonuniform geological body can be divided into several smaller homogeneous geological body combinations. Therefore, the calculation of magnetic anomalies caused by homogeneous polyhedra is more practical and important than other types of methods, especially for high-resolution gravity measurements [8]. With the rapid development of computing facilities, software technology, and numerical methods, the solvability of applied mathematics and geophysical numerical simulation is greatly enhanced. Due to the complex geometry of the underground structure, any practical inversion problem needs to solve the forward problem similar to the underground structure. Therefore, many numerical methods have been introduced into the field of geophysics, such as finite difference, finite element, integral equations, volume integrals, boundary integrals, and hybrid methods [9].
Not only magnetic exploration, but also the Delaunay mesh can be used in other geophysical methods. In 2012, Liu used Delaunay grids for numerical simulation of seismic waves and calculated them using the finite difference method [10]. Li and Zhang used two-dimensional and three-dimensional Delaunay techniques for gravity modeling [11, 12].
After the three-dimensional geological body is decomposed into a finite number of tetrahedrons, the simulated magnetic field is calculated according to the sum of the surface responses of all tetrahedrons.
The advantages of this method in the modeling of complex irregular regions are as follows:(1)Nonuniform grid computing which can be used to accurately describe the geological body, and the calculation results are accurate.(2)The boundary of the model can be performed according to the needs of the model. The engraving and local area are mesh encrypted, while the mesh density of other areas can be loosened according to the characteristics, thereby saving the calculation time of numerical simulation.
In this paper, a three-dimensional geological body is divided into tetrahedrons by means of the finite element method, and the gravity numerical simulation method of the three-dimensional terrain with the gravity center volume method is proposed and realized. The algorithm realizes the modeling of complex structures with terrain and irregular density using unstructured tetrahedral meshes. Finally, numerical examples and field examples are used to verify the results.
2. Mesh Generation
2.1. Basic Concepts
Euclidean space has a geometry consisting of n points that do not coincide with each other. For each point , there is always a continuous space in the vicinity:
We call a unit of the point . units of all points form the graph, as shown in Figure 2, which is in the two-dimensional situation. Points are connected with adjacent relationships in the graph. The vertices of any triangle have circumscribed circles, and the circumscribed circle does not contain other vertices. For the three-dimensional case, the circumscribed circle attribute is extended to the three-dimensional case and is called the circumscribed ball. The larger the angle of the minimum inner angle of the Delaunay triangle, the better the quality of the triangle unit. When dealing with a tetrahedron of a three-dimensional Delaunay mesh, it is much more difficult than a two-dimensional mesh. The inner angle of a tetrahedron does not necessarily perform this characteristic well according to the distribution of spatial point sets. Therefore, a properly distributed three-dimensional point set is required.
2.2. Grid Quality Optimization
Jonathan Richard Shewchuk in 2014 proposed that the Delaunay subdivision is a technique for generating triangular meshes and is suitable for interpolation, finite element, and finite volume methods. The triangle angle in its model should not be too small or too large. The initial mesh elements generated by the Delaunay algorithm have good geometry, but local regions also have poorly shaped cells. For example, when the input surface mesh contains sharp corners or elongated triangles, poor quality units are usually produced near these areas. If these poor quality units are not effectively improved, it will inevitably affect the numerical simulation. In addition, inappropriate cell size and interior point generation strategies can result in poor quality transition units in areas of the dense mesh transitions. Grid optimization is usually a necessary step after mesh generation and before numerical simulation. It plays a vital role in ensuring the reliability and accuracy of numerical simulation. A good way to evaluate the mesh quality is to calculate the ratio of the circumscribed circle radius to the shortest side of the unit. The quotient of the circumscribed circle radius and the shortest side length of the triangle is a measure of the natural improvement of the Delaunay refinement algorithm. The smaller the ratio, the better the quality. We first need to set a value. When the ratio of the radius of the circumcircle of the triangle or tetrahedron to the length of the short side is greater than this value (see Figure 3(a)), the grid updates it (see Figure 3(b)). When the longest side is at the boundary of the model, the above method cannot be used (see Figure 3(c)). When this happens, we use other methods to update the grid, and select the midpoint of the longest side of the triangle t as the new insertion endpoint. Then the midpoint is connected to the other endpoints around to complete the triangle mesh update (see Figure 3). The 3D model uses the face center of the largest area of the tetrahedron [13].
(a)
(b)
(c)
(d)
Before the susceptibility model is meshed, the mesh can be encrypted in some areas according to actual conditions, such as the boundary area of the model that is difficult to describe or the place where the abnormal body is small. By appropriately increasing the number of vertices and updating the mesh, the effect of mesh local encryption is achieved [14]. The encrypted mesh can describe the region more finely, and the accuracy of the magnetic calculation is also guaranteed. In the background area, the number of vertices is appropriately reduced, and the area can be sparsely meshed, which can reduce the memory usage of the computer and save the numerical simulation calculation time (see Figure 4).
3. Forward Calculation
After decomposing the three-dimensional susceptibility model into a series of tetrahedrons, we calculate the response of each tetrahedron to the surface. It is considered that the physical parameters of the unit tetrahedron are uniform. The sum of the calculated values of all tetrahedrons is the surface response of the anomaly:
In equation (2), ∆g represents the vertical component of the gravity anomaly at a certain point on the surface of the geological body q, G is a gravitational constant, is the target geological body density, are the coordinates of stations on the surface, is the vertical distance between the target weight and the surface measurement point, and is the coordinate of the center of gravity (see Figure 5).
Therefore, the surface gravity response of n geological body can be expressed asin which
The expression of the tetrahedron is
Equation (9) can be obtained from equation (8) as
Equation (9) is expanded to obtain the following equation:
According to the coordinate expression of the vector quantity product and the definition of the vector quantity product, the following equation is obtained:
According to the cosine theorem, the following equation is obtained:where p is the modulus of OA, q is the modulus of OB, and r is the modulus of OC. Equation (13) is obtained as follows:
By introducing equations (13) and (12) into equation (10), the following equation is obtained:
4. Synthetic Data
In order to test the method in this paper, we chose to use the sphere model for calculation. The field parameters were total intensity of 40,000 NT, and the magnetic susceptibility of the sphere is 0.01 and that of the background is 0. The sphere has a radius of 10 meters, and the center of the sphere is located at a height of −50 meters. The station is on an undulating surface with an average height of 0 (see Figure 6).
(a)
(b)
(c)
After calculating the magnetic anomaly by using the gravity center volume method proposed in this paper, the results of calculation need to be checked, otherwise the reliability of the magnetic anomaly field obtained by this method cannot be guaranteed. Since the sphere is a regular body, the magnetic anomaly field calculated by the integral method in this paper can be considered as an analytical solution. At the same time, the model is characterized by different number of point sets. Then, mesh generation is performed. The parameters of the two mesh generation are shown in Table 1.
We compared the forward results of the tetrahedral mesh generation model with the analytical solutions; it can be seen that the gravity field of refined mesh calculation is closer to the analytical solution, and the error of the rough mesh is relatively large, but it still meets the requirements of numerical simulation. Therefore, when calculating the actual data, we can use the method of local point set encryption to improve the accuracy of the magnetic anomaly field. This not only satisfies the description of local geological bodies, but also saves a lot of time for gravity numerical calculation (see Figure 7).
(a)
(b)
(c)
(d)
(e)
At the same time, the anomalous curves under different magnetic declination angles are compared, and the fluctuation surface of the line is y = z = 0 (see Figure 8). Thus, the proposed algorithm is effective. The more tetrahedrons there are, the closer the results are to the analytical solution (Figure 9).
(a)
(b)
(c)
(d)
(a)
(b)
The number of tetrahedrons in refined grids is more than ten times larger than that in rough grids, and the time of numerical simulation is also tens of times longer. After the model is discretized by rough grids, its internal angle range is 72.65–132.86, and after the refined grids are discretized, the internal angle range is 72.56–138.89, among them, the maximum value of the internal angle after the refined grids are discretized is larger. The reason is that after the refined grids are depicted, more tetrahedrons will be generated and the quality of the refined grids will be worse. There is a greater chance of decency. However, from the histogram, it can be seen that the precision mesh model has a larger proportion of small angle tetrahedrons than the rough mesh model. For the gravity numerical simulation method proposed in this paper, the smaller the angle of the tetrahedron is, the closer the tetrahedron is to the equilateral tetrahedron and the more accurate the numerical results are. Therefore, although there are some tetrahedrons with larger inner angle in the refined mesh discretization model, it still has better description accuracy compared with the rough mesh discretization model. Therefore, increasing the density of the tetrahedral mesh appropriately is of great significance to the calculation of magnetic anomalies (see Table 1 and Figure 7).
5. Real Data
5.1. Geological Background
The Jinchuan ore-bearing ultrabasic intrusive rock mass is located at the northeastern end of the Longshoushan uplift belt on the southwestern edge, and the mining area is located between 38°N and 39°N. The deep fault (F1) in the northern part of the Longshoushan uplift is adjacent to the tidal fault basin of the Alashan block, while the deep fault in the south is adjacent to the corridor transition zone of the Early Paleozoic fold belt in the North Qilian. The Jinchuan mining area is called the III, I, II, and IV mining areas from west to east. The size, shape, occurrence, and ore storage of rock masses in different mining areas are different. The main rock-forming minerals are olivine, clinopyroxene, orthopyroxene, and a small amount of plagioclase and amphibole. The olivine is mainly composed of noble olivine and a small amount of forsterite. The main rock types are pure peridotite, lherzolite-bearing, lherzolite, plagioclase lherzolite, and olivine lherzolite (see Figure 10).
There are many section models in the four mining areas of the Jinchuan deposit. The depth and lithology information of the section model are determined by drilling data, and the three-dimensional geological model of the Jinchuan deposit can be obtained by spatial interpolation of the two-dimensional section model. The 3D magnetic susceptibility model of the Jinchuan deposit can be obtained by adding the magnetic susceptibility parameter to the geological model of the Jinchuan deposit.
5.2. Modeling of Magnetic Susceptibility
Firstly, the magnetic susceptibility model of the Jinchuan Copper-Nickel Mine No. 2 is established (see Figure 11(a)). According to the information of the section map of the mining area, the mineralization zone and the rich and poor orebodies are buried in the surrounding rock, the section is controlled by drilling, and its accuracy can be guaranteed. At the same time, there are other cross sections and drilling information in this area. Although most of the ore bodies and mineralization zones in other cross sections still extend beyond the drilling control area and lack some data information, they can still be used as a reference standard for modeling. Therefore, the modeling of this cross section is scientific. In order to better depict the boundary of rock orebody and the fine structure of rock orebody and increase the number of vertices at corresponding positions, the magnetic susceptibility model with uneven grid density can be obtained. The magnetic susceptibility model of the Jinchuan mining area is established based on known data (see Figure 11(b)).
(a)
(b)
The green regional density is about 0.7SI, and the main component is basic peridotite. The orebody has a magnetic susceptibility of 0.7 and a geomagnetic field strength of 55,000 NT. Magnetic declination is 56°, and magnetic dip is −1.7°. These parameters are mainly obtained from field acquisition experiments. Finally, the magnetic forward simulation of the mining area is carried out using the method proposed in this paper. The magnetic anomaly is calculated by the method presented in this paper (see Figure 12).
At the same time, we used measured data to compare with numerical simulation (see Figure 13). The measured data are collected in the field of the Jinchuan deposit. The blue line in the figure is the calculated anomaly curve, and the red line is the measured data. It can be seen that the calculation result obtained by using the algorithm of this paper is close to the real data, and the effect is good. Considering the influence of residual magnetism, the calculation error can be considered well.
(a)
(b)
6. Conclusion
Delaunay mesh generation technology is introduced into the mesh generation of gravity numerical simulation in order to effectively characterize the model boundary and refine local features. The gravity center coordinate and volume formula of tetrahedron are used to simulate a three-dimensional magnetic anomaly field. The following conclusions can be drawn from the model test and the analysis of calculation results:(1)Delaunay mesh is an efficient mesh generation algorithm for discrete target geological bodies. This algorithm can not only consider the complex structure of the target body, but also control the grid resolution of the geological model. It can achieve the effect of local encryption, save computational memory, and improve computational efficiency. It has strong adaptability.(2)Delaunay grid generation method can optimize and control the grid according to the need in the process of its generation, so as to improve the quality of grid generation. The results of synthetic data show that a large proportion of tetrahedron angles meet the requirements of the refined mesh model. It provides a guarantee for the accuracy of magnetic numerical simulation.(3)In practical application, the Delaunay mesh generation method can effectively divide complex practical models. At the same time, the abnormal curves calculated by numerical simulation correspond well with the actual abnormal curves. Therefore, it is an effective numerical simulation method to use the Delaunay mesh method to divide the complex geological body model and then calculate the gravity anomaly response using the center of gravity and volume algorithm.
Data Availability
The data used to support the findings of this study have not been made available because the project to which the data belong is still in progress.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Key Research and Development Program of China (2016YFC0600505) and the National Natural Science Foundation of China (41874136).