CN115201927A - A method for combined gravity and magnetic 3D inversion based on cross gradient constraints - Google Patents

Abstract

The invention belongs to the technical field of geophysics applications, and relates to a method for combined gravity and magnetic three-dimensional inversion based on cross gradient constraints. The underground half-space is divided into several regular cuboids in the form of a spatial array, and the initial physical parameters in each cuboid are randomly set, wherein the initial physical parameters are density parameters and magnetization parameters; Set up observation points on the ground; carry out 3D gravity and magnetic anomaly forward calculation with density parameters and magnetization parameters according to the divided cuboid, and obtain the predicted ground observation data; calculate the error according to the actual measured data of the observation point and the predicted ground observation data; Finally, a joint iterative inversion is performed using the cross gradient function to reduce the error to a reasonable range, and an ideal model is obtained. Compared with the traditional inversion, the invention reduces the multi-solution of the gravity and magnetic inversion to a certain extent, improves the resolution of the gravity and magnetic inversion, reduces the corresponding calculation, and improves the efficiency of the inversion.

Description

A method for combined gravity and magnetic 3D inversion based on cross gradient constraints

technical field

The invention belongs to the technical field of geophysics applications, and in particular relates to a method for combined gravity and magnetic three-dimensional inversion based on cross gradient constraints.

Background technique

Geophysical responses are produced by the differences in the physical properties of the subsurface medium. Geophysical inversion is a discipline that infers the distribution of physical properties of subsurface media based on the observed geophysical responses, such as the distribution of physical properties such as density, magnetic susceptibility or resistivity. . This process will inevitably face two basic problems of geophysical exploration, namely the inherent multiple solutions of inversion and the one-sidedness of single-method inversion. And only by combining the guidance of geological theoretical knowledge and taking the road of comprehensive geophysics development, can these two basic problems be well solved. Therefore, on the basis of the single method inversion of geophysics, the joint inversion of geophysics has been gradually developed.

Within the same survey area in the same work area, we can collect data from a variety of geophysical methods, such as gravity anomalies, magnetic anomalies, and resistivity information. Although there seems to be no correlation between the data corresponding to these different geophysical methods, the subsurface medium corresponding to these measured data is unified. Therefore, we can consider comprehensively using these collected geophysical responses to infer the subsurface geology- Geophysical model, including physical parameter values, burial depth, shape and scale of the model. Different physical parameter values satisfy a certain internal correspondence, which can improve the non-uniqueness of inversion and the accuracy of inversion interpretation to a certain extent. Therefore, joint inversion is an important direction of future geophysics development, and it has attracted more and more attention of geophysicists.

With the continuous exploitation of surface and near-surface resources, the depth of geophysical exploration is also increasing, and the problems faced are becoming more and more complex. A single method of geophysical inversion is no longer capable of meeting the current reality. It is very necessary and urgent to study and apply joint geophysical inversion.

At present, the most common inversion methods are sequential inversion and simultaneous joint inversion. Both sequential inversion and simultaneous joint inversion need to use empirical formulas between physical properties as a bridge to connect different types of data. Inaccurate empirical relationships can lead joint inversion astray. Considering the complexity and variability of the tectonic evolution of the earth's materials, the petrophysical properties show different characteristics in different tectonic areas, and it is inappropriate to use the same empirical relationship in a broad study area. Therefore, how to construct the coupling relationship between the parameters of different geophysical models and how to improve the accuracy and efficiency of joint inversion are the key points of joint inversion research.

SUMMARY OF THE INVENTION

Aiming at the technical problems existing in the above-mentioned existing joint inversion, the present invention proposes a method of gravity-magnetic three-dimensional joint inversion based on cross gradient constraint with reasonable design, simple method, wide applicability and high inversion resolution.

In order to achieve the above object, the technical solution adopted in the present invention is that the present invention provides a method for a three-dimensional joint inversion of gravity and magnetism based on cross gradient constraints, comprising the following steps:

