Method, device, equipment and medium for quickly calculating guided wave dispersion curve of viscoelastic anisotropic shale reservoir

Abstract

The invention discloses a method, a device, equipment and a medium for rapidly calculating a viscoelastic anisotropic shale reservoir guided wave dispersion curve. The method comprises the following steps: obtaining parameters of a shale reservoir viscoelastic anisotropic medium; based on a Kelvin-Voigt viscoelastic medium model, calculating a frequency domain stiffness matrix through shale reservoir viscoelastic anisotropy medium parameters; based on the frequency domain stiffness matrix, constructing a guided wave frequency dispersion equation by adopting a Thomson-Haskell transfer matrix method; the guided wave frequency dispersion equation is solved and searched by adopting a fast recursion method, a guided wave phase velocity frequency dispersion curve and an attenuation coefficient curve in the viscoelastic medium are obtained, and the guided wave phase velocity frequency dispersion curve and the attenuation coefficient curve are used for calculating sensitive functions relative to different parameters. Therefore, the viscoelasticity and anisotropy of the medium are combined, errors generated in description of the shale reservoir can be effectively reduced, and therefore the accuracy of reservoir physical property parameter inversion based on the guided wave frequency dispersion property is improved.

Description

Technical Field

[0001] The present invention relates to a method, apparatus, device and medium for quickly calculating the dispersion curve of guided waves in a viscoelastic anisotropic shale reservoir, belonging to the technical field of geological exploration, and particularly to the technical fields of seismic wave forward modeling and inversion. Background Art

[0002] In recent years, unconventional oil and gas have played a key role in clean energy and are regarded as important alternatives to conventional oil and gas resources. Hydraulic fracturing is a key technology for the development of unconventional oil and gas resources such as shale gas. The hydraulic fracturing transformation of shale oil and gas reservoirs requires real-time monitoring and evaluation of its transformation efficiency and safety. The physical properties of shale reservoirs are crucial for geomechanical modeling, fracture imaging, and optimization of hydraulic fracturing design. Shale layers are usually thin, with a thickness ranging from a few meters to dozens of meters. Traditional surface reflection seismic methods have limitations in imaging resolution and are difficult to accurately detect the physical properties of thin reservoirs.

[0003] In a shale reservoir, low-velocity shale layers are located between high-velocity surrounding rocks, forming a special waveguide structure. The longitudinal waves and transverse waves excited by hydraulic fracturing perforation or microseismic sources are confined to propagate within the low-velocity shale layers, continuously reflecting, superposing, and interfering with each other at the top and bottom interfaces of the waveguide structure, thus forming guided waves. In the vertical plane, the P-SV type guided wave formed by the interference of P waves and SV waves is also called the Rayleigh type guided wave because of its similar physical properties to Rayleigh surface waves; in the horizontal plane, the SH type guided wave formed by the interference of SH waves is similar in physical properties to Love surface waves and is also called the Love type guided wave. Guided waves have obvious dispersion phenomena. As a kinematic characteristic, dispersion is closely related to the rock physical parameters of the low-velocity shale layer and the thickness of the shale reservoir. Therefore, guided waves can be used for high-precision imaging of shale reservoir structures. In addition, guided waves have a wide range of applications in the field of geophysics. For example, the guided waves (also called channel waves) propagating in low-velocity coal seams are often used to detect small geological structures in coal seams; the low-velocity fault gouge in the fault zone and the high-velocity surrounding rocks form a waveguide structure, making the seismic waves propagating inside form a special seismic phase - fault trapped wave (or fracture zone guided wave), which can be used for fracture zone structure imaging.

[0004] Shales have different elastic properties in different directions, thus exhibiting obvious anisotropy, which will affect the dispersion characteristics of guided waves propagating in them. At the same time, shales exhibit viscoelasticity, and the energy of guided waves will be dissipated during propagation, resulting in the attenuation of the amplitude. The assumption of a completely elastic isotropic medium in traditional methods is an overly simplified assumption, which cannot accurately describe the propagation characteristics of guided waves in actual formations. Whether in elastic media or viscoelastic media, guided waves have phase velocity dispersion, and the anisotropy of the medium will significantly affect the phase velocity dispersion curve of guided waves. At the same time, the viscoelasticity of the medium will also have a certain impact on the phase velocity dispersion curve of guided waves. In viscoelastic media, there is another dispersion relationship between frequency and attenuation coefficient, which can be expressed as the attenuation dispersion curve of guided waves. The viscoelasticity of the medium will significantly affect the attenuation dispersion curve of guided waves. Therefore, to accurately and comprehensively describe the properties of guided waves, it is necessary to consider viscoelasticity and anisotropy. However, current research on the dispersion properties of guided waves mainly focuses on isotropic completely elastic media. Therefore, how to forward model the dispersion curves of guided waves in viscoelastic anisotropic media to provide a more reliable scientific basis for the exploration and development of shale gas is an urgent problem to be solved at present. Summary of the Invention

[0005] The present invention aims to solve at least one of the technical problems existing in the prior art. For this purpose, in view of the above problems, the object of the present invention is to provide a method, device, equipment and medium for quickly calculating the dispersion curves of guided waves in viscoelastic anisotropic shale reservoirs, which can improve the accuracy of inverting reservoir physical property parameters based on the dispersion properties of guided waves.

[0006] To achieve the above object of the invention, the technical solution adopted by the present invention is as follows:

[0007] In a first aspect, the present invention provides a method for quickly calculating the dispersion curves of guided waves in viscoelastic anisotropic shale reservoirs, including:

