CN113917560B - Three-dimensional heavy magnetic electric shock multi-parameter collaborative inversion method - Google Patents

Abstract

The invention discloses a three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method, comprising: S1, acquiring seismic and electromagnetic data of an area to be measured; S2, performing magnetotelluric two-dimensional and three-dimensional inversion and apparent density respectively according to the seismic and electromagnetic data 3D quantitative inversion of apparent magnetization, gravity, magnetoelectric multi-parameter co-imaging, and gravity interface joint inversion; S3, based on magnetotelluric 2D and 3D inversion, apparent density and apparent magnetization 3D quantitative inversion and gravity, magnetoelectric multi-parameter The result of collaborative imaging is to carry out joint collaborative imaging of gravity, magnetic, electric and seismic physical parameters; S4, combined with the joint collaborative imaging results of gravity, magnetic, and electric seismic physical parameters and the joint inversion results of gravity interface, to construct imaging results for comprehensive geological interpretation, and realize three-dimensional Multi-parameter collaborative inversion of gravity, magneto, electric and seismic.

Description

A 3D Gravity, Magnetism, Electric Seismic Multi-parameter Cooperative Inversion Method

technical field

The invention belongs to the technical field of geological structure imaging, and in particular relates to a three-dimensional gravity, magnetism, electric-seismic multi-parameter collaborative inversion method.

Background technique

Seismic exploration has always been the main means of oil and gas exploration and development, mainly because: ①Seismic imaging technology can be used to obtain imaging of complex underground geological structures, especially reverse time migration (RTM) and full waveform inversion (FWI) can provide higher resolution (2) Seismic reservoir inversion technology can be used to extract reservoir physical parameter information, such as porosity, fluid saturation, permeability, etc., and then predict reservoir fluid type and evaluate reservoir oil and gas. However, under complex geological conditions, such as igneous rock areas, carbonate rock reef development areas, and overthrust fault zones, the serious seismic wave scattering effect and shielding effect lead to weak illumination energy of the underlying strata, and only seismic data can be used to infer the underground structure Multiple solutions exist. In terms of reservoir description and monitoring, when the velocity difference is not large, such as formation water and oil, and reservoirs with different gas saturations, the use of seismic data for reservoir inversion has a certain degree of ambiguity. Because the change of oil (gas) saturation has little effect on seismic wave velocity, seismic data alone cannot predict reservoir fluid type and fluid saturation well.

The magnetotelluric (MT) method is particularly important for solving oil and gas exploration problems under complex geological conditions (such as igneous rock underlying formation imaging, sub-salt imaging, etc.), mainly because: ① electromagnetic signals attenuate when passing through high-resistivity bodies (such as igneous rocks) Very weak, ②Electromagnetic signal is sensitive to low-resistivity sedimentary strata. Therefore, the collected electromagnetic signals contain resistivity information of the underlying sedimentary formation. However, since the MT method is mainly based on low-frequency induced electromagnetic fields, it cannot image relatively thin oil and gas reservoirs well. The basic equations controlling the propagation of seismic waves and electromagnetic waves are the wave equation and the diffusion equation respectively. Although the propagation of these two waves has attenuation and scattering properties, the degrees of attenuation and scattering are different, that is, the attenuation and scattering of seismic waves are relatively weak. .

Contents of the invention

In view of the above-mentioned shortcomings in the prior art, the three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method provided by the present invention solves the problem of multiple solutions and low resolution limitations of a single method in the existing inversion method.

In order to achieve the purpose of the above invention, the technical solution adopted in the present invention is: a three-dimensional gravity, magneto, electric and seismic multi-parameter collaborative inversion method, comprising the following steps:

S1. Obtain seismic and electromagnetic data of the area to be measured;

S2. Carry out 2D and 3D magnetotelluric inversion, 3D quantitative inversion of apparent density and apparent magnetization, multi-parameter co-imaging of gravity, magnetoelectricity and joint inversion of gravity interface according to seismic and electromagnetic data;

S3. According to the results of 2D and 3D inversion of magnetotellurics, 3D quantitative inversion of apparent density and apparent magnetization, and multi-parameter collaborative imaging of gravity, magnetoelectricity, joint collaborative imaging of gravity, magnetic, electric and seismic physical parameters;

S4. Combining the combined imaging results of gravity, magnetic, electric and seismic physical properties and the joint inversion results of gravity and seismic interfaces, constructing imaging results for comprehensive geological interpretation, and realizing three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion.