a. Divide the underground half-space into several regular cuboids in the form of a spatial array, and randomly set the initial physical property parameters in each cuboid, wherein the initial physical property parameters are the density parameter and the magnetization parameter;

b. Set up observation points on the ground at the vertical top of each cuboid;

c. Carry out 3D gravity and magnetic anomaly forward calculation with density parameters and magnetization parameters according to the divided cuboid, and obtain the predicted ground observation data;

d. Carry out inversion calculation according to the data actually measured at the observation point and the predicted ground observation data;

e. Finally, use the cross gradient function for joint inversion to obtain the optimal solution, where the cross gradient function is:

where m ⁽¹⁾ and m ⁽²⁾ are the density model parameters and the magnetization model parameters, respectively, k ₍₁₎ is the difference between the maximum and minimum values of the density model parameter m ⁽¹⁾ , and k ₍₂₎ is the magnetization The difference between the maximum and minimum values of the model parameter m ⁽²⁾ .

Preferably, in the step d, the iterative process of the inversion is divided into an outer loop and an inner loop, wherein the outer loop forward calculates the model fitting data, and selects an appropriate regularization factor according to the size of each initial value 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.

进入子循环。Preferably, in the step d, a larger damping coefficient α is set at the beginning of the outer loop iteration, and an appropriate calculation α can be selected according to the initial value of each item in the objective function.

Enter the sub-loop.

计算得到新模型，通过多次迭代直到迭代次数达到最大或得到模型变化趋于平稳或交叉梯度函数值增加，在迭代过程中，模型和交叉梯度函数收敛，得到结构相似的模型。Preferably, the inner loop is calculated from the outer loop

A new model is obtained by calculation, and iterates many times until the number of iterations reaches the maximum or the model changes tend to be stable or the value of the cross gradient function increases. During the iteration process, the model and the cross gradient function converge, and a model with a similar structure is obtained.

Compared with the prior art, the advantages and positive effects of the present invention are,

1. The present invention provides a method for three-dimensional joint inversion of gravity and magnetism based on cross gradient constraints. For the first time, the three-dimensional gravity and magnetism is used for joint inversion under the condition of cross gradient constraints of density parameters and magnetization parameters. Compared with traditional inversion , which reduces the multi-solution of gravity and magnetic inversion to a certain extent, improves the resolution of gravity and magnetic inversion, reduces the corresponding calculation, and improves the efficiency of inversion.

Description of drawings

In order to illustrate the technical solutions of the embodiments of the present invention more clearly, the following briefly introduces the accompanying drawings used in the description of the embodiments. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained from these drawings without any creative effort.

Figure 1 is a plan view of the Bouguer gravity (△Gb) isoline in the survey area;

Figure 2 is a plan view of the △T contour of the survey area;

Figure 3 is the result of joint inversion of residual density (0 isosurface);

Figure 4 is the result of joint inversion of magnetization (0.6A/m isosurface);

Figure 5 shows the result of joint inversion of residual density (0 isosurface);

Figure 6 is the result of joint inversion of magnetization (0.6A/m isosurface);

Figure 7(a) is the residual density profile of the 2D gravity-magnetic joint inversion result;

Fig. 7(b) is the magnetization profile of the 2D gravity-magnetic joint inversion result;

Fig. 8 is (a) the residual density result of the combined gravity and magnetic inversion of the Y=10000m section;

Fig. 8 is (b) Y=10000m section gravity-magnetic combined inversion magnetization result graph;

Fig. 9 is (a) the residual density result of the combined gravity and magnetic inversion of the Y=15000m section;

Figure 9 is (b) a graph of the result of combined gravity-magnetic inversion magnetization in the Y=15000m section;

Fig. 10 is (a) the residual density result of Y=20000m profile joint gravity and magnetic inversion;

