Fast finite-difference solution method and system for seismoelectric wave field under different saturations

Abstract

The application discloses a quick finite difference solution method and system of seismic-electric wave field under different saturations, belongs to the technical field of seismic wave field numerical simulation in geophysical exploration, and comprises the following steps: setting basic parameters of a model, determining water saturation under different depths according to a Van Genuchten model, applying the water saturation to expressions for sound field and electric field calculation, respectively introducing independent time iteration steps, iteratively calculating the numerical values of the sound field and the electric field, extending the sound field calculation result to the electric field time node, obtaining a calculation result, and verifying the calculation result. The application effectively reduces the calculation amount, improves the accuracy and integrity of the seismic-electric wave field calculation, and makes the simulation result closer to the actual physical process.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation technology of seismic wavefield in geophysical exploration, specifically involving a fast finite difference solution method and system for seismoelectric wavefield under different saturation levels. Background Technology

[0002] In geophysical exploration, seismoelectric effects have attracted much attention because they can reflect the properties of the medium. Seismoelectric effects are a kind of electrodynamic coupling phenomenon related to the electrical double layer in fluid-containing porous media. When seismic waves propagate in porous soil or rock, the fluid in the pores will generate a pressure gradient, which will cause charged ions to move and generate electrical signals. Seismoelectric wave fields are generally divided into two categories: (1) current sources accompanying elastic wave propagation—seismic fields; (2) interface electromagnetic waves induced between two layers with different physicochemical properties. Since seismoelectric signals can reflect the elastic and electrical properties of the medium, they can improve the detection of reservoir characteristics, well logging, environmental monitoring and earthquake prediction, and therefore have attracted much attention in geophysical research. However, actual strata are not always saturated with a single fluid, so it is particularly important to carry out seismoelectric research under partially saturated conditions.

[0003] Existing techniques for modeling partially saturated strata assume a uniform saturation distribution, neglecting the actual impact of water cut distribution. Furthermore, the finite difference algorithm used suffers from low computational efficiency due to its fixed time step. Therefore, a new method is urgently needed. Summary of the Invention

[0004] The purpose of this invention is to provide a fast finite difference solution method and system for seismoelectric wave fields under different saturation levels. By setting different time steps for the sound field and electric field, the computational load is effectively reduced. Furthermore, by using linear prediction technology, the current source density calculated at the sound field time point is extended to a denser electric field time point, thus achieving time decoupling between seismic waves and electromagnetic waves during the iteration process. This improves the accuracy and completeness of the seismoelectric wave field calculation, making the simulation results closer to the actual physical process.

[0005] To achieve the above objectives, this invention provides a fast finite difference solution method for seismoelectric wave fields under different saturation levels, comprising the following steps:

[0006] S1. Set the basic parameters of the model;

[0007] S2. Determine the water saturation at different depths based on the Van Genuchten model;

[0008] S3. Apply the water saturation in S2 to the expressions for calculating the sound field and electric field. By introducing independent time iteration steps, iteratively calculate the values ​​of the sound field and electric field.

[0009] S4. Extend the sound field calculation results in S3 to the electric field time node to obtain the calculation results, and verify the calculation results.

[0010] Preferably, S1 includes the porosity of the porous medium, the density of the solid particles, the water density, the air density, the residual water content of the porous medium, and the saturated water content.

[0011] Preferably, in S2, the water saturation at different depths is calculated using the relationship between matrix suction and depth, and the water content formula, and then used in the calculation of the sound field and electric field. The calculation formula is as follows:

[0012] h(z)=ρ w g(z-z0);

[0013]

[0014] Where z is the depth, z0 is the inlet depth, h(z) is the matrix suction at depth z, θ(h) is the water content, and θ r θ represents the residual water content of the porous medium. s S is the saturated water content, α, h, and m are parameters, and S w (h) represents the water saturation level.

[0015] Preferably, in S3, the electro-electro-coupled equations of the fluid-filled porous medium are decomposed into an acoustic field calculation part and an electric field calculation part.

[0016] Preferably, the sound field calculation part includes setting the time step of the sound field to Δt. s The staggered grid finite difference method is used to update the sound pressure and velocity components.

[0017] Preferably, the electric field calculation part includes setting the time step of the electric field to Δt. EM The coupling between the acoustic field and the electric field is achieved using a linear prediction method.

