Three-dimensional joint inversion method of seismic electromagnetic based on mixed norm

Abstract

The present application relates to the field of oil and gas exploration, specifically to a mixed norm based three-dimensional joint inversion method of seismic electromagnetic, comprising the following steps: firstly, calculating the electromagnetic response of abnormal body at the observation point based on the magnetotelluric method forward calculation, profiling and numbering the underground abnormal body, and solving by using the forward formula and interpolation basis function; then, establishing a mixed norm joint inversion objective function, integrating the magnetotelluric data, and adding a cross gradient term to obtain the final objective function, and introducing a stratum velocity and resistivity conversion formula constraint; finally, using the LBFGS algorithm for iterative optimization. The present application is based on the L1, L2 and Lp mixed norm inversion system and the joint inversion of seismic and electromagnetic cross gradient, which fully gives play to the advantages of the two kinds of data, improves the inversion resolution, focuses the real position of the abnormal body, improves the imaging effect, reduces the acquisition cost, reduces the data processing machine time, and provides a direction for the application of multi-data acquisition data inversion interpretation.

Description

Technical Field

[0001] This invention relates to the field of oil and gas exploration, specifically to a three-dimensional joint inversion method for seismic electromagnetics based on the hybrid norm. Background Technology

[0002] With the continuous exploration of oil and gas resources, shallow and easily accessible oilfields and deposits are becoming increasingly scarce. Future oil and gas exploration will focus on deep and complex areas, thus placing higher demands and challenges on existing geophysical exploration methods, data acquisition, data processing, forward and inverse modeling, and geological interpretation. Exploration of deep and complex oil and gas reservoirs and mineral resources often requires high-quality imaging of complex carbonate fractured-cavity reservoirs and precise characterization of small geological bodies. The need for reacquiring seismic data is increasing, while simultaneously reducing exploration costs. Therefore, there is an urgent need to utilize both old and new seismic data to maximize the extraction of information.

[0003] 3D seismic exploration offers high resolution and is well-suited for characterizing complex oil and gas reservoirs, but its exploration costs are enormous. 3D electromagnetic exploration, on the other hand, boasts significant advantages in terms of low cost and can accurately characterize the morphology of subsurface anomalies. Existing processing methods generally employ single seismic or electromagnetic data inversion techniques, which suffer from strong ambiguity and produce 3D geological models with relatively low resolution.

[0004] Therefore, in order to meet the requirements of low-cost and high-precision development of oil and gas reservoirs, there is an urgent need for a three-dimensional joint inversion method based on the hybrid norm of seismic electromagnetics. Summary of the Invention

[0005] To avoid the aforementioned problems in existing technologies, the present invention aims to provide a three-dimensional joint inversion method for seismic and electromagnetic data based on a hybrid norm. It adopts a hybrid norm inversion system based on L1, L2, and Lp, and employs a joint inversion method using seismic and electromagnetic cross-gradients to fully leverage the advantages of both types of data, improve the resolution of the inversion, and make the anomalies more focused on their true locations.

[0006] The invention provides the following technical solution: a three-dimensional joint inversion method for seismic electromagnetic inversion based on the hybrid norm, comprising the following steps:

[0007] S1: At the observation point, calculate the electromagnetic response of an anomaly at any location in space using the magnetotelluric method forward modeling.

[0008] S2: Establish a joint inversion objective function based on the hybrid norm. Combine the electromagnetic response data in S1 to obtain a joint magnetotelluric inversion objective function. Add a cross gradient term to obtain the final objective function.

[0009] S3: Use the LBFGS algorithm (Limited Memory Quasi-Newton Method) to iteratively optimize the process in step S2, and obtain the continuously updated joint inversion results.

[0010] The present invention is further configured such that, in step S1, the randomly distributed underground anomalies are divided into M×N×L upright hexahedrons and numbered sequentially.

[0011] The present invention further specifies that the forward modeling formula for the magnetotelluric method is:

[0012]

[0013] Where E is the electric field, ω is the angular frequency, ε is the permittivity, σ is the conductivity, and i is the imaginary unit. ▽ is the Hubble density operator, μ is the permeability, and J... s For the magnetotelluric method, J represents the source current density. s Equal to 0;

[0014] The formula for the interpolation basis function is defined as follows:

[0015]

[0016] in, For edge electric field, Let be the vector interpolation basis function, and e be the element number. By discretizing the electric field using the vector interpolation basis function in the electromagnetic forward modeling formula and using the interpolation basis function as a weighting term, the discrete equation can be obtained as follows:

[0017] CE+iωBE+iωS f =0

[0018] Where C is the stiffness matrix, iωB is the mass matrix, and iωS is the mass matrix. f Given a discrete source term matrix, the electromagnetic response at any location in space is obtained by solving the discrete equations and applying the interpolation basis function formula.

[0019] The present invention is further configured such that step S2 specifically comprises:

[0020] S21: Establish a joint inversion objective function for mixed norms, wherein the mixed norms include L1, L2 and Lp mixed norms;

[0021] S22: The electromagnetic data obtained by the magnetotelluric method in step S1 is incorporated into the hybrid norm joint inversion objective function to obtain the joint magnetotelluric inversion objective function.

[0022] S23: By adding a cross-gradient term to the objective function based on the joint magnetotelluric inversion, the final objective function is obtained:

[0023] Φ=Φ d_UAV +Φ d_MT+λΦ m +λ CG Φ CG

[0024] Its gradient formula is:

[0025]

[0026] Where, Φ d_UAV Let Φ be the objective function for the electromagnetic data of the UAV. m For the model objective function based on the L1, L2, and Lp mixed norm, Φ d_MT Let λ be the magnetotelluric objective function, λ be the regularization factor, and Φ be the regularization factor. CG For the cross gradient term, λ CG Cross gradient term weights.

[0027] The present invention is further configured such that step S21 specifically involves establishing the objective function:

[0028] Φ=Φ d_UAV +λΦ m

[0029] in, The objective function for the data; Φ m The objective function for the model is based on the L1, L2, and Lp mixed norms; For UAV observation data; F is the forward operator; * is the complex conjugate; λ is the regularization factor; W d is the data covariance matrix; m is the model parameter for resistivity; T is the matrix transpose; F(m) UAV For UAV electromagnetic forward modeling operators.

[0030] Construct the objective function Φ of the model based on the L1, L2 and Lp mixed norms. m :

[0031]

[0032] Where, Φ ms Φ mf and Φ mp The objective functions for the smoothness constraint (L2 norm), focus constraint (L1 norm), and Lp norm models are W, respectively. m Here is the model covariance matrix; α∈[0,1] are the weight coefficients between the two; m and m ref These are the model parameters and prior model of resistivity, respectively;

[0033] Transform the model parameters to obtain the gradient of the objective function;

[0034]

[0035] The gradient of the objective function established in step S21 is then:

[0036]

[0037] in, The objective function is the data after parameter transformation. This is the objective function of the model after parameter transformation.

[0038] The present invention is further configured such that, in step S22, the magnetotelluric method is incorporated into the joint inversion to improve the inversion depth, and the objective function for the joint magnetotelluric inversion is:

[0039]

[0040] Where, Φ MT The objective function for obtaining electromagnetic data using the magnetotelluric method is denoted as . W represents the gradient of the objective function for magnetotelluric data. d Let F(m) be the data covariance matrix. MT It is a magnetotelluric forward modeling operator.

[0041] Under parameter transformation, the gradient formula for the objective function based on joint magnetotelluric inversion is obtained as follows:

[0042]

[0043] The present invention is further configured such that step S23 specifically involves, in order to further improve the utilization of prior information, adding a cross-gradient term to the objective function based on the joint magnetotelluric inversion:

[0044] Φ=Φ d_UAV +Φ d_MT +λΦ m +λ CG Φ CG

[0045] Where, λ CG Cross gradient term weights, cross gradient term Φ CG It can be represented as

[0046]

[0047] Where, m v For externally inputted known earthquake inversion information, the formation velocity v is introduced, and the earthquake inversion information m is... v Convert to resistivity;

[0048] The cross gradient terms are obtained in the x, y, and z directions respectively:

[0049]

[0050] v x ,v y ,v z Let be the formation velocities in the x, y, and z directions, respectively, which can be represented by vectors as:

[0051]

[0052] Where k, l, z are the directions of separation based on the central difference of the mesh, respectively.

[0053]

[0054] Right now,

[0055] v x =P x_CG m x

[0056] v y =P y_CG m y

[0057] v z =P z_CG m z

[0058] Based on the above formula, we obtain

[0059] Φ CG =(P x_CG m) T (P x_CG m x )+(P y_CG m y ) T (P y_CG m y )+(P z_CG m z ) T (P z_CG m z )

[0060] In the formula, P x_CG P y_CG P z_CG Let m be the coordinates of observation point P in the x, y, and z directions, respectively. x m y m z These are the model parameters for the resistivity of the observation point P in the x, y, and z directions, respectively.

[0061] The gradient of the final objective function is then given by

[0062]

[0063] The present invention is further configured such that, in step S23, by introducing a conversion formula between formation velocity and resistivity, the seismic inversion information m is... v Converting to resistivity for structural constraints:

[0064]

[0065] Where v is the formation velocity, ρ is the formation resistivity, and ρ w Let be the resistivity of pore water, and k be an empirical constant.

[0066] The present invention is further configured such that step S3 specifically involves using the LBFGS algorithm to approximate the product of the inverse of the Hessian matrix of the function and the gradient, and constructing an approximate inverse of the Hessian matrix by storing historical information from the most recent iterations, thereby realizing the iterative update process.

[0067] In summary, the beneficial effects of the above-mentioned technical solution of the present invention are as follows:

[0068] This invention is based on a hybrid norm inversion system of L1, L2, and Lp, and employs a joint inversion method using seismic and electromagnetic cross-gradients to fully leverage the advantages of both types of data. This significantly improves the inversion resolution, making anomalies more accurately located. Therefore, the fusion of seismic and electromagnetic multi-data processing not only improves the imaging of small geological bodies such as fissures and cavities but also reduces acquisition costs and processing time. It is an economical, effective, and feasible method that points the way for future large-scale multi-data acquisition inversion interpretation and application. Attached Figure Description

[0069] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0070] Figure 1 This is a schematic diagram of the underground three-dimensional grid partitioning.

[0071] Figure 2 This is a comparison chart showing the location of anomalies in the real model, the results of traditional electromagnetic 3D inversion, and the inversion results of this invention.

[0072] Figure 3 This is a flowchart of a three-dimensional joint inversion method for seismic electromagnetics based on the hybrid norm. Detailed Implementation

[0073] To enable those skilled in the art to better understand the technical solutions of the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings. Based on the embodiments of the present invention, other similar embodiments obtained by those skilled in the art without creative effort should all fall within the scope of protection of the present invention.

[0074] The present invention will be further described below with reference to the accompanying drawings and preferred embodiments.

[0075] Example 1:

[0076] like Figures 1-3 As shown in the preferred embodiment of the present invention, the seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm includes the following steps:

[0077] S1: At the observation point, calculate the electromagnetic response of an anomaly at any location in space using the magnetotelluric method forward modeling.

[0078] like Figure 1 As shown, the randomly distributed underground anomalies are divided into M×N×L upright hexahedrons and numbered sequentially.

[0079] The forward modeling formula for magnetotelluric method is:

[0080]