[0008] Obtaining the viscoelastic anisotropic medium parameters of the shale reservoir;

[0009] Based on the Kelvin-Voigt viscoelastic medium model, inputting the viscoelastic anisotropic medium parameters of the shale reservoir to calculate the stiffness matrix in the frequency domain;

[0010] Based on the stiffness matrix in the frequency domain, using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation;

[0011] Solving and finding the roots of the guided wave dispersion equation to obtain the phase velocity dispersion curve and attenuation coefficient curve of guided waves in the viscoelastic medium, wherein the phase velocity dispersion curve and attenuation coefficient curve of guided waves are used to calculate the sensitivity functions for different parameters.

[0012] In some possible embodiments, the viscoelastic anisotropic medium parameters of the shale reservoir include density ρ, vertical P-wave velocity V p0 and V s0 ; Thomsen anisotropic parameters ε, γ, and δ; vertical P-wave and S-wave quality factors Q p0 and Q s0 ; quality factor anisotropic parameters ε Q and δ Q and reservoir thickness h.

[0013] In some possible embodiments, each element in the stiffness matrix is as follows:

[0014]

[0015] The real part of the element c ij in the stiffness matrix is used to represent the anisotropy of velocity, and the calculation formula is:

[0016]

[0017] In the formula, other elements are all zero except for the above;

[0018] The imaginary part Q ij of the element c ij in the stiffness matrix is used to represent the anisotropy of the attenuation coefficient, and the calculation formula is:

[0019] Q 33 = Q p0 ;

[0020] Q 55 = Q s0 ;

[0021]

[0022] In the formula, the imaginary part Q ij of the element c ij in the stiffness matrix is similar to the real part ij of the stiffness matrix element c in terms of symmetry, and other elements are all zero.

[0023] In some possible embodiments, based on the frequency-domain stiffness matrix, the Thomsen-Haskell transfer matrix method is used to construct the guided-wave dispersion equation, including: substituting the plane-wave solution into the elastic wave equation to obtain the matrix expression of the guided-wave solution function; obtaining the transfer matrix from the continuity of displacement and stress at each interface; applying the boundary condition that the displacement field at infinity is zero to the upper and lower semi-infinite spaces; and constructing the guided-wave dispersion equation based on the matrix expression, transfer matrix, and boundary condition using the Thomsen-Haskell transfer matrix method.

[0024] In some possible embodiments, the guided-wave dispersion curve is obtained by solving the following dispersion equation:

[0025]

[0026] In the formula, ω, c r and α are respectively the angular frequency, phase velocity and attenuation coefficient of the guided wave, U T represents the boundary condition of the guided wave in the upper half-infinite space, and T1T2…T n-1 T n represents the transfer matrix of the intermediate layer, and V represents the boundary condition of the lower half-infinite space. Among them, the guided wave propagating in the layered medium needs to satisfy the boundary conditions: in the upper half-infinite space, the boundary condition is that the displacement field is zero at infinity; in the intermediate low-velocity layer, the displacement field and stress field are continuous at the interfaces of different layers; in the lower half-infinite space, the boundary condition is that the displacement field is zero at infinity.

[0027] In some possible embodiments, the fast recursive method is used to solve and find the roots of the guided-wave dispersion equation to obtain the guided-wave phase velocity dispersion curve and attenuation coefficient curve in the viscoelastic medium, that is, to traverse all values in the three directions of phase velocity, frequency and attenuation coefficient to find the points that satisfy the dispersion equation equal to zero. The specific process is as follows:

[0028] Fix the frequency ω0, and transform the root finding in the three-dimensional space of phase velocity c, frequency ω0 and attenuation coefficient α into root finding in the two-dimensional plane of phase velocity c and attenuation coefficient α;

[0029] Traverse all values of phase velocity c and attenuation coefficient α in the plane. At this time, the roots of the dispersion equation satisfy the local minimum, complete the root finding of all modes at frequency ω0, and perform the root finding of the next frequency ω0+Δω until the root finding of all frequencies is completed. The recursive method is used to set the initial step sizes Δc and Δα, and then according to the position of the zero point, a new root finding range [c i -Δc, c i +Δc] and [α i -Δα, α i +Δα] is given, and the new step sizes Δc / n and Δα / n are used to continue to refine the root finding range. The number of recursions satisfies the set conditions, and the above process is iterated until the required root finding accuracy is satisfied;

[0030] After solving the zeros at all frequencies, the relationships between the guided-wave phase velocity and frequency and between the guided-wave attenuation coefficient and frequency in the viscoelastic medium are solved from the dispersion equation, that is, the guided-wave phase velocity dispersion curve and attenuation coefficient curve in the viscoelastic medium are obtained.

[0031] In some possible embodiments, the sensitivity functions with respect to different parameters are calculated using the finite difference method, and the calculation formulas are respectively:

[0032]

[0033] Among them, is a sensitivity function of the phase velocity to different parameters, S α is a sensitivity function of the attenuation coefficient to different parameters, f is the frequency, m is the model parameter, and Δm is the perturbation of the model parameter. Among them, the calculation formula for specific parameters is:

[0034]

[0035] Among them, is a sensitivity function of the guided wave phase velocity relative to the shear wave velocity, S ε is a sensitivity function of the guided wave phase velocity relative to the anisotropy parameter ε, S δ is a sensitivity function of the guided wave phase velocity relative to the anisotropy parameter δ, c is the guided wave phase velocity, ΔV s0 , Δε and Δδ represent model perturbations; is a sensitivity function of the guided wave attenuation coefficient relative to the shear wave quality factor, is a sensitivity function of the guided wave attenuation coefficient relative to the longitudinal wave quality factor, and is a sensitivity function of the guided wave attenuation coefficient relative to the anisotropy parameter of the quality factor, α is the guided wave attenuation coefficient, ΔQ s0 , ΔQ p0 , Δε Q and Δδ Q represent model parameter perturbations.