[0018] Preferably, the linear prediction method involves extending the current source density calculated from the sound field to the electric field time node, as shown in the following formula:

[0019]

[0020] Where, j x,z (x,z) represents the current density at position (x,z), where x and z are the abscissa and ordinate in a two-dimensional coordinate system, respectively. k is the sequence number of the sound field time step, and w is the number from 1 to Δt. s / Δt EM integers, (k)Δt s and (k+1)Δt swΔt represents the current time and the time corresponding to the next sound field time step, respectively. EM This represents the electric field offset relative to a certain sound field time step.

[0021] A fast finite difference solution system for seismoelectric wave fields under different saturations, including an initialization module for setting the basic parameters of the model;

[0022] A water saturation distribution calculation module at different formation depths is connected to the initialization module and is used to determine the water saturation at different depths.

[0023] The sound field and electric field calculation module is connected to the water saturation distribution calculation module at the depth of the formation. It is used to apply the water saturation to the expression for sound field and electric field calculation, and to iteratively calculate the values ​​of sound field and electric field by introducing independent time iteration steps respectively.

[0024] The linear prediction module, connected to the sound field and electric field calculation modules, is used to extend the sound field calculation results to the electric field time nodes to obtain the calculation results and verify the calculation results.

[0025] Therefore, the present invention employs the above-mentioned fast finite difference solution method and system for seismoelectric wave fields under different saturations. Compared with the prior art, the present invention has the following significant advantages:

[0026] (1) The present invention sets different time steps for the acoustic field and the electric field respectively, thereby effectively reducing the amount of computation and greatly improving the computational efficiency when dealing with large-scale and complex partially saturated strata models;

[0027] (2) By using linear prediction, this invention extends the current source density calculated at the sound field time point to a more dense electric field time point, realizing the time decoupling of seismic waves and electromagnetic waves in the iterative steps, improving the accuracy and completeness of the seismoelectric wave field calculation, and making the simulation results more consistent with the actual physical process.

[0028] (3) This invention applies water content distribution to unsaturated seismoelectric theory and proposes a nonlinear saturation model, which can more realistically reflect the actual situation of partially saturated strata.

[0029] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0030] Figure 1 The diagram shows the electro-electric effect data simulation of the fast finite difference solution method for the electro-electric wave field under different saturations of this invention, and the system. Figure 1 (a) in the figure shows the simulation diagram of (1) seismic field, (2) interface electromagnetic wave and (3) direct electromagnetic wave of this method; Figure 1(b) in the figure is a comparison between the proposed method and the traditional seismoelectric finite difference method. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. All other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Unless otherwise defined, the technical or scientific terms used in the present invention should have the ordinary meaning understood by those skilled in the art.

[0032] Example 1

[0033] The fast finite difference solution method for the seismoelectric wave field under different saturation levels of the present invention includes the following steps:

[0034] S1. Based on the pore acoustic mechanism, the propagation of seismic waves under the acoustic approximation is described by establishing equations using pressure disturbance and fluid displacement, as follows:

[0035]

[0036] Where P is the sound pressure, u is the displacement, and K is the displacement. u It is the bulk modulus of undrained water, K b It is the bulk modulus of the solid framework, G b It is the shear modulus. K s and K f These are the bulk moduli of the solid matrix and the fluid, respectively. ρ=(1-φ)ρ s +φρ f For the equivalent density, ρ s and ρ f These represent the mass densities of the solid framework and the pore fluid, respectively. F is the force source, φ is the porosity, and T is the time.

[0037] Calculate the source current generated by the flow of pore water relative to the solid phase, based on the effective charge Q per unit pore volume. v Calculate the current density j using Darcy flow w. s electric potential Solve the Maxwell equations under the quasi-static approximation; solve the simultaneous formulas to obtain the relationship between the formula and water saturation, and the calculation method is as follows:

[0038]

[0039] lg Q v =-9.2349-0.8219lgκ0;

[0040]

[0041] ρ f (S w ) = S w ρ w +(1-S w )ρ g ;

[0042]

[0043] Where σ is the conductivity of the medium, K u For the undrained bulk modulus, j s For current density, Let κ be the electric potential, η be the permeability, and η be the electrical potential. f Where is the fluid viscosity, B is the Skempton coefficient, and S is the viscosity of the fluid. w The values ​​represent the water saturation at different depths.

