Method for 3D joint inversion of gravity and magnetic data based on cross-gradient constraint

Abstract

The application belongs to the technical field of geophysical application, and relates to a method of gravity-magnetic three-dimensional joint inversion based on cross gradient constraint. A subsurface half space is divided into a plurality of regular cuboids in a spatial array mode, and initial physical property parameters in each cuboid are randomly set, wherein the initial physical property parameters are density parameters and magnetization strength parameters; an observation point is arranged vertically above each cuboid on the ground; three-dimensional gravity-magnetic anomaly forward calculation is carried out according to the divided cuboids and the density parameters and the magnetization strength parameters, so that predicted ground observation data are obtained; errors are calculated according to actually measured data of the observation points and the predicted ground observation data; finally, joint iterative inversion is carried out by using a cross gradient function, so that the errors are reduced to a reasonable range, and an ideal model is obtained. Compared with a traditional inversion, the application reduces the multi-solution of gravity-magnetic inversion and improves the resolution of gravity-magnetic inversion to a certain extent, reduces corresponding calculation, and improves the inversion efficiency.

Description

Technical Field

[0001] This invention belongs to the field of geophysics application technology, and in particular relates to a method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints. Background Technology

[0002] Geophysical responses arise from differences in the physical properties of subsurface media. Geophysical inversion is a discipline that infers the distribution of subsurface media properties, such as density, magnetic susceptibility, or resistivity, based on observed geophysical responses. This process inevitably faces two fundamental problems in geophysical exploration: the inherent ambiguity of inversion itself and the limitations of single-method inversion. Only by combining geological theoretical knowledge with a comprehensive approach to geophysical development can these two fundamental problems be effectively resolved. Therefore, based on single-method geophysical inversion, combined geophysical inversion has gradually developed.

[0003] Within the same survey area of ​​the same work zone, we can collect data from multiple geophysical methods, such as gravity anomalies, magnetic anomalies, and resistivity information. Although the data from these different geophysical methods may seem unrelated, the subsurface medium corresponding to these measured data is consistent. Therefore, we can consider comprehensively utilizing these collected geophysical responses to infer subsurface geological-geophysical models, including physical property parameters, model burial depth, shape, and size. Different physical property parameters satisfy certain inherent correspondences, which can, to some extent, improve the non-uniqueness of the inversion and enhance the accuracy of inversion interpretation. Therefore, joint inversion is an important direction for future geophysical development and is increasingly attracting the attention of geophysicists.

[0004] As surface and near-surface resources are continuously depleted, the depth of geophysical exploration is also increasing, and the problems faced are becoming more and more complex. Geophysical inversion using a single method is no longer sufficient to address the current practical problems, making the research and application of joint geophysical inversion essential and urgent.

[0005] Currently, the most common inversion methods are sequential inversion and simultaneous joint inversion. Both sequential and simultaneous joint inversion require empirical formulas between physical properties to link different types of data; however, inaccurate empirical relationships can lead to misinterpretations in joint inversion. Considering the complex and variable nature of Earth's tectonic evolution, rock properties exhibit different characteristics in different tectonic regions, making it inappropriate to use the same empirical relationships across a broad study area. Therefore, how to construct coupling relationships between parameters of different geophysical models and how to improve the accuracy and efficiency of joint inversion are currently key aspects of joint inversion research. Summary of the Invention

[0006] To address the technical problems existing in the above-mentioned joint inversion methods, this invention proposes a reasonable, simple, widely applicable, and high-resolution three-dimensional joint inversion method for gravity and magnetic fields based on cross-gradient constraints.

[0007] To achieve the above objectives, the technical solution adopted by this invention is as follows: This invention provides a method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints, comprising the following steps:

[0008] a. Divide the underground semi-space into several regular cuboids in a spatial array, and randomly set the initial physical property parameters in each cuboid, wherein the initial physical property parameters are density parameters and magnetization intensity parameters.

[0009] b. Establish observation points on the ground directly above each cuboid;

[0010] c. Based on the divided cuboids, perform three-dimensional gravity and magnetic anomaly forward modeling calculations using density and magnetization parameters to obtain predicted ground observation data;

[0011] d. Perform inversion calculations based on the actual measured data from the observation points and the predicted ground observation data;

[0012] e. Finally, the optimal solution is obtained by joint inversion using the cross gradient function, where the cross gradient function is:

[0013]

[0014] Where, m (1) m (2) These are the density model parameters and the magnetization model parameters, respectively, k (1) For density model parameters m (1) The difference between the maximum and minimum values, k (2) For the magnetization model parameter m (2) The difference between the maximum and minimum values.

