A meshless gravity and magnetic forward modeling method based on spatial equivalence

Abstract

The present invention relates to the field of gravity and magnetic forward modeling, specifically a gridless gravity and magnetic forward modeling method based on spatial equivalence, comprising the following steps: setting observation points and model positions; setting the positions of discrete nodes and background grids; setting discrete nodes within the model as relevant physical properties; calculating the local kernel matrix and right-hand side terms according to the forward modeling formula and the equivalence principle; calculating the preprocessing matrix and weighting it; and implementing the forward modeling of the gravity and magnetic field according to pseudocode. The present invention improves the gridless method based on the spatial equivalence method, realizes a highly efficient forward modeling of gravity and magnetic data, and provides an effective processing method for efficiently processing gravity and magnetic data in large areas and undulating terrain. It can reduce memory usage by dozens of times and increase computational efficiency by several times, and the finer the subdivision, the more significant the improvement effect.

Description

Technical Field

[0001] The present invention relates to the field of gravity and magnetic forward modeling, and in particular to a gridless gravity and magnetic forward modeling method based on spatial equivalence. Background Art

[0002] Mineral resources are a crucial material foundation for social development and national economic construction. With the growth of the national economy, the demand for mineral resources is increasing. However, while mineral resources in plain areas have been largely explored, shallower resources are becoming increasingly depleted. Currently, prospecting targets are shifting from thickly covered plains to shallower areas with complex topography and geological structures, such as those in mid- and high-altitude mountains, and from shallow depths of less than 500 meters to deeper areas. This not only poses challenges to geophysical exploration equipment, but also places higher demands on the accuracy and efficiency of geophysical data processing technology.

[0003] Forward modeling is the foundation of physical property inversion using gravity and magnetic data. To perform inversion calculations, the forward response relationship between the distribution of subsurface physical properties and the data at the observation points must first be established based on the corresponding mathematical model. Inversion calculations generally employ an iterative approach, so the forward response of the subsurface physical properties must be continuously calculated during the inversion process to determine the residual between it and the actual data. Therefore, the accuracy and efficiency of the forward modeling are crucial determinants of the accuracy and efficiency of the inversion calculations.

[0004] Currently commonly used geophysical property inversion methods typically discretize the subsurface space of the study area into a series of cells and assume constant physical properties within each cell. Therefore, the resolution of the inversion results depends not only on the data resolution but also, to a certain extent, on the size and distribution of the cells. The mainstream approach now involves direct integration of the physical equations satisfied by the gravity and magnetic fields in the spatial or wavenumber domain. In this case, when the number of observation points is 100 × 100 and the number of cells is 100 × 100 × 100, the kernel matrix memory usage is approximately 74.5 Gb, and the kernel matrix generation time (i.e., the forward modeling process) is approximately 9 hours.

[0005] Forward modeling in gravity and magnetic property inversion can also be accomplished by solving the partial differential equations satisfied by the gravity and magnetic fields. Although this method is slightly less accurate than the exact analytical formula method, it is more flexible in implementing forward modeling of complex models, occupies less memory, and is more efficient during inversion calculations. Therefore, it has also been widely studied. Corresponding methods include: finite element method, finite volume method, finite difference method, and meshless method. Although this method can reduce memory usage, it is relatively inefficient. Summary of the Invention

[0006] The algorithm proposed in this invention can greatly reduce memory usage and improve computational efficiency based on traditional gridless methods. The present invention provides the following technical solutions:

[0007] A gridless gravity and magnetic forward modeling method based on spatial equivalence includes the following steps:

[0008] The first step is to set the observation point and model position;

[0009] The second step is to set the positions of discrete nodes and background grid;

[0010] The third step is to set the discrete nodes in the model as related physical properties;

[0011] The fourth step is to calculate the local kernel matrix and the right-hand side term according to the forward modeling formula and the equivalent principle. The local kernel matrix of gravity and magnetic forward modeling is:

[0012] , for gravity forward modeling, the right-hand side term is in the form of:

[0013] ,