[0044] S2. Based on the Van Genuchten model, determine the volumetric water content under different substrate suction conditions. The calculation formula is as follows:

[0045]

[0046] Where θ(h) is the moisture content, θ r θ represents the residual water content of the porous medium. s Let α, h, and m be the saturated water content, and α, h, and m be parameters. Combining this with the relationship between matric suction and depth in aquifer-bearing strata, the relationship is expressed as follows:

[0047] h(z)=ρ w g(z-z0);

[0048] Where z is the depth, z0 is the inlet depth, and h(z) is the matrix suction at depth z. The formula is:

[0049]

[0050] The water saturation S at different depths in different saturation zones was calculated using a combination of methods. w (h). Substituting the calculated water saturation into the sound field and electric field expressions, the seismoelectric wave field under different saturation distributions can be described.

[0051] S3. Given that the propagation speed of electromagnetic waves is much greater than that of seismic waves, a variable time step strategy is adopted. Define Δt. sThe time step for iterating the sound field is determined by placing the sound pressure P at integer grid points and the horizontal and vertical velocity components at half grid points, based on the principle of staggered grid finite difference. This yields a 2N-order spatial accuracy seismic wave differential iteration scheme, used to update the sound pressure and velocity components. The calculation formula is as follows:

[0052]

[0053] Among them, v x Let v be the velocity component in the x-direction. z Let Δx be the velocity component in the z-direction, Δz be the spatial step size along the x-axis, and Δz be the spatial step size along the z-axis. Let Δt be the pressure value at different spatial locations at time step k. EM The time step of the electric field is used to obtain the electric field component j. s , At time t = (l + 1 / 2)Δt EM The expression for is used to calculate electric field-related quantities, and the calculation formula is as follows:

[0054]

[0055] Among them, Q v Let κ be the charge density, η be the permeability, and η be the charge density. f Let be the dynamic viscosity of the fluid, and B be the Skempton coefficient, used to describe the effect of pore pressure variation on effective stress. The spatial difference coefficients, representing the 2Nth order of precision, can be calculated using Taylor expansion. For example, when 2N = 4,

[0056] S4. Using a linear prediction method, the current source density calculated at the sound field time points is extended to more dense electric field time points, achieving time decoupling between seismic waves and electromagnetic waves in the iterative steps. The expression is:

[0057]

[0058] Where, j x,z (x,z) represents the current density at position (x,z), where x and z are the abscissa and ordinate in a two-dimensional coordinate system, respectively. k is the sequence number of the sound field time step, and w is the number from 1 to Δt. s / Δt EM integers, (k)Δt s and (k+1)Δt s wΔt represents the current time and the time corresponding to the next sound field time step, respectively. EMThis represents the electric field offset relative to a certain acoustic field time step. The current density can be calculated iteratively using the seismic wave time step, thus achieving time decoupling between seismic waves and electromagnetic waves in the iteration steps, thereby improving the algorithm efficiency for calculating the seismoelectric wave field.

[0059] In practical implementation, the basic parameters of the model are first determined, as shown in the table below:

[0060] Table 1. Porous Media and Fluid Parameters

[0061]

[0062]

[0063] S1. Based on the data in Table 1, determine the basic parameters of the model, including the porosity of the porous medium, the density of the solid particles, the water density, the air density, the residual water content of the porous medium, and the saturated water content. Calculate the equivalent density ρ using the formula:

[0064] ρ=(1-φ)ρ s +φρ f ;

[0065] And based on the water saturation S at different depths w To adjust the equivalent density ρ.

[0066] To calculate the water saturation distribution within a depth range of 0 to 70 meters, it is first necessary to use the given formula for the relationship between matrix suction and depth:

[0067] h(z)=ρ w g(z-z0);

[0068] Determine the matrix suction h(z) at each depth z, and assume the inlet depth z0 = 50.

[0069] S2. Next, the water content θ(h) at a specific depth is calculated using the Van Genuchten model formula. The calculation formula is as follows:

[0070]

[0071] Where, the parameter α = 0.01, n = 2, and

[0072] Then, according to the definition of water saturation:

[0073]