[0015] Preferably, in step d, the inversion iterative process is divided into an outer loop and an inner loop. The outer loop performs forward modeling to fit the data and selects an appropriate regularization factor based on the magnitude of the initial values ​​of each term in the objective function. The inner loop modifies the model to obtain an ideal model corresponding to the regularization factor, calculates the error, and iterates multiple times by reducing the regularization factor until the number of iterations reaches the maximum or the error reaches the ideal value to achieve structural similarity.

[0016] Preferably, in step d, a relatively large damping coefficient α is set at the beginning of the outer loop iteration. An appropriate damping coefficient α can be selected based on the magnitude of the initial values ​​of each term in the objective function. Enter the sub-loop.

[0017] Preferably, the inner loop is calculated using the outer loop. A new model is calculated and iterated multiple times until the number of iterations reaches its maximum, the model changes tend to stabilize, or the value of the cross gradient function increases. During the iteration process, the model and the cross gradient function converge, resulting in a structurally similar model.

[0018] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0019] 1. This invention provides a method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints. For the first time, it performs joint inversion of three-dimensional gravity and magnetic fields under the condition of cross-gradient constraints on density parameters and magnetization parameters. Compared with traditional inversion, it reduces the ambiguity of gravity and magnetic inversion to a certain extent, improves the resolution of gravity and magnetic inversion, reduces the corresponding calculations, and improves the efficiency of inversion. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 A plan view of the Bouguer gravity (ΔGb) contour lines in the survey area;

[0022] Figure 2 This is a plan view of the contour lines of the △T area in the survey region;

[0023] Figure 3 The image shows the results of the joint inversion of residual density (0 isosurface).

[0024] Figure 4 The image shows the results of the joint inversion of magnetization (0.6 A / m isosurface).

[0025] Figure 5 The image shows the results of the joint inversion of residual density (0 isosurface).

[0026] Figure 6 The image shows the results of the joint inversion of magnetization (0.6 A / m isosurface).

[0027] Figure 7(a) shows the residual density profile of the two-dimensional gravity and magnetic joint inversion results;

[0028] Figure 7(b) shows the magnetization profile of the two-dimensional gravity and magnetic joint inversion results;

[0029] Figure 8(a) shows the residual density results obtained from the combined gravity and magnetic inversion of the Y = 10000m profile;

[0030] Figure 8(b) shows the magnetization results obtained by combined gravity and magnetic inversion of the Y = 10000m profile;

[0031] Figure 9(a) shows the residual density results obtained from the combined gravity and magnetic inversion of the Y = 15000m profile;

[0032] Figure 9(b) shows the results of the combined gravity and magnetic inversion magnetization inversion of the Y = 15000m profile;

[0033] Figure 10(a) shows the residual density results obtained from the combined gravity and magnetic inversion of the Y = 20000m profile;

[0034] Figure 10(b) shows the magnetization results obtained by combined gravity and magnetic inversion of the Y = 20000m profile;

[0035] Figure 11(a) shows the residual density results obtained from the combined gravity and magnetic inversion of the Y = 25000m profile;

[0036] Figure 11(b) shows the magnetization results obtained by combined gravity and magnetic inversion of the Y = 25000m profile;

[0037] Figure 12(a) shows the residual density results obtained from the combined gravity and magnetic inversion of the Y = 30000m profile;

[0038] Figure 12(b) shows the magnetization results obtained by combined gravity and magnetic inversion of the Y = 30000m profile;

[0039] Figure 13 A separate inversion plot of the cross gradient values ​​for Y = 20000m;

[0040] Figure 14 The plot shows the cross gradient values ​​obtained by joint inversion for Y = 20000m. Detailed Implementation

[0041] To better understand the above-mentioned objectives, features, and advantages of the present invention, the present invention will be further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0042] Numerous specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways than those described herein, and therefore the invention is not limited to the specific embodiments disclosed in the following specification.

[0043] Example 1: This example provides a method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints.

[0044] First, the underground semi-space is divided into several regular cuboids in a spatial array, and initial physical property parameters are randomly set for each cuboid, including density and magnetization parameters. The calculation process of gravity anomaly is summarized into the following two main steps: (1) Calculate the gravitational potential generated at the observation point by the remaining mass of a specific geological body relative to the surrounding rock according to Newton's law of universal gravitation; (2) Calculate the derivative of the above gravitational potential along the vertical direction to obtain the gravity anomaly generated by the geological body at the corresponding observation point. Therefore, in this embodiment, an observation point is set on the ground directly above each cuboid.

