CN115201927A

CN115201927A - Method for heavy magnetic three-dimensional joint inversion based on cross gradient constraint

Info

Publication number: CN115201927A
Application number: CN202210876606.6A
Authority: CN
Inventors: 张扬; 杜利明; 梅贞华; 胡创业; 陈琦; 雷会东; 刘升台; 冯娜; 姚晓杰; 伊善强
Original assignee: CHINA METALLURGICAL GEOLOGY BUREAU GEOLOGICAL EXPLORATION INSTITUTE OF SHANDONG ZHENGYUAN
Current assignee: CHINA METALLURGICAL GEOLOGY BUREAU GEOLOGICAL EXPLORATION INSTITUTE OF SHANDONG ZHENGYUAN
Priority date: 2022-07-25
Filing date: 2022-07-25
Publication date: 2022-10-18

Abstract

The invention belongs to the technical field of geophysical application, and relates to a gravity-magnetic three-dimensional joint inversion method based on cross gradient constraint. Dividing the underground half space into a plurality of regular cuboids in a space array mode, and randomly setting initial physical parameters in each cuboid, wherein the initial physical parameters are density parameters and magnetization intensity parameters; setting an observation point on the ground right above each cuboid; carrying out three-dimensional gravity magnetic anomaly forward modeling calculation according to the divided cuboids by using density parameters and magnetization parameters to obtain predicted ground observation data; calculating errors according to actually measured data of the observation points and predicted ground observation data; and finally, performing joint iterative inversion by using a cross gradient function, so that the error is reduced to a reasonable range, and obtaining an ideal model. Compared with the traditional inversion, the method reduces the multi-solution property of the heavy magnetic inversion to a certain extent, improves the resolution of the heavy magnetic inversion, reduces corresponding calculation and improves the inversion efficiency.

Description

Method for heavy magnetic three-dimensional joint inversion based on cross gradient constraint

Technical Field

The invention belongs to the technical field of geophysical application, and particularly relates to a method for heavy magnetic three-dimensional joint inversion based on cross gradient constraint.

Background

Geophysical response is generated by differences in the physical properties of the subsurface medium, and geophysical inversion is a discipline for inferring the physical distribution of the subsurface medium from the observed geophysical response, such as the distribution of physical properties such as density, magnetic susceptibility or resistivity. This process will necessarily face two fundamental problems of geophysical exploration, namely the inherent ambiguity of inversion itself and the one-sidedness of single-method inversion. However, the two basic problems can be well solved only by combining the guidance of geological theoretical knowledge and walking a road integrating the development of geophysical science. Therefore, on the basis of the inversion of a single geophysical method, the combined geophysical inversion is gradually developed.

In the same measuring area range of the same work area, data of various geophysical methods can be collected, such as gravity anomaly, magnetic anomaly, resistivity information and the like. Although there appears to be no correlation between the data corresponding to these different geophysical methods, the subsurface medium corresponding to these measured data is uniform, and therefore, we can consider using these collected geophysical responses in combination to estimate a geological-geophysical model of the subsurface, including the values of the physical parameters, the buried depth, shape, size of scale, and the like of the model. Different physical property parameter values meet a certain internal corresponding relation, so that the non-uniqueness of inversion can be improved to a certain extent, and the accuracy of inversion interpretation is improved. Therefore, joint inversion is an important direction for future geophysical development and is also attracting more and more attention of geophysical workers.

With the continuous exploitation of earth surface and near-earth surface resources, the depth of geophysical exploration is continuously increased at present, the problems are more and more complicated, the geophysical inversion of a single method cannot be competent for the actual problems at present, and the research and the application of the geophysical joint inversion are necessary and urgent.

At present, the most common inversion modes in the prior art are sequential inversion and synchronous joint inversion, and the sequential inversion and the synchronous joint inversion both need to use an empirical formula between physical properties as a bridge to link different types of data, and an inaccurate empirical relationship can bring the joint inversion into the way. Considering the complex variability of the evolution of the earth's material structure, petrophysical properties exhibit different characteristics in different tectonic zones, and it is not appropriate to use the same empirical relationships in a wide area of investigation. Therefore, how to construct the coupling relationship between different geophysical model parameters and how to improve the accuracy and efficiency of joint inversion are the key points of the current joint inversion research.

Disclosure of Invention

Aiming at the technical problems of the conventional joint inversion, the invention provides the heavy magnetic three-dimensional joint inversion method based on the cross gradient constraint, which has the advantages of reasonable design, simple method, wide applicability and high inversion resolution.

