An inversion method for arbitrary variable density interfaces

Abstract

The present invention discloses an inversion method for arbitrary variable density interfaces, and specifically relates to the field of gravity exploration technology. Based on the method of three-dimensional prism partitioning, the present invention constructs a three-dimensional density interface model in a three-dimensional space grid, and by performing Taylor series expansion on the density interface forward function and retaining the first-order term, the linear function of the density interface model forward is approximately obtained and brought into the regularization objective function, and the density interface inversion iterative equation is obtained by minimizing the objective function, and the density interface inversion model is solved using a linear conjugate gradient algorithm. By iteratively correcting the density interface multiple times, the density interface and density model of the study area are inverted. The present invention realizes the accurate acquisition of density interface distribution under complex geological structures, and based on the fast algorithm of three-dimensional gravity inversion, realizes batch rapid inversion of gravity-density interfaces, laying a foundation for the exploration and development of complex geological structures.

Description

Technical Field

[0001] The present invention belongs to the technical field of gravity exploration, and in particular relates to an inversion method for arbitrary variable density interfaces. Background Art

[0002] The gravity-density interface corresponds to geological structures such as basement undulation, spatial distribution of layers, tectonic uplift, and Moho surface, and has important theoretical and practical significance in studying regional structures, delineating oil and gas prospective areas, and crust-mantle structure.

[0003] Frequency-domain inversion methods emerged in the 1970s and rapidly developed and applied due to their computational speed. Spatial-domain density interface inversion methods emerged earlier, dating back to the 1950s. Subsequently, they have been widely studied and used, resulting in the development of a variety of inversion methods. Typical methods include empirical formulas, direct iteration methods, ridge regression methods, regularization methods, compressed mass surface methods, series methods, and spline function methods. Direct iteration methods, ridge regression methods, and regularization methods are the most widely studied and used methods, while research on other methods is relatively limited.

[0004] Whether in frequency domain forward modeling or spatial domain forward modeling, the current density interface modeling method basically uses a prism partitioning model to construct a density interface model. This is done by partitioning the spatial domain horizontally and then using the buried depth of the top or bottom surface of the prism to represent the depth of the density interface in the vertical direction. When performing variable density forward modeling, this partitioning method requires that the density change be summarized as a function, and then the forward modeling formula is derived based on the function. This results in the known density change having to be expressed by a finite function, for example, a linear density function that only changes with depth, a quadratic polynomial density function, a cubic polynomial density function, a parabolic density function, a hyperbolic density function, an exponential density function, or a two-dimensional density change function based on a polynomial or a three-dimensional density change function.

[0005] The above function is generally smooth and can well fit the density distribution of simple geological structures. However, due to the influence of tectonic action, large changes in stratum morphology, missing strata, and other sudden changes in density, actual geological structures are often complex, and the above function cannot accurately express the density distribution of geological structures. Therefore, for most geological structures, the existing density variation law cannot accurately express the actual density distribution characteristics of geological structures. It is urgent to propose an inversion method for arbitrary variable density interfaces to accurately obtain the density distribution of complex geological structures. Summary of the Invention

[0006] In order to solve the problem that the existing density variation function is difficult to accurately describe the density distribution characteristics of complex geological structures, the present invention proposes an inversion method for arbitrary variable density interfaces. The present invention realizes the accurate acquisition of density interface distribution in complex geological structures, and performs rapid calculation based on three-dimensional gravity inversion, which is conducive to the rapid inversion of gravity-density interfaces in batches.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] An inversion method for arbitrary variable density interfaces comprises the following steps:

[0009] Step 1: Obtain density attribute information of the study area, construct a three-dimensional density interface model in a three-dimensional space grid, and use three-dimensional prism forward modeling to determine the gravity field of the three-dimensional density interface model of the study area;

[0010] Step 2: Based on the forward gravity field of the three-dimensional density interface model of the study area, determine the objective function of the density interface regularization inversion and establish the density interface inversion calculation model;

[0011] Step 3: Perform density interface inversion based on the density interface inversion calculation model to obtain the density interface and density model of the study area.

[0012] Preferably, in step 1, a three-dimensional prism model is used to represent the density interface of the study area, and the density attribute information of the area above and below the density interface is determined based on the density attribute information of the study area.

[0013] Preferably, in step 1, the gravity field is determined based on the three-dimensional prism model gravity forward modeling as:

[0014] g=g(h) (1)

[0015] Where g is the surface gravity value and h is the density interface depth.

[0016] Preferably, in step 2, the objective function of the density interface regularization inversion is determined according to the gravity field of the study area, as shown in formula (2):

[0017]

[0018] Where, is the objective function of the density interface regularization inversion, W d is the gravity data weighting matrix, g obs is the gravity observation value, g(h) is the gravity field function, λ is the regularization parameter, W h is the density interface weighting matrix, h pre is the density interface reference data;

[0019] Since the gravity field function g(h) is a nonlinear function, the Taylor series expansion of the gravity field function g(h) is performed and the first-order term is retained to obtain:

[0020]

[0021] Among them, the Jacobian matrix A is used for discrete three-dimensional grid gravity inversion i Medium element A jk The calculation formula is:

[0022]

[0023] Where i is the number of gravity inversion calculations, h is i is the density interface depth corresponding to the i-th gravity inversion calculation, A i is the density interface depth h i The corresponding Jacobian matrix, g j is the jth observation point, h k is the kth interface depth point, G jk is the interface depth point h k The three-dimensional prism model is located at the observation point g j Forward coefficient, m k is the density value of the three-dimensional prism model, Δh k is the longitudinal length of the three-dimensional prism model;

[0024] The gravity field function g(h) after Taylor series expansion is brought into the objective function of the density interface regularization inversion, and the density interface depth h is derived to obtain:

[0025]

[0026] Where T is the transposed matrix;

[0027] The value of the objective function of the density interface regularization inversion is set to zero, and the density interface inversion iterative formula is obtained as follows:

[0028]

[0029] Add formula (6) twice get:

[0030]

[0031] By inverting the matrix on the left side of formula (7), the density interface inversion calculation model is obtained, as shown in formula (8):

[0032]

[0033] Preferably, in step 3, a linear conjugate gradient algorithm is used to solve the density interface inversion calculation model, and the density interface and density model of the study area are obtained by inversion through multiple iterative corrections of the density interface, which specifically includes the following steps:

[0034] Step 3.1, inversion parameter initialization;

[0035] Input gravity observation data g obs , select the density interface reference data h pre , set the initial value h0 of the density interface depth and the density interface weighting matrix W h , gravity data weighting matrix W d and the maximum number of corrections Nmax, the preset correction accuracy and the density distribution between known interfaces in the density interface inversion calculation model;

[0036] Step 3.2: Use discrete 3D grids for gravity inversion, and calculate the density interface depth h i and the density interface depth h i The corresponding density model calculates the density interface depth h i The corresponding Jacobian matrix A i ;

[0037] Step 3.3, solve the density interface inversion calculation model based on the linear conjugate gradient method to obtain the corrected density interface depth h i+1 , determine the density interface depth h i+1 The corresponding density interface;

[0038] Step 3.4, calculate the correction value Δh of the density interface depth. If the correction value Δh is less than the preset correction accuracy or the number of corrections has reached the preset maximum number of corrections, proceed to step 3.5. Otherwise, update the number of corrections i and return to step 3.2 for correction.

[0039] Step 3.5, output density interface and density model.

[0040] Preferably, in step 3.3, solving the density interface inversion calculation model based on the linear conjugate gradient method to obtain the density interface corresponding to the corrected density interface depth specifically includes the following steps:

[0041] Step 3.3.1, iterative calculation parameter initialization;

[0042] Preset the gradient accuracy value ε of the density interface inversion calculation model, set the linear conjugate gradient calculation parameters Q and b, Set the initial gradient r0 of the density interface inversion calculation model, the initial search direction p0 of the iterative calculation, and the initial value m0 of the iterative update, r0 = b-Qm0, p0 = r0, m0 = 0.001;

[0043] Step 3.3.2: Calculate the step size of the iterative search direction based on the gradient of the density interface inversion calculation model and the iterative search direction, as shown in formula (9):

[0044]

[0045] Where, α k is the step size of the iterative search direction during the k-th iteration calculation, r k is the gradient value of the density interface inversion calculation model during the kth iteration calculation, p k The search direction calculated for the kth iteration;

[0046] Step 3.3.3, update the iterative calculation results and the gradient of the density interface inversion calculation model:

[0047] m k+1 =m k +α k p k (10)

[0048] r k+1 =r k -α k Qp k (11)