[0081] Where E is the electric field, ω is the angular frequency, ε is the permittivity, σ is the conductivity, i is the imaginary unit, ▽ is the Hamiltonian operator, μ is the permeability, and J is the magnetic field. s For the magnetotelluric method, J represents the source current density. s Equal to 0;

[0082] The formula for the interpolation basis function is defined as follows:

[0083]

[0084] in, For edge electric field, Let be the vector interpolation basis function, and e be the element number. By discretizing the electric field using the vector interpolation basis function in the forward modeling formula of the magnetotelluric method, and using the interpolation basis function as a weighting term, the discrete equation can be obtained as follows:

[0085] CE+iωBE+iωS f =0

[0086] Where C is the stiffness matrix, iωB is the mass matrix, and iωS is the mass matrix. f Given a discrete source term matrix, the electromagnetic response at any location in space is obtained by solving the discrete equations and applying the interpolation basis function formula.

[0087] S2: Establish a joint inversion objective function based on the hybrid norm. Combine the electromagnetic response data in S1 to obtain a joint magnetotelluric inversion objective function. Add a cross gradient term to obtain the final objective function.

[0088] S21: Establish a joint inversion objective function for mixed norms, wherein the mixed norms include L1, L2 and Lp mixed norms;

[0089] Establish the objective function:

[0090] Φ=Φ d_UAV +λΦ m

[0091] in, The objective function for the data; Φ m The objective function for the model is based on the L1, L2, and Lp mixed norms; For UAV observation data; F is the forward operator; * is the complex conjugate; λ is the regularization factor; W d Let be the data covariance matrix; m be the model parameters for resistivity; and T be the matrix transpose; T(m) UAV Forward modeling operator for UAV electromagnetic data.

[0092] Construct the objective function Φ of the model based on the L1, L2 and Lp mixed norms. m :

[0093]

[0094] Where, Φ ms Φ mf and Φ mp The objective functions for the smoothness constraint (L2 norm), focus constraint (L1 norm), and Lp norm models are W, respectively. m Here is the model covariance matrix; α∈[0,1] are the weight coefficients between the two; m and m ref These are the model parameters and prior model of resistivity, respectively;

[0095] Transform the model parameters to obtain the gradient of the objective function;

[0096]

[0097] The gradient of the objective function established in step S21 is then:

[0098]

[0099] in, The objective function is the data after parameter transformation. This is the objective function of the model after parameter transformation.

[0100] S22: The electromagnetic data obtained by the magnetotelluric method in step S1 is incorporated into the hybrid norm joint inversion objective function to obtain the joint magnetotelluric inversion objective function.

[0101] In joint inversion with magnetotelluric methods, incorporating magnetotelluric methods improves the inversion depth. Therefore, the objective function for joint magnetotelluric inversion is:

[0102]

[0103] Where, Φ MT The objective function for obtaining electromagnetic data using the magnetotelluric method is denoted as . Let W be the gradient of the magnetotelluric objective function. d Let F(m) be the data covariance matrix. MT It is a magnetotelluric forward modeling operator.

[0104] Under parameter transformation, the gradient formula for the objective function based on joint magnetotelluric inversion is obtained as follows:

[0105]

[0106] S23: By adding a cross-gradient term to the objective function based on the joint magnetotelluric inversion, the final objective function is obtained:

[0107] To further improve the utilization of prior information, a cross-gradient term is added to the objective function based on the joint magnetotelluric inversion:

[0108] Φ=Φ d_UAV +Φ d_MT +λΦ m +λ CG Φ CG

[0109] Where, λ CG Cross gradient term weights, Φ d_UAV Let Φ be the objective function of the data. m For the model objective function based on the L1, L2, and Lp mixed norm, Φ d_MT Let λ be the magnetotelluric objective function, λ be the regularization factor, and Φ be the regularization factor. CG This is the cross gradient term.

[0110] Cross gradient term Φ CG It can be represented as

[0111]