Figure 10 is (b) the result of the combined gravity and magnetic inversion of the magnetization at the Y=20000m section;

Fig. 11 is (a) the residual density result of the joint gravity-magnetic inversion of the Y=25000m section;

Fig. 11 is (b) the result of the combined gravity and magnetic inversion of the magnetization in the Y=25000m section;

Figure 12 is (a) the residual density result of Y=30000m profile joint gravity-magnetic inversion;

Fig. 12 is (b) the result of combined gravity-magnetic inversion magnetization in the Y=30000m section;

Fig. 13 is a graph of the cross gradient value of Y=20000m inversion alone;

FIG. 14 is a cross-gradient value map of the joint inversion of Y=20000m.

Detailed ways

In order to more clearly understand the above objects, features and advantages of the present invention, the present invention will be further described below with reference to the accompanying drawings and embodiments. It should be noted that the embodiments of the present application and the features in the embodiments may be combined with each other in the case of no conflict.

Many specific details are set forth in the following description to facilitate a full understanding of the present invention, however, the present invention may also be implemented in other ways than those described herein, and therefore, the present invention is not limited to the specific embodiments of the following disclosure. limit.

Embodiment 1, this embodiment provides a method for 3D joint inversion of gravity and magnetism based on cross gradient constraints

First, the underground half-space is divided into several regular cuboids in the form of a spatial array, and the initial physical parameters in each cuboid are randomly set, wherein the initial physical parameters are density parameters and magnetization parameters. The calculation process of gravitational anomaly can be summarized into the following two main steps: (1) Calculate the gravitational potential at the observation point generated by the residual mass of a specific geological body relative to the surrounding rock according to Newton’s universal gravitational formula; (2) Calculate the above-mentioned gravitational potential along the The derivative of the vertical direction is used to obtain the gravity anomaly generated by the geological body for the corresponding observation point. To this end, in this embodiment, an observation point is set up on the ground at the vertical top of each cuboid.

Specifically, it is assumed that the surface observation plane is divided into a P×Q grid, the underground space is divided into L layers, and the grid of each layer is M×N. At this time, the number of integrations to be completed to calculate a forward modeling is: P×Q×M×N×L, which undoubtedly brings great difficulties to the calculation. Therefore, it is assumed here that the number of observation points on the ground is M×N, and each observation point is located vertically directly above the cuboid unit arranged according to the regular grid. According to the relative relationship between the observation point and the model unit, the amount of calculation can be greatly reduced. .

Suppose there are M observation points on the observation section, and there are M rectangular units corresponding to each layer of the underground model. For the first observation point, there are M relative distances r1, r2, ..., rM relative to a certain underground layer of model grid; for the second observation point, it is relative to a certain underground layer of model grid There are two relative distances (a pair) of the lattice that are the same, that is, r1=r3, so there are r1, r2, ..., rM, a total of M-1 relative distances; and so on, the Mth observation point and The existing relative distances of the first observation point are exactly the same, and the existing relative distances are r1, r2, ..., rM. Based on the above analysis, it can be seen that for the M models in the first layer of this profile, each observation point needs to calculate M distances, and M observation points need to calculate M2 distances. However, the distances of the subsequent M-1 measuring points are all different combinations of the calculated distances of the first measuring point. In fact, we only need to calculate the relative distance of M existences. For the three-dimensional model, there is a similar regularity, and only one dimension needs to be added on the basis of the above analysis.

According to the divided cuboid, the forward calculation of the three-dimensional gravity and magnetic anomaly is calculated by the density parameter and the magnetization parameter, and the predicted ground observation data is obtained.

Let the volume of the geological body be v and the density difference be σ. Let (ξ, η, ζ) be the coordinates of any point in the geological body. At this time, the volume element of this point can be expressed as dv=dξdηdζ, and its residual mass is written as dm=σdv=σdξddζ. The geological body is at the observation point A(x , y, z) at the gravitational potential is expressed as:

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