[0036] In a second aspect, the present invention also provides a device for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir, including:

[0037] A parameter acquisition unit configured to obtain the viscoelastic anisotropic medium parameters of the shale reservoir;

[0038] A stiffness matrix calculation unit configured to calculate the stiffness matrix in the frequency domain based on the Kelvin-Voigt viscoelastic medium model through the viscoelastic anisotropic medium parameters of the shale reservoir;

[0039] An equation construction unit configured to construct a guided wave dispersion equation by using the Thomsen-Haskell transfer matrix method based on the stiffness matrix in the frequency domain;

[0040] An equation solving unit configured to solve and find the root of the guided wave dispersion equation to obtain the guided wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium, where the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity functions relative to different parameters.

[0041] In a third aspect, the present invention provides an electronic device, including: at least one processor; and a memory communicatively connected to the processor; wherein, the memory stores instructions executable by the processor, and when the instructions are executed by the processor, the processor is enabled to execute the method as described above.

[0042] In a fourth aspect, the present invention provides a computer-readable storage medium storing one or more programs, the one or more programs including computer instructions for causing a computer to execute the method as described above.

[0043] Due to the above technical solutions adopted by the present invention, it has the following characteristics:

[0044] 1. Compared with the traditional viscoelastic isotropic forward modeling method and elastic anisotropic forward modeling, the present invention can more accurately reflect reservoir information based on the viscoelastic anisotropic theory. The forward modeling theory based on a completely elastic medium cannot obtain the parameter quality factor that can reflect wave amplitude attenuation, and at the same time, the forward modeling theory based on an isotropic medium cannot accurately reflect the velocity information of the medium. Therefore, by combining the viscoelasticity and anisotropy of the medium, the error generated in the description of shale reservoirs can be effectively reduced, thereby improving the accuracy of inverting reservoir physical parameters based on the dispersion properties of guided waves, and further contributing to a deeper understanding of the mechanical properties, stress state, and fracture distribution of the reservoir, thus improving the accuracy of reservoir evaluation and the rationality of development decisions.

[0045] 2. The present invention uses a fast recursive method to solve and find the roots of the guided wave dispersion equation, which can find the roots and solve more efficiently and stably. Even in the high-frequency and high-order parts, there will be no missed roots and jumping orders. The idea of recursion is introduced in the process of finding the roots, ensuring the computational efficiency of finding the roots.

[0046] 3. The present invention can calculate the sensitivity functions of the Rayleigh-type guided wave phase velocity dispersion curve to the vertical shear wave velocity, Thomsen anisotropic parameters ε and δ through numerical solutions; the Rayleigh-type guided wave attenuation coefficient curve to the vertical P-wave and S-wave velocity quality factors, quality factor anisotropic parameters ε Q and δ Q The sensitivity functions of, through the sensitivity analysis of the above parameters, provide support for parameter inversion in viscoelastic anisotropic media.

[0047] In summary, the present invention can be widely applied in geological exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] By reading the following detailed description of the preferred embodiments, various other advantages and benefits will become clear to those of ordinary skill in the art. The drawings are only for the purpose of showing the preferred embodiments and are not considered to be a limitation of the present invention. Throughout the drawings, the same reference numerals are used to represent the same components. In the drawings:

[0049] Figure 1 Flow chart of the method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to an embodiment of the present invention;

[0050] Figure 2 Schematic diagram of the solution process of the guided wave dispersion equation by the fast recursive method according to an embodiment of the present invention;

[0051] Figure 3 Forward modeling of the dispersion curve of the guided wave phase velocity in the viscoelastic anisotropic medium model according to an embodiment of the present invention;

[0052] Figure 4 Forward modeling of the guided wave attenuation coefficient curve in the viscoelastic anisotropic medium model according to an embodiment of the present invention;

[0053] Figure 5 Sensitivity analysis of the guided wave phase velocity to the vertical shear wave velocity according to an embodiment of the present invention;

[0054] Figure 6 Sensitivity analysis of the guided wave phase velocity to the Thomsen anisotropic parameter ε according to an embodiment of the present invention;

[0055] Figure 7 Sensitivity analysis of the guided wave phase velocity to the Thomsen anisotropic parameter δ according to an embodiment of the present invention;

[0056] Figure 8 Sensitivity analysis of the guided wave attenuation coefficient to the quality factor of the vertical longitudinal wave velocity according to an embodiment of the present invention;

[0057] Figure 9 Sensitivity analysis of the guided wave attenuation coefficient to the quality factor of the vertical shear wave velocity according to an embodiment of the present invention;

[0058] Figure 10 Sensitivity analysis of the guided wave attenuation coefficient to the quality factor anisotropic parameter ε Q of the present invention;

[0059] Figure 11 Sensitivity analysis of the guided wave attenuation coefficient to the quality factor anisotropic parameter δ Q of the present invention;

[0060] Figure 12 Structural diagram of the electronic device according to an embodiment of the present invention. Detailed implementation manners

[0061] It should be understood that the terms used herein are for the purpose of describing particular example embodiments only and are not intended to be limiting. Unless the context clearly dictates otherwise, the singular forms "a", "an", and "the" as used herein may also include the plural forms. The terms "comprising", "including", "containing", and "having" are inclusive and thus specify the presence of stated features, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, steps, operations, elements, components, and / or combinations thereof. The method steps, processes, and operations described herein are not to be construed as necessarily requiring their performance in the particular order described or illustrated, unless the order of performance is explicitly stated. It should also be understood that additional or alternative steps may be used.

