CN114200541B - Three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraint - Google Patents

Abstract

A 3D gravity-magnetism joint inversion method based on cosine dot product gradient constraints, including gravity, gravity gradient and magnetic method 3D forward modeling fast algorithm based on geometric grid method; 3D cosine dot product gradient function is derived based on the law of cosines, based on A three-dimensional joint inversion algorithm of the truncated Gauss-Newton method, which can not only avoid storing the Jacobian matrix, but also reduce the number of forward calculations when the conjugate gradient method is used to solve linear equations. The full tensor gravity gradient data is added to the traditional gravity and magnetic method data, which can further improve the resolution ability of the gravity-magnetic joint inversion; In the gradient and magnetic method data, the advantages of the three aspects are fully utilized to effectively improve the calculation efficiency, memory consumption and resolution ability of the 3D gravity-magnetic joint inversion.

Description

A 3D Gravity-Magnetic Joint Inversion Method Based on Cosine Dot Product Gradient Constraint

technical field

The invention relates to the processing technical field of a joint inversion of three-dimensional gravity and magnetic data, and particularly proposes a high-precision and high-efficiency three-dimensional gravity and magnetic joint inversion method based on cosine dot product gradient constraints.

Background technique

As the exploration target gets deeper and the exploration environment becomes more and more complex, it is difficult to accurately infer the attribute information of the underground target body based on a single geophysical data. Joint inversion is one of the effective methods for comprehensive interpretation of multiple geophysical methods, and it is an effective means to reduce inversion multi-solutions, improve the utilization of multiple data, and improve the inconsistency of inversion model structures. The cross-gradient joint inversion does not depend on the prior petrophysical relationship, and does not need to make any assumptions about the initial model. It is a flexible and effective joint inversion method, which has been applied to two or more geophysical methods. . However, the cross-gradient function directly cross-multiplies model parameters of different orders of magnitude and units, only considering whether the gradient directions are parallel, ignoring the influence of different model parameter gradient amplitudes on the cross-gradient constraint weights of different regions, and the need to calculate more in three-dimensional cases. components increase the computational complexity. Therefore, how to construct the coupling relationship between parameters of different geophysical models, and how to improve the accuracy and efficiency of joint inversion are the key points of joint inversion research.

Contents of the invention

The purpose of the present invention is to provide a high-precision and high-efficiency three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraints, so as to solve the problems of slow calculation speed and large memory consumption of the existing gravity-magnetic joint inversion.

As shown in Figure 1, the present invention comprises the following steps:

1) Input the observed data of gravity, gravity gradient and magnetic method, input the initial model and reference model of density and magnetization, expected fitting difference and maximum number of iterations;

2) Construct the 3D gravity-magnetic joint inversion objective functions Φ _g (m _g ) and Φ _κ (m _κ ) based on cosine dot product gradient constraints;

3) Calculate the regularization factor and the structural constraint weight factor;

4) Calculate the gravity and magnetic data covariance matrix, model integral sensitivity matrix, focus weighting matrix, and calculate the cosine dot product gradient function and its Jacobian matrix;

5) Using the truncated Gauss-Newton method to optimize and solve the objective function of the three-dimensional gravity-magnetic joint inversion constrained by the gradient of the cosine dot product, and obtain the update amount of the density and magnetization model;

6) Calculate the forward modeling response of the density and magnetization update model, and solve the fitting difference of each iteration of the gravity and magnetic method. If each method meets the expected fitting difference or reaches the maximum number of iterations, the joint inversion iteration ends.

Beneficial effects of the present invention:

1. The cosine dot product gradient function of the present invention does not need to calculate multiple components in a three-dimensional case, which simplifies the complexity of calculating traditional cross gradient functions and improves the calculation efficiency of structural coupling. At the same time, it avoids the calculation impact caused by the magnitude and dimension differences of model parameters, and is more suitable for large-scale three-dimensional multi-parameter joint inversion calculations.

2. The truncated Gauss-Newton method of the present invention is a combination of the Gauss-Newton method and the conjugate gradient method. When this method is applied to the joint inversion, it can avoid storing the Jacobian matrix, and can greatly reduce the number of forward calculations, which can effectively The computational efficiency and memory consumption of 3D gravity-magnetic joint inversion can be greatly improved.

3. The present invention adds full tensor gravity gradient data to the traditional gravity and magnetic method data, which can further improve the resolving power of gravity-magnetic joint inversion.

Description of drawings

Figure 1 is a flow chart of the three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraints;