Divide the underground half-space into upright cuboids, and set the physical parameters in each cube to be constant constants. The discretized integral solution of the integral formula is given, and the specific form is as follows:

The underground medium is divided into M cuboids, and N observation points are arranged on the observation surface. At this time, the gravity anomaly generated by the jth cuboid at the ith observation point can be expressed as:

Δg=G _ij σ _j

where Gij is a known quantity determined according to 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.

The relationship between the magnetic potential and the gravitational potential can be expressed by the Poisson formula:

where U is the magnetic potential, G is the gravitational constant, M is the magnetization, σ is the density difference, and V is the gravitational potential. Available:

In the formula, I is the geomagnetic inclination angle, A' is the azimuth angle of the survey line, and

M _x =M·cosI·sinA; M _y =M·cosI·cosA; M=M·sinI

Substitute Za, Hax, and Hay into the expression of geomagnetic anomaly ΔT, we can calculate:

ΔT=H _ax cosIcosA'+H _ay cosIsinA'+Z _a sinI

Then, the inversion calculation is performed according to the data actually measured at the observation point and the predicted ground observation data.

The cross gradient function has the following properties

(1) When the change directions of the two physical parameters (density parameter and magnetization parameter) of the joint inversion are parallel, or one of the physical parameters remains unchanged, the cross gradient function of the two is equal to zero;

(2) When the change directions of the two physical parameters (density parameter and magnetization parameter) of the joint inversion are not parallel, the cross gradient function of the two is not equal to zero.

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

In the two-dimensional case, the cross gradient vector only contains non-zero components in the strike direction (y-axis direction), and the forward difference approximation is used in the numerical calculation:

In 3D space:

After discretizing the physical property model, the cross gradient also has its corresponding discrete form. When using a cuboid (rectangular in the two-dimensional case) grid, it can be discretized using difference:

In the formula, i, j and k are the unit numbers of the three-dimensional model along the coordinate direction, and x, y and z represent the physical property difference of the physical unit in the three-coordinate direction. The cross gradient vector can be written as:

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

where:

constraint:

τ=τ(m ⁽¹⁾ ,m ⁽²⁾ )=0

为先验模型，C_m、C_d分别为模型协方差矩阵和数据协方差矩阵，g为正演算子，α正则化项的权重因子(阻尼系数)，

σ为数据观测误差，为方便进行目标函数最优化计算，计算目标函数均取L2范数。Where m ⁽¹⁾ and m ⁽²⁾ are different geophysical models, and d ⁽¹⁾ and d ⁽²⁾ are the corresponding geophysical data

is the prior model, C _m and C _d are the model covariance matrix and the data covariance matrix, respectively, g is the forward operator, the weight factor (damping coefficient) of the α regularization term,

σ is the data observation error. In order to facilitate the optimization calculation of the objective function, the calculated objective function takes the L2 norm.

To solve the objective function formula of the joint inversion, it needs to be linearized. The Taylor series expansion method is used to linearize the forward function and the cross gradient function. To find the solution of the function with constraints, the Lagrange multiplication can be used. Subfad:

W=Φ(m)+2Λ ^T [t(m ₀ )+B(mm ₀ )]

where Λ is the Lagrangian operator, and B is the derivative of the cross gradient function.

The model update amount and Lagrangian operator can be obtained by taking the partial derivative of the model parameters equal to zero:

Δm ⁽ⁱ⁾ =N ^(i)-1 (n ⁽ⁱ⁾ -B ^(i)T Λ)

in:

The corresponding Lagrange multipliers:

Model update expression:

对模型的一阶导，扩展为：B is

The first derivative of the model is extended to:

为无交叉梯度约束的正则化最小二乘模型更新：

Update for a regularized least-squares model without cross-gradient constraints:

When the forward calculus operator is linear,

