CN110398782B - A Joint Regularization Inversion Method for Gravity Data and Gravity Gradient Data - Google Patents

Abstract

The invention relates to a joint regularization inversion method for gravity data and gravity gradient data, comprising the following steps: step 1: acquiring gravity data, performing inversion calculation on the gravity data, and obtaining an inversion result; The inversion result m takes the absolute value and normalizes it. For the zero value in the result, it is set to a minimum value, and the processed result is used as a weighting matrix; Step 3: Obtain the gravity gradient data, and compare the gravity gradient data with all The above weighting matrix is simultaneously applied to the inversion calculation, and the inversion result of the joint inversion of the gravity gradient data and the weighting matrix is obtained. The invention effectively utilizes the geological information contained in the gravity data, constructs a weighted matrix by using the gravity data, applies the weighted matrix to the gravity gradient inversion calculation, improves the resolution of the inversion result, and makes the bottom boundary of the inversion result clearer and more accurate.

Description

A Joint Regularization Inversion Method for Gravity Data and Gravity Gradient Data

technical field

The invention relates to the technical field of potential field data inversion methods, in particular to a joint regularization inversion method of gravity data and gravity gradient data.

Background technique

The gravity data is obtained by measuring the vertical component of the potential field, and the gravity gradient data is obtained by measuring the change of the gravity field in three directions. If the Cartesian coordinate system is taken as an example, the vertical downward direction is expressed by the positive direction of the Z axis, and the gravity The data are the first-order partial derivatives of the gravity field in the Z-axis direction, and the gravity gradient data are the second-order partial derivatives of the gravity field in the three directions of the X, Y, and Z axes. Comparing the gravity data and the gravity gradient data in the frequency domain, it can be found that the gravity data contains more low-frequency information, and the gravity gradient data contains more high-frequency information. Therefore, the simultaneous application of gravity and gravity gradient data to the inversion calculation can achieve information complementarity and improve the inversion effect.

Since the potential field lacks the resolution in the depth direction, it is necessary to introduce a depth weighting function in the inversion of gravity data and gravity gradient data. The setting of the depth weighting function is related to the data decay rate. The decay rate of gravity data with distance is proportional to r-2, and the decay rate of gravity gradient with distance is proportional to r-3, where r represents the distance between two mass blocks. However, in most joint inversion methods of gravity and gravity gradient, the same depth weighting function is used for gravity and gravity gradient data to offset the decay rate, and the difference in decay rate with distance is not considered. At the same time, in the common joint inversion of gravity and gravity gradient data, the gravity data and gravity gradient data are generally directly formed into a matrix and participate in the inversion calculation. Because the magnitude of gravity data is lower than that of gravity gradient data, and in the same matrix In this method, the geological information contained in the gravity data and gravity gradient data cannot be extracted well, resulting in a low vertical resolution of the inversion results.

SUMMARY OF THE INVENTION

In view of the deficiencies of the prior art, the purpose of the present invention is to provide a joint regularization inversion method of gravity data and gravity gradient data, which can solve the problem that the vertical resolution of the inversion results in the joint inversion of gravity data and gravity gradient data is not high. question.

The technical scheme for realizing the object of the present invention is: a joint regularization inversion method for gravity data and gravity gradient data, comprising the following steps:

Step 1: Acquire gravity data, perform inversion calculation on the gravity data, and obtain the inversion result. The objective function φ _d of inverting the gravity data is formula (1):

表示二-范数的平方运算，d表示观测数据，A表示灵敏度矩阵，m表示模型参数，A_w＝AW^-1，m_w＝Wm，W表示深度加权函数，是一个对角矩阵；in,

Represents the square operation of the two-norm, d represents the observation data, A represents the sensitivity matrix, m represents the model parameter, A _w =AW ^-1 , m _w =Wm, and W represents the depth weighting function, which is a diagonal matrix;

Step 2: Take the absolute value of the inversion result m obtained in Step 1, and perform normalization processing. For the zero value in the result, set it as a minimum value, and the processed result is used as a weighting matrix;

Step 3: Obtain the gravity gradient data, apply the gravity gradient data and the weighting matrix to the inversion calculation at the same time, and obtain the inversion result of the joint inversion of the gravity gradient data and the weighting matrix.

Further, in the step 1, the specific process of inverting the gravity data includes:

(1) The number of iterations is represented by i, the maximum number of iterations is set to n, the initial value of i is 0, i=0, 1, 2, ..., n, determine the value of the initial model m ₀ , and obtain formula (2) :

m _w0 =Wm ₀ ------(2)

(2) Calculate the initial gradient I ₀ of the objective function according to formula (3):