In order to achieve the above object, the present invention adopts a technical solution that the present invention provides a method for gravity-magnetic three-dimensional joint inversion based on cross gradient constraint, comprising the following steps:

a. dividing the underground half space into a plurality of regular cuboids in a space array mode, and randomly setting initial physical parameters in each cuboid, wherein the initial physical parameters are density parameters and magnetization intensity parameters;

b. setting an observation point on the ground right above each cuboid;

c. carrying out three-dimensional gravity magnetic anomaly forward modeling calculation according to the divided cuboids by using density parameters and magnetization parameters to obtain predicted ground observation data;

d. performing inversion calculation according to the actually measured data of the observation points and the predicted ground observation data;

e. and finally, performing joint inversion by using a cross gradient function to obtain an optimal solution, wherein the cross gradient function is as follows:

wherein m is ⁽¹⁾ 、m ⁽²⁾ Respectively density model parameter and magnetization model parameter, k ₍₁₎ As a density model parameter m ⁽¹⁾ Is the difference between the maximum and minimum values of (c), k ₍₂₎ For the magnetization model parameter m ⁽²⁾ The difference between the maximum value and the minimum value of (c).

Preferably, in the step d, an iterative process of inversion is divided into an outer loop and an inner loop, wherein the outer loop forward modeling calculates model fitting data, and a proper regularization factor is selected according to the size of each initial value in the target function; and modifying the model by the inner loop to obtain an ideal model corresponding to the regularization factor, calculating an error, and iterating for multiple times by reducing the regularization factor until the iteration times reach the maximum or the error reaches an ideal value so as to achieve structural similarity.

Preferably, in the step d, a larger damping coefficient alpha is set at the beginning of the outer loop iteration, and a proper calculation alpha can be selected according to the size of each initial value in the objective function to calculate

And entering a sub-cycle.

Preferably, the inner loop is calculated from the outer loop

Calculating to obtain a new model, and repeating the iteration for multiple times until the iteration number reaches the maximum or obtaining the modelThe change tends to be stable or the value of the cross gradient function increases, and in the iterative process, the model and the cross gradient function converge to obtain a model with a similar structure.

Compared with the prior art, the invention has the advantages and positive effects that,

1. the invention provides a gravity-magnetic three-dimensional joint inversion method based on cross gradient constraint, which is characterized in that joint inversion is carried out by three-dimensional gravity under the condition of the cross gradient constraint of density parameters and magnetization parameters for the first time, compared with the traditional inversion, the multi-solution of gravity-magnetic inversion is reduced to a certain extent, the resolution of gravity-magnetic inversion is improved, corresponding calculation is reduced, and the inversion efficiency is improved.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings needed to be used in the description of the embodiments are briefly introduced below, and it is obvious that the drawings in the following description are some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings based on these drawings without inventive labor.

FIG. 1 is a plan view of a region under test with a Booth gravity (Δ Gb) contour;

FIG. 2 is a plan view of a contour line of a measurement area Δ T;

FIG. 3 is a graph of the results of the joint inversion of residual densities (0 iso-surface);

FIG. 4 is a graph of the results of joint inversion of magnetization (0.6A/m iso-surface);

FIG. 5 is a graph of the results of the joint inversion of residual densities (0 iso-surface);

FIG. 6 is a graph of the results of joint inversion of magnetization (0.6A/m iso-surface);

FIG. 7 (a) is a two-dimensional gravity-magnetic joint inversion result residual density profile;

FIG. 7 (b) is a magnetization profile of the result of two-dimensional gravity-magnetic joint inversion;

fig. 8 is a graph of (a) residual density results of Y =10000m section gravity-magnetic joint inversion;

fig. 8 is (b) a graph of magnetization results of a Y =10000m section gravity-magnetic joint inversion;

fig. 9 is a graph of (a) residual density results of Y =15000m section gravity-magnetic joint inversion;

fig. 9 is (b) a graph of magnetization results of Y =15000m section remagnetism joint inversion;

fig. 10 is a (a) Y =20000m section gravity-magnetic joint inversion residual density result diagram;

fig. 10 is a (b) Y =20000m section-gravity-magnetic joint inversion magnetization result graph;

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

fig. 11 is a graph of (b) magnetization results of the Y =25000m section gravity-magnetic joint inversion;

fig. 12 is a (a) Y =30000m section gravity-magnetic joint inversion residual density result graph;

fig. 12 is (b) a graph of magnetization results of Y =30000m section gravity-magnetic joint inversion;

fig. 13 is a graph of Y =20000m individual inverted cross-gradient values;