[0062] Although the terms first, second, third, etc. may be used herein to describe multiple elements, components, regions, layers, and / or sections, these elements, components, regions, layers, and / or sections should not be limited by these terms. These terms may be used only to distinguish one element, component, region, layer, or section from another. Unless the context clearly indicates otherwise, terms such as "first", "second", and other numerical terms when used herein do not imply an order or sequence. Thus, a first element, component, region, layer, or section discussed below may be referred to as a second element, component, region, layer, or section without departing from the teachings of the example embodiments.

[0063] For ease of description, spatial relative relationship terms may be used herein to describe the relationship of one element or feature shown in the figures to another element or feature, such as "inner", "outer", "inside", "outside", "below", "above", etc. Such spatial relative relationship terms are intended to include different orientations of the device in use or operation in addition to the orientation depicted in the figures.

[0064] Since the current research on the dispersion properties of guided waves mainly focuses on isotropic and fully elastic media, it is unable to accurately and comprehensively describe the properties of guided waves. The method, device, equipment, and medium for quickly calculating the dispersion curves of guided waves in viscoelastic anisotropic (vertical transversely isotropic: VTI) shale reservoirs provided by the present invention include: obtaining the viscoelastic anisotropic medium parameters of the shale reservoir; based on the Kelvin-Voigt viscoelastic medium model, inputting the viscoelastic anisotropic medium parameters of the shale reservoir to calculate the stiffness matrix in the frequency domain; based on the stiffness matrix in the frequency domain, using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation; solving and finding the roots of the guided wave dispersion equation to obtain the phase velocity dispersion curve and attenuation coefficient curve of the guided wave in the viscoelastic medium, where the phase velocity dispersion curve and attenuation coefficient curve of the guided wave are used to calculate the sensitivity functions for different parameters. Therefore, the present invention derives the analytical expression of the guided wave dispersion curve in viscoelastic VTI media based on the Kelvin-Voigt viscoelastic medium model and the Thomsen-Haskell transfer matrix method, calculates the sensitivity functions of the guided wave phase velocity dispersion curve relative to the vertical shear wave velocity, Thomsen parameters ε and δ, and calculates the sensitivity functions of the guided wave attenuation coefficient curve relative to the quality factor of the vertical shear wave velocity, the quality factor of the vertical longitudinal wave velocity, and the attenuation anisotropy parameter ε Q and δ Q to improve the accuracy of inverting reservoir physical properties parameters based on the dispersion properties of guided waves.

[0065] The exemplary embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present invention are shown in the drawings, it should be understood that the present invention can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present invention can be more thoroughly understood and the scope of the present invention can be fully conveyed to those skilled in the art.

[0066] Example 1: As Figure 1 shown, the method for quickly calculating the dispersion curves of guided waves in viscoelastic anisotropic shale reservoirs provided in this embodiment includes:

[0067] S1. Obtain the viscoelastic anisotropic medium parameters of the shale reservoir from comprehensive information such as logging and geology.

[0068] In this embodiment, the viscoelastic anisotropic medium parameters of the shale reservoir include density ρ, vertical longitudinal and transverse wave velocities V p0 and V s0 ; Thomsen anisotropy parameters ε and δ; vertical longitudinal and transverse wave quality factors Q p0 and Q s0 ; quality factor anisotropy parameter ε qand δ Q and the reservoir thickness h.

[0069] S2. Based on the Kelvin-Voigt viscoelastic medium model, calculate the stiffness matrix C in the frequency domain through the viscoelastic anisotropic medium parameters of the shale reservoir.

[0070] In this embodiment, the calculation process of the stiffness matrix C in the frequency domain is as follows:

[0071] For the viscoelastic VTI medium, in the frequency domain, the constitutive equation can be expressed as:

[0072] σ(ω) = C(ω)∈(ω);

[0073] where σ = (σ xx , σ yy , σ zz , σ yz , σ xz , σ xy ) T represents the stress vector, ∈ represents the strain vector, C represents the stiffness matrix of the viscoelastic VTI medium, ω is the angular frequency, x, y, and z represent the coordinate axes, and u, v, and w respectively represent the displacement fields along the x, y, and z directions.

[0074] In the viscoelastic medium, the elements in the stiffness matrix C can be expressed in the form of complex numbers in the frequency domain. Using the Kelvin-Voigt viscoelastic medium model, each element in the stiffness matrix can be expressed in the following format:

[0075]

[0076] The real part of the stiffness matrix element c ij can be calculated from the vertical P-wave and S-wave velocities V and V p0 and the Thomson anisotropic parameters ε, δ, and γ, and is used to represent the velocity anisotropy, as shown below: S0

[0077]

[0078]

[0079] In the formula, all other elements are zero except for the above.

[0080] The imaginary part Q ij of the stiffness matrix element c ij can be calculated from the vertical P-wave and S-wave quality factors Q p0 and Q s0 and the quality factor anisotropic parameter ε Q , δ​Q and γ Q It is calculated and used to represent the anisotropy of the attenuation coefficient, as follows:

[0081] Q 33 = Q p0 ;

[0082] Q 55 = Q s0 ;

[0083]

[0084] In the formula, the elements in the Q matrix and the stiffness matrix C are similar in symmetry.

[0085] S3. Based on the above frequency-domain stiffness matrix C, the Thomsen-Haskell transfer matrix method is used to construct the guided-wave dispersion equation.

