Inversion method of mineralization based on seismoelectric effect

Abstract

The invention discloses an inversion method for mineralization based on seismoelectric effect, comprising: (1) constructing an initial model for detecting geological structure; (2) obtaining observed seismic electromagnetic field data; (3) using the Pride equation group to perform forward calculation on the model to obtain predicted seismic electromagnetic field data; when the first iterative calculation is performed, the model is the initial model, and when the other iterative calculations are performed, the model is the updated model; (4) constructing the objective function of one-dimensional regularized inversion of seismoelectric effect according to the observed seismic electromagnetic field data, the predicted seismic electromagnetic field data and the model roughness; (5) taking the objective function to be minimized and its derivative being 0 as the goal, using the Gauss-Newton method to perform iterative inversion calculation of the update amount of the model parameters to obtain the updated model; (6) iteratively repeating steps (3)-(5) until the iteration is terminated to obtain the final inversion model. The method can detect underground mineralization using seismoelectric data.

Description

Technical Field

[0001] The invention belongs to the field of geological exploration, and in particular relates to an inversion method of mineralization based on seismoelectric effect. Background Art

[0002] In natural porous media such as soil and rock, the pore fluid (usually water) usually contains charged ions to form a pore electrolyte solution. The porous strata are generally electrically neutral, but because the solid phase skeleton selectively adsorbs certain ions in the pore solution, there is a net residual charge in the originally neutral solution, forming a double electric layer near the pore wall. Due to the substitution and bond breaking effects, the surface of the solid phase matrix is ​​negatively charged, and some cations in the solution are adsorbed to the surface of the mineral to form an adsorption layer. At the same time, the ions in the solution are redistributed to form a diffusion layer. The adsorption layer and the diffusion layer together constitute a double electric layer. When elastic waves propagate in porous media, the pore fluid moves relative to the solid phase skeleton, and the net residual charged ions in the pore electrolyte solution move with it to form an electric current and an electromagnetic field. Therefore, in porous media, not only are the solid phase motion and fluid motion coupled, but the elastic wave field and the electromagnetic field are also coupled. Frenkel first studied the coupling between seepage movement and electromagnetic field in theory. Based on his work, Biot established a complete set of acoustic theory of fluid-saturated porous media. Pride established the control equations of macroscopic seismic-electric coupling based on the work of Frenkel and Biot. At present, the research on seismic-electric effect is mainly focused on field and laboratory observations, and some scholars have conducted numerical simulation research on it, including different mechanisms such as electrokinetic effect, geodynamo effect, piezoelectric effect, and piezomagnetic effect. (Pride & Haartsen, 1996; Haartsen & Pride, 1997; Hu & Wang, 1999; Garambois & Dietrich, 2002; Guan & Hu, 2008; Gao & Hu, 2010; Hu & Gao, 2011).

[0003] The reservoirs of oil and gas resources are mainly porous media containing fluids. Single seismic exploration cannot distinguish the properties of underground reservoir fluids, and single electrical exploration cannot accurately distinguish the location of underground reservoirs. Seismoelectric exploration combines the advantages of seismic exploration and electrical exploration, and can more accurately analyze the unknown formation parameters of the corresponding formations. Therefore, in geophysical exploration, for porous media containing fluids, seismoelectric exploration methods based on the electrokinetic effect of the double electric layer are sensitive to the parameters of porous media, making them a potential means of detecting underground media, and will have broad application prospects in the fields of oil, gas and water resource exploration.

[0004] However, at present, few scholars have conducted inversion research using seismoelectric signals. Guan et al. (2013) proposed a method to invert seismoelectric logging permeability in fluid-saturated porous formations. Macchioli-Grande et al. (2020) studied the possibility of using Bayesian inversion to infer glacial medium parameters from a set of noisy synthetic data. However, these inversion studies are all based on the assumption of uniform half-space, and there are still certain defects in imaging underground structures. Summary of the invention

[0005] In view of the above, an object of the present invention is to provide a method for inversion of mineralization based on seismoelectric effect, which can realize the detection of underground mineralization parameters using seismoelectric data.

