CN112199859A - Method for joint inversion of gravity gradient data - Google Patents

Abstract

The invention provides a method for jointly inverting gravity gradient data, relates to the technical field of geophysical inversion, and aims to solve the problems of insufficient spatial resolution and the like of the existing gravity gradient data inversion, a gradient depth weighting function, physical property constraint based on pair/index transformation and the like are introduced on the basis of focusing inversion, 6 components in full tensor gradient data are combined, and a nonlinear conjugate gradient method is used for realizing three-dimensional density inversion. Meanwhile, in order to improve the visualization effect and operability of inversion, a software platform with visualization function is developed for the proposed inversion method based on Python language, PyQt, Matplotlib and other tool packages, and the method effectively improves the distinguishing capability of adjacent geologic bodies and the imaging capability of the longitudinal space of the target body; the software platform and the development method have the advantages of easy use, practicability and the like.

Description

Method for joint inversion of gravity gradient data

Technical Field

The invention relates to the technical field of geophysical inversion, in particular to a method for jointly inverting gravity gradient data.

Background

The density distribution of the underground space can be obtained by utilizing gravity exploration data to realize three-dimensional inversion, which is a common interpretation method in geological exploration. Compared with gravity anomaly data, the signal to noise ratio of the gravity gradient data is higher, particularly the full-tensor gradient data contains more geological information, and the joint inversion carried out by using a plurality of gravity gradient data is beneficial to improving the defect of poor longitudinal spatial resolution of the current position field data and further improving the distinguishing capability of adjacent geologic bodies.