[0086] In this embodiment, the plane-wave solution is substituted into the elastic wave equation to obtain the matrix expression of the guided-wave solution function; the transfer matrix is obtained from the continuity of displacements and stresses at each interface; the boundary condition of zero displacement field at infinity is applied to the upper and lower semi-infinite spaces; the Thomsen-Haskell transfer matrix method is used to construct the guided-wave dispersion equation based on the matrix expression, the transfer matrix, and the boundary condition.

[0087] Furthermore, the elastic wave equation is composed of the motion equation and the constitutive equation. The motion equation in the frequency domain can be expressed in the following form:

[0088]

[0089]

[0090] The constitutive equation in the frequency domain can be expressed in the following form:

[0091]

[0092] In the shale reservoir, the Rayleigh-type guided wave is formed by the interference of P waves and SV waves. Only the waves propagating along the x-axis need to be considered because the vertical transversely isotropic medium is symmetric about the z-axis. Under this condition:

[0093] v = 0

[0094] σ xy = σ yy = σ yz = 0.

[0095] Substituting the above conditions into the constitutive equation and the motion equation, the following wave equation is obtained:

[0096]

[0097] In the case of only considering Rayleigh-type guided waves (x-z plane), the velocity anisotropy parameter γ and the attenuation anisotropy parameter γ Q are ignored. The solutions of the displacement field and the stress field in the frequency domain in the above elastic wave equation are set as:

[0098] w(ω) = y1(z)e -ikx 2πδ(ω′ - ω);

[0099] σ zz (ω) = y2(z)e -ikx 2πδ(ω′ - ω);

[0100] u(ω) = iy3(z)e -ikx 2πδ(ω′ - ω);

[0101] σ zx (ω) = iy4(z)e -ikx 2πδ(ω′ - ω).

[0102] where δ is the Dirac function, and y = (y1, y2, y3, y4) T is the function to be solved. The guided wave number is extended to the complex domain in the viscoelastic medium as follows:

[0103]

[0104] where k, ω, c r and α are the wave number, angular frequency, phase velocity and attenuation coefficient of the guided wave respectively.

[0105] The function y(z) can be expressed in the form of the product of the matrix D(z) and the amplitude vector F:

[0106] y(z) = D(z)F;

[0107] where:

[0108]

[0109] where F1, F2, F3, F4 represent amplitudes, and the elements of the matrix D(z) are calculated by the following formulas:

[0110] C1(z) = coshυ1z;

[0111] C3(z) = coshυ3z;

[0112] S1(z) = sinhu1z;

[0113] S3(z) = sinhυ3z;

[0114] d1 = c 33 υ1 - c 13 kε1;

[0115] d2 = c 33 υ3 - c 13 kε3;

[0116] d3 = c 55 (k + υ1ε1);

[0117] d4 = c 55 (k + υ3ε3).

[0118] Among them, υ1, υ 2、 υ3, υ4 are eigenvalues, representing the vertical wavenumbers of the up and down longitudinal and transverse waves, and are given by the following formula:

[0119]

[0120] M1 = c 55 (ρω 2 - c 55 k 2 ) + c 33 (ρω 2 - c 11 k 2 ) + (c 13 + c 55 ) 2 k 2 ;

[0121] M2 = (ρω 2 - c 55 k 2 )(ρω 2 - c 11 k 2 );

[0122] Among them, ε1 and ε3 are given by the following formula:

[0123]

[0124] The guided wave propagating in the layered medium needs to satisfy the boundary conditions: in the upper half - infinite space, the boundary condition is that the displacement field is zero at infinity; in the middle low - velocity layer, the displacement field and the stress field are continuous at the interfaces of different layers; in the lower half - infinite space, the boundary condition is that the displacement field is zero at infinity.

[0125] The guided - wave dispersion curve is obtained by solving the following dispersion equation:

[0126] D(ω, c r , α r ) = det|U TT1T2…T n-1 T n V|=0;

[0127] where U T represents the boundary condition of the guided wave in the upper half-infinite space, and T1T2…T n-1 T n represents the transfer matrix of the intermediate layer, and V represents the boundary condition of the lower half-infinite space:

[0128]

[0129] d5 = c 55 (υ1ε1 - υ3ε3);

[0130] d6 = c 33 (υ3ε1 - υ1ε3).

[0131] S4. Use the fast recursive method to solve and find the roots of the guided wave dispersion equation.

[0132] In this embodiment, it can be seen from the dispersion equation that in the viscoelastic anisotropic medium, the dispersion curve of the guided wave can be regarded as existing in a three-dimensional space (frequency, phase velocity, attenuation coefficient), and a new root-finding method needs to be used to solve the dispersion curve of the guided wave in the viscoelastic medium.

[0133] Furthermore, the method adopted by the present invention is the fast recursive dispersion equation root-finding algorithm, which traverses all values in the three directions of phase velocity, frequency, and attenuation coefficient to find the points that satisfy the dispersion equation equal to zero, and improves the efficiency and accuracy through the recursive algorithm. The specific calculation method is as follows:

[0134] The first step is to fix the frequency ω0, and transform the root-finding in the three-dimensional space of phase velocity c, frequency ω0, and attenuation coefficient α into root-finding in the two-dimensional plane of phase velocity c and attenuation coefficient α;