[0006] In order to achieve the above-mentioned invention object, the present invention provides the following technical solutions:

[0007] A method for inverting mineralization based on seismoelectric effect comprises the following steps:

[0008] (1) constructing an initial model based on inversion geology, where the model parameter is mineralization;

[0009] (2) Obtaining observed seismic electromagnetic field data that can be used for inversion;

[0010] (3) The Pride equations are used to perform forward calculations on the model to obtain predicted seismic electromagnetic field data as forward modeling results; when the first iterative calculation is performed, the model is an initial model, and when the other iterative calculations are performed, the model is an updated model;

[0011] (4) Construct the objective function of one-dimensional regularized inversion of seismoelectric effect based on observed seismic electromagnetic field data, predicted seismic electromagnetic field data and model roughness;

[0012] (5) Taking the objective function as the minimum, which is equivalent to its derivative being 0, the Gauss-Newton method is used to iteratively invert and calculate the update amount of the model parameters to obtain the updated model;

[0013] (6) Iteratively repeat steps (3)-(5) until the iteration is terminated, and the latest updated model obtained is used as the final inversion model.

[0014] In one embodiment, when solving the objective function, the sensitivity of the observed seismic electromagnetic field data is calculated, and the specific calculation is:

[0015]

[0016] Among them, m j represents the model parameters of the jth model, Δm j represents the update amount of model parameters, Fs (m j ) represents the forward response data of the jth model, i.e., the predicted seismic electromagnetic field data, s represents the index of the forward response data, J sj represents the objective function of the j-th model constructed with respect to the s-th forward response data;

[0017] The observed seismic electromagnetic field data are selected according to the sensitivity of the observed seismic electromagnetic field data for calculating the update amount of the model parameters.

[0018] In one embodiment, the method of screening the observed seismic electromagnetic field data for calculation of the update amount of the model parameters based on the sensitivity of the observed seismic electromagnetic field data includes: setting the noise error corresponding to the observed seismic electromagnetic field data exceeding the sensitivity threshold to a maximum value so that these data do not play a role in the inversion process.

[0019] In one embodiment, during the iterative inversion calculation process, upper and lower limits of conductivity are imposed during the inversion process to reduce false abnormal data, specifically including:

[0020] The conductivity conversion function based on logarithmic parameters changes the model parameters into the conversion domain, and imposes conductivity upper and lower limits on the inversion results to reduce the multiple solutions of the inversion. That is, the model parameter m is converted into k Transform to the model parameter h in the transformation domain k , and then h k Use the following formula to transform into m k , where a k is the lower limit, the value range is less than 0.001, b k is the upper limit, and its value range is greater than or equal to 1.

[0021] h k =log(m k -a k )-log(b k -m k ),a k <m k <b k ,k=1,2,...,M.

[0022]

[0023] In one embodiment, during the iterative inversion calculation process, the objective function is solved to obtain the fitting difference between the predicted seismic electromagnetic field data and the observed seismic electromagnetic field data, and the fitting difference being less than a set threshold is used as the inversion iteration termination condition. When the inversion iteration termination condition is not met, the model is updated according to the update amount of the model parameters obtained by solving the objective function.

[0024] Compared with the prior art, the present invention has the following beneficial effects:

[0025] In the inversion method of mineralization based on seismoelectric effect provided by an embodiment of the present invention, on the basis of constructing an initial model for detecting geological structure and obtaining observed seismic electromagnetic field data, the objective function of one-dimensional regularized inversion of seismoelectric effect is constructed according to the observed seismic electromagnetic field data, the predicted seismic electromagnetic field data obtained by forward modeling and the model roughness, and the Gauss-Newton method is used to perform iterative inversion to calculate the update amount of model parameters to obtain an updated model. In this way, the final inversion model that conforms to the detected geological structure can be obtained, which greatly improves the efficiency and accuracy of seismoelectric inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0027] Figure 1 is a flow chart of a method for inverting mineralization based on seismoelectric effect provided by an embodiment of the present invention;