[0014] in, G is the gravitational constant; For the number ( i , j , k )'s background grid's length, width, and height; 、 For the number ( i , j , k )'s background grid interpolation function and β Directional derivative, Yes The general representation of β ∈{ x , y , z}; Yes The general representation of β ∈{ x , y , z}; For the number ( i , j , k )’s background grid supports the density of nodes within the domain; W and W E is the Gauss-Legendre integral weight coefficient, W is the 8th-order identity matrix, W E is the 4th-order identity matrix, ;

[0015] , , They are x , y , z Total number of directional background grids; n x , n y , n z The study areas x , y , z Total number of directional background grids; n xe , n ye , n ze The expansion area x , y , z The total number of directional background grids. For magnetic forward modeling, the right-hand side term is:

[0016] ,

[0017] in, is the magnetic permeability in vacuum; For the number ( i , j , k )’s background grid supports the magnetization of nodes within the domain, Yes , , , , The general expression of For the number ( i , j , k )’s second-order derivative of the interpolation function of the background grid; β 1, β 2∈{ x , y , z}; cos( φ Mβ1 ) is the value of cos( φ Mx ), cos( φ My ), cos( φ Mz ) is a general representation of cos( φ Mβ1 ) is the magnetization intensity in x , y , z Direction cosines of direction; cos( φTβ2 ) is the value of cos( φ Tx ), cos( φ Ty ), cos( φ Tz ) is a general representation of cos( φ Tβ2 ) is the Earth's magnetic field x , y , z Direction cosines of direction;

[0018] The fifth step is to calculate the preprocessing matrix and weight it. The calculation formula is: , , M 1, M 2 are the left and right preprocessing matrices respectively; D is the kernel matrix; x is the quantity to be solved; b is the right-hand term, diag() means extracting the diagonal elements of the matrix;

[0019] The sixth step is to realize the forward modeling of gravity and magnetic field based on pseudo code.

[0020] As a further solution of the present invention: discrete nodes include regularly distributed nodes and irregularly distributed nodes, and irregular nodes are not set in the expanded area of ​​the background grid and the outermost layer of the study area, and regular nodes in the study area are evenly arranged.

[0021] As a further solution of the present invention: the background grid is a regular cube, and the corner points of the background grid are regularly distributed nodes.

[0022] As a further solution of the present invention, the relationship between the local kernel matrix in the region where the support domain morphology of the background grid and the relative position relationship of the adjacent nodes are consistent is:

[0023] ,in, K (( a , b , c ), ( d , e , f )) indicates that the number in the matrix is ​​( a , b , c ) and the node numbered ( d , e , f )'s node relationship value; 1≤ a , d ≤ N x , 1≤b , e ≤ N x , 1≤ c , f ≤ N x ; N x = n x +2 n xe , N y = n y +2 n ye , N z = n z +2 n ze They are x , y , z Total number of directional background grids, n x , n y , n z The study areas x , y , z Total number of directional background grids; n xe , n ye , n ze The expansion area x , y , z The total number of directional background grids.

[0024] As a further solution of the present invention: when the number of the point pair ( a' , b' , c' )and( d' , e' , f' ) satisfies the following relationship:

[0025] ,

[0026] at this time: ,in, K E (( a , b , c ), ( u ,v , w )) indicates that the number in the matrix is ​​( a , b , c ) and the node numbered ( u , v , w )'s relationship value; s is the Gauss-Legendre integration order, and ceil() means rounding up.

[0027] As a further solution of the present invention: the preprocessing matrix includes M K and M KE Two parts,

[0028] ;

[0029] .

[0030] As a further solution of the present invention: the pseudo code is as follows:

[0031] 1) ;

[0032] 2) , ;

[0033] 3) , , ;

[0034] 4) ;

[0035] 5) ;

[0036] 6) ;

[0037] 7) ;

[0038] 8) ;

[0039] 9) ;

[0040] 10) , , ;

[0041] 11) , ;

[0042] 12) ;

[0043] 13) , where the subscript k Indicates the k iterations, D is the kernel matrix, x is the gravity magnetic field, b The right-hand term.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] This paper improves upon the meshless method based on a spatial equivalence approach, achieving a highly efficient forward modeling of gravity and magnetic data. This approach provides an effective method for efficiently processing gravity and magnetic data over large areas and in undulating terrain. It can reduce memory usage by tens of times and increase computational efficiency several times, with the improvement becoming more pronounced the finer the meshing. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 Schematic diagram of node layout.