The physical property range constraints of the depth weighting and the transfer function method are introduced. The depth weighting adopts the approximate function proposed by Li and Oldenburg. The joint depth weighting matrix is the combination of the single gravity and magnetic method depth weighting diagonal matrix:

The physical parameter conversion function adopts the formula proposed by Commer and its inverse transformation, then the conversion joint model parameters are:

Its partial derivative matrix is:

Introducing the structural coupling factor β, the iterative equation of joint inversion becomes:

In the formula, the expressions of matrix N and vector n are:

In the formula, g is the transformed joint forward operator. The vector t satisfies the following system of linear equations:

The structural coupling factor β is introduced into the above formula. Using this formula for calculation can avoid the time-consuming singular value decomposition calculation and greatly improve the solution efficiency. It is worth mentioning that the choice of coupling factor is independent of other parameter choices. When it is close to infinity, the structural coupling disturbance term is approximately zero, and the inversion is close to a single inversion. When the coupling factor is large enough, it is equivalent to no damping, and the structural disturbance information is maximized.

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

进入子循环，由子循环得到与阻尼系数对应的理想的模型，计算misfit，通过减小阻尼系数α多次迭代，直到迭代次数达到最大或misfit达到理想值。(1) Main loop: Set a larger damping coefficient α at the beginning of the iteration, and select a suitable calculation α according to the initial value of each item in the objective function. calculate

Enter the sub-loop, obtain the ideal model corresponding to the damping coefficient from the sub-loop, calculate the misfit, and iterate many times by reducing the damping coefficient α until the number of iterations reaches the maximum or the misfit reaches the ideal value.

计算得到新模型，通过多次迭代直到迭代次数达到最大或得到模型变化趋于平稳或交叉梯度函数值增加，在迭代过程中，模型和交叉梯度函数收敛(计算中交叉梯度函数的约束式用线性近似)，得到结构相似的模型。Sub-loop: use the value calculated in the main loop

The new model is obtained by calculation, and iterates for many times until the number of iterations reaches the maximum or the model change tends to be stable 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 in the calculation is linear. approximation) to obtain a model with a similar structure.

Relative change of adjacent models:

Finally, when the empirical relationship between physical properties is 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 found:

Add k to prevent the instability of numerical calculation. In general, k is the expected change amplitude of m. Before the joint inversion, separate inversions are performed, and then the maximum and minimum values of each model parameter of the individual inversion results are obtained. The difference of the values is the corresponding k value. For the combined gravity and magnetic inversion, m ⁽¹⁾ and m ⁽²⁾ are the density model parameters and the magnetization model parameters, respectively.

Measured data inversion

Figure 1 is a plan view of the Bouguer gravity (△Gb) isoline in the survey area. The distance between the survey lines and the measuring points is 500m. As shown in the figure, the data reflects the Bouguer gravity anomaly in a certain area of Jiaodong Peninsula. The gravity in the survey area is The NE trend is dominant in the stepped zone, and the lithology on both sides has obvious density difference, which is consistent with the strike of the Jiaojia fault zone in the survey area. Figure 2 is a plan view of the △T contour of the survey area, which is the result after polarizing. The distance between the survey lines is 100m, and the distance between the survey points is 500m. As shown in the figure, the NE-NEE direction has an obvious positive anomaly, which is different from the gradient of gravity. The shape of the belt is similar. During the joint inversion of the cross-gradient method, the operation of the cross-gradient derivative Bk is involved. Bk is a matrix with the size of (3×Nm)×Nm, and Nm is the number of model parameters, so the calculation of tk requires huge memory, so a relatively rough grid is adopted, and the size of the model parameters is small. The grid is divided into 35 × 30 × 20, and the vertical and horizontal coordinates of 18,000 to 43,000 m are taken to carry out separate gravity and magnetic inversion and joint inversion under the same grid. Then, three-dimensional gravity and magnetic inversion and joint inversion are carried out to obtain five inversion result profiles of ordinate Y=1000m, ordinate Y=15000m, ordinate Y=20000m, ordinate Y=25000m, and ordinate Y=30000m. Compare the joint inversion with the individual inversion, and take one of the profiles Y=20000m to make the cross gradient value of the inversion result of this profile.