[0045] Specifically, assuming the surface observation plane is divided into P×Q grids, and the underground space is divided into L layers, with each layer having an M×N grid, then the number of integrations required for one forward modeling operation would be P×Q×M×N×L, which undoubtedly poses a significant computational challenge. Therefore, we assume here that there are M×N observation points on the ground, each located directly above a rectangular cell arranged in a regular grid. Based on the relative relationship between the observation points and the model cells, the computational load can be greatly reduced.

[0046] Suppose there are M observation points on the observation profile, and each layer of the underground model corresponds to M rectangular elements. For the first observation point, there are M possible relative distances (r1, r2, ..., rM) with respect to a certain layer of the underground model mesh. For the second observation point, two (a pair) of its relative distances to a certain layer of the underground model mesh are the same, i.e., r1 = r3, so there are a total of M-1 possible relative distances (r1, r2, ..., rM). Similarly, the relative distances between the Mth observation point and the first observation point are exactly the same, and their relative distances are r1, r2, ..., rM. Based on the above analysis, for the M models in the first layer of this profile, each observation point needs to have its distance calculated M times, and the distances for the M observation points need to be calculated M2 times. However, the distances for the subsequent M-1 observation points are simply different combinations of the distances calculated for the first observation point. In fact, we only need to calculate the M possible relative distances. A similar pattern exists for the 3D model; we only need to add one dimension to the above analysis.

[0047] Based on the division of the cuboid, the three-dimensional gravity and magnetic anomaly forward modeling is performed using density and magnetization parameters to obtain the predicted ground observation data.

[0048] Let the volume of the geological body be v, and the density difference be σ. Let (ξ,η,ζ) be the coordinates of any point within the geological body. The volume element at that point can then be expressed as dv = dξdηdζ, and its remaining mass is written as dm = σdv = σdξddζ. The gravitational potential of the geological body at the observation point A(x,y,z) is expressed as:

[0049]

[0050] Where G is the gravitational constant, r represents the distance from any point (ξ,η,ζ) in the geological body to the observation point, and the gravity anomaly can be written as:

[0051]

[0052] The underground semi-space is divided into upright cuboids, with the physical properties within each cube set to constants. A discretized integral solution for the integral formula is given, in the following form:

[0053]

[0054] The subsurface medium is divided into M cuboids, and N observation points are set up on the observation surface. The gravity anomaly generated by the j-th cuboid at the i-th observation point can be expressed as:

[0055] Δg=G ij σ j

[0056] In the formula, Gij is a known quantity determined based on the relative position of the j-th cuboid at the i-th observation point, and σj is the residual density of the j-th cuboid.

[0057] The relationship between magnetic potential and gravitational potential can be expressed using Poisson's formula:

[0058]

[0059] In the formula, U is the magnetic potential, G is the gravitational constant, M is the magnetization, σ is the density difference, and V is the gravitational potential. We can obtain:

[0060]

[0061] In the formula, I is the geomagnetic inclination angle, A' is the azimuth angle of the survey line, and additionally...

[0062] M x =M·cosI·sinA; M y =M·cosI·cosA; M=M·sinI

[0063] Substituting Za, Hax, and Hay into the expression for geomagnetic anomaly ΔT, we can calculate:

[0064] ΔT=H ax cosIcosA'+H ay cosIsinA'+Z a sinI

[0065] Then, inversion calculations are performed based on the actual measured data from the observation points and the predicted ground observation data.

[0066] The cross gradient function has the following properties

[0067] (1) When the two physical property parameters (density parameter and magnetization parameter) of the joint inversion change in parallel directions, or when one of the physical property parameters remains unchanged, the cross gradient function of the two is equal to zero.

[0068] (2) When the two physical property parameters (density parameter and magnetization parameter) in the joint inversion change in different directions, the cross gradient function of the two is not equal to zero.

[0069] These properties indicate that joint inversion based on cross-gradient coupling improves the reliability of inversion results by enhancing the consistency of the structure reflected by different physical property parameters.

[0070] In the two-dimensional case, the cross gradient vector contains non-zero components only in the directional direction (y-axis direction), and is approximated by forward difference in numerical calculations:

[0071]

[0072] In three-dimensional space:

[0073]

[0074] After discretizing the physical property model, the cross gradient also has a corresponding discretized form. When using a cuboid mesh (rectangular in two-dimensional cases), it can be discretized using finite difference:

[0075]

[0076] In the formula, i, j, and k are the element numbers of the 3D model along the coordinate directions, and x, y, and z represent the differences in physical properties of the elements in the three coordinate directions. The cross gradient vector can be written as:

[0077]