[0135] The second step is to traverse all values of phase velocity c and attenuation coefficient α in the plane. At this time, the roots of the dispersion equation satisfy the local minimum, complete the root-finding of all modes at frequency ω0, and perform the root-finding of the next frequency ω0 + Δω until the root-finding of all frequencies is completed. To improve the efficiency and accuracy of root-finding, the present invention adopts a recursive method. The initial step sizes Δc and Δα can be set to larger values, and then new root-finding ranges [c i -Δc, c i +Δc] and [α i -Δα, α i +Δα] are given according to the position of the zero point, and the new step sizes Δc / n and Δα / n are used to continue to accurately find the root-finding range. The recursive times are about 10 times, and the above process is iterated until the required root-finding accuracy is satisfied.

[0136] After solving for the zeros at all frequencies, the relationship between the phase velocity of guided waves and frequency and the relationship between the attenuation coefficient of guided waves and frequency in a viscoelastic medium can be obtained from the dispersion equation.

[0137] S5. Calculate the sensitivity functions of the guided wave phase velocity dispersion curve and attenuation coefficient curve based on finite differences with respect to different parameters.

[0138] In this embodiment, by means of the phase velocity dispersion curve and attenuation coefficient curve of guided waves in a viscoelastic anisotropic medium, the sensitivity of the phase velocity and attenuation coefficient of guided waves to different parameters in a viscoelastic medium can be further analyzed. The present invention uses the finite difference method to approximate the partial derivatives of the guided wave phase velocity dispersion curve and attenuation coefficient curve with respect to different parameters, that is, the sensitivity functions:

[0139]

[0140] where is the sensitivity function of the phase velocity to different parameters, s α is the sensitivity function of the attenuation coefficient to different parameters, f is the frequency, m is the model parameter, and Δm is the perturbation of the model parameter.

[0141] Furthermore, through sensitivity analysis of all parameters affecting the guided wave phase velocity dispersion and attenuation coefficient curves, the parameters that have obvious effects on the two curves are obtained. Among them, the parameters sensitive to the guided wave phase velocity dispersion curve are the vertical shear wave velocity, Thomsen anisotropic parameters ε and δ; the parameters sensitive to the guided wave attenuation coefficient curve are the vertical P-wave velocity quality factor, vertical S-wave velocity quality factor, quality factor anisotropic parameter ε Q and δ Q . By analyzing the sensitivity of different parameters, a theoretical basis can be further provided for the inversion of these parameters.

[0142] The following details the application of the method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir of the present invention through specific embodiments.

[0143] The forward modeling method for the Rayleigh-type guided wave phase velocity dispersion curve and quality factor dispersion curve of a viscoelastic anisotropic shale reservoir provided in this embodiment, where the model parameters for forward modeling are shown in Table 1:

[0144]

[0145] Based on the model parameters in Table 1, using the Kelvin-Voigt viscoelastic medium model, calculate the stiffness matrix C in the frequency domain through the viscoelastic anisotropic medium parameters of the shale reservoir, construct the guided wave dispersion equation using the Thomsen-Haskell transfer matrix method, solve and find the roots of the guided wave dispersion equation by the fast recursive method, and the obtained guided wave phase velocity dispersion curve is asFigure 3 as shown, and the attenuation coefficient curve of the guided wave is as Figure 4 shown. For Figure 3 the fundamental order of the guided wave phase velocity dispersion curve in

[0146]

[0147]

[0148] wherein, is the sensitivity function of the guided wave phase velocity relative to the shear wave velocity, S ε is the sensitivity function of the guided wave phase velocity relative to the anisotropic parameter ε, S δ is the sensitivity function of the guided wave phase velocity relative to the anisotropic parameter δ, c is the guided wave phase velocity, f is the frequency, ΔV s0 , Δε and Δδ represent model perturbations, which are used to approximate the partial derivatives. The sensitivity analysis regarding the vertical shear wave velocity is as Figure 5 shown, and the sensitivity analysis regarding the Thomsen parameters ε and δ is as Figure 6 and Figure 7 shown. The guided wave phase velocity is most sensitive to the shear wave velocity structure V s0 , and has a certain sensitivity to the anisotropic parameters ε and δ. Therefore, it can be used for the inversion of the anisotropic parameters of shale reservoirs. For Figure 4 the fundamental order of the guided wave attenuation coefficient dispersion curve in Q analyze the sensitivity of the dispersion curve to the vertical shear wave velocity quality factor, the vertical compressional wave velocity quality factor, the Thomsen parameter ε Q describing attenuation and δ

[0149]

[0150] wherein, is the sensitivity function of the guided wave attenuation coefficient relative to the shear wave quality factor, is the sensitivity function of the guided wave attenuation coefficient relative to the compressional wave quality factor, and are the sensitivity functions of the guided wave attenuation coefficient relative to the anisotropic parameters of the quality factor, α is the guided wave attenuation coefficient, f is the frequency, ΔQ s0 , ΔQ p0 , Δε Q and Δδ Q represent model parameter perturbations. The sensitivity analysis regarding the vertical compressional wave velocity quality factor is as Figure 8 shown, and the sensitivity analysis regarding the vertical shear wave velocity quality factor is as Figure 9As shown, the quality factor anisotropy parameters ε Q and δ Q sensitivity analysis is as Figure 10 and Figure 11 shown. The guided wave attenuation coefficient is most sensitive to the quality factor parameters in the low-velocity layer. The guided wave attenuation coefficient curve and the sensitivity function can provide a basis for the study and inversion of the viscoelastic properties of shale reservoirs.