[0049] Where m k+1 is the updated iterative calculation result, m k is the result of the kth iteration, r k+1 is the gradient value of the density interface inversion calculation model after the update;

[0050] Step 3.3.3, if the updated density interface inversion calculation model gradient value is the second norm || r k+1 ||2 is less than the preset gradient accuracy value ε, then the result m calculated according to the kth iteration k , determine the corrected density interface depth h i+1 , as shown in formula (12):

[0051] h i+1 =m k +h pre (12)

[0052] Calculate the corrected density interface depth h i+1 Then, proceed to step 3.3.5;

[0053] If the updated density interface inversion calculation model gradient value of the second norm || r k+1 If ||2 is not less than the preset gradient accuracy value ε, proceed to step 3.3.4;

[0054] Step 3.3.4: Update the search direction and number of iterations based on the updated density interface inversion calculation model gradient value, and then return to step 3.3.2 to continue the iterative calculation.

[0055] The search direction of the iterative calculation after the update is:

[0056]

[0057] Where p k+1 Compute the search direction for the updated iteration;

[0058] Step 3.3.5, end the iterative calculation and output the corrected density interface depth h i+1 and the density interface depth h i+1 The corresponding density interface.

[0059] The beneficial technical effects brought about by the present invention are:

[0060] The present invention proposes an inversion method for arbitrary variable density interfaces. Unlike traditional density interface inversion methods, the method of the present invention must summarize the interface density changes into a function and derive a forward formula, and then perform density interface inversion based on the forward derivation. The method of the present invention is based on a three-dimensional space grid, and divides the three-dimensional density interface model into multiple three-dimensional prisms within the three-dimensional space grid. By assigning each three-dimensional prism a separate density attribute, the density change of the interface within the three-dimensional space grid is regarded as a spatial distribution based on the function, and the spatial distribution of the interpolation can be constrained based on the known density attribute information, thereby realizing the inversion of arbitrary variable density interfaces. The method is suitable for situations where density interfaces are interlaced, which is conducive to accurately obtaining the density interface distribution under complex geological structures, and provides a basis for the exploration and development of complex geological structures.

[0061] The present invention is based on a three-dimensional gravity rapid inversion method, which improves the computational efficiency of three-dimensional density interface inversion, has a fast inversion speed and low computational cost, is conducive to batch rapid inversion of gravity-density interfaces, and is suitable for promotion and application in the process of complex geological structure interpretation. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 This is a schematic diagram of the density interface represented by the three-dimensional prism model according to the present invention.

[0063] Figure 2 Schematic diagram of the Jacobian matrix in the present invention.

[0064] Figure 3 Schematic diagram of each layer interface in the three-dimensional density interface model of this embodiment.

[0065] Figure 4 Schematic diagram of each density model in the three-dimensional density interface model of this embodiment.

[0066] Figure 5 is the density distribution at the density interface in this embodiment. Figure 5 (a) is the density distribution at interface 1, Figure 5 (b) shows the density distribution at interface 2.

[0067] Figure 6 This is a schematic diagram of the gravity field determined by gravity forward modeling in this embodiment.

[0068] Figure 7 This is the iterative inversion result of this embodiment. Figure 7 (a) is the gravity forward field after the first correction of the three-dimensional density interface model. Figure 7 (b) is the gravity forward field after the first correction of the three-dimensional density interface model. Figure 7 (c) is the density interface after the first correction. Figure 7 (d) is the difference between the density interface and the theoretical interface after the first correction. Figure 7 (e) is the gravity forward field after the third correction of the three-dimensional density interface model. Figure 7 (f) is the gravity forward field after the third correction of the three-dimensional density interface model. Figure 7 (g) is the density interface after the third correction. Figure 7 (h) is the difference between the density interface and the theoretical interface after the third correction. Figure 7 (i) is the gravity forward field after the 10th correction of the three-dimensional density interface model. Figure 7 (j) is the gravity forward field after the 10th correction of the three-dimensional density interface model. Figure 7 (k) is the density interface after the 10th correction. Figure 7 (l) is the difference between the density interface and the theoretical interface after the 10th correction. Figure 7 (m) is the gravity forward field after the 17th correction of the three-dimensional density interface model. Figure 7 (n) is the gravity forward field after the 17th correction of the three-dimensional density interface model. Figure 7 The middle (o) is the density interface after the 17th correction. Figure 7 (p) is the difference between the density interface and the theoretical interface after the 17th correction. DETAILED DESCRIPTION