[0047] Figure 2 A schematic diagram of the background grid.

[0048] Figure 3 It is a cubic spline weight function.

[0049] Figure 4 Schematic diagram of the node and background grid layout of this method.

[0050] Figure 5 This is the first node distribution diagram of this method.

[0051] Figure 6 This is a schematic diagram of the second node distribution of this method.

[0052] Figure 7 This is the model location diagram.

[0053] Figure 8 is the forward modeling result of gravity and magnetic field.

[0054] Figure 9 This is the technical process of the present invention. DETAILED DESCRIPTION

[0055] The following will provide a clear and complete description of the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of them. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0056] The element-free Galerkin method is based on the calculus of variations to solve differential equations. For any variational problem:

[0057] (1),

[0058] in, I , F , f , g is a universal function; r is the independent variable; V is the integration space; Γ is the boundary of the integration space; α is a constant.

[0059] See Figure 1 In the element-free Galerkin method, by V A series of scattered points are arranged in the problem domain, and the function value of each node is calculated, which is the solution of the differential equation in the entire problem domain.

[0060] See Figure 2 In order to realize the integration of the system equation in the problem domain, it is necessary to construct a set of background grids covering the entire problem domain, and convert the continuous integral in the entire problem domain into the sum of the integrals in each background grid. The integration of the system equation in a background grid is usually implemented using a certain numerical integration method, and the integral node value is obtained by interpolation of the calculation nodes within a certain range. Therefore, the final integral form only contains the information of the calculation nodes, and the background grid is only used to assist the integration. The interpolation method adopted in the present invention is the Moving least square method (MLS), the integration method is the Gauss-Legendre integration method, and the essential boundary condition adopts the Lagrange multiplier method.

[0061] MLS interpolation is based on the least squares method and uses the nodes in a certain area to construct an interpolation function, and then calculates the value of the point to be interpolated. u ( r ) in the solution domain Ω in n nodes r i ( i = 1, 2, ... , n )( n is a positive integer) the function value is u i = u ( r i ), then the function value u ( r 0) At the calculation point r Neighborhood of 0 Ω0 can be approximately expressed as:

[0062] (2),

[0063] in, r = [ r 1, r 2, ... , r n ] is the calculation point r Neighborhood of 0 Ω Coordinates of each point within 0; p ( r )=[ p 1( r ), p 2( r ), ... , p m ( r )] is the basis function; m is the number of basis function terms; a ( r 0) is the function coefficient; u p ( r 0, r ) is the interpolation function at the calculation point r The calculated value is 0.

[0064] Basis functions need m Subcomplete ( m is a positive integer), and the following conditions must be met:

[0065] (3),

[0066] in, C k ( Ω ) indicates that there are up to k The function space of the continuous derivative of the order. The present invention adopts the following form:

[0067] (4).

[0068] The function coefficients are calculated by constructing the L2 norm objective function and taking the minimum value. First, the following objective function is constructed:

[0069] (5),

[0070] in, ; ; is a weight function, which makes the farther nodes no longer affect the interpolation node, and the farther the point is from the interpolation point, the smaller the effect it has on it.

[0071] When , we can get:

[0072] (6).

[0073] but:

[0074] (7).

[0075] The gradient at this point is:

[0076] (8),

[0077] in: β , β 1, β 2∈{ x , y , z}; ; ; From the calculation process of the MLS interpolation method, we can see that the smoothness and continuity of the interpolation method depend on the smoothness of the basis function and the weight function. The basis function is obviously smooth, and the weight function needs to meet the following conditions:

[0078] 1. The weight function value is greater than 0 in the solution domain and equal to 0 outside the solution domain;

[0079] 2. From the calculation point r It decreases monotonically from 0 to the boundary of the solution domain;

[0080] 3. The weight function is a smooth function.

[0081] There are many weight functions that meet the above conditions, such as Gaussian weight function, exponential function, cubic spline function, and quartic spline function. The present invention adopts cubic spline function as the weight function.