[0151] Example 2: The above Example 1 provides a method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir. Correspondingly, this example provides a device for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir. The device provided in this example can implement the method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir in Example 1, and the device can be implemented in a software, hardware, or a combination of software and hardware manner. For the convenience of description, when describing this example, various units are described separately according to their functions. Of course, in implementation, the functions of each unit can be implemented in the same or multiple software and / or hardware. For example, the device may include integrated or separate functional modules or functional units to execute the corresponding steps in the methods of Example 1. Since the device of this example is basically similar to the method embodiment, the description process of this example is relatively simple, and the relevant parts can refer to the partial description of Example 1. The embodiments of the device for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir provided by the present invention are only illustrative.

[0152] Specifically, the device for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir provided in this example includes:

[0153] A parameter acquisition unit configured to obtain the viscoelastic anisotropic medium parameters of the shale reservoir;

[0154] A stiffness matrix calculation unit configured to calculate the stiffness matrix in the frequency domain based on the Kelvin-Voigt viscoelastic medium model through the viscoelastic anisotropic medium parameters of the shale reservoir;

[0155] An equation construction unit configured to construct a guided wave dispersion equation based on the stiffness matrix in the frequency domain using the Thomsen-Haskell transfer matrix method;

[0156] An equation solving unit configured to solve and find the roots of the guided wave dispersion equation to obtain the guided wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium, where the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function for different parameters.

[0157] Embodiment 3: This embodiment provides an electronic device corresponding to the method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir provided in Embodiment 1. The electronic device can be an electronic device for a client, such as a mobile phone, a laptop computer, a tablet computer, a desktop computer, etc., to execute the method of Embodiment 1.

[0158] As Figure 12 shown, the electronic device includes a processor, a memory, a communication interface, and a bus. The processor, the memory, and the communication interface are connected through the bus to complete communication with each other. A computer program that can run on the processor is stored in the memory. When the processor runs the computer program, it executes the method of Embodiment 1. The implementation principle and technical effects are similar to those of Embodiment 1 and will not be elaborated here. Those skilled in the art can understand that Figure 12 the structure shown in

[0159] is only a block diagram of some structures related to the solution of this application, and does not constitute a limitation on the computing device to which the solution of this application is applied. The specific computing device may include more or fewer components than those shown in the figure, or combine certain components, or have a different component layout.

[0160] In a preferred embodiment, the above-mentioned logical instructions in the memory can be implemented in the form of a software functional unit and, when sold or used as an independent product, can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a part of this technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to enable a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in various embodiments of this application. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), and optical discs that can store program codes.

[0160] In a preferred embodiment, the processor can be a general-purpose processor of various types such as a central processing unit (CPU) and a digital signal processor (DSP), which is not limited here.

[0161] Embodiment 4: This embodiment provides a computer-readable storage medium storing one or more programs. The one or more programs include computer instructions that, when executed by a computer, cause the computer to execute the method provided in Embodiment 1 above.

[0162] In a preferred embodiment, the computer-readable storage medium may be a tangible device that holds and stores instructions executed by a computer, such as, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any combination of the foregoing. The computer-readable storage medium stores computer program instructions that cause a computer to execute the method provided in the first embodiment above.

[0163] This application is described with reference to the flowcharts and / or block diagrams of methods, apparatus (devices), and computer program products according to embodiments of the present application. It should be understood that each flow and / or block in the flowchart and / or block diagram, and the combination of flows and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, such that the instructions executed by the processor of the computer or other programmable data processing device generate means for implementing the functions specified in Figure 1 one or more of the flows Figure 1 or blocks or combinations of blocks.

[0164] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a particular manner, such that the instructions stored in the computer-readable memory generate a manufacture including instruction means that implement the functions specified in Figure 1 one or more of the flows Figure 1 or blocks or combinations of blocks.

[0165] These computer program instructions can also be loaded onto a computer or other programmable data processing device, such that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process, so that the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in Figure 1 one or more of the flows Figure 1 or blocks or combinations of blocks.

[0166] Each embodiment in this specification is described in a progressive manner. For the same or similar parts among the embodiments, reference can be made to each other, and the key point of each embodiment is to illustrate the differences from other embodiments. In the description of this specification, the descriptions with reference to terms such as "a preferred embodiment", "furthermore", "specifically", "in this embodiment", etc. mean that the specific features, structures, materials or characteristics described in connection with the embodiment or example are included in at least one embodiment or example of the embodiments of this specification. In this specification, the schematic expressions of the above terms do not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in a suitable manner in any one or more embodiments or examples. In addition, without contradiction, those skilled in the art can combine and combine the different embodiments or examples described in this specification and the features of different embodiments or examples.

[0167] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for quickly calculating the dispersion curve of guided waves in a viscoelastic anisotropic shale reservoir, characterized in that, Including: Obtaining the viscoelastic anisotropic medium parameters of the shale reservoir; Based on the Kelvin-Voigt viscoelastic medium model, inputting the viscoelastic anisotropic medium parameters of the shale reservoir to calculate the stiffness matrix in the frequency domain; Based on the stiffness matrix in the frequency domain, using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation; Solving and finding the roots of the guided wave dispersion equation to obtain the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve in the viscoelastic medium. Among them, the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve are used to calculate the sensitivity functions for different parameters.

2. The method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 1, characterized in that The viscoelastic anisotropic medium parameters of shale reservoirs include density ρ, vertical P-wave velocity V p0 and V s0 ; Thomsen anisotropic parameters ε, γ and δ; vertical P-wave and S-wave quality factors Q p0 and Q s0 ; quality factor anisotropic parameters ε Q and δ Q and reservoir thickness h.