[0069] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0070] Example 1

[0071] The present invention proposes an inversion method for arbitrary variable density interfaces, which specifically includes the following steps:

[0072] Step 1: Based on the seismic data and well logging data of the study area, the density attribute information of the study area is obtained, and a three-dimensional density interface model is constructed in the three-dimensional space grid. The density attributes of the area above and below the density interface are set according to the density attribute information of the study area, and the gravity field of the three-dimensional density interface model of the study area is determined by using three-dimensional prism forward modeling.

[0073] Since the density attribute in the three-dimensional prism model is known information and the depth of the density interface is unknown, the gravity forward modeling is a function of the density interface. The three-dimensional prism forward modeling is used to determine the gravity field of the study area, such as Figure 1 As shown, the gravity field is determined as shown in formula (1):

[0074] g=g(h) (1)

[0075] Where g is the surface gravity value and h is the density interface depth.

[0076] Step 2: Based on the gravity field of the three-dimensional density interface model in the study area, determine the objective function of the density interface regularization inversion and establish a density interface inversion calculation model.

[0077] According to the gravity field of the study area, the objective function of the density interface regularization inversion is determined as shown in formula (2):

[0078]

[0079] Where, is the objective function of the density interface regularization inversion, W d is the gravity data weighting matrix, g obs is the gravity observation value, g(h) is the gravity field function, λ is the regularization parameter, W h is the density interface weighting matrix, h pre It is the reference data of density interface.

[0080] Since the gravity field function g(h) is a nonlinear function and cannot be directly expanded into a linear function, the gravity field function g(h) is expanded using Taylor series and the first-order term is retained to obtain an approximate linear function:

[0081]

[0082] Among them, the Jacobian matrix A is used for discrete three-dimensional grid gravity inversion i Medium element A jk The calculation formula is:

[0083]

[0084] Where i is the number of gravity inversion calculations, h is iis the density interface depth corresponding to the i-th gravity inversion calculation, A i is the density interface depth h i The corresponding Jacobian matrix, g j is the jth observation point, h k is the kth interface depth point, G jk is the interface depth point h k The three-dimensional prism model is located at the observation point g j The forward coefficients of Figure 2 As shown, the forward coefficient G jk is the forward modeled value of unit density, m k is the density value of the three-dimensional prism model, Δh k is the longitudinal length of the three-dimensional prism model.

[0085] The gravity field function g(h) after Taylor series expansion is brought into the objective function of the density interface regularization inversion, and the density interface depth h is derived to obtain:

[0086]

[0087] Where T is the transposed matrix.

[0088] The value of the objective function of the density interface regularization inversion is set to zero, and the density interface inversion iterative formula is obtained as follows:

[0089]

[0090] Add formula (6) twice After finishing, we can get:

[0091]

[0092] By inverting the matrix on the left side of formula (7), the density interface inversion calculation model is obtained, as shown in formula (8):

[0093]

[0094] Step 3: Perform density interface inversion based on the density interface inversion calculation model. Use the linear conjugate gradient algorithm to solve the density interface inversion calculation model. By performing multiple iterative corrections on the density interface, the density interface and density model of the study area are obtained by inversion. Specifically, the following steps are included:

[0095] Step 3.1, inversion parameter initialization;

[0096] Input gravity observation data g obs , select the density interface reference data h pre , set the initial value h0 of the density interface depth and the density interface weighting matrix Wh , gravity data weighting matrix W d and the maximum number of corrections Nmax, the preset correction accuracy and the density distribution between known interfaces in the density interface inversion calculation model.

[0097] Step 3.2: Use discrete 3D grids for gravity inversion, and calculate the density interface depth h i and the density interface depth h i The corresponding density model calculates the density interface depth h i The corresponding Jacobian matrix A i .

[0098] Step 3.3, solve the density interface inversion calculation model based on the linear conjugate gradient method to obtain the corrected density interface depth h i+1 , determine the density interface depth h i+1 The corresponding density interface.

[0099] Step 3.3 specifically includes the following steps:

[0100] Step 3.3.1, iterative calculation parameter initialization;