[0082] (9),

[0083] in, ; d 0 is the radius of the solution domain. Figure 3 is a function form, Figure 3 The function form in can obviously meet the above requirements for the weight function; ; .

[0084] The Gauss quadrature formula is an interpolation quadrature formula constructed using the Lagrange interpolation method.

[0085] For a function f ( x )exist[ a ,b ]Integral in the interval , for a given weight function ρ ( x ), , and its Gauss quadrature formula is:

[0086] (10),

[0087] Among them, the weight coefficient W k It can be obtained by the basis function of the interpolation function:

[0088] (11).

[0089] For the Gauss-Legendre quadrature formula, the interpolation nodes and weight functions in the Gauss quadrature formula can be obtained from the Legendre polynomial. The non-normalized Legendre polynomial can be obtained by the following recursive formula:

[0090] (12).

[0091] Legendre polynomials P n ( x ) = 0 root { x pk}∈[-1, 1], k = 1, 2, ... , n . Then the interpolation node coordinates in the Gauss-Legendre quadrature formula { x k} can be expressed as:

[0092] (13).

[0093] Weight coefficient { W k}for:

[0094] (14).

[0095] The final integral formula is:

[0096] (15).

[0097] The Gauss-Legendre quadrature formula for triple integrals is:

[0098] (16).

[0099] Since the approximate function of MLS does not have the properties of Kronecker delta function, that is, N j ( r i ) ≠ δ ij Therefore, it does not have the interpolation property, that is, Therefore, the essential boundary conditions cannot be naturally satisfied in the meshless method using this method, and additional functions need to be added to impose the essential boundary conditions. They are mainly divided into boundary configuration method, Lagrange multiplier method, modified variational principle, penalty function method and finite element coupling method. This paper adopts the Lagrange multiplier method to impose boundary conditions and introduces it into the variational problem described by formula (1):

[0100] (17),

[0101] in, is the Lagrange multiplier vector; for The values ​​at the boundary nodes of the background grid can be approximated by corresponding methods for different data; W E is the weight coefficient matrix calculated based on the Gauss-Legendre integration method at the boundary; N The interpolation function constructed for the MLS interpolation method.

[0102] Then its weak variational form is:

[0103] (18),

[0104] in, K , K E The matrix obtained by discretizing the integral in the variational problem through certain interpolation functions and integration methods has different forms for different differential equations.

[0105] Taking into account and The arbitrariness of , we can get:

[0106] (19).

[0107] According to the calculation principle of the element-free Galerkin method, it is known that during the calculation process of this method, an interpolation function needs to be constructed in the neighborhood of the integration node in each discrete unit, which will consume a lot of computing time. Therefore, the present invention proposes a high-efficiency algorithm based on spatial equivalence.

[0108] 1. Equivalent calculation

[0109] See Figure 4 , Figure 4 The layout of nodes and background grids in the present invention follows the following rules: 1) Discrete nodes are divided into two types: regularly distributed nodes and irregularly distributed nodes, and irregular nodes are not set in the outermost background grid of the expansion area and the study area;

[0110] 2) Regular nodes are evenly distributed within the study area;

[0111] 3) The node spacing in the expanded area increases as it approaches the boundary, and becomes consistent with the node spacing in the study area when it approaches the study area;

[0112] 4) The background grid is a regular cube, and the corner points are regularly distributed nodes;

[0113] 5) Use interpolation equations and integral equations of the same order globally;

[0114] 6) The support domain of an integration node (the neighborhood for constructing the interpolation function of the node) is the background grid where it is located and its adjacent background grids, that is, each time the interpolation function is constructed, the adjacent 64 nodes are used.

[0115] Obviously, when applying the element-free Galerkin method, it is no longer necessary to calculate the interpolation function and integral function in each background grid to obtain the corresponding local kernel matrix. Instead, all local kernel matrices can be obtained by calculating only a small number of local kernel matrices through spatial equivalence relations. Figure 4 It can be found that the local kernel matrix calculated for the background grid that does not contain irregular nodes in the support domain has the following characteristics:

[0116] 1) For a background mesh with a symmetric relationship (about its center or about its central axis), its support domain shape and node distribution within the support domain also have a symmetric relationship, so the local kernel matrix constructed by it also has a symmetric relationship. Since the elements in the kernel matrix constructed based on the element-free Galerkin method are the relationship values ​​between any two points, the above content can be expressed as the following functional relationship:

[0117] (20).

[0118] Among them, the subscript ( i , j , k ) indicates that the matrix is ​​the equation in number ( i , j , k )’s background grid is integrated to obtain the local kernel matrix; (( a , b , c ),( d , e , f )) indicates that the value is the number in the background grid support domain (a , b , c ) and the node numbered ( d , e , f )'s node's relationship value in the matrix; i = 1, 2, ... , N x ; j =1, 2, ... , N y ; k = 1, 2, ... , N z ; N x = n x + 2 n xe ; N y = n y + 2 n ye ; N z = n z + 2 n ze ; a , b , c, d , e , f ∈ {1, 2, 3, 4}. And:

[0119] when i 1 = i hour, a 1 = a , d 1 = d ;when i 1 = N x －i +1, a 1 = 5- a , d 1 = 5- d ;

[0120] when j 1 = j hour, b 1 = b , e 1 = e ;when j 1 = Ny －j +1, b 1 = 5- b , e 1 = 5- e ;

[0121] when k 1 = k hour, c 1 = c , f 1 = f ;when k 1 = N z －k +1, c 1 = 5- c , f 1 = 5- f .

[0122] because K E It is only relevant to the boundary elements, so:

[0123] (twenty one).

[0124] Among them, the subscript β ( i , j , k ) indicates that the matrix is ​​the equation in number ( i , j , k ) in the background grid β The local kernel matrix obtained by the surface integral operation in the direction; β ∈ { x , y , z};(( a , b , c ),( u , v , w )) indicates that the value is the number in the background grid support domain ( a , b , c ) and the node numbered ( u , v , w )'s Gauss-Legendre integration node in the matrix; u , v , w = 1, 2, ... , s ; sis the order of Gauss-Legendre integration. And:

[0125] when i 1 = i hour, a 1 = a , u 1 = u ;when i 1 = N x －i +1, a 1 = 5- a , u 1 = s - u +1;

[0126] when j 1 = j hour, b 1 = b , v 1 = v ;when j 1 = N y －j +1, b 1 = 5- b , v 1 = s - v +1;

[0127] when k 1 = k hour, c 1 = c , w 1 = w ;when k 1 = N z －k +1, c 1 = 5- c , w 1 = s - w +1.

[0128] 2) In areas where the support domain of the background grid and the relative position relationship of adjacent nodes are consistent, the local kernel matrices constructed are completely consistent:

[0129] (twenty two).

[0130] Therefore, the background grid to be processed in this step is composed of N x × Ny × N z Reduced to ( n xe + 1) × ( n ye + 1) × ( n ze + 1) + n e indivual, n e is the number of background grids with irregular nodes in the support domain.

[0131] Similarly, for K E ,have

[0132] (twenty three).

[0133] Obviously, the above equivalent operations significantly reduce the computational complexity of the kernel matrix, but directly combining them into a global kernel matrix still cannot effectively reduce memory usage and improve inversion computational efficiency. Therefore, when combining the individual local kernel matrices into a global kernel matrix, the present invention proposes a sparse storage method based on the equivalent relationship of the local kernel matrices.

[0134] The elements in the global kernel matrix constructed based on the element-free Galerkin method are the comprehensive relationship values ​​between any two nodes, which are obtained by summing the relationship values ​​of these two nodes in all local kernel matrices. Due to the symmetry of the local kernel matrices, there are a large number of identical values ​​in the global kernel matrix, based on which they can be sparsely stored.