Further, in the step S2, the two-dimensional and three-dimensional magnetotelluric inversion includes the following sub-steps:

A1. According to seismic and electromagnetic data, carry out 2D and 3D inversion of magnetotelluric coordinates and dipon inversion in spherical cylindrical coordinates;

A2. According to seismic and electromagnetic data, carry out structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling;

A3. The results of 2D and 3D inversion of magnetotelluric coordinates in spherical coordinates and the results of 2D and 3D magnetotelluric inversion and modeling with structured adaptive finite elements are used as the results of 2D and 3D inversion of magnetotellurics;

Wherein, step A1 includes the following sub-steps:

A11. Based on the acquired electromagnetic data, construct the integral form of the Maxwell equation that the electric field intensity vector and the magnetic field intensity vector satisfy:

∮H·dl=∫∫J·dS=∫∫σE·ds

∮E·dl=∫∫iωμ ₀ ·HdS

In the formula, E and H are the electric field intensity vector and the magnetic field intensity vector respectively, J is the current density, dl is the contour of the closed integral, dS is the area covered by the contour, σ is the medium conductivity, ω is the circular frequency, μ ₀ is the vacuum permeability;

A12. In the spherical coordinate system, along the r, θ and φ directions, use several parallel spherical surfaces to divide the spherical coordinate system into several grid units at different intervals, so as to realize the division of the spherical coordinate system into a staggered grid discretization;

A13. Discretize the integral form of the Maxwell equation in the discretized spherical coordinate system, and then determine the staggered grid finite difference equation, and use it as the result of the two-dimensional and three-dimensional magnetotelluric inversion and dipon inversion in the spherical coordinate system;

Step A2 includes the following sub-steps:

A21. For seismic and electromagnetic data, the magnetotelluric adaptive vector finite element method, the linear equation system solution method based on the auxiliary space preconditioning algorithm, the parallel method based on the grid partition technology, and the adaptive finite element forward modeling method are used in sequence. It is processed by magnetotelluric three-dimensional inversion method;

A22. The results obtained from the 3D magnetotelluric inversion method based on adaptive finite element forward modeling are used as the results of structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling;

Wherein, in the step A21, the magnetotelluric three-dimensional inversion method based on adaptive finite element forward modeling specifically includes the following sub-steps:

T1. According to the data of the area to be measured obtained by the parallel method based on grid partitioning technology, and the area to be inverted determined according to it;

T2. Discretize the area to be inverted, and use the discretized grid as the inversion grid;

Among them, the inversion grid does not change during the inversion process;

T3. In each iteration of the inversion process, set a forward modeling grid for each frequency;

Among them, the initial value of the forward modeling grid is the same as that of the inversion grid;

T4. In the iteration of the inversion process, use the posterior error estimate to perform several times of adaptive grid refinement on the forward modeling grid, and calculate the partial derivative on the final forward modeling grid;

T5. Map the partial derivative calculated on the final forward modeling grid back to the corresponding inversion grid to realize the three-dimensional magnetotelluric inversion.

Further, in the step S2, the three-dimensional quantitative inversion of apparent density and apparent magnetization includes the following sub-steps:

B1. Use the regularized descending technique to obtain the gravity and magnetic anomaly data in the whole space of the area to be tested;

B2. Based on the drilling constraint technology, change point constraints or line constraints into volume constraints, and use probabilistic imaging technology to search for the grid interval where the gravity and magnetic anomaly data volume exists, and realize the dimensionality reduction of the solution space;

B3. For the inversion objective function of the solution space after dimensionality reduction, a multi-method inversion system is used to perform three-dimensional quantitative inversion of apparent density and apparent magnetization;

Among them, the methods in the multi-method inversion system include the golden section inversion method based on correlation search, the 3D gravity inversion method based on fast proximal objective function optimization, and the petrophysical relationship constraint inversion method based on machine learning and multivariate geostatistics. method.

Further, the constrained inversion method of petrophysical relationships based on machine learning and multivariate geostatistics in the step B3 includes the following sub-steps:

R1. With the logging data and seismic frame as constraints, establish a cross variogram reflecting the spatial statistical petrophysical relationship;

R2. Cluster the velocity and density distributions obtained by separate seismic inversion and gravity inversion respectively, and reassign the velocity and density values to the clustering results according to the statistical results of rock sample data and well logging data, and obtain new The speed and density model of ;