[0074] Calculate the water saturation S at the corresponding depth. w(h). The entire process involves starting with determining the matrix suction at depth z, converting it into water content using the Van Genuchten formula, and finally calculating the water saturation to describe the soil moisture content at different depths.

[0075] S3. When performing separate calculations of the acoustic and electric fields in seismoelectric theory, the acoustic field is calculated first. The given bulk modulus K of the porous medium is used. b (Pa)=1.1×10 10 Pa and the skeleton shear modulus G of the porous medium b (Pa)=1.05×10 10 Pa calculates the undrained bulk modulus K. u This calculation is based on water saturation. Assume the power source F = F0sin(ωt), where F0 = 100 N and frequency f = 10 Hz. The calculation is performed using the staggered-grid finite difference method with a time step Δt. s =1×10 -5 The sound pressure P and velocity components are updated iteratively.

[0076] For the electric field component, the conductivity σ and charge density Q are calculated based on different degrees of saturation. v Fluid viscosity η f The current density j is determined by combining the permeability κ of the porous medium with other parameters. s and electric potential The distribution of the electric field. The iterative calculation of the electric field uses a smaller time step Δt. EM =1×10 -8 To ensure accuracy, the entire process starts with fundamental physical parameters and performs precise numerical simulations of the acoustic and electric fields, aiming to meticulously analyze the interaction between seismic waves and electromagnetic fields.

[0077] Using the sound field time step Δt respectively s =1×10 -5 s and electric field time step Δt EM =1×10 -8 The system performs independent iterative calculations. Within each time step, the sound pressure P, velocity components, and electric field components are updated according to the corresponding differential iterative scheme.

[0078] S4. The current source density in the sound field calculation is extended to a more dense electric field time node using a linear prediction method. The calculation formula is as follows:

[0079]

[0080] Where, j x,z(x,z) represents the current density at position (x,z), where x and z are the abscissa and ordinate in a two-dimensional coordinate system, respectively. k is the sequence number of the acoustic field time step, and w is the sequence number of the electric field time step. (k)Δt s and (k+1)Δt s wΔt represents the current time and the time corresponding to the next sound field time step, respectively. EM This is the electric field offset relative to a certain sound field time step, used to interpolate and calculate the current density at intermediate moments.

[0081] like Figure 1 As shown, this method can effectively simulate all three types of seismic electromagnetic wave fields. Figure 1 Figure (b) shows the comparison between the proposed method and the traditional seismoelectric finite difference method. It can be seen that the results of the two methods are in high agreement, thus verifying the correctness of the proposed method.

[0082] Therefore, this invention employs the aforementioned fast finite difference solution method and system for seismoelectric wave fields under different saturations. By setting different time steps for the sound field and electric field, the computational load is effectively reduced. Furthermore, by using linear prediction techniques, the current source density calculated at the sound field time point is extended to a denser electric field time point, thus achieving time decoupling between seismic waves and electromagnetic waves during the iteration process. This improves the accuracy and completeness of the seismoelectric wave field calculation, making the simulation results closer to the actual physical process.