[0135] For two nodes that do not coexist with irregular nodes in any support domain, their numbers are recorded as ( a , b , c )and( d , e , f ) (Note that here a ~ f is used for global numbering of nodes, and is used for local numbering within the support domain in equations (21) and (22), and satisfies:

[0136] ,

[0137] Its relationship value in the overall kernel matrix is:

[0138] (24), among which, K (( a ,b , c ), ( d , e , f )) indicates that the number in the matrix is ​​( a , b , c ) and the node numbered ( d , e , f )'s node relationship value; 1≤ a , d ≤ N x , 1≤ b , e ≤ N x , 1≤ c , f ≤ N x ; N x = n x +2 n xe , N y = n y +2 n ye , N z = n z +2 n ze They are x , y , z Total number of directional background grids, n x , n y , n z The study areas x , y , z Total number of directional background grids; n xe , n ye , n ze The expansion area x , y , z The total number of directional background grids.

[0139] Combining equations (20), (22) and (24), the elements in the matrix have the following equality relationship:

[0140] 1) Symmetrical equality between point pairs: The relationship values ​​of point pairs symmetrical about the central axis or center are equal:

[0141] (25).

[0142] 2) Translation equality relationship: For pairs of points with equal distances and the same node distribution within the support domain containing the pair, the relationship values ​​between the pairs are equal:

[0143] (26).

[0144] 3) Symmetric equality relationship between points:

[0145] (27).

[0146] See Figure 5 , based on the above symmetric relationship, calculate Figure 5 The relationship values ​​of the point pairs in the two blue boxes in the kernel matrix can be obtained by calculating the relationship values ​​of the two blue boxes (if the two boxes overlap, the point pairs in the overlapping area need to calculate the relationship values ​​including the influence of irregular nodes, and also need to calculate the relationship values ​​excluding the influence of irregular nodes, so that their values ​​can be used for their equivalent point pairs).

[0147] The relationship value between any pair of points that are not related to the irregular nodes can be obtained according to the following relationship:

[0148] First, the number of the point pairs with known relationship values ​​( a' , b' , c' )and( d' , e' , f' )satisfy:

[0149] ,

[0150] Then we have:

[0151] (28).

[0152] See Figure 6 , for the matrix K E It needs to be calculated Figure 6 The relationship value of the point pair in the blue box (here it is not two discrete nodes, but a discrete node and an integral node) can be obtained according to the following relationship:

[0153] First, discrete nodes with known relationship values ​​( a' , b' , c') and Gauss-Legendre integration nodes ( u' , v' , w' ) number satisfies:

[0154] ,

[0155] Then we have:

[0156] (29).

[0157] 2. Solve the calculation

[0158] The forward modeling process based on the element-free Galerkin method is the process of solving the linear equation system (19). The previous article introduced the sparse storage method for the kernel matrix. For this large sparse linear equation system, this paper adopts the conjugate gradient method and certain preprocessing methods to achieve efficient solution calculation.

[0159] The preprocessing method is a method that makes the properties of the matrix after preprocessing better, that is, the distribution of eigenvalues ​​is more concentrated, thereby improving the convergence speed of the iterative solution method. Its general form is:

[0160] (30).

[0161] in, M 1, M 2 are the left and right preprocessing matrices respectively; D is the kernel matrix; x is the quantity to be solved; b The right-hand term.

[0162] In the present invention, the core matrix D It is a symmetric matrix and has been stored sparsely. If the generation efficiency of the preprocessing matrix is ​​too slow or the symmetry and sparsity of the kernel matrix are destroyed, it is obviously not worth the loss. Therefore, this paper uses the following form to construct the preprocessing matrix:

[0163] (31).

[0164] Among them, diag() means extracting the diagonal elements of the matrix. Therefore, the two preprocessing matrices are the same and are diagonal matrices. Obviously,

[0165] (32).

[0166] That is, the matrix after preprocessing is still a symmetric matrix. At the same time, since the preprocessing matrix M 1 , M 2 by the kernel matrix A Generated by the sparse form of the kernel matrix, its calculation formula can be expressed as MK and M KE Two parts:

[0167] (33).

[0168] (34).

[0169] Therefore it has similar characteristics to the kernel matrix:

[0170] (35).

[0171] (36).

[0172] And the elements in the final preprocessed kernel matrix can be expressed as:

[0173] (37).

[0174] By comparing equations (28), (29), (35) and (36), it can be found that the symmetric form of the preprocessed matrix is ​​exactly the same as the symmetric form of the original kernel matrix. Therefore, the sparse storage form of the preprocessed kernel matrix is ​​exactly the same as the sparse storage form of the original kernel matrix, and will not be repeated here.