R3. Using the cross variogram function as the prior function of the inversion objective function, based on the new velocity and density, carry out the gravity inversion constrained by the petrophysical relationship.

Further, in the step S2, performing gravity, magnetoelectric multi-parameter collaborative imaging includes resistivity imaging and multi-parameter imaging;

Among them, the method of resistivity imaging is as follows: using the cross-gradient operator to jointly invert the three-dimensional controllable electromagnetic data and seismic data to obtain resistivity imaging;

The method of multi-parameter imaging is specifically as follows:

Y1. Based on the seismic wave equation and Maxwell equation, the numerical simulation of seismic wave field and electromagnetic wave field is carried out by using the staggered grid finite difference method;

Y2. Based on the data fitting items of electromagnetic data and seismic data obtained by numerical simulation, model regularization items and cross-gradient constraints, construct a joint inversion objective function based on cross-gradient electromagnetic and seismic correspondence;

Y3, using the linear conjugate gradient algorithm to solve the constructed joint inversion objective function;

Y4. Based on the joint inversion objective function solving process, joint electromagnetic and seismic inversion is performed to realize multi-parameter imaging.

Further, in the step Y4, when performing joint electromagnetic and seismic inversion, the methods adopted are as follows: using data fitting based on two adjacent iterative processes, adopting an adaptive regularization method to adjust the regularization factor, and adjusting the regularization factor based on K Means clustering and regression analysis carry out adaptive correction, cross gradient weighting, frequency division multi-scale inversion and multi-grid division and set physical constraint range inversion.

Further, in the step S2, the method for performing the joint inversion of the severe earthquake interface is specifically as follows:

The stochastic inversion method and the cross-gradient inversion method combined with gravity and seismic data are used to invert the subsurface density interface.

Further, in the step S3, the method for performing joint collaborative imaging of gravity, magnetic, electric and seismic physical property parameters is specifically:

Based on the comprehensive response characteristics of gravity, magnetism, electricity and earthquake in the joint collaborative imaging results of gravity, magnetic, electric and seismic physical property parameters and the joint inversion result of gravity interface, the gravity, magnetic, electric and seismic modeling-inversion is carried out, and the inversion results are different. Physical parameters of the geophysical field, and use the physical parameters of different geophysical fields such as gravity, magnetism, electricity and seismic to constrain each other and invert collaborative imaging.

The beneficial effects of the present invention are:

(1) The invention innovatively constructs the collaborative inversion and imaging technology of heavy earthquake, electric shock and gravity, magnetic and electric shock based on multi-parameter modeling and constraints such as density, velocity, resistivity and magnetic susceptibility;

(2) The method of the present invention is applied in the exploration of oil and gas resources and mineral resources, and has achieved outstanding application effects. Based on the internal connection of rock physics, new algorithms and simulation techniques have been introduced to form a joint inversion technology of heavy earthquake electric interface and The multi-parameter collaborative simulation imaging technology of gravity, magnetism, electric and seismic has realized high-precision imaging of various physical parameters (density, resistivity, speed);

(3) Using the method of machine learning, the multi-parameter joint collaborative simulation imaging technology such as density, magnetization, and resistivity has been developed, realizing the leap from a single method to a multi-method and multi-dimensional collaborative interpretation technology. An innovative technology for extracting and enhancing magnetoelectric weak information combination.

Description of drawings

Fig. 1 is a flow chart of the three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method provided by the present invention.

Detailed ways

The specific embodiments of the present invention are described below so that those skilled in the art can understand the present invention, but it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art, as long as various changes Within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are included in the protection list.

As shown in Figure 1, a three-dimensional gravity, magneto, electric and seismic multi-parameter collaborative inversion method includes the following steps:

S1. Obtain seismic and electromagnetic data of the area to be measured;

S2. Carry out 2D and 3D magnetotelluric inversion, 3D quantitative inversion of apparent density and apparent magnetization, multi-parameter co-imaging of gravity, magnetoelectricity and joint inversion of gravity interface according to seismic and electromagnetic data;

S3. According to the results of 2D and 3D inversion of magnetotellurics, 3D quantitative inversion of apparent density and apparent magnetization, and multi-parameter collaborative imaging of gravity, magnetoelectricity, joint collaborative imaging of gravity, magnetic, electric and seismic physical parameters;