[0028] Figure 2 1 is a comparison diagram of the inversion model result provided by the embodiment of the present invention and the real model. The positions of the middle abnormal layers are respectively located at (a) 2-4 km; (b) 4-6 km; (c) 6-8 km; (d) 11-13 km. The blue dotted line is the depth of the earthquake source.

[0029] Figure 3 1 is a comparison diagram of observed seismic electromagnetic field data, predicted seismic electromagnetic field data and residuals provided by an embodiment of the present invention, wherein the positions of the middle abnormal layers are respectively located at (a) 2-4 km; (b) 4-6 km; (c) 6-8 km; (d) 11-13 km, and the residuals are magnified 10 times;

[0030] Figure 4 1 is the covariance matrix of the inversion model result provided by the embodiment of the present invention, and the positions of the intermediate abnormal layers are respectively located at (a) 2-4 km; (b) 4-6 km; (c) 6-8 km; and (d) 11-13 km. DETAILED DESCRIPTION

[0031] To make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific implementation methods described herein are only used to explain the present invention and do not limit the scope of protection of the present invention.

[0032] Figure 1 FIG. 1 is a flow chart of a method for inverting mineralization based on seismoelectric effect provided by an embodiment of the present invention. Figure 1 As shown, the inversion method of mineralization based on seismoelectric effect provided in the embodiment includes the following steps:

[0033] Step 1: Build an initial model for detecting geological structures.

[0034] When constructing the initial model, set the geometric model information, transmission data information and control parameters. The geometric model information includes elastic properties, electrical properties, etc. The specific parameters are shown in Table 1. The transmission data information includes the number and location of the source and receiving points. The control parameters include the maximum number of iterations and the regularization factor λ. The model parameter that needs to be updated during the iteration process is the mineralization.

[0035] Table 1

[0036]

[0037]

[0038] Step 2: Obtain observed seismic electromagnetic field data that can be used for inversion.

[0039] In the embodiment, the observed seismic electromagnetic field data used for inversion may come from existing measured data. When it is measured data, it needs to be pre-processed for denoising before being used.

[0040] Step 3: Use the Pride equations to perform forward modeling on the model to obtain predicted seismic electromagnetic field data as forward modeling results.

[0041] In the embodiment, when the Pride equations are used to perform forward calculation on the model, the Pride equations are expressed as:

[0042]

[0043]

[0044]

[0045]

[0046]

[0047]

[0048] in, represents the Laplace operator, σ(ω) represents the dynamic conductivity, ω represents the circular frequency, u is the solid phase displacement, w=φ(uu f) is the seepage displacement, u f is the average fluid displacement, φ is the porosity, P is the pore fluid pressure, τ is the stress tensor, I is the unit tensor, ρ f is the pore fluid density, ρ is the equivalent density of the porous medium, η is the pore fluid viscosity, κ(ω) is the dynamic permeability, L(ω) represents the electrokinetic coupling coefficient, H, C, M and G are four separate elastic moduli, H is the magnetic field, E is the electric field, f and F are the body force densities acting on the fluid phase and the entire porous medium, respectively, ε and μ are the dielectric constant and magnetic permeability of the porous medium;

[0049] The model parameters and center frequency are introduced into the Pride equations, and the forward modeling results are obtained by calculation and solution. Among them, when calculating the solid phase displacement u, seepage displacement w and average fluid displacement u f The center frequency is required when . The model parameters are all the parameters shown in Table 1.

[0050] Since step 3 is an iterative calculation process, in the first iteration, the model is the initial model, and in other iterative calculations except the first iteration, the model is the updated model.

[0051] Step 4: construct the objective function of one-dimensional regularized inversion of seismoelectric effect based on observed seismic electromagnetic field data, predicted seismic electromagnetic field data and model roughness.

[0052] In the embodiment, an objective function based on the L2 norm is used to regularize the inversion process to optimize the inversion update model. Specifically, the objective function of the regularized inversion is constructed according to the observed seismic electromagnetic field data, the predicted seismic electromagnetic field data obtained by forward modeling, and the model roughness:

[0053]

[0054] in is the data fitting term, is the model roughness constraint term, β is the regularization factor, and the specific forms of u and k are:

[0055] u=W d (d obs -d prd ), (8)