[0101] Preset the gradient accuracy value ε of the density interface inversion calculation model, set the linear conjugate gradient calculation parameters Q and b, Set the initial gradient r0 of the density interface inversion calculation model, the initial search direction p0 of the iterative calculation, and the initial value m0 of the iterative update, r0 = b-Qm0, p0 = r0, m0 = 0.001.

[0102] Step 3.3.2: Calculate the step size of the iterative search direction based on the gradient of the density interface inversion calculation model and the iterative search direction, as shown in formula (9):

[0103]

[0104] Where, α k is the step size of the iterative search direction during the k-th iteration calculation, r k is the gradient value of the density interface inversion calculation model during the kth iteration calculation, p k Search direction calculated for the kth iteration.

[0105] Step 3.3.3, update the iterative calculation results and the gradient of the density interface inversion calculation model:

[0106] m k+1 =m k +α k p k (10)

[0107] r k+1 =rk -α k Qp k (11)

[0108] Where m k+1 is the updated iterative calculation result, m k is the result of the kth iteration, r k+1 The gradient value of the inversion calculation model for the updated density interface.

[0109] Step 3.3.3, if the updated density interface inversion calculation model gradient value is the second norm || r k+1 ||2 is less than the preset gradient accuracy value ε, then the result m calculated according to the kth iteration k , determine the corrected density interface depth h i+1 , as shown in formula (12):

[0110] h i+1 =m k +h pre (12)

[0111] Calculate the corrected density interface depth h i+1 Then proceed to step 3.3.5.

[0112] If the updated density interface inversion calculation model gradient value of the second norm || r k+1 If ||2 is not less than the preset gradient accuracy value ε, proceed to step 3.3.4.

[0113] In step 3.3.4, based on the updated density interface inversion calculation model gradient value, the search direction and number of iterative calculations are updated. The iterative calculation search direction is updated as shown in formula (13):

[0114]

[0115] Where p k+1 Compute the search direction for the updated iteration.

[0116] After updating the search direction and the number of iterative calculations k, return to step 3.3.2 to continue the iterative calculation.

[0117] Step 3.3.5, end the iterative calculation and output the corrected density interface depth h i+1 and the density interface depth h i+1 The corresponding density interface.

[0118] Step 3.4, calculate the correction amount Δh of the density interface depth. When the correction amount Δh of the density interface depth is less than the preset correction accuracy or the number of corrections has reached the preset maximum number of corrections, enter step 3.5; otherwise, update the number of corrections i and continue to return to step 3.2 for correction.

[0119] Step 3.5, output density interface and density model.

[0120] Example 2

[0121] In order to verify the feasibility of the inversion method for arbitrary variable density interfaces proposed in the present invention and its ability to handle complex models, this embodiment adopts the inversion method for arbitrary variable density interfaces described in Example 1 to perform density interface inversion, which specifically includes the following steps:

[0122] Step 1: Obtain the density attribute information of the study area and establish a three-dimensional density interface model controlled by a double interface, such as Figure 3 and Figure 4 As shown, the density interface in the 3D density interface model is interlaced. The 3D density interface model is divided into 91×81×50 3D prism models with dimensions of 100m×100m×50m. In this embodiment, the origin coordinates of the 3D prism models in the 3D space grid are (50, 50, 0), and the vertical coordinate is set to be negative downward.

[0123] In this embodiment, the three-dimensional density interface model is composed of three layers of density models. The spatial range of the first layer of density model is controlled by the surface interface and interface 1. The surface interface is the plane where sur0_h(x, y) = 0, and the interface depth sur1_h(x, y) of interface 1 is:

[0124] sur1_h(x,y)=-(18+(sin(x / 2000×π)+sin(y / 1000×π))×3)×40 (14)

[0125] The spatial distribution of the density ρ1(x,y,z) of the first layer model is:

[0126] ρ1(x,y,z)=-(x+y-8700) / 85000+exp(z / 200) / 1500-0.65 (15)

[0127] Where (x, y, z) is the spatial position of the point in the three-dimensional density interface model.

[0128] The spatial range of the second-layer density model is controlled by interface 1 and interface 2. The interface depth sur2_h(x,y) of interface 2 is:

[0129]

[0130] The spatial distribution of the density ρ1(x,y,z) of the second layer model is:

[0131] ρ2(x,y,z)=(x+y-8700) / 85000+exp(z / 333) / 4000-0.45 (17)

[0132] The spatial range of the third-layer density model is controlled by interface 2 and the bottom interface of the model space. The density of the third-layer model is ρ3(x, y, z)=0.