S4. Combining the combined imaging results of gravity, magnetic, electric and seismic physical properties and the joint inversion results of gravity and seismic interfaces, constructing imaging results for comprehensive geological interpretation, and realizing three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion.

In step S2 of this embodiment, performing 2D and 3D inversion of magnetotelluric elements includes the following sub-steps:

A1. According to seismic and electromagnetic data, carry out 2D and 3D inversion of magnetotelluric coordinates and dipon inversion in spherical cylindrical coordinates;

A2. According to seismic and electromagnetic data, carry out structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling;

A3. The results of 2D and 3D inversion of magnetotelluric coordinates in spherical coordinates and the results of 2D and 3D magnetotelluric inversion and modeling with structured adaptive finite elements are used as the results of 2D and 3D inversion of magnetotellurics;

Wherein, step A1 includes the following sub-steps:

A11. Based on the acquired electromagnetic data, construct the integral form of the Maxwell equation that the electric field intensity vector and the magnetic field intensity vector satisfy:

∮H·dl=∫∫J·dS=∫∫σE·ds

∮E·dl=∫∫iωμ ₀ ·HdS

In the formula, E and H are the electric field intensity vector and the magnetic field intensity vector respectively, J is the current density, dl is the contour of the closed integral, dS is the area covered by the contour, σ is the medium conductivity, ω is the circular frequency, μ ₀ is the vacuum permeability;

A12. In the spherical coordinate system, along the r, θ and φ directions, use several parallel spherical surfaces to divide the spherical coordinate system into several grid units at different intervals, so as to realize the division of the spherical coordinate system into a staggered grid discretization;

A13. Discretize the integral form of the Maxwell equation in the discretized spherical coordinate system, and then determine the staggered grid finite difference equation, and use it as the result of the two-dimensional and three-dimensional magnetotelluric inversion and dipon inversion in the spherical coordinate system;

Specifically, in the frequency range of magnetotelluric research (about 10 ⁵ ～10 ^-5 Hz), since the frequency is relatively lower and satisfies the quasi-steady electromagnetic field, the displacement current can be ignored, and the magnetic permeability of the underground medium is taken as the vacuum permeability, Construct the corresponding integral form of the Maxwell equation; at the same time, because the discretization of the differential form of the Max equation system in the spherical coordinate system is relatively complicated, the discretization of the integral form of the Max equation system results in a linear equation, so a grid section is required before the discretization point.

The above step A2 includes the following sub-steps:

A21. For seismic and electromagnetic data, the magnetotelluric adaptive vector finite element method, the linear equation system solution method based on the auxiliary space preconditioning algorithm, the parallel method based on the grid partition technology, and the adaptive finite element forward modeling method are used in sequence. It is processed by magnetotelluric three-dimensional inversion method;

A22. The results obtained from the 3D magnetotelluric inversion method based on adaptive finite element forward modeling are used as the results of structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling;

Specifically, in step A21, in the 3D magnetotelluric inversion method based on adaptive finite element forward modeling, the model established by 3D forward modeling is used, and the model established by forward inversion is firstly passed through the magnetotelluric method adaptive vector finite element method Carry out adaptive refinement on the grid. Adaptive grid refinement can automatically refine the area with large error according to the posterior error estimation. Compared with the global grid refinement, the error decrease speed is significantly faster than the global grid refinement; The linear equation solution method of the spatial preconditioning algorithm solves the objective function of the inversion; the parallel method of the grid partitioning technology distributes computing tasks to multiple processes according to the grid partitioning, and each process uses MPI to communicate during the computing process to keep in sync. The forward calculation algorithm is used to realize the inversion, and the grid separation of the forward and inversion is carried out, which solves the problem of inversion accuracy.

Based on this, in the above step A21, the 3D magnetotelluric inversion method based on adaptive finite element forward modeling specifically includes the following sub-steps:

T1. According to the data of the area to be measured obtained by the parallel method based on grid partitioning technology, and the area to be inverted determined according to it;

T2. Discretize the area to be inverted, and use the discretized grid as the inversion grid;

Among them, the inversion grid does not change during the inversion process;

T3. In each iteration of the inversion process, set a forward modeling grid for each frequency;

Among them, the initial value of the forward modeling grid is the same as that of the inversion grid;

T4. In the iteration of the inversion process, use the posterior error estimate to perform several times of adaptive grid refinement on the forward modeling grid, and calculate the partial derivative on the final forward modeling grid;