[0056] k=W m (mm ref ). (9)

[0057] In the above formula, d obs is the observed data vector, d prd is the electromagnetic response calculated by model m, m refis the reference model. If the reference model is given initially, the reference model remains unchanged. If it is not given, the first iteration is the initial model, and each subsequent iteration is an updated model. d is a diagonal matrix whose elements are the inverse of the noise in the observed response. and The general form of can be expressed as:

[0058]

[0059] Where x represents u or k in the above formula. The objective function of the nth iteration can be expressed as

[0060]

[0061] and

[0062] u=W d (d obs -d n-1 -Jδm), (12)

[0063] k=W m (m n-1 +δm-m ref ), (13)

[0064] Where, d n-1 is the model m obtained from the previous iteration n-1 The calculated response vector, δm = m n -m n-1 , J is the sensitivity matrix.

[0065] Step 5, with the goal of minimizing the objective function, which is equivalent to its derivative being 0, the Gauss-Newton method is used to iteratively invert and calculate the update amount of the model parameters to obtain an updated model.

[0066] In the embodiment, when solving the objective function of the regularized inversion, the sensitivity matrix (Jacobian matrix) is approximately calculated by the first-order forward difference method, that is, the sensitivity of the observed seismic electromagnetic field data is calculated. The specific calculation is:

[0067]

[0068] Among them, m j represents the model parameters of the jth model, Δm j represents the update amount of model parameters, F s (m j ) represents the forward response data of the jth model, i.e., the predicted seismic electromagnetic field data, s represents the index of the forward response data, J sj Represents the objective function of the j-th model constructed with respect to the s-th forward response data.

[0069] In the embodiment, the sensitivity is used as a screening benchmark for the observed seismic electromagnetic field data, and the observed seismic electromagnetic field data is screened for calculation of the update amount of the model parameters according to the sensitivity of the observed seismic electromagnetic field data. Specifically, the noise error corresponding to the observed seismic electromagnetic field data exceeding the sensitivity threshold is set to be extremely large, so that these observed seismic electromagnetic field data do not play a role in the inversion process.

[0070] In the embodiment, it is also necessary to calculate the update model roughness, that is, when the update amount m of the model parameter of the current iteration is obtained, according to the reference model m 0 , calculate the model roughness constraint term φ according to formulas (7) and (9) m (k) to update the objective function for the next inversion.

[0071] In the embodiment, in order to make the inversion result closer to the actual situation, in the iterative inversion calculation process, upper and lower limits of conductivity are imposed on the inversion process to reduce false abnormal data, specifically including: changing the model parameters to the conversion domain (called logarithmic transformation) based on the conductivity conversion function of the logarithmic parameter, and imposing upper and lower limits of conductivity on the inversion result to reduce the multiple solutions of the inversion.

[0072] That is, use formula (15) to transform the model parameter m k Transform to the model parameter h in the transformation domain k , and then h k Using formula (16) to transform into m k , where a k is the lower limit, the value range is less than 0.001, preferably 1e-5, b k is the upper limit, the value range is greater than or equal to 1, preferably 1,

[0073] h k =log(m k -a k )-log(b k -m k ),a k <m k <b k ,k=1,2,...,M. (15)

[0074]

[0075] Step 6, iteratively repeat steps 3 to 5 until the iteration is terminated, and the latest updated model is used as the final inversion model.

[0076] In the embodiment, each time the iterative calculation is performed, steps 3 to 5 are continued according to the updated model updated in the previous iteration until the iteration is terminated, the inversion termination condition is met or the maximum number of iterations is reached, and the final inversion model that conforms to the detected geological structure is obtained, and the forward response is output at the same time.

[0077] Specifically, in the iterative inversion calculation process, the objective function is solved to obtain the fitting difference between the predicted seismic electromagnetic field data and the observed seismic electromagnetic field data, and the fitting difference being less than the set threshold is used as the inversion iteration termination condition. When the inversion iteration termination condition is not met, the model is updated according to the update amount of the model parameters obtained by solving the objective function.