fig. 14 is a graph of Y =20000m jointly inverted cross-gradient values.

Detailed Description

In order that the above objects, features and advantages of the present invention can be more clearly understood, the present invention will be further described with reference to the accompanying drawings and examples. It should be noted that the embodiments and features of the embodiments of the present application may be combined with each other without conflict.

In the following description, numerous specific details are set forth in order to provide a thorough understanding of the present invention, however, the present invention may be practiced otherwise than as specifically described herein and, therefore, the present invention is not limited to the specific embodiments disclosed in the following description.

Embodiment 1, the embodiment provides a gravity-magnetic three-dimensional joint inversion method based on cross gradient constraint

Firstly, dividing an underground half space into a plurality of regular cuboids in a space array mode, and randomly setting initial physical parameters in each cuboid, wherein the initial physical parameters are density parameters and magnetization intensity parameters. The calculation process of the gravity anomaly is summarized as the following two main steps: (1) Calculating the gravitational potential of the specific geologic body relative to the residual mass of the surrounding rock at the observation point according to a Newton's universal gravitation formula; (2) And calculating the derivative of the gravitational potential along the vertical direction to obtain the gravity anomaly of the geologic body to the corresponding observation point. For this reason, an observation point is set up on the ground right above the vertical direction of each rectangular parallelepiped in the present embodiment.

Specifically, assuming that the earth surface observation plane is divided by P × Q grids, the underground space is divided into L layers, and the grid of each layer is M × N, and the number of integration times required to be completed by calculating one forward course is P × Q × M × N × L, which undoubtedly brings great difficulty to calculation. Therefore, assuming that the number of observation points on the ground is M × N, each observation point is located right above the rectangular units arranged according to the regular grid, and the calculation amount can be greatly reduced according to the relative relationship between the observation points and the model units.

And M observation points exist on the observation section, and M rectangular units correspondingly exist in each layer of the underground model. For a first observation point, the relative distances of r1, r2, rM and rM exist relative to a layer of model grids in the underground; for the second observation point, two (one pair) of relative distances to a certain layer of model mesh in the underground are the same, i.e. r1= r3, so that r1, r 2. By analogy, the relative distance of the mth observation point to the first observation point is identical, and the relative distance of the mth observation point to the first observation point is r1, r 2. From the above analysis, it can be seen that for the M models of the first layer of the section, 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 different combinations of the distances calculated for the first measuring point. In fact, we only need to calculate the M existing relative distances. Similar regularity exists for the three-dimensional model, and only one dimension needs to be added on the basis of the analysis.

And calculating three-dimensional gravity magnetic anomaly forward modeling calculation according to the density parameters and the magnetization parameters of the divided cuboids to obtain predicted ground observation data.

Let v be the volume of the geologic body and σ be the density difference. Let (ξ, η, ζ) be the coordinates of any point in the geologic volume, where the volume element of that point can be represented as dv = d ξ d η d ζ, with its residual mass written as dm = σ dv = σ d ξ dd ζ, the gravitational potential of the geologic volume at observation point a (x, y, z) is represented as:

where G is a universal gravitational constant, r represents a distance from any point (ξ, η, ζ) in the geologic body to an observation point, and the gravity anomaly can be written as:

the underground half space is divided into vertical cuboids, and the physical property parameters in each cuboid are constant. A discretized integral solution of the integral formula is given, in particular in the form:

the underground medium is divided into M cuboids, N observation points are distributed on an observation surface, and 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 from the relative position of the jth cuboid at the ith observation point, and σ j is the residual density of the jth cuboid.

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

in the formula, U is a magnetic potential, G is a universal gravitation constant, M is magnetization intensity, sigma is density difference, and V is a gravitational potential. The following can be obtained:

wherein 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

substituting Za, hax and Hay into the geomagnetic anomaly delta T expression to obtain:

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

and then, performing inversion calculation according to the data actually measured by the observation point and the predicted ground observation data.

The cross-gradient function has the following properties

(1) When the change directions of two physical parameters (density parameter and magnetization parameter) of the joint inversion are parallel, or one of the physical parameters is kept unchanged, the cross gradient function of the two physical parameters 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 physical parameters is not equal to zero.

These properties indicate that: the joint inversion based on the cross-gradient coupling achieves the purpose of improving the reliability of the inversion result by enhancing the consistency of the structures reflected by different physical parameters.

In the two-dimensional case, where the cross-gradient vector contains non-zero components only in the direction of the trend (y-direction), the numerical calculation is approximated by the forward difference:

in three-dimensional space:

after discretizing the physical model, the cross-gradient also has its corresponding discrete form. When a cuboid (rectangular in the case of two dimensions) grid is employed, it can be discretized using differential pairs:

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

combining the density, the magnetization intensity and the speed parameter, the gravity magnetic-shock observation data, the cross gradient vector, a first derivative matrix of the cross gradient function and the like into a unified vector, and considering the target function:

in the formula:

and (3) constraint:

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

wherein m is ⁽¹⁾ 、m ⁽²⁾ For different geophysical models, d ⁽¹⁾ 、d ⁽²⁾ Corresponding geophysical data

Being a prior model, C _m 、C _d Respectively, a model covariance matrix and a data covariance matrix, g is a forward operator, a weight factor (damping coefficient) of an alpha regularization term,

and sigma is a data observation error, and L2 norms are taken for calculating the objective function for conveniently carrying out the optimization calculation of the objective function.

In order to solve the target function formula of the joint inversion, the target function formula needs to be linearized, a forward function and a cross gradient function are linearized by using a taylor series expansion method, a solution of a function with constraint conditions is solved, and the solution can be obtained by adopting a lagrange multiplier method:

W＝Φ(m)+2Λ ^T [t(m ₀ )+B(m-m ₀ )]

where Λ is the lagrange operator and B is the derivative of the cross-gradient function.

And (3) solving the model parameter to obtain model updating quantity and Lagrange operator with the partial derivative equal to zero:

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

wherein:

the corresponding lagrange multiplier:

model update expression:

b is

For a first derivative of the model, the extension is:

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

when the forward operator is linear in nature,

and (3) introducing depth weighting and physical property range constraint of a conversion function method, wherein the depth weighting adopts an approximate function proposed by Li and Oldenburg, and a combined depth weighting matrix is a combination of a single gravity and magnetic method depth weighting diagonal matrix:

and the physical property parameter conversion function adopts the inverse conversion of a formula provided by the Commer, and then the conversion joint model parameters are:

the partial derivative matrix is:

introducing a structural coupling factor beta, and changing a joint inversion iteration equation into:

in the formula, the expressions of the matrix N and the vector N are respectively:

wherein g is the transformed joint syndrome. The vector t satisfies the following linear system of equations:

structural coupling factor beta is introduced into the formula, and the formula is used for calculation, so that time-consuming singular value decomposition calculation can be avoided, and the solving efficiency is greatly improved. It is worth mentioning that the selection of the coupling factor is independent of the other parameter selections. When the inversion is close to infinite small, the structural coupling disturbance term is approximate to zero, and the inversion approaches to single inversion; when the coupling factor is large enough, the structure disturbance information is maximally reserved equivalently without damping.

The iterative process of inversion is divided into an inner loop and an outer loop, the outer loop is a main loop, the main loop enables the model to be fitted with data, the inner loop is a sub-loop, and the sub-loop modifies the model to achieve structural similarity.

(1) Main circulation: and (3) setting a larger damping coefficient alpha at the beginning of iteration, and selecting a proper calculation alpha according to the size of each initial value in the objective function. Computing

And entering a sub-loop, obtaining an ideal model corresponding to the damping coefficient by the sub-loop, calculating misfit, and iterating for multiple times by reducing the damping coefficient alpha until the iteration times reach the maximum or the misfit reaches an ideal value.

And (3) sub-circulation: calculated in the main loop

And calculating to obtain a new model, and iterating for multiple times until the iteration times reach the maximum or the obtained model changes to be stable or the cross gradient function value is increased, wherein in the iteration process, the model and the cross gradient function are converged (a constraint formula of the cross gradient function in calculation is linearly approximated) to obtain a model with a similar structure.

Relative variation of neighboring models:

finally, when the empirical relationship between physical properties is unknown or ambiguous, the method is very suitable for introducing cross gradient constraints to carry out joint inversion. As long as the model parameters are coupled, the cross-gradient function can be solved:

k is added to prevent instability of numerical calculation, in general, k is an expected change amplitude of m, separate inversion is respectively carried out before joint inversion, then the difference between the maximum value and the minimum value of each model parameter of a separate inversion result is a corresponding k value, and for gravity and magnetic joint inversion, m is a maximum value and a minimum value of each model parameter of the separate inversion result are used as k values ⁽¹⁾ 、m ⁽²⁾ Respectively, density model parameters and magnetization model parameters.

Inversion of measured data