T5. Map the partial derivative calculated on the final forward modeling grid back to the corresponding inversion grid to realize the three-dimensional magnetotelluric inversion.

In step S2 of this embodiment, the three-dimensional quantitative inversion of apparent density and apparent magnetization includes the following sub-steps:

B1. Use the regularized descending technique to obtain the gravity and magnetic anomaly data in the whole space of the area to be tested;

In this step, the downward continuation of the potential field is a typical ill-posed problem. Due to errors in actual observations and some interference information, high-frequency oscillations will occur when the potential field is extended down to the vicinity of the field source, often covering up useful information . The regularization method is one of the most stable algorithms, which can solve the instability in the extension;

B2. Based on the drilling constraint technology, change point constraints or line constraints into volume constraints, and use probabilistic imaging technology to search for the grid interval where the gravity and magnetic anomaly data volume exists, and realize the dimensionality reduction of the solution space;

In this step, according to the known prior information such as drilling, the plane control constraints composed of these plane control points can be established, and the corresponding physical parameter constraint arrays can be established from the control points to form the constraint conditions, so as to realize the control of the physical property inversion model. Constrained modeling and constrained inversion. The density distribution of geological bodies has a certain continuity and is correlated with the spatial distribution of gravity anomalies, and drilling data can be combined with gravity anomalies and regularized normalized extension results to carry out well constraints and constraint expansion. Considering that on the basis of knowing a certain drilling constraint, it is judged whether the known drilling constraint has encountered a geological source body. If it has encountered a geological source body, combined with the distribution of residual gravity anomalies and the regularized descending results, the quantitative interval constraint range is given, and the constraint To expand, a point constrained by the well is expanded to a body, and the grid unit that the well passes through is a fixed value, and other expanded areas are further iteratively modified with the inversion iteration process to achieve dimensionality reduction in the solution space.

B3. For the inversion objective function of the solution space after dimensionality reduction, a multi-method inversion system is used to perform three-dimensional quantitative inversion of apparent density and apparent magnetization;

The methods in the above multi-method inversion system include the golden section inversion method based on correlation search, the 3D gravity inversion method based on fast proximal objective function optimization, and the petrophysical relationship constrained inversion method based on machine learning and multivariate geostatistics. ; The above inversion method overcomes the difficult problem of inversion multi-solution, and improves the reliability of apparent density and apparent magnetization inversion.

Among them, the golden section method is a classic algorithm in optimization calculation, which has the characteristics of simple structure and high global optimization accuracy. The basic idea of this method is to gradually narrow the search range based on the principle of symmetry and proportional contraction.

The mathematical theory basis of the 3D gravity inversion method based on the fast proximal objective function optimization is the Lipschitz continuum in the mathematical analysis.

The constrained inversion method of petrophysical relations based on machine learning and multivariate geostatistics uses machine learning and multivariate geostatistics to form an algorithm process, which can better reflect the spatial statistical petrophysical relations; based on this, in the above step R3 Constrained inversion method of petrophysical relationships based on machine learning and multivariate geostatistics

Include the following sub-steps:

R1. With the logging data and seismic frame as constraints, establish a cross variogram reflecting the spatial statistical petrophysical relationship;

R2. Cluster the velocity and density distributions obtained by separate seismic inversion and gravity inversion respectively, and reassign the velocity and density values to the clustering results according to the statistical results of rock sample data and well logging data, and obtain new The speed and density model of ;

R3. Using the cross variogram function as the prior function of the inversion objective function, based on the new velocity and density, carry out the gravity inversion constrained by the petrophysical relationship.

Compared with the traditional inversion method, the model obtained by the above inversion method can better fit the gravity field and petrophysical relationship, and is more conducive to geological interpretation.

The gravity, magnetoelectric multi-parameter collaborative imaging in step S2 of this embodiment includes resistivity imaging and multi-parameter imaging;

Among them, the method of resistivity imaging is as follows: using the cross-gradient operator to jointly invert the three-dimensional controllable electromagnetic data and seismic data to obtain resistivity imaging; the three-dimensional controllable electromagnetic data and seismic data contain complementary information , the joint inversion of controllable electromagnetic and seismic data can obtain more reliable imaging of subsurface structures.

The method of multi-parameter imaging is specifically as follows:

Y1. Based on the seismic wave equation and Maxwell equation, the numerical simulation of seismic wave field and electromagnetic wave field is carried out by using the staggered grid finite difference method;