[0175] The pseudo code of the conjugate gradient algorithm combined with this preprocessing method is as follows:

[0176] 1) ;

[0177] 2) , ;

[0178] 3) , , ;

[0179] 4) ;

[0180] 5) ;

[0181] 6) ;

[0182] 7) ;

[0183] 8) ;

[0184] 9) ;

[0185] 10) , , ;

[0186] 11) , ;

[0187] 12) ;

[0188] 13) .

[0189] 3. Gravity and magnetic forward modeling formula

[0190] The variational problem corresponding to gravity forward modeling is:

[0191] (38).

[0192] in, g is the gravitational field; ρ is the density; G ≈ 6.67 × 10 -11 N·m 2 kg -2 is the universal gravitational constant. Based on this, the local kernel matrix and right-hand side term of the gravity forward model can be constructed:

[0193] (39).

[0194] in, For the number ( i , j , k )'s background grid's length, width, and height; For the number ( i , j , k )’s background grid interpolation function and its derivatives; For the number ( i , j , k )'s background grid supports the density of nodes within the domain.

[0195] The variational problem corresponding to magnetic forward modeling is:

[0196] (40).

[0197] in, T is the magnetic field; M is the magnetization intensity; μ 0 ≈ 4π×10 -7 T·m / A is the magnetic permeability in vacuum; cos( φ Mx ),cos( φ My ), cos( φ Mz ) is the magnetization intensity M Direction cosines; cos( φ Tx ), cos( φ Ty ), cos( φ Tz ) are the direction cosines of the Earth's magnetic field.

[0198] Based on this, the local kernel matrix and right-hand side term of the magnetic forward model can be constructed:

[0199] (41).

[0200] in, d x,ijk 、 d y,ijk 、 d z,ijk For the number ( i , j , k )'s background grid's length, width, and height; N ijk 、 N β,ijk 、 N β1β2,ijk For the number ( i , j , k )’s background grid interpolation function and its derivatives; M ijk For the number ( i , j , k )’s background mesh supports the magnetization of nodes within the domain.

[0201] The specific implementation of the present invention is described in detail below with reference to specific embodiments.

[0202] See Figure 7-Figure 9 , calculated using the method of the present invention Figure 7 The gravity and magnetic fields of the model are shown. At the same time, the calculation time of gravity forward modeling under different segmentation conditions is statistically analyzed and given (since the gravity and magnetic forward modeling is similar, the improvement ratio using this method is consistent). The model parameters are shown in Table 1.

[0203] Table 1 Model parameters

[0204] Serial number 1 2 Location (km) (1.0,1.0,-0.6) (0.5,0.5,0.5) Length, width and height (km) (2.1, 1.0, -0.8) (0.7, 0.7, 0.7) <![CDATA[Density (g / cm 3 )]]> 0.5 0.5 Magnetization intensity (A / m) 0.5 0.5 Inclination and declination (60°,10°) (60°,10°)

[0205] The forward modeling effect of gravity and magnetic field is as follows Figure 8 As shown in Table 2, the calculation time and memory usage under different segmentation conditions are shown in Table 2.

[0206] Table 2 Computation time and memory usage

[0207] Number of nodes / number of mesh elements 20×20×10 40×40×20 80×80×40 Number of observation points 20×20 40×40 80×80 This method calculates the time 20.7s 76.0s 848s The memory usage of this method 31.5Mb 69.7Mb 611Mb Integration method calculation time 5.2s 165.9s 5308.4s Integration method memory usage 12.2Mb 309.6Mb 12500Mb Computation time ratio 0.25 2.18 6.25 Memory usage ratio 0.39 4.44 20.45

[0208] The data in Table 2 show that this method can reduce memory usage by several times and improve computational efficiency, and the finer the subdivision, the more significant the improvement in efficiency.

