Unstructured grid joint inversion of gravity and magnetic data based on equivalent surface source and wave number domain

Abstract

The present application is suitable for the field of geophysical exploration technology, and provides a non-structured grid joint inversion method for gravity and magnetic based on equivalent surface source and wave number domain, comprising the following steps: using surface source equivalence to optimize the calculation of gravity and magnetic kernel function matrix; the wave number domain method calculates the kernel function matrix of gravity and magnetic by combining the characteristics of the equivalent surface of non-structured tetrahedral subdivision. In practical application, the present application can more quickly calculate the gravity and magnetic kernel function matrix under non-structured tetrahedral grid subdivision, thereby improving the calculation efficiency of the joint inversion of gravity and magnetic under non-structured tetrahedral grid subdivision, making the joint inversion of gravity and magnetic under non-structured tetrahedral grid subdivision with higher inversion resolution and field source description capability have better development in the fields of large-area mineral resource exploration and fine structure exploration. Through simulation experiment verification, the efficiency improvement of kernel function calculation in the joint inversion of gravity and magnetic under non-structured tetrahedral grid subdivision can reach more than 90%.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration technology, and particularly relates to an unstructured grid gravity and magnetic joint inversion method based on equivalent surface source and wavenumber domain. Background Technology

[0002] Joint inversion of gravity and magnetic potential field data is a method for inverting the distribution of subsurface physical properties using observed gravity and magnetic data. It requires dividing the subsurface space to be inverted into closely spaced grid cells. A kernel function matrix is ​​established using a forward modeling formula for each grid cell at the surface observation point, and then a system of linear equations is established using the measured gravity and magnetic data at the observation point, the kernel function matrix, and the unknown physical properties at each grid cell. Solving this system of linear equations yields the physical property values ​​at the subsurface grid cells. In joint inversion of gravity and magnetic potential field data, two commonly used grid subdivision methods are structured hexahedral grids and unstructured tetrahedral grids. Unstructured tetrahedral grid subdivision can better characterize the irregular boundaries of undulating terrain and subsurface geological bodies and has higher inversion resolution. Therefore, to obtain better inversion results, current joint inversion of gravity and magnetic potential field data tends to favor unstructured tetrahedral grid subdivision.

[0003] The most time-consuming part of the joint inversion calculation of gravity and magnetic potential fields under unstructured tetrahedral mesh is the calculation of its kernel function matrix. Current kernel function matrix calculations utilize Gauss's theorem, decomposing the volume integral in the forward modeling formula of each tetrahedral mesh element into the sum of the surface integrals of each tetrahedral face. However, calculating these surface integrals requires multiple techniques, including coordinate rotation and Green's theorem, resulting in complex formulas and long computation times. Furthermore, its application presupposes a regular arrangement of mesh elements and that the observation point is located at the horizontal center of the top of the surface mesh element, imposing strict requirements on mesh generation and observation point location. Current rapid calculation methods rely on the structure and symmetry of the mesh generation, making them suitable for structured hexahedral mesh generation. Unstructured tetrahedral mesh generation lacks structure and symmetry, thus preventing the use of existing methods to improve the efficiency of joint inversion of gravity and magnetic potential field data under unstructured tetrahedral mesh generation. This severely limits the application of joint inversion methods for gravity and magnetic potential field data under unstructured tetrahedral mesh generation in large-scale resource exploration and its ability to finely characterize subsurface geological bodies. To address this, we propose a joint inversion method for gravity and magnetic fields based on an equivalent surface source and an unstructured grid in the wavenumber domain. Summary of the Invention

[0004] The purpose of this invention is to provide a joint inversion method for gravity and magnetic fields based on an equivalent surface source and the wavenumber domain, which aims to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution:

[0006] The unstructured grid gravity and magnetic field joint inversion method based on equivalent surface source and wavenumber domain includes the following steps:

[0007] The kernel function matrices of gravity and magnetism are optimized by utilizing the equivalence of surface sources; the wavenumber domain method is used to calculate the kernel function matrices of gravity and magnetism by combining the characteristics of equivalent surfaces of unstructured tetrahedral partitioning.

[0008] Furthermore, the kernel function matrix of gravity has the following matrix form:

[0009] Formula 1:

[0010] Formula 2:

[0011] In Equation 1, each element of the matrix is ​​calculated using the forward modeling formula for the l-th tetrahedron of unit density at the k-th observation point, i.e., Equation 2; in Equation 2, F k,i For the forward modeling value of the triangular face, Let f(obs(k),tetra) be the angle between the outward normal of the i-th face of the tetrahedral element and the z-axis. l (i) is the forward gravity formula for the triangular face of the i-th face of the l-th tetrahedron with unit density at the k-th observation point when the outward normal is aligned with the z-axis. Using the F values ​​of each adjacent tetrahedron in Equations 1 and 2... k,i The characteristic of opposites is represented by -F at the k-th observation point corresponding to the i-th face of the l-th tetrahedron. k,i Replace F in the adjacent tetrahedron of the l-th tetrahedron k,i Reduce F by nearly half k,i Calculation process.

[0012] Furthermore, a wavenumber domain calculation method is used to reduce the time complexity of calculating the kernel function matrix, therefore F k,i The two-dimensional plane data composed of all points is calculated using the wavenumber domain formula, as follows:

[0013] Formula 3:

[0014] In Equation 3, F is F k,i The data of the two-dimensional plane composed of all points, where G is the gravitational constant, G = 6.67 × 10⁻⁶. - 11 m 3 / (kg·s 2 ), z0 is the height of the observation surface, k x and k y Let x and y be the wavenumbers in the x and y directions, respectively, and k be the vector of wavenumbers in the wavenumber domain, with the absolute value of k being k. x ,k yThe square root of the sum of squares, j is the node number of the integration node used to calculate the definite integral using Gaussian numerical integration, J is the total number of Gaussian integration nodes used, γ is the component of the outward normal per unit length in the z-direction, and w j Let be the weight of the j-th node in the Gaussian integral, and i be the imaginary unit. The coordinates of the Gaussian nodes are in the spatial domain.

[0015] Furthermore, according to the Poisson equation, the formula for calculating the kernel function matrix of magnetism is composed of the derivative form of the formula for calculating the kernel function of gravity, and has the following relationship:

[0016] Formula 4:

[0017] Formula 5:

[0018] Formula 6:

[0019] Formula 7:

[0020] In equations 4-7, μ0 is the permeability in vacuum, which is 4π × 10⁻⁶. -7 H / m, δ and i are the magnetization declination and magnetization tilt, respectively; I0 is the tilt of the geomagnetic field; D0 is the angle between magnetic north and true north; V is the gravitational potential, which is the... The integral in the z-direction, and the integral and differentiation of V in the x, y, and z directions in the wavenumber domain, can all be achieved by the following formula:

[0021] Formula 8:

[0022] Formula 9:

[0023] Compared with the prior art, the beneficial effects of the present invention are:

[0024] In practical applications, this invention enables faster calculation of the gravity and magnetic kernel function matrix under unstructured tetrahedral mesh partitioning, thereby improving the computational efficiency of joint gravity and magnetic inversion under unstructured tetrahedral mesh partitioning. This allows for better development of joint gravity and magnetic inversion under unstructured tetrahedral mesh partitioning, which offers higher inversion resolution and field source characterization capabilities, in fields such as large-area mineral resource exploration and fine-structure exploration. The algorithm provided by this invention improves upon the shortcomings of existing methods, and simulation experiments verify that the efficiency improvement in kernel function calculation during joint gravity and magnetic inversion under unstructured tetrahedral mesh partitioning can reach over 90%. Attached Figure Description

[0025] Figure 1 This is a schematic diagram of the unstructured tetrahedral mesh partitioning and the common face of adjacent tetrahedrons in this invention.

[0026] Figure 2 The graph shows a comparison of kernel function computation time between the method of this invention and the traditional method under different numbers of partitioning units.

[0027] Figure 3 The following are examples of the inversion results of gravity and magnetic anomalies of three prisms in Embodiment 1 of the present invention: (a) shows the gravity forward modeling results and model distribution of the three prism models; (b) shows the magnetic forward modeling results and model distribution of the three prism models; (c) shows the density results of the joint gravity and magnetic inversion under unstructured mesh partitioning using the method of the present invention; and (d) shows the magnetization results of the joint gravity and magnetic inversion under unstructured mesh partitioning using the method of the present invention. Detailed Implementation

[0028] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0029] The specific implementation of the present invention will be described in detail below with reference to specific embodiments.

[0030] An embodiment of the present invention provides a joint inversion method for gravity and magnetic fields based on an equivalent surface source and the wavenumber domain using an unstructured grid, comprising the following steps:

[0031] The kernel function matrices of gravity and magnetism are optimized by utilizing the equivalence of surface sources; the wavenumber domain method is used to calculate the kernel function matrices of gravity and magnetism by combining the characteristics of equivalent surfaces of unstructured tetrahedral partitioning.

[0032] As a preferred embodiment of the present invention, the objective function formula for obtaining the density and magnetization vectors by the joint gravity and magnetic inversion method under unstructured mesh partitioning is as follows:

[0033]