FIG. 1 is a plane view of a grid distribution gravity (delta Gb) contour line of a measurement area, the distance between measurement lines and the interval between measurement points are both 500m, as shown in the figure, data reflects abnormal grid distribution gravity of a certain area of the Jiaodong peninsula, a gravity step belt in the measurement area takes the NE trend as the main direction, lithology on two sides has obvious density difference, and the trend corresponds to the trend condition of a Jiajia fracture zone in the measurement area. FIG. 2 is a plan view of the contour lines of the measurement area Δ T, as a result of polarization, with a measurement line spacing of 100m and a measurement point spacing of 500m, as shown in the figure, the NE-NEE direction has a distinct positive anomaly, similar to the morphology of the gravity step band. In the process of carrying out the cross-gradient method joint inversion, the operation of the cross-gradient derivative Bk is involved. Bk is a matrix with the size of (3 multiplied by Nm) multiplied by Nm, nm is the number of model parameters, so that the calculation of tk needs huge memory, a coarser grid subdivision is adopted, the size of the model parameters is smaller, for the measurement area, the grid subdivision is 35 multiplied by 30 multiplied by 20, and areas with the ordinate of 5000-35000 m and the abscissa of 18000-43000 m are respectively subjected to independent gravity-magnetic inversion and joint inversion under the same subdivision grid. And then carrying out three-dimensional gravity magnetic individual inversion and joint inversion, and making five inversion result sections with the ordinate Y =1000m, the ordinate Y =15000m, the ordinate Y =20000m, the ordinate Y =25000m and the ordinate Y =30000 m. And (3) comparing the joint inversion with the single inversion, taking one section Y =20000m, and making a cross gradient value of the section inversion result.

By analyzing fig. 3 to 12, the form of the result of the independent inversion is disordered, the abnormal distribution is irregular, and the deviation between the abnormal amplitude and the actual situation is large, so that the required information is difficult to analyze. In the joint inversion diagram, the residual density result shows that the north and the west of the measurement area are the positive residual density anomaly, the south and the east are the negative residual density anomaly, the trend of the density residual density 0 isosurface is approximately NNE-NE, and the change rule approximately accords with the trend of the Jiaojia fracture zone in the area. The contour profile of magnetization of 0.6A/m is approximately NE, and a depth of 1000m or less indicates a positive magnetization anomaly.

Comparing fig. 13 and fig. 14, it can be seen that the joint inversion result has better structural similarity than the individual inversion result, ensuring that the fitting difference is small enough and also satisfying the constraint of the cross gradient value. The theoretically obtained model not only meets the data constraint, but also has better structural similarity between the models, and reduces the multiple solution of inversion.

The above description is only a preferred embodiment of the present invention, and not intended to limit the present invention in other forms, and any person skilled in the art may apply the above modifications or changes to the equivalent embodiments with equivalent changes, without departing from the technical spirit of the present invention, and any simple modification, equivalent change and change made to the above embodiments according to the technical spirit of the present invention still belong to the protection scope of the technical spirit of the present invention.

Claims

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

b. setting an observation point on the ground right above each cuboid;

wherein m is ⁽¹⁾ 、m ⁽²⁾ Respectively density model parameter and magnetization model parameter, k ₍₁₎ As a density model parameter m ⁽¹⁾ Maximum sum ofDifference of minimum values, k ₍₂₎ As a parameter m of the magnetization model ⁽²⁾ The difference between the maximum value and the minimum value of (d).

2. The cross-gradient-constraint-based gravity-magnetic three-dimensional joint inversion method according to claim 1, wherein in the step d, an iterative process of inversion is divided into an outer loop and an inner loop, wherein the outer loop forward models model fitting data, and an appropriate regularization factor is selected according to the size of each initial value in a target function; and modifying the model by the inner loop to obtain an ideal model corresponding to the regularization factor, calculating an error, and iterating for multiple times by reducing the regularization factor until the iteration times reach the maximum or the error reaches an ideal value so as to achieve structural similarity.

3. The cross-gradient-constraint-based gravity-magnetic three-dimensional joint inversion method according to claim 2, wherein in the step d, a larger damping coefficient alpha is set at the beginning of outer loop iteration, and a proper calculation alpha can be selected according to the initial values of each item in the objective function to calculate

And entering a sub-loop.

4. The cross-gradient constraint-based gravity-magnetic three-dimensional joint inversion method according to claim 3, wherein the inner loop is obtained by calculation in an outer loop

And calculating to obtain a new model, and iterating for multiple times until the iteration times reach the maximum or the obtained model changes to be stable or the cross gradient function value is increased, wherein in the iteration process, the model and the cross gradient function are converged to obtain a model with a similar structure.