However, since the inversion has multi-solution, constraints or a priori information need to be introduced in the actual calculation to improve the uniqueness of the solution. The method disclosed in document [1] is to provide a range of physical property distribution and to impose the upper limit or the lower limit of the range if a result is not within the range (wherein document [1] is Li Y.,2001,3-D investment of Gravity Gradiometer data.2001 SEG annular meeting. society of application geographities.). However, this method destroys the conjugacy of the known search process, i.e., each iteration forces a change in the calculated solution. Meanwhile, the gravity gradient data has the property of attenuation along with the distance, so that the inversion result may have the problems of lower longitudinal resolution and the like. Some researchers have introduced depth weighting functions in documents [2-4] (where document [2] Li Y, Oldenburg D W.3-D inversion of the magnetic data [ J ]. Geophysics,1996,61(2): 394-. Therefore, it is necessary to use a constraint method for maintaining the conjugation and to use an effective depth weighting function.

Disclosure of Invention

Aiming at the defects of the prior art, the invention provides a method for jointly inverting gravity gradient data, which comprises the following steps:

step 1: constructing a target function phi (m) of three-dimensional physical inversion according to the acquired gravity gradient data;

in the formula, phi_d、φ_mRespectively representing a data non-fitting function and a model objective function; λ is a regularization parameter; g represents a sensitivity matrix, d represents gravity gradient data; m represents the inverted model parameter vector, m_jRepresenting the jth model parameter in a vector M, M representing the length of the vector M; σ represents a positive number that avoids generating singular values;

step 2: introducing the constraint range of the model parameter through the formula (2) to construct a model parameter m_jIs converted into a variable x_jLet j equal to 1,2, …, M, obtain the conversion variable corresponding to each model parameter by using formula (3), and represent all conversion variables by conversion vector x;

in the formula, a and b represent model parameters m_jLower limit value and upper limit value of (1);

and step 3: calculating a depth weighting function f (z);

where ω is a constant, r is a proportionality coefficient, dz and z are the depths of the underground divided cubes and the maximum depth of the inversion region, respectively, and z is the maximum depth of the inversion region_cRepresenting a depth coefficient based on prior geological information;

and 4, step 4: and (3) solving the target function phi (m) by using a nonlinear conjugate gradient method in combination with the depth weighting function f (z) to obtain the optimal inversion model parameters.

The step 4 comprises the following steps:

step 4.1: setting a maximum number of iterations k_maxInitializing a model parameter vector, and calculating by formula (3) to obtain an initial conversion vector x₀Calculating an initial gradient g using a depth weighting function f (z)₀；

In the formula, λ₀Representing an initial regularization parameter, λ₀＝0；

Step 4.2: when k equals 1 in the 1 st iteration, let λ₁＝λ₀Obtaining regularization parameters of 1 st iteration and initializing search direction vectors p₀＝g₀；

Step 4.3: calculating the search step length alpha in the k iteration by using an Armijo search method_kWherein k is 1,2_max；

Step 4.4: using the formula x_k＝x_k-1+α_kp_kComputing the transformation vector x at the kth iteration_k；

Step 4.5: convert vector x_kEach element x in (1)_jCalculating by using the formula (2) to obtain a corresponding model parameter m_jUsing all model parameters as vectors m_kExpressed by the formula d_k＝Gm_kCalculating fitting value d at k iteration_k；

Step 4.6: according to d_kD calculating the mean square error rms of the kth iteration_k；

Step 4.7: (ii) determination of rms_kWhether it is not greater than the threshold value epsilon,if yes, stopping iteration; if not, entering next iteration;

step 4.8: when the 2 nd iteration k is 2, the formula phi is firstly used_d＝||Gm-d||²Calculate phi_dUsing the formula

Calculate phi_mUsing the formula λ₂＝φ_d/φ_mCalculating a regularization parameter of the 2 nd iteration;

when k is more than or equal to 3, directly using formula lambda_k＝λ_k-12 updating regularization parameter lambda of the kth iteration_k；

Step 4.9: updating the gradient g using equation (6)_k；

Step 4.10: calculating an intermediate variable beta using equation (7)_k；

Step 4.11: updating the search direction p with equation (8)_k；

p_k＝g_k+β_kp_k-1 (8)

Step 4.12: let k be 1,2_maxWhen k is<k_maxRepeating the step 4.3 to the step 4.11 to carry out iterative calculation; mean square error rms at kth iteration_kStopping iteration and outputting the vector m when the value is less than or equal to epsilon_kThen vector m_kAll elements in (1) are the optimal inversion model parameters.

The invention has the beneficial effects that:

the invention provides a method for jointly inverting gravity gradient data, which effectively utilizes prior geological information such as physical properties, depth range and the like, and further improves the result precision compared with the existing method; the combination of the multi-component gravity gradient data improves the resolution ratio of underground space three-dimensional density imaging, better distinguishes adjacent geologic bodies and determines the buried depth of the top and bottom surfaces of a target body.

Drawings

FIG. 1 is a flow chart of a method of joint inversion of gravity gradient data in the present invention;

FIG. 2 is a software development flow diagram of a method for joint inversion of gravity gradient data according to the present invention;

FIG. 3 is a diagram of theoretical model positions in the present invention;

FIG. 4 is a graph of inversion results in accordance with the present invention;

FIG. 5 is a graph of measured data of the Wentton salt dome region in the present invention, in which (a) shows gravity gradient data V_xxComponent, graph (b) represents gravity gradient data V_xyComponent, graph (c) shows gravity gradient data V_xzComponent, plot (d) shows gravity gradient data V_yyComponent, graph (e) shows gravity gradient data V_yzComponent (d), graph (f) shows gravity gradient data V_zzA component;

fig. 6 is a diagram showing inversion results of measured data in the venturi region according to the present invention.

Detailed Description

The invention is further described with reference to the following figures and specific examples.

As shown in fig. 1, a method for jointly inverting gravity gradient data can accurately calculate three-dimensional density distribution, and effectively improve spatial resolution of a result, including the following steps:

step 1: constructing a target function phi (m) of three-dimensional physical inversion according to the acquired gravity gradient data, wherein the gravity gradient data can be obtained by constructing a theoretical model for simulation, or acquiring actually measured data by using a measuring instrument;

in the formula, phi_d、φ_mRespectively representing a data non-fitting function and a model objective function; λ is a regularization parameter; g denotes the sensitivity matrix and d denotes the gravity gradientData; m represents the inverted model parameter vector, m_jRepresenting the jth model parameter in a vector M, M representing the length of the vector M; σ represents a positive number that avoids generating singular values;

in the present embodiment, the acquired gravity gradient data is first subjected to conventional gridding processing, and V in the full-tensor gravity gradient data is selected_xx、V_xy、V_xz、V_yy、V_yz、V_zzThe six components are subjected to joint inversion.

Step 2: introducing the constraint range of the model parameter through the formula (2) to construct a model parameter m_jIs converted into a variable x_jLet j equal to 1,2, …, M, obtain the conversion variable corresponding to each model parameter by using formula (3), and represent all conversion variables by conversion vector x;

in the formula, a and b represent model parameters m_jLower limit value and upper limit value of (1);

and step 3: calculating a depth weighting function f (z);

where ω is an empirical constant, r is a proportionality coefficient, dz and z are the depths of the underground divided cubes and the maximum depth of the inversion region, respectively, and z is the maximum depth of the inversion region_cRepresenting a depth coefficient based on prior geological information;

and 4, step 4: combining with a depth weighting function f (z), solving an objective function phi (m) by using a nonlinear conjugate gradient method to obtain an optimal inversion model parameter, wherein the optimal inversion model parameter comprises the following steps:

step 4.1: setting a maximum number of iterations k_maxInitializing the model parameter vector and passingCalculating to obtain an initial conversion vector x by formula (3)₀Calculating an initial gradient g using a depth weighting function f (z)₀；

In the formula, λ₀Representing an initial regularization parameter, λ₀＝0；

Step 4.2: when k equals 1 in the 1 st iteration, let λ₁＝λ₀Obtaining regularization parameters of 1 st iteration and initializing search direction vectors p₀＝g₀；

Step 4.3: calculating the search step length alpha in the k iteration by using an Armijo search method_kWherein k is 1,2_max；

Step 4.4: using the formula x_k＝x_k-1+α_kp_kComputing the transformation vector x at the kth iteration_k；

Step 4.5: convert vector x_kEach element x in (1)_jCalculating by using the formula (2) to obtain a corresponding model parameter m_jUsing all model parameters as vectors m_kExpressed by the formula d_k＝Gm_kCalculating fitting value d at k iteration_k；

Step 4.6: according to d_kD calculating the mean square error rms of the kth iteration_k；

Step 4.7: (ii) determination of rms_kIf not, stopping iteration; if not, entering next iteration;

step 4.8: when the 2 nd iteration k is 2, the formula phi is firstly used_d＝||Gm-d||²Calculate phi_dUsing the formula

Calculate phi_mUsing the formula λ₂＝φ_d/φ_mCalculating a regularization parameter of the 2 nd iteration;

when k is more than or equal to 3, directly using formula lambda_k＝λ_k-12 updating regularization parameter lambda of the kth iteration_k；

Step 4.9: updating the gradient g using equation (6)_k；

Step 4.10: calculating an intermediate variable beta using equation (7)_k；

Step 4.11: updating the search direction p with equation (8)_k；

p_k＝g_k+β_kp_k-1 (8)

Step 4.12: let k be 1,2_maxWhen k is<k_maxRepeating the step 4.3 to the step 4.11 to carry out iterative calculation; mean square error rms at kth iteration_kStopping iteration and outputting the vector m when the value is less than or equal to epsilon_kThen vector m_kAll elements in the data are the optimal inversion model parameters; when k is more than or equal to k_maxWhen it is time, the operation is stopped.

In order to improve the utility of the present application, a software development method based on a gravity gradient data joint inversion method is provided below, as shown in fig. 2, including the following steps:

1) building a development environment on a Windows system by using tool packages such as Python language, PyQt, Matplotlib and the like through tools such as Eric and Microsoft Visual Studio and the like;