[0133] Extract the density distribution above the density interface from the model, such as Figure 5 As shown, from Figure 5 As can be seen from the figure, the density distribution above both interfaces fluctuates, with an overall trend of low in the southwest and high in the northeast. Because Interface 2 penetrates Interface 1, some blank areas appear on Interface 1. The density in the area where Interface 2 penetrates Interface 1 experiences a sudden change, resulting in a complex overall density change.

[0134] The gravity field of the study area is determined by forward modeling of a three-dimensional prism. The forward modeled gravity field of the study area is as follows: Figure 6 As shown, the gravity measurement point is located on the surface interface and coincides with the horizontal center point of the three-dimensional prism. In order to better conform to actual gravity measurement and reduce the influence of model boundary effects, in this embodiment, each boundary of the gravity forward field is reduced to 10 measurement points before being used for inversion calculation.

[0135] Step 2: Based on the forward gravity field of the study area, determine the objective function of the density interface regularization inversion and establish the density interface inversion calculation model as follows:

[0136]

[0137] Step 3: Perform density interface inversion based on the density interface inversion calculation model, use the linear conjugate gradient algorithm to solve the density interface inversion calculation model, and obtain the density interface and density model of the study area by performing multiple iterative corrections on the density interface.

[0138] In the density interface inversion process of this embodiment, interface 2 in the three-dimensional density interface model is set as the target interface. The initial value of the density interface depth of interface 2 is set to h0 = -1383, that is, the density interface depth of interface 2 is set to the average value of the theoretical interface depth. For the theoretical model experiment, since the density distribution, fixed interface and other factors are known and accurate, the inversion result of the target interface is unique. Therefore, in this embodiment, there is no need to set a reference density interface. The density interface weighting matrix W h Set as depth weighting matrix, set gravity data weighting matrix W d =1, maximum number of corrections Nmax = 17.

[0139] According to the inversion process, after 17 interface corrections, the inversion result of interface 2 in the three-dimensional density interface model is obtained, as shown in Figure 7 As shown in the figure, the iterative inversion process shows that as the number of corrections increases, the interface inversion error and the fitted gravity field error of the corrected model gradually decrease, and the inversion results tend to stabilize. The final density inversion interface error is generally small, with larger errors occurring only at the boundaries of the inversion area. In some areas, the maximum error is 168 meters. In the center of the inversion area, the inversion error is generally less than 50 meters, which is less than the longitudinal length of a single three-dimensional prism.

[0140] The above experiments verify that the density interface inversion calculation model constructed by the present invention is capable of inverting any variable density interface. It also further proves that when the known information is rich and accurate, the three-dimensional density interface inversion method of the present invention can achieve excellent inversion effect.

[0141] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. An inversion method for arbitrary variable density interfaces, characterized in that: The specific steps include: Step 1: Obtain density attribute information of the study area, construct a three-dimensional density interface model in a three-dimensional space grid, and use three-dimensional prism forward modeling to determine the gravity field of the three-dimensional density interface model of the study area; Step 2: Based on the forward gravity field of the three-dimensional density interface model of the study area, determine the objective function of the density interface regularization inversion and establish the density interface inversion calculation model; Step 3: Perform density interface inversion based on the density interface inversion calculation model to obtain the density interface and density model of the study area; In step 1, a three-dimensional prism model is used to represent the density interface of the study area, and density attribute information of the area above and below the density interface is determined based on the density attribute information of the study area; In step 1, the gravity field is determined based on the three-dimensional prism model gravity forward modeling as: g=g(h) (1) Where g is the surface gravity value, h is the density interface depth; In step 2, the objective function of the density interface regularization inversion is determined according to the gravity field of the study area, as shown in formula (2): Where, is the objective function of the density interface regularization inversion, W d is the gravity data weighting matrix, g obs is the gravity observation value, g(h) is the gravity field function, λ is the regularization parameter, W h is the density interface weighting matrix, h pre is the density interface reference data; Since the gravity field function g(h) is a nonlinear function, the Taylor series expansion of the gravity field function g(h) is performed and the first-order term is retained to obtain: Among them, the Jacobian matrix A is used for discrete three-dimensional grid gravity inversion i Medium element A jk The calculation formula is: Where i is the number of gravity inversion calculations, h is i is the density interface depth corresponding to the i-th gravity inversion calculation, A i is the density interface depth h i The corresponding Jacobian matrix, g j is the jth observation point, h k is the kth interface depth point, G jk is the interface depth point h k The three-dimensional prism model is located at the observation point g j Forward coefficient, m k is the density value of the three-dimensional prism model, Δh k is the longitudinal length of the three-dimensional prism model; The gravity field function g(h) after Taylor series expansion is brought into the objective function of the density interface regularization inversion, and the density interface depth h is derived to obtain: Where T is the transposed matrix; The value of the objective function of the density interface regularization inversion is set to zero, and the density interface inversion iterative formula is obtained as follows: Add formula (6) twice get: By inverting the matrix on the left side of formula (7), the density interface inversion calculation model is obtained, as shown in formula (8):