Through the analysis of Fig. 3 to Fig. 12, the results of individual inversion are cluttered, with irregular distribution of anomalies, and the magnitude of the anomaly deviates greatly from the actual situation, so it is difficult to analyze the required information. In the joint inversion diagram, the residual density results show that the northwest of the survey area is a positive residual density anomaly, and the southeast is a negative residual density anomaly. The direction of the belt is roughly the same. The magnetization 0.6A/m isosurface trend is roughly NE, and the magnetization is positive abnormality below the depth of 1000m.

Comparing Fig. 13 and Fig. 14, it can be seen that the joint inversion results have better structural similarity than the individual inversion results, which ensures that the fitting difference is small enough and also satisfies the constraint of the cross gradient value. The theoretically obtained model not only satisfies the data constraints, but also has good structural similarity between the models, which reduces the multi-solution of the inversion.

The above are only preferred embodiments of the present invention, and are not intended to limit the present invention in other forms. Any person skilled in the art may use the technical content disclosed above to make changes or modifications to equivalent changes. The embodiments are applied to other fields, but any simple modifications, equivalent changes and modifications made to the above embodiments according to the technical essence of the present invention still belong to the protection scope of the technical solutions of the present invention without departing from the content of the technical solutions of the present invention.

Claims (4)

1. a method based on the combined gravity and magnetic three-dimensional inversion of cross gradient constraint, is characterized in that, comprises the following steps: a. Divide the underground half-space into several regular cuboids in the form of a spatial array, and randomly set the initial physical property parameters in each cuboid, wherein the initial physical property parameters are the density parameter and the magnetization parameter; b. Set up observation points on the ground at the vertical top of each cuboid; c. Carry out 3D gravity and magnetic anomaly forward calculation with density parameters and magnetization parameters according to the divided cuboid, and obtain the predicted ground observation data; d. Carry out inversion calculation according to the data actually measured at the observation point and the predicted ground observation data; e. Finally, use the cross gradient function for joint inversion to obtain the optimal solution, where the cross gradient function is:

where m ⁽¹⁾ and m ⁽²⁾ are the density model parameters and the magnetization model parameters, respectively, k ₍₁₎ is the difference between the maximum and minimum values of the density model parameter m ⁽¹⁾ , and k ₍₂₎ is the magnetization The difference between the maximum and minimum values of the model parameter m ⁽²⁾ . 2. The method for combined gravity and magnetic three-dimensional inversion based on cross gradient constraints according to claim 1, wherein in the step d, the iterative process of the inversion is divided into an outer loop and an inner loop, wherein the outer loop The forward calculation model fits the data, and the appropriate regularization factor is selected according to the size of the initial values of the objective function; the inner loop modifies the model to obtain an ideal model corresponding to the regularization factor, calculates the error, and reduces the regularization factor multiple times. Iterate until the number of iterations reaches a maximum or the error reaches an ideal value to achieve structural similarity. 3. The method for combined gravity and magnetic three-dimensional inversion based on cross gradient constraints according to claim 2, wherein in the step d, the outer loop iteration starts to set a larger damping coefficient α, which can be determined according to the objective function. Select the appropriate calculation α for the size of each initial value, and calculate

Enter the sub-loop. 4 . The method for combined gravity and magnetic three-dimensional inversion based on cross gradient constraints according to claim 3 , wherein the inner loop is obtained by calculating in the outer loop. 5 .

A new model is obtained by calculation, and iterates many times until the number of iterations reaches the maximum or the model changes tend to be stable or the value of the cross gradient function increases. During the iteration process, the model and the cross gradient function converge, and a model with a similar structure is obtained.

Images