I ₀ =A _w (A _w m _w0 -d)------(3)

For the first iterative calculation, the search direction is d ₀ =-I ₀ ;

(3) Calculate the i-th search step size k _{(i) according to formula (4)} :

Among them, d _i represents the observation data of the ith iteration, and I _i represents the gradient of the objective function corresponding to the ith iteration;

(7) Calculate the result mi of the _ith iteration according to formula (5), and calculate the intermediate variable r according to formula (6):

Among them, _mw(i-1) represents the result of the i-1th iteration, and _mw(i) is an intermediate variable;

r=Am _i -d------(6)

小于等于预设值或者迭代次数达到最大迭代次数n时，终止迭代，否则继续按步骤(6)进行；(8) When

When it is less than or equal to the preset value or the number of iterations reaches the maximum number of iterations n, the iteration is terminated, otherwise, proceed to step (6);

设定新的搜索方向为d_i+1＝-I_i+1+β_i+1d_i，并跳转至步骤(3)处理。(6) Calculate the gradient I _i+1 of the objective function corresponding to the i+1th iteration, and calculate the intermediate parameters

Set the new search direction as d _i+1 =-I _i+1 +β _i+1 d _i , and jump to the processing of step (3).

Further, in the step 3, the gravity gradient data and the weighting matrix are simultaneously applied to the inversion calculation, and constraints are added through the regularization method, and the gravity gradient data and the weighting matrix are jointly inverted. The objective function φ is formula (7):

φ＝φ _d +αφ _s ------(7)

Among them, φ _s represents the stability function, which is calculated according to formula (8):

φ _s = α _s ∫∫∫ _V m ² d _v +α _x ∫∫∫ _V m _x ² d _v +α _y ∫∫∫ _V m _y ² d _v +α _z ∫ ∫∫ _V m _z ² d _v − -----(8)

Among them, α _s , α _x , α _y , and α _z represent the weights of each component, which are constants, m _x , my _y , and m _z represent the first-order partial derivatives of the model parameter m in the three directions of x, y, and z, and V represents the integral area, and the weighting function W is introduced into the objective function φ, so the objective function φ of formula (7) becomes formula (9):

Formula (9) is the objective function that the gravity gradient data and the weighting matrix are applied to the inversion calculation. The conjugate gradient algorithm is used to solve the formula (9), and the inversion result is obtained. The inversion result is the gravity data and gravity gradient data. Inversion results after joint regularization inversion.

The beneficial effects of the invention are as follows: the invention effectively utilizes the geological information contained in the gravity data, constructs a weighted matrix by using the gravity data, applies the weighted matrix to the gravity gradient inversion calculation, improves the resolution of the inversion results, and the inversion results Bottom border is sharper and more accurate.

Description of drawings

Fig. 1 is the schematic flow chart of the present invention;

Figure 2 shows the results of joint inversion of gravity data and gravity gradient data using the traditional method;

FIG. 3 is a joint inversion result of gravity data and gravity gradient data according to the present invention.

specific implementation

Below, in conjunction with the accompanying drawings and specific embodiments, the present invention is further described:

As shown in Figures 1 to 3, a joint regularization inversion method for gravity data and gravity gradient data includes the following steps:

Step 1: Acquire gravity data, perform inversion calculation on gravity data, and its objective function φ _d is as in formula (1), that is, the inversion of gravity data is calculated according to formula (1):

表示二-范数的平方运算，d表示观测数据，是一个向量，大小为n₁×1，n₁表示观测点个数，A表示灵敏度矩阵，可以采用现有的通用计算公式得到，其大小为n₁×n₂，n₂表示测区划分的小长方体个数，A_ij为A中的元素，表示编号为i的小长方体在编号为j的测点引起的重力或重力梯度值，密度为1g/cm³，m表示模型参数，为向量，也即是本步骤通过公式(1)需要得到的反演结果，φ_d表示拟合函数。in,

Represents the square operation of the two-norm, d represents the observation data, is a vector, the size is n ₁ × 1, n ₁ represents the number of observation points, A represents the sensitivity matrix, which can be obtained by using the existing general calculation formula, its size is n ₁ ×n ₂ , n ₂ represents the number of small cuboids divided by the survey area, A _ij is the element in A, and represents the gravity or gravity gradient value caused by the small cuboid number i at the measuring point number j, the density is 1g/cm ³ , m represents the model parameter, and is a vector, that is, the inversion result that needs to be obtained by formula (1) in this step, and φ _d represents the fitting function.

Since the decay rate of gravity data with distance is proportional to r-2, and the decay rate of gravity gradient with distance is proportional to r-3, the depth weighting function W needs to be introduced to offset this effect, and formula (2) is obtained:

Among them, A _w =AW ^-1 , m _w =Wm, W is a diagonal matrix, the conjugate gradient algorithm is used to solve the formula (2), and m is obtained, and the specific process is as follows:

(1) The number of iterations is represented by i, the maximum number of iterations is set to n, the initial value of i is 0, that is, i=0, 1, 2, ..., n, and the value of the initial model m ₀ is determined. In the case of verification information, generally take the value of 0, so as to obtain formula (3):

m _w0 =Wm ₀ ------(3)

(2) Calculate the initial gradient I ₀ of the objective function according to formula (4):

I ₀ =A _w (A _w m _wi -d)------(4)

For the first iterative calculation, the search direction is d ₀ =-I ₀ ;

(3) Calculate the i-th search step k _{(i) according to formula (5)} :

Among them, d _i represents the observation data of the ith iteration, and I _i represents the gradient of the objective function corresponding to the ith iteration;

(4) Calculate the result mi of the _ith iteration, as in formula (6), if m _i is less than or equal to the preset value, it means that the value of the element in A conforms to the geological information, and the intermediate value is calculated according to formula (7). variable r:

Among them, _mw(i-1) represents the result of the i-1th iteration, and _mw(i) is an intermediate variable;

r=Am _i -d------(7)

小于等于预设的极小值(例如10^-6或更小)或者迭代次数达到最大迭代次数n时，终止迭代，否则继续按步骤(6)(5) When

If it is less than or equal to the preset minimum value (for example, 10 ^-6 or less) or the number of iterations reaches the maximum number of iterations n, terminate the iteration, otherwise continue to step (6)

conduct:

设定新的搜索方向为d_i+1＝-I_i+1+β_i+1d_i，并跳转至步骤(3)处理。(6) Calculate the gradient I _i+1 of the objective function corresponding to the i+1th iteration, and calculate the intermediate parameters

Set the new search direction as d _i+1 =-I _i+1 +β _i+1 d _i , and jump to the processing of step (3).

Step 2: Take the absolute value of the inversion result m obtained in Step 1, and perform normalization processing. For the zero value in the result, set it to a minimum value (for example, ^10-6 or less), and the processed result is used as weighting matrix;

Step 3: Obtain gravity gradient data, apply the gravity gradient data and the weighting matrix obtained in step 2 to the inversion calculation at the same time, and obtain the inversion result.

The inversion calculation in this step adopts the smooth inversion algorithm to calculate, and the specific process is as follows:

Constraints are added through the regularization method, and the objective function φ is formula (8):

φ＝φ _d +αφ _s ------(8)

Among them, φ _s represents the stability function, which is calculated according to formula (9):

φ _s = α _s ∫∫∫ _V m ² d _v +α _x ∫∫∫ _V m _x ² d _v +α _y ∫∫∫ _V m _y ² d _v +α _z ∫ ∫∫ _V m _z ² d _v − -----(9)

Among them, α _s , α _x , α _y , and α _z represent the weights of each component, which are coefficients and constants, which can be preset or given according to empirical values, and m _x , my _y , and m _z represent the model parameters m in x, The first-order partial derivatives in the three directions of y and z, V represents the integral area, and φ _s can also be represented in matrix form, such as formula (10):

Wi = α _i D _i , _i = x, y, z, where Wi represents the difference operator, and D _i is the matrix calculated according to the finite difference operator in different directions _, and the weight function W is introduced into the objective function, so , the objective function φ of formula (8) becomes formula (11):

Formula (11) is the objective function that the gravity gradient data and the weighting matrix are applied to the inversion calculation. Similarly, the conjugate gradient algorithm is used to solve the formula (11), and the inversion result is obtained, and the inversion result is also the present invention. The final inversion result is the inversion result after the joint regularization inversion of the gravity data and the gravity gradient data. The specific process is the same as that of step 1 to solve the objective function, and the specific process will not be repeated.

It can be seen from the above description that the present invention effectively utilizes the geological information contained in the gravity data, constructs a weighted matrix by using the gravity data, applies the weighted matrix to the inversion calculation of the gravity gradient, improves the resolution of the inversion result, and the bottom boundary of the inversion result is clearer and more accurate.