3. The method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 2, wherein Each element in the stiffness matrix is: The element c in the stiffness matrix ij The real part of is used to represent the anisotropy of velocity, and the calculation formula is: In the formula, other elements are all zero except this; Element c in the stiffness matrix ij The imaginary part Q ij , which is used to represent the anisotropy of the attenuation coefficient, and the calculation formula is as follows: Q 33 = Q p0 ; Q 55 = Q s0 ; In the formula, the imaginary part Q of the element c in the stiffness matrix ij is similar to the symmetry of the real part of the element c in the stiffness matrix ij and the other elements are all zero. ij The real part ​ 4. The method for rapidly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 1, characterized in that Based on the stiffness matrix in the frequency domain, using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation, including: substituting the plane wave solution into the elastic wave equation to obtain the matrix expression of the guided wave solution function; obtaining the transfer matrix from the continuity of displacements and stresses at each interface; applying the boundary condition that the displacement field is zero at infinity to the upper and lower semi-infinite spaces; using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation based on the matrix expression, the transfer matrix and the boundary condition.

5. The method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 4, wherein The guided wave dispersion curve is obtained by solving the following dispersion equation: D(ω, c r , α r ) = det|U T T1T2...T n-1 T n V| = 0; where ω, c r and α are the angular frequency, phase velocity and attenuation coefficient of the guided wave, respectively, and U T represents the boundary condition of the guided wave in the upper half-infinite space, and T1T2...T n-1 T n represents the transfer matrix of the intermediate layer, and V represents the boundary condition of the lower half-infinite space. Among them, the guided wave propagating in the layered medium needs to satisfy the boundary conditions: in the upper half-infinite space, the boundary condition is that the displacement field is zero at infinity; in the intermediate low-velocity layer, the displacement field and stress field are continuous at the interfaces of different layers; in the lower half-infinite space, the boundary condition is that the displacement field is zero at infinity.

6. The method for quickly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 5, characterized in that Using the fast recursive method to solve and find the roots of the guided wave dispersion equation to obtain the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve in the viscoelastic medium, that is, traversing all values in the three directions of the phase velocity, frequency and attenuation coefficient to find the points that satisfy the dispersion equation equal to zero. The specific process is: Fixing the frequency ω0, transforming the root finding in the three-dimensional space of the phase velocity c, frequency ω0 and attenuation coefficient α into the root finding in the two-dimensional plane of the phase velocity c and the attenuation coefficient α; Traverse all values of the phase velocity c and attenuation coefficient α in the plane. At this time, the roots of the dispersion equation satisfy the local minimum. Complete the root finding for all modes at the frequency ω0, and then perform the root finding for the next frequency ω0 + Δω until the root finding for all frequencies is completed. Use a recursive method to set the initial step sizes Δc and Δα, and then give a new root finding range [c i -Δc, c i +Δc] and [α i -Δα, α i +Δα] according to the position of the zero point. Continue to refine the root finding range with the new step sizes Δc / n and Δα / n. The number of recursions satisfies the set conditions, and iterate the above process until the required root finding accuracy is met; After solving all the zeros at different frequencies, solving the relationship between the guided wave phase velocity and frequency and the relationship between the guided wave attenuation coefficient and frequency from the dispersion equation, that is, obtaining the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve in the viscoelastic medium.

7. The method for rapidly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to claim 2, characterized in that, Calculating the sensitivity functions for different parameters is based on the finite difference method, and the calculation formulas are respectively: Among them, is a sensitivity function of the phase velocity to different parameters, S α is a sensitivity function of the attenuation coefficient to different parameters, f is the frequency, m is the model parameter, and Δm is the perturbation of the model parameter. Among them, the calculation formula for specific parameters is as follows: Among them, is a sensitive function of the guided wave phase velocity relative to the shear wave velocity, S ε is a sensitive function of the guided wave phase velocity relative to the anisotropy parameter ε, S δ is a sensitive function of the guided wave phase velocity relative to the anisotropy parameter δ, c is the guided wave phase velocity, ΔV s0 , Δε and Δδ represent model perturbations; is a sensitive function of the guided wave attenuation coefficient relative to the shear wave quality factor, is a sensitive function of the guided wave attenuation coefficient relative to the longitudinal wave quality factor, and are sensitive functions of the guided wave attenuation coefficient relative to the anisotropy parameter of the quality factor, α is the guided wave attenuation coefficient, ΔQ s0 , ΔQ p0 , Δε Q and Δδ Q represent model parameter perturbations.

8. An apparatus for quickly calculating the dispersion curve of guided waves in a viscoelastic anisotropic shale reservoir, characterized in that, Including: A parameter acquisition unit configured to obtain the viscoelastic anisotropic medium parameters of the shale reservoir; A stiffness matrix calculation unit configured to calculate the stiffness matrix in the frequency domain based on the Kelvin-Voigt viscoelastic medium model through the viscoelastic anisotropic medium parameters of the shale reservoir; An equation construction unit configured to construct the guided wave dispersion equation based on the stiffness matrix in the frequency domain using the Thomsen-Haskell transfer matrix method; An equation solving unit configured to solve and find the roots of the guided wave dispersion equation to obtain the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve in the viscoelastic medium. Among them, the dispersion curve of the guided wave phase velocity and the attenuation coefficient curve are used to calculate the sensitivity functions for different parameters.

9. An electronic device, characterized in that, Including: At least one processor; And a memory communicatively connected to the processor; wherein, the memory stores instructions executable by the processor, and the instructions are executed by the processor so that the processor can execute the method according to any one of claims 1-7.

10. A computer-readable storage medium storing one or more programs, characterized in that, The one or more programs include computer instructions for causing a computer to perform the method according to any one of claims 1-7.

Images