[0112] Where, m v For externally inputted known earthquake inversion information, the formation velocity v is introduced, and the earthquake inversion information m is... v Convert to resistivity.

[0113] It should be noted that if the formation velocity is used directly, the large difference in magnitude between resistivity and formation velocity will result in an unbalanced inversion effect of the two types of information. Therefore, by introducing a conversion formula between formation velocity and resistivity, the inversion information of the earthquake can be converted into resistivity for structural constraints, which can effectively solve the above problem.

[0114] In geophysical exploration, formation velocity and resistivity are two key parameters. Formation velocity is usually obtained through seismic wave velocity measurements, while resistivity is obtained through electromagnetic methods. Understanding the relationship between these two is of great significance for joint inversion.

[0115] The relationship between formation velocity and resistivity is mainly affected by the following factors:

[0116] 1) Rock type: Different types of rocks have different velocity and resistivity characteristics.

[0117] 2) Porosity and water content: Porosity and water content affect the resistivity of the formation and the propagation speed of seismic waves.

[0118] 3) Mineral composition: Differences in mineral composition can lead to changes in formation velocity and resistivity.

[0119] The Archie formula is one of the most commonly used formulas for describing formation resistivity, and it is used to calculate the resistivity of saturated porous media.

[0120]

[0121] Where: a and m are empirical constants related to rock type and structure, φ is porosity, and S w It represents water saturation, where n is a constant related to saturation, and ρ is the saturation level. w It is the resistivity of pore water.

[0122] Based on extensive experimental and measured data, researchers have proposed several empirical formulas to describe the relationship between formation velocity and resistivity. These formulas are typically extracted from experimental data using methods such as regression analysis.

[0123] The Faust formula is an empirical formula used in seismic exploration to describe the relationship between formation velocity and resistivity. Proposed by Lawrence Y. Faust in 1951, this formula is based on extensive measured data and empirical regression analysis, aiming to provide a simple method for estimating formation velocity. The following is a detailed description of the principle of the Faust formula. The Faust formula describes the relationship between formation velocity and resistivity, and its typical form is:

[0124]

[0125] Where v is the formation velocity, ρ is the formation resistivity, and ρ w Let k be the resistivity of pore water, and k be an empirical constant. The value of the empirical constant k needs to be calibrated using measured data. Different geological conditions and measurement areas may result in different k values. Therefore, in practical applications, specific field calibration is usually required to ensure the accuracy of the formula. (1) Collect measured data: Collect a large amount of formation velocity and resistivity data through field measurements in a specific area. (2) Regression analysis: Use statistical analysis methods to perform regression analysis on the collected data to determine the specific value of k.

[0126] Faust's formula connects the relationship between formation velocity and resistivity based on the fundamental principles of seismic and electromagnetic wave propagation. Its core idea can be understood through the following physical concepts:

[0127] 1. Relationship between formation velocity and porosity: The propagation velocity of seismic waves is closely related to the porosity of rocks. The higher the porosity, the greater the water content in the rock, and the lower the propagation velocity of seismic waves, because the presence of water increases the flexibility of the medium and slows down the seismic wave velocity.

[0128] 2. Relationship between resistivity and porosity: Resistivity is also related to the porosity of rocks. The higher the porosity, the greater the water content in the rock, and the lower the resistivity is usually, because water has high electrical conductivity, which increases the conductivity of the rock.

[0129] 3. Establishment of the empirical formula: Through extensive experimental data, Faust determined the relationship between formation velocity and resistivity using regression analysis, and proposed the aforementioned empirical formula for formation velocity V. This formula assumes the resistivity ρ of pore water... w The resistivity ρ of the formation is relatively stable, while the resistivity ρ of the formation depends on the porosity and water content of the rock.

[0130] The cross gradient terms are obtained in the x, y, and z directions respectively:

[0131]

[0132] v x ,v y ,v z Let be the formation velocities in the x, y, and z directions, respectively, which can be represented by vectors as:

[0133]

[0134] Where k, l, z are the directions of separation based on the central difference of the mesh, respectively.

[0135]

[0136] Right now,

[0137] v x =P x_CG m x

[0138] v y =P y_CG m y

[0139] v z =P z_CG m z

[0140] Based on the above formula, we obtain

[0141] Φ CG =(P x_CG m) T (P x_CG m x )+(P y_CG m y ) T (P y_CG m y )+(P z_CG m z ) T (P z_CG m z )

[0142] In the formula, P x_cG P y_cG P z_CG Let m be the coordinates of observation point P in the x, y, and z directions, respectively. x m y m z These are the model parameters for the resistivity of the observation point P in the x, y, and z directions, respectively.

[0143] The final gradient of the objective function is then:

[0144]

[0145] S3: The process of step S2 is iteratively optimized using the LBFGS algorithm (Limited Memory Quasi-Newton method). The product of the inverse of the Hessian matrix of the function and the gradient is approximated by the LBFGS algorithm. By storing historical information from the most recent iterations, an approximate inverse of the Hessian matrix is ​​constructed to realize the iterative update process.

[0146] This embodiment uses a real model as an example, and applies both traditional electromagnetic three-dimensional inversion and the method described in this invention. The results are as follows: Figure 2 As shown; where Figure 2Column A in the diagram represents the positions of three different anomalies in the real model within three-dimensional space. Column B represents the results of the three anomalies corresponding to column A based on traditional electromagnetic three-dimensional inversion. Column C represents the results of the three anomalies corresponding to column A based on the method of this invention. It can be seen that compared to traditional electromagnetic three-dimensional inversion, this invention significantly improves the inversion resolution, making the anomalies more focused on their true locations.

[0147] The above description is merely a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should also be considered within the scope of protection of the present invention.

Claims

1. A three-dimensional joint inversion method for seismic electromagnetic systems based on the hybrid norm, characterized in that, Includes the following steps: S1: At the observation point, calculate the electromagnetic response of an anomaly at any location in space using the magnetotelluric method forward modeling. S2: Establish a joint inversion objective function based on the hybrid norm. Combine the electromagnetic response data in S1 to obtain a joint magnetotelluric inversion objective function. Add a cross gradient term to obtain the final objective function. S3: Use the LBFGS algorithm to iteratively optimize the process in step S2 to obtain the continuously updated joint inversion results.

2. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 1, characterized in that, In step S1, the randomly distributed underground anomalies are divided into M×N×L upright hexahedrons and numbered sequentially; the coordinates of the observation point P(x, y, z) are consistent with the grid division of the upright hexahedrons.

3. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 2, characterized in that, The forward modeling formula for magnetotelluric method is: Where E is the electric field, ω is the angular frequency, ε is the permittivity, σ is the conductivity, and i is the imaginary unit. Here, μ is the density operator, and J is the permeability. s For the magnetotelluric method, J represents the source current density. s Equal to 0; The formula for the interpolation basis function is defined as follows: in, For edge electric field, Using vector interpolation basis functions, the electric field can be discretized using these basis functions in the electromagnetic forward modeling formula. With the interpolation basis functions as weighting terms, the discrete equations are obtained as follows: CE+iωBE+iωS f =0 Where C is the stiffness matrix, iωB is the mass matrix, and iωS is the mass matrix. f Given a discrete source term matrix, the electromagnetic response at any location in space is obtained by solving the discrete equations and applying the interpolation basis function formula.

4. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 2, characterized in that, Step S2 specifically involves: S21: Establish a joint inversion objective function for mixed norms, wherein the mixed norms include L1, L2 and Lp mixed norms; S22: The electromagnetic response data obtained by the magnetotelluric method in step S1 is incorporated into the hybrid norm joint inversion objective function to obtain the joint magnetotelluric inversion objective function. S23: By adding a cross-gradient term to the objective function based on the joint magnetotelluric inversion, the final objective function is obtained: F=F d_UAV +F d_MT +λΦ m +λ CG F CG Its gradient formula is: Where, Φ d_UAV Let Φ be the objective function for the electromagnetic data of the UAV. m For the model objective function based on the L1, L2, and Lp mixed norm, Φ d_MT Let λ be the magnetotelluric objective function, λ be the regularization factor, and Φ be the regularization factor. CG For the cross gradient term, λ CG The weights of the cross gradient terms, where ▽ is the Habitat density operator. The objective function is the data after parameter transformation. The magnetotelluric objective function after parameter transformation. The objective function of the model after parameter transformation. This is the cross gradient term after parameter transformation.

5. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 4, characterized in that, Step S21 specifically involves establishing the objective function: F=F d_UAV +λΦ m in, The objective function for the data; Φ m The objective function for the model is based on the L1, L2, and Lp mixture norms; For UAV observation data; F is the forward operator; * is the complex conjugate; λ is the regularization factor; W d is the data covariance matrix; m is the model parameter for resistivity; T is the matrix transpose; F(m) UAV For the electromagnetic forward modeling operator of unmanned aerial vehicles; Construct the objective function Φ of the model based on the L1, L2 and Lp mixed norms. m : Where, Φ ms Φ mf and Φ mp The objective functions for the L2 norm, L1 norm, and Lp norm models are respectively, W. m Here is the covariance matrix of the resistivity model; α∈[0,1] are the weighting coefficients; m and m ref These are the model parameters and prior model of resistivity, respectively; Transform the model parameters: The gradient of the objective function established in step S21 is then: in, The objective function is the data after parameter transformation. This is the objective function of the model after parameter transformation.

6. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 5, characterized in that, Step S22 specifically involves integrating the magnetotelluric method into the joint inversion process to improve the inversion depth. Therefore, the objective function for the joint magnetotelluric inversion is: Where, Φ MT The objective function for obtaining electromagnetic data using the magnetotelluric method is denoted as . W represents the gradient of the objective function for magnetotelluric data. d Let F(m) be the data covariance matrix. MT For magnetotelluric forward modeling operators; Under parameter transformation, the gradient formula for the objective function based on joint magnetotelluric inversion is obtained as follows:

7. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 6, characterized in that, In step S23, λ CG Cross gradient term weights, cross gradient term Φ CG It can be represented as: Where, m v For externally inputted known earthquake inversion information, the formation velocity v is introduced, and the earthquake inversion information m is... v Convert to resistivity; The cross gradient terms are obtained in the x, y, and z directions respectively: v x ,v y ,v z Let be the formation velocities in the x, y, and z directions, respectively, which can be represented by vectors as: Where k, l, z are the directions of separation based on the central difference of the mesh, respectively. Right now, v x =P x_CG m x v y =P y_CG m y v z =P z_CG m z therefore, Φ CG =(P x_CG m) T (P x_CG m x )+(P y_CG m y ) T (P y_CG m y )+(P z_CG m z ) T (P z_CG m z ) In the formula, Let m be the coordinates of observation point P in the x, y, and z directions, respectively. x m y m z These are the model parameters for the resistivity of the observation point P in the x, y, and z directions, respectively.

8. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 7, characterized in that, In step S23, by introducing the conversion formula between formation velocity and resistivity, the seismic inversion information m is... v Converting to resistivity for structural constraints: Where v is the formation velocity, ρ is the formation resistivity, and ρ w Let be the resistivity of pore water, and k be an empirical constant.

9. The seismic electromagnetic three-dimensional joint inversion method based on the hybrid norm according to claim 1, characterized in that, Step S3 specifically involves using the LBFGS algorithm to approximate the product of the inverse of the Hessian matrix of the function and the gradient. By storing historical information from the most recent iterations, an approximate inverse of the Hessian matrix is ​​constructed, thus realizing the iterative update process.

Images