As shown in Figure 2 and Figure 3, it is the result of inversion of the actual gravity and gravity gradient data of the Vinton Salt Dome in the United States. The left side of Figure 2 and Figure 3 is a horizontal section, and the abscissa and ordinate represent the east-west and north-south directions. Position coordinates, the middle column and the right column in Figure 2 and Figure 3 are vertical sections, the ordinate represents the depth, and the abscissa represents the horizontal position. Fig. 2 shows the gravity data gz and four different combinations of gravity gradient components (g _zz , g _xx |g _yy |g _zz , g _xx |gxy|g _yy |g _zz , g _xx |g _xy |g _xz ) |g _yy |g _yz |g _zz ) joint inversion results, and Fig. 3 shows the results obtained by using the method of the present invention for the corresponding component combination.

Compared with the existing method for joint inversion of gravity data and gravity gradient data, the inversion result of the present invention is clearer in the deep boundary, the convergence of the inversion result is better, and the resolution of the inversion result is improved.

The embodiment disclosed in this specification is only an illustration of the unilateral feature of the present invention, and the protection scope of the present invention is not limited to this embodiment, and any other functionally equivalent embodiments fall within the protection scope of the present invention. For those skilled in the art, various other corresponding changes and deformations can be made according to the technical solutions and concepts described above, and all these changes and deformations should fall within the protection scope of the claims of the present invention.

Claims (3)

1. a joint regularization inversion method of gravity data and gravity gradient data, is characterized in that, comprises the steps: Step 1: Acquire gravity data, perform inversion calculation on the gravity data, and obtain the inversion result. The objective function φ _d of inverting the gravity data is formula (1):

in,

Represents the square operation of the two-norm, d represents the observation data, A represents the sensitivity matrix, m represents the model parameter, A _w =AW ^-1 , m _w =Wm, and W represents the depth weighting function, which is a diagonal matrix; Step 2: Take the absolute value of the inversion result m obtained in Step 1, and perform normalization processing. For the zero value in the result, set it as a minimum value, and the processed result is used as a weighting matrix; Step 3: Obtain the gravity gradient data, apply the gravity gradient data and the weighting matrix to the inversion calculation at the same time, and obtain the inversion result of the joint inversion of the gravity gradient data and the weighting matrix. 2. The joint regularization inversion method of gravity data and gravity gradient data according to claim 1, is characterized in that, in described step 1, the concrete process that gravity data is carried out inversion calculation comprises: (1) The number of iterations is represented by i, the maximum number of iterations is set to n, the initial value of i is 0, i=0, 1, 2, ..., n, determine the value of the initial model m ₀ , and obtain formula (2) : m _w0 =Wm ₀ ------(2) (2) Calculate the initial gradient I ₀ of the objective function according to formula (3): I ₀ =A _w (A _w m _w0 -d)------(3) For the first iterative calculation, the search direction is d ₀ =-I ₀ ; (3) Calculate the i-th search step size k _{(i) according to formula (4)} :

Among them, d _i represents the observation data of the ith iteration, and I _i represents the gradient of the objective function corresponding to the ith iteration; (4) Calculate the result mi of the _ith iteration according to formula (5), and calculate the intermediate variable r according to formula (6):

Among them, _mw(i-1) represents the result of the i-1th iteration, and _mw(i) is an intermediate variable; r=Am _i -d------(6) (5) When

When it is less than or equal to the preset value or the number of iterations reaches the maximum number of iterations n, the iteration is terminated, otherwise, proceed to step (6); (6) Calculate the gradient I _i+1 of the objective function corresponding to the i+1th iteration, and calculate the intermediate parameters

Set the new search direction as d _i+1 =-I _i+1 +β _i+1 d _i , and jump to the processing of step (3). 3. The gravity data and gravity gradient data joint regularization inversion method according to claim 1, is characterized in that, in described step 3, to gravity gradient data and described weighting matrix are simultaneously applied to inversion calculation, by The regularization method adds constraints, and the objective function φ of the joint inversion of the gravity gradient data and the weighting matrix is formula (7): φ＝φ _d +αφ _s ------(7) Among them, φ _s represents the stability function, which is calculated according to formula (8): φ _s = α _s ∫∫∫ _V m ² d _v +α _x ∫∫∫ _V m _x ² d _v +α _y ∫∫∫ _V m _y ² d _v +α _z ∫ ∫∫ _V m _z ² d _v − -----(8) Among them, α _s , α _x , α _y , and α _z represent the weights of each component, which are constants, m _x , my _y , and m _z represent the first-order partial derivatives of the model parameter m in the three directions of x, y, and z, and V represents the integral area, and the weighting function W is introduced into the objective function φ, so the objective function φ of formula (7) becomes formula (9):

Formula (9) is the objective function that the gravity gradient data and the weighting matrix are applied to the inversion calculation. The conjugate gradient algorithm is used to solve the formula (9), and the inversion result is obtained. The inversion result is the gravity data and gravity gradient data. Inversion results after joint regularization inversion.

Images