Figure 2 is a schematic diagram of three-dimensional underground grid division;

Fig. 3 is the schematic diagram of density and magnetization theoretical model;

Figure 4 is a schematic diagram of a three-dimensional gravity and magnetic inversion density model;

Fig. 5 is a schematic diagram of a three-dimensional gravity and magnetic inversion magnetization model.

Detailed ways

As shown in Figure 1, a 3D gravity-magnetic joint inversion method based on cosine dot product gradient constraints includes the following steps:

1) Input the observed data of gravity, gravity gradient and magnetic method, input the initial model and reference model of density and magnetization, expected fitting difference and maximum number of iterations;

2) Construct the 3D gravity-magnetic joint inversion objective functions Φ _g (m _g ) and Φ _κ (m _κ ) based on cosine dot product gradient constraints;

3) Calculate the regularization factor and the structural constraint weight factor;

4) Calculate the gravity and magnetic data covariance matrix, model integral sensitivity matrix, focus weighting matrix, and calculate the cosine dot product gradient function and its Jacobian matrix;

5) Using the truncated Gauss-Newton method to optimize and solve the objective function of the three-dimensional gravity-magnetic joint inversion constrained by the gradient of the cosine dot product, and obtain the update amount of the density and magnetization model;

6) Calculate the forward modeling response of the density and magnetization update model, and solve the fitting difference of each iteration of the gravity and magnetic method. If each method meets the expected fitting difference or reaches the maximum number of iterations, the joint inversion iteration ends.

In step 1), the observation data of gravity, gravity gradient and magnetic method are input, and the observation data are obtained by establishing a theoretical model of density and magnetization to calculate the forward response.

In step 2), the objective function of the 3D gravity-magnetic joint inversion based on the cosine dot product gradient constraint is constructed as follows:

Among them, C _dg is the data covariance matrix of gravity and gravity gradient observation data d _g ; C _dκ is the data covariance matrix of magnetic observation data d _κ ; C _mg and C _mκ are the density m _g and magnetization m _κ respectively Model covariance matrix, m _gref and m _κref are reference models, set to the model obtained in the last iteration; λ _g , λ _κ , β _g and β _g are weighting factors; f _g (m _g ) and f _κ (m _κ ) are respectively the forward modeling response of gravity and gravity gradient and magnetic method, and t _d (m _g ,m _κ ) is the cosine dot product gradient function. The expression is as follows:

和

为密度和磁化强度模型梯度，|·|为范数。Among them, ε is a very small constant,

and

is the density and magnetization model gradient, |·| is the norm.

In step 3), the regularization factor and structural constraint weight factor are calculated.

Set initial regularization factors λ _g =0.1; λ _κ =0.001. With the increase of iterations, the model constraints may increase. In order to ensure that the convergence of the objective function is the smallest in the whole area, the regularization factor needs to be corrected at this time. The expression is as follows:

In order to ensure the iterative convergence of the joint inversion, we will reduce the regularization factor proportionally when the fitting difference increases or the fitting difference changes slowly:

Among them, k is the number of iterations, q is the scaling factor, set q=0.6;

The calculation expression of structural constraint weight factor is as follows:

where C _g and C _κ are weighting coefficients for density and magnetization structural constraints.

In step 4), calculate the gravity and magnetic data covariance matrix, model integral sensitivity matrix, focus weighting matrix, and calculate the cosine dot product gradient function and its Jacobian matrix;

The data covariance matrix is the inverse diagonal matrix of the noise error of gravity, gravity gradient and magnetic method data. When there is no noise interference, it is expressed as a unit diagonal matrix. The calculation expression of the focusing matrix is as follows:

W _s =(m ² +ε ² ) ^(-1/2) (9)

The model integral sensitivity matrix expression is as follows:

Among them, A is the Jacobian matrix.

The calculation expression of the cosine dot product gradient function is as follows:

Among them, η _g and η _κ are normalization operators, the second subscripts c, b, d and r of the model parameters m _g and m _κ represent the units at different positions of the three-dimensional grid respectively, and Δx, Δy, Δz represent is the distance between grid cells in different directions, as shown in Figure 2.

The Jacobian matrix expression of the cosine dot product gradient function is as follows:

Among them, B _g and B _κ are the partial derivatives of the cosine dot product gradient constraint function with respect to density and magnetization, respectively, namely the Jacobian matrix.

In step 5), the truncated Gauss-Newton method is used to optimize and solve the objective function of the three-dimensional gravity-magnetic joint inversion constrained by the gradient of the cosine dot product, and the updated quantities of the density and magnetization models are obtained as follows:

A_g和A_κ为重磁正演响应雅克比矩阵。Among them, k is the number of iterations,

A _g and A _κ are Jacobian matrices of gravity and magnetic forward modeling responses.

Formulas (13) and (14) are not solved directly, but the Gauss-Newton iterative linear equation is approximated by using the conjugate gradient method, and the accuracy of the solution is controlled by a forcing term. The expression of the forcing term is as follows:

In step 6), the forward modeling response of the density and magnetization update model is calculated, and the fitting difference of each iteration of the gravity and magnetic method is solved;

By analyzing the positional relationship between the underground model unit and the observation point, it can be found that there is an equivalence between the relative positional relationship between the model units and the observation point on the same layer. Using this equivalence relationship, only the first model of each layer Calculate the kernel function of the geometric grid to obtain the kernel function matrix of the unit, and the kernel function matrix of other model units in this layer can be obtained by searching the equivalence of the geometric grid, so as to quickly calculate the gravity and magnetic forward modeling response; The expression for calculating the difference is as follows:

Among them, d is the observation data, f is the forward modeling response, and N is the number of observation points. If each method meets the expected poor fit or reaches the maximum number of iterations, the joint inversion iteration ends.

The three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraints provided by the present invention is verified below.

As shown in Figure 3, a double anomalous body combination model is established; the model includes two anomalous bodies with the same size and different buried depths, the size of the two anomalous bodies is 200m×200m×200m, and the residual density is 1g/cm ³ and a magnetization of 1A/m. The spatial scale x, y, z of the underground model is 1000m×1000m×600m, the underground model is divided into 20×20×12 closely arranged regular cubic units, each unit size is 50m×50m×50m, and the residual density of the background is 0g/cm ³ , the background magnetization is 0.001A/m, the number of ground observation points is 20×20, and the observation points are located in the center of the grid.

First, the effectiveness of the cosine dot product-gradient joint inversion algorithm is verified by comparing the results of gravity and magnetic inversion alone and joint inversion. Then, the gravity gradient data is added to the gravity-magnetic joint inversion to carry out the cosine dot product gradient joint inversion calculation to verify the resolution ability and calculation efficiency of the algorithm. All the above methods are carried out under the framework of the truncated Gauss-Newton method, and the initial model is set as the residual density model of the uniform half space 0g/cm ³ and the magnetization model of 0.001A/m.

The separate inversion results of gravity and magnetic data are shown in Fig. 4(a,b,c) and Fig. 5(a,b,c). It can be found that the gravity inversion density model alone cannot restore the geometric profile and physical parameter values of the deep anomaly, while the magnetization model can better restore the geometric shape and physical parameter values of the anomalous body. When the cosine dot product gradient structure constraints are added, the joint inversion results of gravity and magnetic data (Fig. 4d, e, f and Fig. 5d, e, f) can effectively improve the resolution of the density model, especially the deep anomaly of the density model The geometry of the anomalous body is obviously restored, and the physical parameter values of the anomalous body are also close to the real values. Compared with the single inversion, the cosine dot product gradient joint inversion obtains a density and magnetization model with better structural consistency. . Finally, the gravity gradient data is added to carry out the joint inversion calculation of the three-dimensional cosine dot product gradient constraints of gravity, gravity gradient and magnetic method data. The density and magnetization models are shown in (4g,h,l) and (5g,h, As shown in l), due to the increase in the amount of data, the density model has been further improved, and it is more consistent with the real density model, especially due to the constraints of the cosine dot product gradient structure, a model with higher consistency with the magnetization structure has been obtained; at the same time , the calculation time of this example is about 217.5s, and the maximum memory consumption is about 0.117GB. Compared with the traditional gravity-magnetic three-dimensional joint inversion algorithm, the algorithm proposed by the present invention is a high-efficiency, high-resolution joint inversion algorithm.

Claims (3)

1. A three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraint is characterized by comprising the following steps: the method comprises the following steps:

1) Inputting observation data of gravity, gravity gradient and a magnetic method, inputting an initial model and a reference model of density and magnetization intensity, and expecting a fitting difference and a maximum iteration number;

2) Construction of three-dimensional gravity-magnetic joint inversion target function phi based on cosine dot product gradient constraint _g (m _g ) And phi _κ (m _κ )；

3) Calculating a regularization factor and a structural constraint weight factor;