2. The inversion method for arbitrary variable density interface according to claim 1, characterized in that: In step 3, the linear conjugate gradient algorithm is used to solve the density interface inversion calculation model. By performing multiple iterative corrections on the density interface, the density interface and density model of the study area are obtained by inversion, which specifically includes the following steps: Step 3.1, inversion parameter initialization; Input gravity observation data g obs , select the density interface reference data h pre , set the initial value h0 of the density interface depth and the density interface weighting matrix W h , gravity data weighting matrix W d and the maximum number of corrections Nmax, the preset correction accuracy and the density distribution between known interfaces in the density interface inversion calculation model; Step 3.2: Use discrete 3D grids for gravity inversion, and calculate the density interface depth h i and the density interface depth h i The corresponding density model calculates the density interface depth h i The corresponding Jacobian matrix A i ; Step 3.3, solve the density interface inversion calculation model based on the linear conjugate gradient method to obtain the corrected density interface depth h i+1 , determine the density interface depth h i+1 The corresponding density interface; Step 3.4, calculate the correction value Δh of the density interface depth. If the correction value Δh is less than the preset correction accuracy or the number of corrections has reached the preset maximum number of corrections, proceed to step 3.

5. Otherwise, update the number of corrections i and return to step 3.2 for correction. Step 3.5, output density interface and density model.

3. The inversion method for arbitrary variable density interface according to claim 2, characterized in that: In step 3.3, the density interface inversion calculation model is solved based on the linear conjugate gradient method to obtain the density interface corresponding to the corrected density interface depth, which specifically includes the following steps: Step 3.3.1, iterative calculation parameter initialization; Preset the gradient accuracy value ε of the density interface inversion calculation model, set the linear conjugate gradient calculation parameters Q and b, Set the initial gradient r0 of the density interface inversion calculation model, the initial search direction p0 of the iterative calculation, and the initial value m0 of the iterative update, r0 = b-Qm0, p0 = r0, m0 = 0.001; Step 3.3.2: Calculate the step size of the iterative search direction based on the gradient of the density interface inversion calculation model and the iterative search direction, as shown in formula (9): Where, α k is the step size of the iterative search direction during the k-th iteration calculation, r k is the gradient value of the density interface inversion calculation model during the kth iteration calculation, p k The search direction calculated for the kth iteration; Step 3.3.3, update the iterative calculation results and the gradient of the density interface inversion calculation model: m k+1 =m k +α k p k (10) r k+1 =r k -a k Qp k (11) Where m k+1 is the updated iterative calculation result, m k is the result of the kth iteration, r k+1 is the gradient value of the density interface inversion calculation model after the update; Step 3.3.3, if the updated density interface inversion calculation model gradient value is the second norm || r k+1 ||2 is less than the preset gradient accuracy value ε, then the result m calculated according to the kth iteration k , determine the corrected density interface depth h i+1 , as shown in formula (12): h i+1 =m k +h pre (12) Calculate the corrected density interface depth h i+1 Then, proceed to step 3.3.5; If the updated density interface inversion calculation model gradient value of the second norm || r k+1 If ||2 is not less than the preset gradient accuracy value ε, proceed to step 3.3.4; Step 3.3.4: Update the search direction and number of iterations based on the updated density interface inversion calculation model gradient value, and then return to step 3.3.2 to continue the iterative calculation. The search direction of the iterative calculation after the update is: Where p k+1 Compute the search direction for the updated iteration; Step 3.3.5, end the iterative calculation and output the corrected density interface depth h i+1 and the density interface depth h i+1 The corresponding density interface.

Images