2) aiming at the data format characteristics of input data and output results, modules such as a graphical interface, data management, drawing and the like of a software platform are designed by utilizing PyQt and Matplotlib;

3) writing the inversion code using C language, generating dll file by Microsoft Visual Studio, wherein the pseudo code programmed is as follows:

wherein k is_maxIs the maximum number of iterations.

4) Self-defining an init function, adding all variables (input, output, intermediate variables and the like) involved in calculation into the init function, and using the init function for transferring parameters; the custom class "class custom class", is transmitted into the algorithm function in dll;

5) designing a module (comprising functions of calculating, drawing, storing data and the like) and a graphical interface thereof aiming at inverted input parameters; connecting the graphical interface with the functional function through a signal slot mechanism in the PyQt;

6) calling a self-defined init function and class custom class in the main function to realize an inversion function;

7) after the compilation runs, the software is tested.

By using the software written above, the space imaging capability and noise immunity of the inversion method are tested by using a bicube combined model containing 5% of Gaussian noise, the position and inversion result of the theoretical model are respectively shown in fig. 3 and 4, and the inversion parameters are set as follows: lambda [ alpha ]₀0, density constraint value a-0.1, b-0.5, depth weighting function ω -0.001, r-20, z_c325, dz 1000. The inversion result is compared with the actual position of the model, the inversion method disclosed by the invention favorably distinguishes the adjacent cubes and calculates the top and bottom depths of the two cubes, the numerical range of the result is consistent with that of the model, and the effectiveness and the noise immunity of the method are proved.