[0078] The density, magnetization, and velocity parameters, gravity and magnetic seismic observation data, cross gradient vectors, and the first-order derivative matrix of the cross gradient function are combined into a unified vector, considering the objective function:

[0079]

[0080] In the formula:

[0081]

[0082] constraint:

[0083] τ=τ(m (1),m (2) ) = 0

[0084] Where m (1) m (2) For different geophysical models, d (1) d (2) Corresponding geophysical data As a priori model, C m C d Here, are the model covariance matrix and the data covariance matrix, respectively; g is the forward modeling operator; and α is the weighting factor (damping coefficient) of the regularization term. σ represents the data observation error. To facilitate the optimization calculation of the objective function, the L2 norm is used for all objective functions.

[0085] To solve the objective function of the joint inversion, it needs to be linearized. This can be done using the Taylor series expansion method to linearize the forward function and the cross-gradient function. The solution for the constrained function can then be obtained using the Lagrange multiplier method.

[0086] W=Φ(m)+2Λ T [t(m0)+B(m-m0)]

[0087] In the formula, Λ is the Lagrange operator and B is the derivative of the cross gradient function.

[0088] Taking the partial derivatives of the model parameters and finding them to be zero yields the model update parameters and the Lagrange operator:

[0089] Δm (i) =N (i)-1 (n (i) -B (i)T Λ)

[0090] in:

[0091]

[0092]

[0093] The corresponding Lagrange multipliers:

[0094]

[0095] Model update expression:

[0096]

[0097] B is The first derivative of the model is extended as follows:

[0098]

[0099] Update the regularized least squares model without cross-gradient constraints:

[0100]

[0101] When the forward operator is linear

[0102] Introducing depth weighting and the property range constraint of the transformation function method, the depth weighting adopts the approximation function proposed by Li and Oldenburg, and the joint depth weighting matrix is ​​a combination of the depth weighting diagonal matrices of the single gravity and magnetic method:

[0103]

[0104] The property parameter transformation function adopts the inverse transformation of the formula proposed by Commer, and the transformed joint model parameters are as follows:

[0105]

[0106] Its partial derivative matrix is:

[0107]

[0108] Introducing the structural coupling factor β, the joint inversion iterative equation becomes:

[0109]

[0110] In the formula, the expressions for matrix N and vector n are respectively:

[0111]

[0112] In the formula, g is the transformed joint forward operator. Vector t satisfies the following system of linear equations:

[0113]

[0114] The above formula introduces a structural coupling factor β. Using this formula for calculation avoids time-consuming singular value decomposition calculations, greatly improving solution efficiency. It is worth noting that the choice of the coupling factor is independent and unrelated to the choice of other parameters. When approaching infinity, the structural coupling perturbation term is approximately zero, and the inversion approaches a single inversion; while when the coupling factor is sufficiently large, it is equivalent to no damping, maximizing the preservation of structural perturbation information.

[0115] The inversion iterative process is divided into two loops: the outer loop, also known as the main loop, which makes the model fit the data, and the inner loop, also known as the sub-loop, which modifies the model to achieve structural similarity.

[0116] (1) Main loop: At the beginning of the iteration, a relatively large damping coefficient α is set. An appropriate damping coefficient α can be selected based on the magnitude of the initial values ​​of each term in the objective function. Calculation Enter the sub-loop, obtain the ideal model corresponding to the damping coefficient from the sub-loop, calculate misfit, and iterate multiple times by decreasing the damping coefficient α until the number of iterations reaches the maximum or misfit reaches the ideal value.

[0117]

[0118] Sub-loop: Calculated using the values ​​from the main loop The new model is calculated and iterated multiple times until the number of iterations reaches its maximum, the model changes tend to stabilize, or the value of the cross gradient function increases. During the iteration process, the model and the cross gradient function converge (the constraint formula of the cross gradient function is approximated by linearity in the calculation), resulting in a structurally similar model.

[0119] Relative changes between adjacent models:

[0120]

[0121] Finally, when the empirical relationships between physical properties are unknown or unclear, it is very suitable to introduce cross-gradient constraints for joint inversion. As long as the model parameters are coupled, the cross-gradient function can be obtained:

[0122]

[0123] To prevent numerical instability, k is added. Generally, k is the expected variation amplitude of m. Before joint inversion, separate inversions are performed. The difference between the maximum and minimum values ​​of each model parameter in the separate inversion results is the corresponding k value. For joint gravity and magnetic inversion, m... (1) m (2) These are the density model parameters and the magnetization model parameters, respectively.

[0124] Inversion of measured data