Y2. Based on the data fitting items of electromagnetic data and seismic data obtained by numerical simulation, model regularization items and cross-gradient constraints, construct a joint inversion objective function based on cross-gradient electromagnetic and seismic correspondence;

Y3, using the linear conjugate gradient algorithm to solve the constructed joint inversion objective function;

In this step, the conjugate gradient method (Conjugate Gradient) is a method between the steepest descent method and the Newton method. It only needs to use the first-order derivative information, but it overcomes the shortcomings of the steepest descent method. In order to overcome the shortcomings of Newton's method that needs to store and calculate the Hesse matrix and find its inversion, the conjugate gradient method is not only one of the most useful methods for solving large linear equations, but also one of the most effective algorithms for solving large nonlinear optimization. Among various optimization algorithms, the conjugate gradient method is a very important one. Its advantage is that it requires less memory, has step convergence, high stability, and does not require any external parameters.

Y4. Based on the joint inversion objective function solving process, joint electromagnetic and seismic inversion is performed to realize multi-parameter imaging.

In the above step Y4, when performing joint electromagnetic and seismic inversion, the methods adopted are as follows: using data fitting based on two adjacent iterative processes, using adaptive regularization method to adjust the regularization factor, based on K-means clustering and Regression analysis performs adaptive correction, cross-gradient weighting, multi-scale inversion with frequency division, multi-grid division, and inversion with set physical constraint range.

Specifically, the adaptive regularization method refers to a method of adaptively adjusting the value of the regularization factor according to certain criteria during the inversion process; based on the K-value clustering algorithm to determine the boundaries of geological bodies or formations, the clustering Algorithm is a classic clustering method, which is simple, fast and automatic, and has been widely used in image processing, pattern recognition and other fields. After clustering the resistivity and velocity models using the K-means clustering technique, the resistivity or velocity values were adjusted using regression analysis techniques. Therefore, the complete adaptive correction method is:

(1) Using K-means clustering technology to divide the resistivity and velocity models into different regions and extract the anomaly regions of interest;

(2) In the resistivity and velocity model, corresponding to the anomalous body area of the same geological body, the relationship between the resistivity and velocity in the anomalous body area is obtained by using the regression analysis technique;

(3) For the same geological body in the resistivity and velocity model, use the relationship between resistivity and velocity obtained from the logarithmic regression analysis, and use a physical property with high reliability to correct another physical property, The relative resolution of electromagnetic (resistivity) and seismic (velocity) for a certain geological body depends on the actual situation.

When using the frequency-division multi-scale method, firstly, the controllable source electromagnetic data of all frequencies are divided into three frequency bands from low frequency to high frequency: low frequency band, middle frequency band, and high frequency band, and then use the electromagnetic data of these three frequency bands to The controllable source electromagnetic and seismic joint inversion is carried out at each level, and the electromagnetic data of each scale has a level-by-level inclusion relationship in the joint inversion.

In this example, the specific method of frequency-division multi-scale inversion method is given: first, a resistivity and velocity model is obtained by joint inversion using low-frequency electromagnetic data and seismic data, and then the joint inversion of low-frequency electromagnetic and seismic data Obtain the resistivity and velocity model as the initial model for the joint inversion of low- and mid-frequency electromagnetic data and seismic data. The resistivity and velocity models obtained by inverting low- and medium-frequency electromagnetic and seismic data are used as the initial model for joint inversion of low, medium and high-frequency (ie, full-frequency) electromagnetic and seismic data, until the termination condition is met and the iterative process stops.

The basic idea of the multi-grid inversion scheme is to divide the model area to be inverted into a series of grids of different scales (thickness), and use coarser grids to invert low-frequency controllable source electromagnetic data. With the addition of high-frequency electromagnetic data, the subdivision grid used is getting finer and finer. The key to the multi-grid inversion scheme is the conversion from coarse grid to fine grid. The model parameters of the grid are directly transmitted to the finer grid it contains, as the initial model for inversion on the finer grid. This method has the characteristics of small calculation amount and easy implementation.

When performing cross-gradient weighting: in the joint inversion process, different weight factors are given to the corresponding grids according to the distance to the section where the seismic line is located, so that the resistivity and velocity within a certain distance from the section where the seismic line is located have a strong correlation. In this way, the imaging accuracy of resistivity data volumes in areas without seismic lines in this range can be improved to a certain extent.