[0034] Where A ρ and A M These are the kernel function matrices for gravity and the nucleus, respectively. M and ρ are the magnetization and density corresponding to each subdivision unit, d1 and d2 are the gravity and magnetic data at each observation point, and u1, u2, a1, and a2 are regularization factors used to balance the weights of multiple L2 norms; their values ​​are determined by the L-curve method. and The cross-gradient operator for density and magnetization is used to constrain the structural consistency of density and magnetization during inversion. W is the depth weighting function, a fundamental weighting function commonly needed in potential field inversion, used to balance the weights of subdivided units at different depths. The kernel function matrices for gravity and magnetism are A. ρ A MThis is the most time-consuming calculation part of the inversion process, and it is also the part that is improved in this invention.

[0035] As a preferred embodiment of the present invention, taking the gravity kernel function matrix as an example, its matrix form is as follows:

[0036]

[0037]

[0038] In equation (2), each element of the matrix is ​​calculated using the forward modeling formula for the l-th tetrahedron of unit density at the k-th observation point, i.e., equation (3); in equation (3), F k,i It is the forward modeling value of the triangle. It is the angle between the outward normal of the i-th face of the tetrahedral element and the z-axis, f(obs(k),tetra l (i) is the forward gravity formula for the triangular face of the l-th tetrahedron with unit density at the k-th observation point when the outward normal is in the same direction as the z-axis.

[0039] In unstructured tetrahedral meshing, most tetrahedra have four adjacent tetrahedra that coplanar with them, and the coplanar forms are as follows: Figure 1 As shown. When tetrahedrons are coplanar, the common face of adjacent tetrahedrons has the same corner coordinates for two adjacent tetrahedrons. In equation (3), f(obs(k),tetra l The calculated values ​​of (i) are the same, while the common surface has opposite directions of the outward normal in the forward formula of two adjacent tetrahedrons, i.e., equation (3). Therefore They are opposites. Therefore, when calculating the kernel function matrix A... ρ In this case, it is not necessary to calculate the forward modeling result of the unit density tetrahedron for each matrix element according to equation (3). F in k,i It can be replaced by the inverse of the forward modeling result of the triangular face of the adjacent tetrahedron.

[0040] After utilizing this equivalent surface source feature, it is only necessary to calculate the forward modeling result of each face under the unstructured tetrahedral mesh, and then use the opposites of each face to replace adjacent coplanar tetrahedra, thereby forming the tetrahedral forward modeling result and establishing the kernel function matrix A. ρ For unstructured tetrahedral meshes, the forward modeling of each triangular face is reduced from two calculations to only one, thus reducing the computational load by nearly half.

[0041] In a preferred embodiment of the present invention, the kernel function forward modeling under conventional unstructured tetrahedral mesh partitioning is calculated using analytical solutions, and its triangular face forward modeling value F k,iThe computation is complex and time-consuming. Replacing it with a wavenumber domain calculation method can effectively reduce the time complexity of the formula, lowering the time complexity of the kernel function matrix calculation from O(k,l) to O(l), thereby improving the efficiency of kernel function calculation in the inversion. This allows the gravity and magnetic joint inversion method under unstructured tetrahedral mesh partitioning to have better application effects in regions with large amounts of gravity and magnetic data. Therefore, F k,i The two-dimensional plane data (F) formed under all point conditions can be calculated using the wavenumber domain formula, as follows:

[0042]

[0043] Where G is the gravitational constant, G = 6.67 × 10⁻⁶ -11 m 3 / (kg·s 2 z0 is the observation surface height, k is the wavenumber vector in the wavenumber domain, and the absolute value of k is the wavenumber in the x and y directions (k x ,k y The square root of the sum of squares of the result, j is the integration node number used to calculate the definite integral using Gaussian numerical integration, J is the total number of Gaussian integration nodes used, γ is the component of the outward normal per unit length in the z-direction, and w j Let be the weight of the j-th node in the Gaussian integral, and i be the imaginary unit. The coordinates of the Gaussian nodes are in the spatial domain.

[0044] Therefore, by combining the wavenumber domain method with the characteristics of the equivalent surface of the unstructured tetrahedral partition, the gravity kernel function matrix can be calculated efficiently, and the same method can be used to obtain the calculation method of the magnetic kernel function matrix.

[0045] They have the following relationship:

[0046]

[0047]

[0048]

[0049]

[0050] Where μ0 is the permeability in vacuum, which is 4π × 10⁻⁶. -7 H / m, G is the gravitational constant, δ and i are the magnetization declination and magnetization tilt, respectively, I0 is the tilt of the geomagnetic field, D0 is the angle between magnetic north and true north, B x B y and B z This represents the three components of the Earth's magnetic field. V is the gravitational potential, which is a matrix element. The integral in the z-direction, and the integral and differentiation of V in the x, y, and z directions in the wavenumber domain, can all be achieved by the following formula:

[0051]

[0052]

[0053] Where, k x and k y denoted as x and y, respectively, and k is the vector of wavenumbers in the wavenumber domain.

[0054] Therefore, an efficient method for calculating the magnetic kernel function matrix under unstructured tetrahedral mesh can be determined by equations (4) to (10).

[0055] In this embodiment of the invention, the computational efficiency of kernel functions is compared by designing different partitioning methods for 160,000 observation points with different numbers of partitioning units. Finally, the computational time optimization effect of the method of this invention is compared with that of traditional analytical solution methods. Figure 2 As shown.

[0056] Example 1: The gravity and magnetic inversion of three underground prisms was calculated using the method of the present invention. The inversion results are as follows: Figure 3 As shown, the positions and boundaries of the three prisms are accurately obtained, and the computation time is reduced by 93.2% compared with the traditional analytical solution for the joint gravity and magnetic inversion under unstructured mesh partitioning.

[0057] The above are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make several modifications and improvements without departing from the concept of the present invention, and these should also be considered within the scope of protection of the present invention. These modifications and improvements will not affect the effectiveness of the implementation of the present invention or the practicality of the patent.

Claims

1. A joint inversion method for gravity and magnetic fields based on equivalent surface sources and wavenumber domain unstructured grids, characterized in that, Includes the following steps: The kernel function matrices for gravity and magnetism are optimized by utilizing the equivalence of surface sources; the wavenumber domain method is used to calculate the kernel function matrices for gravity and magnetism by combining the characteristics of equivalent surfaces obtained by unstructured tetrahedral partitioning. The kernel function matrix of gravity has the following matrix form: Formula 1: ; Formula 2: ; In Equation 1, each element of the matrix is ​​calculated using the forward modeling formula for the l-th tetrahedron of unit density at the k-th observation point, i.e., Equation 2; in Equation 2, For the forward modeling value of the triangular face, Let be the angle between the outward normal of the i-th face of the tetrahedral element and the z-axis. This is the forward gravity formula for the triangular face of the l-th tetrahedron with unit density at the k-th observation point when the outward normal is aligned with the z-axis direction; In each row of the matrix in Equation 1, there are multiple pairs of adjacent tetrahedrons sharing the same triangular face, and the paired elements correspond to the calculation in Equation 2. They have an opposite relationship; this characteristic can be used to pass through the l-th tetrahedron - The forward gravity modeling result for the triangular faces of the adjacent tetrahedrons of the l-th tetrahedron; The time complexity of calculating the kernel function matrix is ​​reduced by employing a wavenumber domain calculation method. The two-dimensional plane data composed of all points is calculated using the wavenumber domain formula, as follows: Formula 3: ; In Equation 3, F is The data of the two-dimensional plane composed of all points, where G is the gravitational constant, G = 6.67 × 10⁻⁶. -11 m 3 / (kg∙s 2 ), z0 is the height of the observation surface, k x and k y Let x and y be the wavenumbers in the x and y directions, respectively, and k be the vector of wavenumbers in the wavenumber domain, with the absolute value of k being k. x , k y The square root of the sum of squares, j is the node number of the integration node used to calculate the definite integral using Gaussian numerical integration, J is the total number of Gaussian integration nodes used, γ is the component of the outward normal per unit length in the z-direction, and w j Let be the weight of the j-th node in the Gaussian integral, and i be the imaginary unit. The coordinates of the Gaussian nodes are in the spatial domain.

2. The unstructured grid gravity and magnetic field joint inversion method based on equivalent surface source and wavenumber domain as described in claim 1, characterized in that, A method for calculating the magnetic kernel function matrix is ​​obtained by combining the wavenumber domain calculation method with the characteristics of the equivalent surface of the unstructured tetrahedral partition. According to the Poisson equation, the formula for calculating the magnetic kernel function matrix is ​​a combination of the derivative forms of the formula for calculating the gravity kernel function, and has the following relationship: Formula 4: ; Formula 5: ; Formula 6: ; Formula 7: ; In equations 4-7, It is the permeability in vacuum, which is H / m, Let i and i be the magnetization declination and magnetization tilt, respectively, and I0 be the tilt angle of the Earth's magnetic field. is the angle between magnetic north and true north; V is the gravitational potential, which is the... The integral in the z-direction, and the integration and differentiation of V in the x, y, and z directions in the wavenumber domain, are achieved by the following formula: Formula 8: ; Formula 9: .

Images