[0083] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A fast finite-difference solution method for seismoelectric wave fields under different saturation levels, characterized in that, Includes the following steps: S1. Set the basic parameters of the model; the basic parameters of the model include the porosity of the porous medium, the density of solid particles, the water density, the air density, the residual water content of the porous medium, and the saturated water content. S2. Based on the Van Genuchten model, combined with the relationship between matrix suction and depth and the water content formula, calculate the water saturation at different depths. The relationship between matrix suction and depth is as follows: ； In the formula, For depth, This refers to the depth of the air intake point. In depth Matrix suction at the location, For the density of water, It is the acceleration due to gravity; The Van Genuchten formula for water content is: ； In the formula, Moisture content, This refers to the residual water content of the porous medium. This is the saturated water content. , and For parameters; Fit feature parameters to the Van Genuchten model; The formula for calculating water saturation is: ； In the formula, Water saturation; S3. Apply the water saturation value from S2 to the expressions for calculating the sound and electric fields. Iterate the calculations of the sound and electric fields by introducing independent time iteration steps. Specifically: The equivalent fluid density, fluid bulk modulus, and fluid viscosity required for sound field calculation are corrected by water saturation. The corrected fluid parameters are then substituted into the sound field wave equations. Simultaneously, the electrical conductivity, permeability, and effective charge per unit pore volume of the medium required for electric field calculation are corrected by water saturation. The corrected electrical parameters are then substituted into the quasi-static Maxwell's equations for the electric field. Equivalent fluid density Fluid bulk modulus and fluid viscosity Represented as: ； ； ； In the formula, The mass density of the pore fluid; The density of water; air density; The bulk modulus of water; The bulk modulus of air; The dynamic viscosity of water; The dynamic viscosity of air; dielectric conductivity Penetration rate and effective charge per unit pore volume Represented as: ； ； ； In the formula, The conductivity of the aqueous phase; Porosity of porous media; For the tortuosity of porous media; This is the base permeability when the medium is fully saturated; This is the saturation correction index; The residual water saturation of the porous medium; The effective charge of the pore volume; Given that the propagation speed of electromagnetic waves is much greater than that of seismic waves, a variable time step strategy is adopted, defining the time iteration step of the sound field. , to sound pressure By placing the components at integer grid points and the horizontal and vertical velocity components at half grid points respectively, a 2N-order spatially accurate seismic wave differential iteration scheme is obtained, which is used to update the sound pressure and velocity components. The sound field iteration calculation formula is as follows: ； ； ； In the formula, for Directional velocity component, for Directional velocity component, for Spatial step size of the axis, for Spatial step size of the axis, In time step Pressure values ​​at different spatial locations; Discrete indexes of two-dimensional computational grids; This represents the traversal index of the difference summation term; This refers to the sequence number of the sound field iteration time step; It is the half-value of the difference order; For 2N-order precision spatial difference weighting coefficients; Equivalent density; The bulk modulus is the undrained bulk modulus. Shear modulus; Define the electric field independent time iteration step size The source components of the electric field current are obtained. , In time The expression for is calculated using the following formula: ； ； ； In the formula, This refers to the independent time step sequence number of the electric field; Integer grid Electric field The dielectric potential of the time layer; For penetration rate; The dynamic viscosity of the fluid; The dielectric conductivity; Skempton coefficient; S4. The results of the sound field calculation in S3 are extended to the electric field time node using a linear prediction method to obtain the calculation results, and the calculation results are verified.

2. The fast finite difference solution method for seismoelectric wave fields under different saturations according to claim 1, characterized in that, In S4, the linear prediction method is as follows: Extending the current source density calculation of the sound field to the electric field time node, the calculation formula is as follows: ； in, For in position Current density at that point and These are the x-coordinate and y-coordinate in a two-dimensional coordinate system, respectively. This is the sequence number of the sound field time step. From 1 to integers, and These represent the current moment and the moment corresponding to the next sound field time step, respectively. This represents the electric field offset relative to a certain sound field time step.

3. A fast finite difference solution system for seismic wave fields under different saturation levels, characterized in that, include: The initialization module is used to set the basic parameters of the model; the basic parameters include the porosity of the porous medium, the density of the solid particles, the water density, the air density, the residual water content of the porous medium, and the saturated water content. The water saturation distribution calculation module at formation depth is connected to the initialization module and is used to calculate the water saturation at different depths based on the VanGenuchten model combined with the relationship between matrix suction and depth and the water content formula. The sound field and electric field calculation module is connected to the water saturation distribution calculation module at the depth of the formation. It is used to correct the equivalent fluid density, fluid bulk modulus, and fluid viscosity required for sound field calculation, and the medium conductivity, permeability, and effective charge per unit pore volume required for electric field calculation, respectively, based on the calculated water saturation. The corrected fluid parameters are substituted into the sound field wave equations and the corrected electrical parameters are substituted into the electric field quasi-static Maxwell equations. Based on the characteristic that the propagation speed of electromagnetic waves is much greater than that of seismic waves, independent sound field time iteration steps and electric field time iteration steps are introduced respectively. The sound pressure, velocity components, medium potential, and current source components are iteratively solved using a 2N-order spatial precision finite difference scheme to obtain the step-by-step iterative values ​​of the sound field and electric field. The linear prediction module, connected to the sound field and electric field calculation modules, is used to extend and match the sound field calculation output to the electric field time node using the linear prediction method, complete the conversion of the seismoelectric wave field coupling result, and verify the final coupling calculation result.

Images