[0209] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A gridless gravity and magnetic forward modeling method based on spatial equivalence, characterized in that: The following steps are involved: The first step is to set the observation point and model position; The second step is to set the positions of discrete nodes and background grid; The third step is to set the discrete nodes in the model as related physical properties; The fourth step is to calculate the local kernel matrix and the right-hand side term according to the forward modeling formula and the equivalent principle. The local kernel matrix of gravity and magnetic forward modeling is: , for gravity forward modeling, the right-hand side term is in the form of: ， in, G is the gravitational constant; For the number ( i , j , k )'s background grid's length, width, and height; 、 For the number ( i , j , k )'s background grid interpolation function and β Directional derivative, Yes The general representation of β ∈{ x , y , z }; Yes The general representation of β ∈{ x , y , z }; For the number ( i , j , k )’s background grid supports the density of nodes within the domain; W and W E is the Gauss-Legendre integral weight coefficient, W is the 8th-order identity matrix, W E is the 4th-order identity matrix, ; , , They are x , y , z Total number of directional background grids; n x , n y , n z The study areas x , y , z Total number of directional background grids; n xe , n ye , n ze The expansion area x , y , z The total number of directional background grids. For magnetic forward modeling, the right-hand side term is: ， in, is the magnetic permeability in vacuum; For the number ( i , j , k )’s background grid supports the magnetization of nodes within the domain, Yes , , , , The general expression of For the number ( i , j , k )’s second-order derivative of the interpolation function of the background grid; β 1, β 2∈{ x , y , z }; cos( φ Mβ1 ) is the value of cos( φ Mx ), cos( φ My ), cos( φ Mz ) is a general representation of cos( φ Mβ1 ) is the magnetization intensity in x , y , z Direction cosines of direction; cos( φ Tβ2 ) is the value of cos( φ Tx ), cos( φ Ty ), cos( φ Tz ) is a general representation of cos( φ Tβ2 ) is the Earth's magnetic field x , y , z Direction cosines of direction; The fifth step is to calculate the preprocessing matrix and weight it. The calculation formula is: , , M 1, M 2 are the left and right preprocessing matrices respectively; D is the kernel matrix; x is the quantity to be solved; b is the right-hand term, diag() means extracting the diagonal elements of the matrix; The sixth step is to realize the forward modeling of gravity and magnetic field based on pseudo code.

2. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 1 is characterized in that: The discrete nodes include regularly distributed nodes and irregularly distributed nodes, and irregular nodes are not set in the expanded area of ​​the background grid and the outermost layer of the research area, and the regular nodes in the research area are evenly arranged.

3. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 2 is characterized in that: The background grid is a regular cube, and the corner points of the background grid are regularly distributed nodes.

4. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 1 is characterized in that: The relationship between the local kernel matrix in the area where the support domain shape of the background grid and the relative position relationship of the adjacent nodes are consistent is: ,in, K (( a , b , c ), ( d , e , f )) indicates that the number in the matrix is ​​( a , b , c ) and the node numbered ( d , e , f )'s node relationship value; 1≤ a , d ≤ N x , 1≤ b , e ≤ N x , 1≤ c , f ≤ N x ; N x = n x +2 n xe , N y = n y +2 n ye , N z = n z +2 n ze They are x , y , z Total number of directional background grids, n x , n y , n z The study areas x , y , z Total number of directional background grids; n xe , n ye , n ze The expansion area x , y , z The total number of directional background grids.

5. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 4 is characterized in that: When the point pair number ( a' , b' , c' )and( d' , e' , f' ) satisfies the following relationship: ,at this time: ,in, K E (( a , b , c ), ( u , v , w )) indicates that the number in the matrix is ​​( a , b , c ) and the node numbered ( u , v , w )'s relationship value; s is the Gauss-Legendre integration order, and ceil() means rounding up.

6. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 1 is characterized in that: The preprocessing matrix includes M K and M KE Two parts, ； 。 7. The gridless gravity and magnetic forward modeling method based on spatial equivalence according to claim 1 is characterized in that: For any system of equations to be solved Dx = b , the pseudo code is as follows: ; ， ； ， ， ； ； ； ； ； ； ； ， ， ； ， ； ； , where the subscript k Indicates the k iterations, D is the kernel matrix, x is the gravity magnetic field, b The right-hand term.

Images