The practical applicability and feasibility of the inversion method and software for testing the measured data in the Wendon salt dome region are adopted, the measured data in the Wendon salt dome region are shown in the graphs (a) to (f) in FIG. 5 (E (Herpot) is the unit of gravity gradient data in the graph), and the inversion results are respectively shown in FIG. 6. The inversion parameters were set as follows: lambda [ alpha ]₀0, density constraint value a-0.1, b-0.5, depth weighting function ω -0.001, r-20, z_c225 and dz 1000. The results show that the density distribution calculated by the inversion method is similar to the existing geological data and other people research results, but the underground cover rock imaging effect is clearer and betterThe generated software has good mapping effect, and the practicability and feasibility of the method are proved.

Claims (2)

1. A method for joint inversion of gravity gradient data is characterized by comprising the following steps:

step 1: constructing a target function phi (m) of three-dimensional physical inversion according to the acquired gravity gradient data;

in the formula, phi_d、φ_mRespectively representing a data non-fitting function and a model objective function; λ is a regularization parameter; g represents a sensitivity matrix, d represents gravity gradient data; m represents the inverted model parameter vector, m_jRepresenting the jth model parameter in a vector M, M representing the length of the vector M; σ represents a positive number that avoids generating singular values;

step 2: introducing the constraint range of the model parameter through the formula (2) to construct a model parameter m_jIs converted into a variable x_jLet j equal to 1,2, …, M, obtain the conversion variable corresponding to each model parameter by using formula (3), and represent all conversion variables by conversion vector x;

in the formula, a and b represent model parameters m_jLower limit value and upper limit value of (1);

and step 3: calculating a depth weighting function f (z);

where ω is a constant, r is a proportionality coefficient, dz and z are the depths of the underground divided cubes and the maximum depth of the inversion region, respectively, and z is the maximum depth of the inversion region_cRepresenting a depth coefficient based on prior geological information;

and 4, step 4: and (3) solving the target function phi (m) by using a nonlinear conjugate gradient method in combination with the depth weighting function f (z) to obtain the optimal inversion model parameters.

2. The method for jointly inverting gravity gradient data according to claim 1, wherein the step 4 comprises:

step 4.1: setting a maximum number of iterations k_maxInitializing a model parameter vector, and calculating by formula (3) to obtain an initial conversion vector x₀Calculating an initial gradient g using a depth weighting function f (z)₀；

In the formula, λ₀Representing an initial regularization parameter, λ₀＝0；

Step 4.2: when k equals 1 in the 1 st iteration, let λ₁＝λ₀Obtaining regularization parameters of 1 st iteration and initializing search direction vectors p₀＝g₀；

Step 4.3: calculating the search step length alpha in the k iteration by using an Armijo search method_kWherein k is 1,2_max；

Step 4.4: using the formula x_k＝x_k-1+α_kp_kComputing the transformation vector x at the kth iteration_k；

Step 4.5: convert vector x_kEach element x in (1)_jCalculating by using the formula (2) to obtain a corresponding model parameter m_jUsing all model parameters as vectors m_kExpressed by the formula d_k＝Gm_kCalculating fitting value d at k iteration_k；

Step (ii) of4.6: according to d_kD calculating the mean square error rms of the kth iteration_k；

Step 4.7: (ii) determination of rms_kIf not, stopping iteration; if not, entering next iteration;

step 4.8: when the 2 nd iteration k is 2, the formula phi is firstly used_d＝||Gm-d||²Calculate phi_dUsing the formula

Calculate phi_mUsing the formula λ₂＝φ_d/φ_mCalculating a regularization parameter of the 2 nd iteration;

when k is more than or equal to 3, directly using formula lambda_k＝λ_k-12 updating regularization parameter lambda of the kth iteration_k；

Step 4.9: updating the gradient g using equation (6)_k；

Step 4.10: calculating an intermediate variable beta using equation (7)_k；

Step 4.11: updating the search direction p with equation (8)_k；

p_k＝g_k+β_kp_k-1 (8)

Step 4.12: let k be 1,2_maxWhen k is<k_maxRepeating the step 4.3 to the step 4.11 to carry out iterative calculation; mean square error rms at kth iteration_kStopping iteration and outputting the vector m when the value is less than or equal to epsilon_kThen vector m_kAll elements in (1) are the optimal inversion model parameters.

Images