[0078] In the embodiment, the termination condition of the inversion iteration is that the fitting error is less than 1, wherein the calculation of the fitting error RMS is:

[0079]

[0080] Among them, Ndata is the number of observed seismic electromagnetic field data, s is the index of observed seismic electromagnetic field data, d pre，s and d obs,s is the forward response (i.e. predicted seismic electromagnetic field data) obtained by forward calculation and the observed seismic electromagnetic field data, W d is the data variance matrix.

[0081] The embodiment also provides an inversion example using the above-mentioned inversion method of mineralization based on seismoelectric effect to verify the correctness of the inversion algorithm.

[0082] Figure 2 The inversion model results provided by the embodiment of the present invention are compared with the real model. The positions of the middle abnormal layers are respectively (a) 2-4 km; (b) 4-6 km; (c) 6-8 km; (d) 11-13 km. The horizontal dotted line is the depth of the earthquake source; the thickness of the abnormal layer is 2 km, the depth of the earthquake source is 10 km, and the earthquake source is a double-couple point source M xz =M zx =M 0 =4.37×10 17 N·m, equivalent to an earthquake with a moment magnitude of 5.7, the receiving point is located at (100km, 5km, 1m), and the source time function is:

[0083]

[0084] Among them, the center frequency f 0 The pulse duration is 4s and the frequency range is 0-1Hz.

[0085] Figure 3The following are comparison diagrams of observed seismic electromagnetic field data, predicted seismic electromagnetic field data and residuals provided by the embodiment of the present invention. The positions of the middle abnormal layers are respectively located at (a) 2-4 km; (b) 4-6 km; (c) 6-8 km; (d) 11-13 km, and the residuals are magnified 10 times. Figure 2 and 3 It can be seen that the inversion model is in good agreement with the true model, and the response curve can also be well fitted. Figure 2 There is a certain gap in (d) because the seismic electromagnetic field data is mainly sensitive to the formation parameters above the focal depth. All these can show that the mineralization inversion algorithm based on the seismoelectric effect of the present invention is reliable and has certain practicality. Figure 4 is the covariance matrix of the inversion model result provided by the embodiment of the present invention. Figure 4 It can be seen that the inversion is stably convergent and the gap between the inversion result and the true model is small, which further proves the correctness of the inversion algorithm.

[0086] The specific implementation methods described above provide a detailed description of the technical solutions and beneficial effects of the present invention. It should be understood that the above is only the most preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, supplements and equivalent substitutions made within the scope of the principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for inversion of mineralization based on seismoelectric effect, characterized in that: The following steps are involved: (1) constructing an initial model for detecting geological structures, wherein the model parameter is the mineralization; (2) Obtaining observed seismic electromagnetic field data that can be used for inversion; (3) The Pride equations are used to perform forward calculations on the model to obtain predicted seismic electromagnetic field data as forward modeling results; when the first iterative calculation is performed, the model is an initial model, and when the other iterative calculations are performed, the model is an updated model; (4) Construct the objective function of one-dimensional regularized inversion of seismoelectric effect based on observed seismic electromagnetic field data, predicted seismic electromagnetic field data and model roughness; (5) Taking the objective function as the minimum, which is equivalent to its derivative being 0, the Gauss-Newton method is used to iteratively invert and calculate the update amount of the model parameters to obtain the updated model; (6) Iteratively repeat steps (3)-(5) until the iteration is terminated, and the latest updated model obtained is used as the final inversion model.

2. The inversion method of mineralization based on seismoelectric effect according to claim 1 is characterized in that: When the Pride equations are used to perform forward calculations on the model, the Pride equations are expressed as: in, represents the Laplace operator, σ(ω) represents the dynamic conductivity, ω represents the circular frequency, u is the solid phase displacement, w=φ(uu f ) is the seepage displacement, u f is the average fluid displacement, φ is the porosity, P is the pore fluid pressure, τ is the stress tensor, I is the unit tensor, ρ f is the pore fluid density, ρ is the equivalent density of the porous medium, η is the pore fluid viscosity, κ(ω) is the dynamic permeability, L(ω) represents the electrokinetic coupling coefficient, H, C, M and G are four separate elastic moduli, H is the magnetic field, E is the electric field, f and F are the body force densities acting on the fluid phase and the entire porous medium, respectively, ε and μ are the dielectric constant and magnetic permeability of the porous medium; The model parameters and center frequency are introduced into the Pride equations, and the forward modeling results are obtained by calculation and solution. Among them, when calculating the solid phase displacement u, seepage displacement w and average fluid displacement u f The center frequency is required.