In step S2 of this embodiment, the method for performing the joint inversion of the severe-seismic interface is specifically as follows:

The stochastic inversion method and the cross-gradient inversion method combined with gravity and seismic data are used to invert the subsurface density interface. Through the above joint inversion method, the accuracy of interface depth inversion is improved, especially the identification ability of deep structures is enhanced.

In step S3 of this embodiment, the method for performing joint collaborative imaging of gravity, magnetic, electric, and seismic physical property parameters is specifically:

Based on the comprehensive response characteristics of gravity, magnetism, electricity and earthquake in the joint collaborative imaging results of gravity, magnetic, electric and seismic physical property parameters and the joint inversion results of gravity interface, the gravity, magnetic, electric and seismic modeling-inversion is carried out, and the inversion is different. Physical parameters of the geophysical field, and use the physical parameters of different geophysical fields such as gravity, magnetism, electricity and seismic to constrain each other and invert collaborative imaging.

Claims (4)

1. A three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method is characterized in that it comprises the following steps: S1. Obtain seismic and electromagnetic data of the area to be measured; S2. Carry out 2D and 3D magnetotelluric inversion, 3D quantitative inversion of apparent density and apparent magnetization, multi-parameter collaborative imaging of gravity, magnetoelectricity, and joint inversion of gravity interface based on seismic and electromagnetic data; S3. According to the results of 2D and 3D inversion of magnetotellurics, 3D quantitative inversion of apparent density and apparent magnetization, and multi-parameter collaborative imaging of gravity, magnetoelectricity, joint collaborative imaging of gravity, magnetic, electric and seismic physical parameters; S4. Combining the physical parameters of gravity, magnetism, electric shock, collaborative imaging results, and joint inversion results of heavy shock interface, construct imaging results for comprehensive geological interpretation, and realize three-dimensional gravity, magnetism, electric shock and multi-parameter collaborative inversion; In the step S2, the two-dimensional and three-dimensional magnetotelluric inversion includes the following sub-steps: A1. According to seismic and electromagnetic data, carry out 2D and 3D inversion of magnetotelluric coordinates and dipon inversion in spherical cylindrical coordinates; A2. According to seismic and electromagnetic data, carry out structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling; A3. The results of 2D and 3D inversion of magnetotelluric coordinates in spherical coordinates and the results of 2D and 3D magnetotelluric inversion and modeling with structured adaptive finite elements are used as the results of 2D and 3D inversion of magnetotellurics; Wherein, step A1 includes the following sub-steps: A11. Based on the acquired electromagnetic data, construct the integral form of the Maxwell equation that the electric field intensity vector and the magnetic field intensity vector satisfy: ∮H·dl=∫∫J·dS=∫∫σE·ds ∮E·dl=∫∫iωμ ₀ ·HdS In the formula, E and H are the electric field intensity vector and the magnetic field intensity vector respectively, J is the current density, dl is the contour of the closed integral, dS is the area covered by the contour, σ is the medium conductivity, ω is the circular frequency, μ ₀ is the vacuum permeability; A12. In the spherical coordinate system, along the r, θ and φ directions, use several parallel spherical surfaces to divide the spherical coordinate system into several grid units at different intervals, so as to realize the division of the spherical coordinate system into a staggered grid discretization; A13. Discretize the integral form of the Maxwell equation in the discretized spherical coordinate system, and then determine the staggered grid finite difference equation, and use it as the result of the two-dimensional and three-dimensional magnetotelluric inversion and dipon inversion in the spherical coordinate system; Step A2 includes the following sub-steps: A21. For seismic and electromagnetic data, the magnetotelluric adaptive vector finite element method, the linear equation system solution method based on the auxiliary space preconditioning algorithm, the parallel method based on the grid partition technology, and the forward modeling based on the adaptive finite element method are used in sequence. It is processed by magnetotelluric three-dimensional inversion method; A22. The results obtained from the 3D magnetotelluric inversion method based on adaptive finite element forward modeling are used as the results of structured adaptive finite element 2D and 3D magnetotelluric inversion and modeling; Wherein, in the step A21, the magnetotelluric three-dimensional inversion method based on adaptive finite element forward modeling specifically includes the following sub-steps: T1. According to the data of the area to be measured obtained by the parallel method based on grid partitioning technology, and the area to be inverted determined according to it; T2. Discretize the area to be inverted, and use the discretized grid as the inversion grid; Among them, the inversion grid does not change during the inversion process; T3. In each iteration of the inversion process, set a forward modeling grid for each frequency; Among them, the initial value of the forward modeling grid is the same as that of the inversion grid; T4. In the iteration of the inversion process, use the posterior error estimate to perform several times of adaptive grid refinement on the forward modeling grid, and calculate the partial derivative on the final forward modeling grid; T5. Map the partial derivative calculated on the final forward modeling grid back to the corresponding inversion grid to realize the three-dimensional magnetotelluric inversion; In the step S2, the three-dimensional quantitative inversion of apparent density and apparent magnetization includes the following sub-steps: B1. Use the regularized descending technique to obtain the gravity and magnetic anomaly data in the whole space of the area to be tested; B2. Based on the drilling constraint technology, change point constraints or line constraints into volume constraints, and use probabilistic imaging technology to search for the grid interval where the gravity and magnetic anomaly data volume exists, and realize the dimensionality reduction of the solution space; B3. For the inversion objective function of the solution space after dimensionality reduction, a multi-method inversion system is used to perform three-dimensional quantitative inversion of apparent density and apparent magnetization; Among them, the methods in the multi-method inversion system include the golden section inversion method based on correlation search, the 3D gravity inversion method based on fast proximal objective function optimization, and the petrophysical relationship constraint inversion method based on machine learning and multivariate geostatistics. method; In the step S2, performing gravity, magnetoelectric multi-parameter collaborative imaging includes resistivity imaging and multi-parameter imaging; Among them, the method of resistivity imaging is as follows: using the cross-gradient operator to jointly invert the three-dimensional controllable electromagnetic data and seismic data to obtain resistivity imaging; The method of multi-parameter imaging is specifically as follows: Y1. Based on the seismic wave equation and Maxwell equation, the numerical simulation of seismic wave field and electromagnetic wave field is carried out by using the staggered grid finite difference method; Y2. Based on the data fitting items of electromagnetic data and seismic data obtained by numerical simulation, model regularization items and cross-gradient constraints, construct a joint inversion objective function based on cross-gradient electromagnetic and seismic correspondence; Y3, using the linear conjugate gradient algorithm to solve the constructed joint inversion objective function; Y4. Based on the joint inversion objective function solving process, perform joint electromagnetic and seismic inversion to realize multi-parameter imaging; in the step S2, the method for joint inversion of the gravity-seismic interface is specifically: The stochastic inversion method and the cross-gradient inversion method combined with gravity and seismic data are used to invert the subsurface density interface. 2. the three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method according to claim 1, is characterized in that, the petrophysical relationship constrained inversion method based on machine learning and multivariate geostatistics in the step B3 comprises the following subdivisions: step: R1. With the logging data and seismic frame as constraints, establish a cross variogram reflecting the spatial statistical petrophysical relationship; R2. Cluster the velocity and density distributions obtained by separate seismic inversion and gravity inversion respectively, and reassign the velocity and density values to the clustering results according to the statistical results of rock sample data and well logging data, and obtain new The speed and density model of ; R3. Using the cross variogram function as the prior function of the inversion objective function, based on the new velocity and density, carry out the gravity inversion constrained by the petrophysical relationship. 3. The three-dimensional gravity, magnetic, electric and seismic multi-parameter collaborative inversion method according to claim 1, characterized in that, in the step Y4, when carrying out the joint inversion of electromagnetic and seismic, the methods adopted are as follows: using phase-based The data fitting of adjacent two iterations adopts adaptive regularization method to adjust regularization factor, adaptive correction based on K-means clustering and regression analysis, cross-gradient weighting, frequency-division multi-scale inversion and multi-grid division And set physical constraint range inversion. 4. The three-dimensional gravity-magnetism-electric-seismic multi-parameter collaborative inversion method according to claim 1, characterized in that, in the step S3, the method for carrying out the joint collaborative imaging of the gravity-magnetic-electric-seismic physical property parameters is specifically: Based on the comprehensive response characteristics of gravity, magnetism, electricity, and earthquake in the joint collaborative imaging results of gravity, magnetic, electric, and seismic physical property parameters and the joint inversion results of gravity interface, the gravity, magnetic, electrical, and seismic modeling-inversion is carried out, and the inversion results are different. Physical parameters of the geophysical field, and use the physical parameters of different geophysical fields such as gravity, magnetism, electricity and seismic to constrain each other and invert collaborative imaging.

Images