[0125] Figure 1 The figure shows a plane map of Bouguer gravity (ΔGb) contour lines in the survey area. The spacing between survey lines and the interval between survey points are both 500m. As shown in the figure, the data reflects the Bouguer gravity anomaly in a certain area of ​​the Jiaodong Peninsula. The gravity gradient zone in the survey area is mainly NE-trending, and there are obvious density differences in the lithology on both sides, which is consistent with the trend of the Jiaojia fault zone in the survey area. Figure 2The figure shows a planar map of the ΔT contour lines in the survey area, after polarization. The survey line spacing is 100m, and the measurement point spacing is 500m. As shown in the figure, there is a significant positive anomaly in the NE-NEE direction, similar to the morphology of the gravity gradient zone. During the joint inversion using the cross-gradient method, the calculation of the cross-gradient derivative Bk is involved. Bk is a matrix of size (3×Nm)×Nm, where Nm is the number of model parameters. Therefore, the calculation of tk requires a large amount of memory. Thus, a coarser grid partitioning is adopted, and the model parameter size is smaller. For this survey area, the grid partitioning is 35×30×20. The region with ordinates from 5000 to 35000m and abscissas from 18000 to 43000m is selected and subjected to separate gravity and magnetic inversions and joint inversions under the same grid partitioning. Then, three-dimensional gravity and magnetic inversions were performed individually and jointly, producing five inversion profiles with ordinates Y=1000m, Y=15000m, Y=20000m, Y=25000m, and Y=30000m. The joint and individual inversions were compared, and the cross gradient value of the inversion results for one of the profiles Y=20000m was calculated.

[0126] Through the Figures 3-1 2. Analysis revealed that the individual inversion results were chaotic, with irregular anomaly distribution and significant deviations in anomaly amplitude from actual values, making it difficult to extract the necessary information. In the joint inversion map, the residual density results showed a positive residual density anomaly in the northwest of the survey area and a negative residual density anomaly in the southeast. The residual density 0 isosurface roughly trended NNE-NE, and its variation pattern largely matched the strike of the Jiaojia fault zone in this region. The magnetization intensity 0.6 A / m isosurface roughly trended NE, exhibiting a positive magnetization intensity anomaly below a depth of 1000 m.

[0127] Compare Figure 13 and Figure 14 The joint inversion results show better structural similarity than the individual inversion results, ensuring a sufficiently small fitting difference while also satisfying the constraints of the cross gradient values. Theoretically, the resulting model satisfies the data constraints, and the good structural similarity between models reduces the ambiguity of the inversion.

[0128] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments for application in other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for 3D joint inversion of gravity and magnetic data based on cross-gradient constraint, characterized in that, Includes the following steps: a. Divide the underground semi-space into several regular cuboids in a spatial array, and randomly set the initial physical property parameters in each cuboid, wherein the initial physical property parameters are density parameters and magnetization intensity parameters; b. Establish observation points on the ground directly above each cuboid; c. Based on the divided cuboids, perform three-dimensional forward modeling calculations of gravity and magnetic anomalies using density and magnetization parameters to obtain predicted ground observation data; d. Perform inversion calculations based on the actual measured data from the observation points and the predicted ground observation data; e. Finally, the optimal solution is obtained by joint inversion using the cross gradient function, where the cross gradient function is: wherein , are density model parameters and magnetization model parameters, respectively, is a density model parameter is the difference between the maximum and minimum values of is the difference between the maximum and minimum values of the magnetization model parameter is the difference between the maximum and minimum values of the magnetization model parameter In step d, the inversion iterative process is divided into an outer loop and an inner loop. The outer loop performs forward modeling to fit the data and selects an appropriate regularization factor based on the magnitude of the initial values ​​of each term in the objective function. The inner loop modifies the model to obtain an ideal model corresponding to the regularization factor, calculates the error, and iterates multiple times by reducing the regularization factor until the number of iterations reaches the maximum or the error reaches the ideal value to achieve structural similarity.

2. The method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints according to claim 1, characterized in that, In step d, the outer loop iteration starts with a relatively large damping coefficient α. An appropriate damping coefficient α is selected based on the initial values ​​of each term in the objective function, and the calculation is performed. Enter the sub-loop, where, This is for updating the regularized least squares model without cross-gradient constraints.

3. The method for three-dimensional joint inversion of gravity and magnetic fields based on cross-gradient constraints according to claim 2, characterized in that, The inner loop is calculated using the outer loop. The new model is calculated and iterated multiple times until the number of iterations reaches its maximum, the model changes tend to stabilize, or the value of the cross gradient function increases. During the iteration process, the model and the cross gradient function converge, resulting in a structurally similar model.

Images