3. The inversion method of mineralization based on seismoelectric effect according to claim 1 is characterized in that: The objective function of the constructed one-dimensional regularized inversion of seismoelectric effect is expressed as: in, is the data fitting term, is the model roughness constraint term, β is the regularization factor, and the specific forms of u and k are: u=W d (d obs -d prd ), (8) k=W m (m-m ref ). (9) In the above formula, d obs is the observed seismic electromagnetic field data vector, d prd is the electromagnetic response calculated by model m, that is, the predicted seismic electromagnetic field data vector, m ref is the reference model, W d is a diagonal matrix whose elements are the inverse of the noise in the observed response, and The general form of can be expressed as: Where x represents u or k in the above formula, and the objective function of the nth iteration can be expressed as and u=W d (d obs -d n-1 -Jδm), (12) k=W m (m n-1 +δm-m ref ), (13) Where, d n-1 is the model m obtained from the previous iteration n-1 The calculated response vector, δm = m n -m n-1 , J is the sensitivity matrix.

4. The inversion method of mineralization based on seismoelectric effect according to claim 1 is characterized in that: When solving the objective function, the sensitivity of the observed seismic electromagnetic field data is calculated. The specific calculation is: Among them, m j represents the model parameters of the jth model, Δm j represents the update amount of model parameters, F s (m j ) represents the forward response data of the jth model, i.e., the predicted seismic electromagnetic field data, s represents the index of the forward response data, J sj represents the objective function of the j-th model constructed with respect to the s-th forward response data; The observed seismic electromagnetic field data are selected according to the sensitivity of the observed seismic electromagnetic field data for calculating the update amount of the model parameters.

5. The inversion method of mineralization based on seismoelectric effect according to claim 4 is characterized in that: The method of screening the observed seismic electromagnetic field data according to the sensitivity of the observed seismic electromagnetic field data for calculating the update amount of the model parameters includes: The noise errors corresponding to the observed seismic electromagnetic field data exceeding the sensitivity threshold are set to be extremely large so that these data do not play a role in the inversion process.

6. The inversion method of mineralization based on seismoelectric effect according to claim 1 is characterized in that: During the iterative inversion calculation process, upper and lower limits of conductivity are imposed to reduce false abnormal data, including: The conductivity conversion function based on logarithmic parameters changes the model parameters into the conversion domain, and imposes conductivity upper and lower limits on the inversion results to reduce the multiple solutions of the inversion. That is, the model parameter m is converted into k Transform to the model parameter h in the transformation domain k , and then h k Using formula (16) to transform into m k , where a k is the lower limit, the value range is less than 0.001, b k is the upper limit, and its value range is greater than or equal to 1. h k =log(m k -a k )-log(b k -m k ),a k <m k <b k ,k=1,2,...,M. (15) 7. The inversion method of mineralization based on seismoelectric effect according to claim 1 is characterized in that: During the iterative inversion calculation process, the objective function is solved to obtain the fitting difference between the predicted seismic electromagnetic field data and the observed seismic electromagnetic field data. The fitting difference being less than the set threshold is used as the inversion iteration termination condition. When the inversion iteration termination condition is not met, the model is updated according to the update amount of the model parameters obtained by solving the objective function.

8. The inversion method of mineralization based on seismoelectric effect according to claim 7 is characterized in that: The inversion iteration termination condition is that the fitting error is less than 1, wherein the calculation of the fitting error RMS is: Among them, Ndata is the number of observed seismic electromagnetic field data, s is the index of observed seismic electromagnetic field data, d pre，s and d obs,s It is the predicted seismic electromagnetic field data and the observed seismic electromagnetic field data, W d is the data variance matrix.

Images