4) Calculating a weight magnetic data covariance matrix, a model integral sensitivity matrix, a focusing weighting matrix, and calculating a cosine dot product gradient function and a Jacobian matrix thereof;

5) Utilizing a truncated gauss-newton method to carry out optimization solution on a three-dimensional gravity-magnetic joint inversion target function constrained by cosine dot product gradient, and obtaining the updating amount of a density and magnetization intensity model;

6) Calculating forward response of the density and magnetization updating model, solving fitting difference of each iteration of the gravity and magnetization method, and finishing joint inversion iteration if each method meets expected fitting difference or reaches maximum iteration times;

in the step 2), a three-dimensional gravity-magnetic joint inversion target function based on cosine dot product gradient constraint is constructed as follows:

wherein, C _dg For gravity and gravity gradient observation data d _g The data covariance matrix of (a); c _dκ For magnetic observation of data d _κ The data covariance matrix of (a); c _mg And C _mκ Respectively is density m _g And magnetization m _κ Model covariance matrix of (1), m _gref And m _κref Setting the model obtained by the last iteration as a reference model; lambda [ alpha ] _g And λ _κ As a regularization factor, beta _g And beta _g Constraining the weight factor for the structure; f. of _g (m _g ) And f _κ (m _κ ) Expressed as the forward response of gravity and gravity gradient and magnetism, respectively, t _d (m _g ,m _κ ) The cross gradient function is a cosine dot product gradient function, a plurality of components do not need to be calculated under the three-dimensional condition, the complexity of calculating the traditional cross gradient function is simplified, and the calculation efficiency of structural coupling is improved; meanwhile, the calculation influence caused by the magnitude order of the model parameters and the dimension difference is avoided, and the method is more suitable for large-data-volume three-dimensional multi-parameter joint inversion calculation; the cosine dot product gradient function is expressed as follows:

where, epsilon is a very small constant,

and

is the density and magnetization model gradient, | · | is a norm;

in step 3), the constraint weight factors of the gravity and magnetic structure are obtained by an adaptive method, and the updating expression of the constraint weight factors of the gravity and magnetic structure in each iteration is as follows:

wherein, C _g And C _κ Weighting coefficients for density and magnetization;

in step 4), the cosine dot product gradient function calculation expression is as follows:

wherein eta is _g And η _κ For normalization operator, model parameter m _g And m _κ The second lower corner marks c, b, d and r respectively represent the units at different positions of the three-dimensional grid, and Δ x, Δ y and Δ z represent the distances between grid units in different directions;

the calculation of the three-dimensional cosine dot product gradient function Jacobian matrix is more complex compared with two-dimensional calculation, each row of non-zero elements of the Jacobian matrix B is changed into eight elements from six elements, and the expression is as follows:

wherein, B _g And B _κ The partial derivatives of the cosine dot product gradient constraint function with respect to density and magnetization, respectively.

2. The three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraint according to claim 1, characterized in that:

in the step 5), solving the joint inversion by using the traditional Gauss-Newton method needs to occupy a large amount of memory space to store the Jacobian matrix; the truncated gauss-newton method is a combination of the gauss-newton method and a conjugate gradient method, each iteration is performed with an external cycle and an internal cycle, the external cycle adopts the gauss-newton method to calculate a coefficient matrix and a right-end term, the internal cycle approximately solves a gauss-newton iteration linear equation by using the conjugate gradient method, and the accuracy of the solution is controlled by a forcing term; the method not only avoids storing the Jacobian matrix, but also can greatly reduce forward calculation times, and can further improve the calculation efficiency and memory consumption when being applied to three-dimensional gravity-magnetic joint inversion with cosine dot product gradient constraint; the coerced expression is as follows:

wherein min is expressed as the minimum value, k is the iteration number, A _g And A _κ Is a gravity magnetic forward response Jacobian matrix,

3. the three-dimensional gravity-magnetic joint inversion method based on cosine dot product gradient constraint according to claim 1, characterized in that:

in the step 6), a geometric grid method is used for carrying out three-dimensional forward calculation on gravity, gravity gradient and a magnetic method, so that forward calculation time can be greatly reduced; meanwhile, the forward modeling needs to be solved for multiple times by truncating the internal cycle of the Gauss-Newton method, so that the calculation efficiency of the three-dimensional gravity-magnetic joint inversion can be improved from the three aspects of the forward modeling, the coupling mode and the inversion by combining the geometric lattice method with the truncating Gauss-Newton method and the cosine dot product gradient structure constraint.

Images