Method and device for fast calculating dispersion curve of waveguide in viscoelastic anisotropic shale reservoir, equipment and medium

Abstract

The disclosed method, device, equipment and medium for quickly calculating the dispersion curve of a viscoelastic anisotropic shale reservoir waveguide comprise: obtaining viscoelastic anisotropic medium parameters of a shale reservoir; calculating a frequency domain stiffness matrix based on a Kelvin-Voigt viscoelastic medium model and the viscoelastic anisotropic medium parameters of the shale reservoir; constructing a waveguide dispersion equation based on the frequency domain stiffness matrix and using a Thomsen-Haskell transfer matrix method; and solving the waveguide dispersion equation by using a fast recursive method to find roots and obtain a waveguide phase velocity dispersion curve and an attenuation coefficient curve in the viscoelastic medium, wherein the waveguide phase velocity dispersion curve and the attenuation coefficient curve are used to calculate a sensitivity function relative to different parameters. Therefore, the viscoelasticity and anisotropy of the medium are combined, which can effectively reduce errors in the description of the shale reservoir, thereby improving the accuracy of reservoir physical property parameter inversion based on waveguide dispersion properties.

Description

TECHNICAL FIELD

[0001] The application relates to a method and device for quickly calculating a dispersion curve of a viscoelastic anisotropic shale reservoir guided wave, an equipment and a medium, and relates to the technical field of geological exploration, in particular to the technical field of seismic wave forward and inversion. BACKGROUND

[0002] In recent years, unconventional oil and gas has played a key role in clean energy and is considered an important alternative 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. Hydraulic fracturing of shale oil and gas reservoirs needs to monitor and evaluate its modification efficiency and safety in real time. The physical properties of shale reservoirs are crucial for geomechanical modeling, fracture imaging, and hydraulic fracturing design optimization. Shale layers are usually thin, with a thickness ranging from a few meters to tens 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] The low-velocity shale layer in the shale reservoir is located between the high-velocity surrounding rock, forming a special waveguide structure. The P and S waves excited by hydraulic fracturing perforation or microseismic sources are limited to propagate within the low-velocity shale layer, constantly reflected, superimposed, and interfered at the top and bottom interfaces of the waveguide structure, thereby forming guided waves. In the vertical plane, P-SV type guided waves are formed by the interference of P and SV waves, and because of the similar physical properties to Rayleigh surface waves, they are also called Rayleigh type guided waves; in the horizontal plane, SH type guided waves are formed by the interference of SH waves, and because of the similar physical properties to Love surface waves, they are also called Love type guided waves. Guided waves have obvious dispersion phenomenon, and dispersion, as a kinematic characteristic, 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 the structure of the shale reservoir. In addition, guided waves have a wide range of applications in the field of geophysics. For example: the guided waves (also known as 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 rock form a waveguide structure, and the seismic waves propagating inside form a special seismic phase-fault trapping wave (or fault zone guided wave) which can be used for fault zone structure imaging.

[0004] The shale has different elastic properties in different directions, thus showing obvious anisotropy, which will affect the dispersion characteristics of the guided waves propagating therein. Meanwhile, the shale shows viscoelasticity, and the energy of the guided waves will be lost during the propagation, thus causing the amplitude to attenuate. The assumption of the fully elastic isotropic medium in the traditional method is a too simplified assumption, which cannot accurately describe the propagation characteristics of the guided waves in the actual formation. In the elastic medium or the viscoelastic medium, the guided waves have phase velocity dispersion, and the anisotropy of the medium will significantly affect the phase velocity dispersion curve of the guided waves, and the viscoelasticity of the medium will also have certain influence on the phase velocity dispersion curve of the guided waves. In the viscoelastic medium, there is another dispersion relationship between the frequency and the attenuation coefficient, which can be expressed as the attenuation dispersion curve of the guided waves, and the viscoelasticity of the medium will significantly affect the attenuation dispersion curve of the guided waves. Therefore, in order to accurately and comprehensively describe the properties of the guided waves, the viscoelasticity and the anisotropy need to be considered. However, the current research on the dispersion properties of the guided waves mainly focuses on the isotropic fully elastic medium. Therefore, how to forward calculate the dispersion curve of the guided waves in the viscoelastic anisotropic medium is a problem to be solved at present, which can provide more reliable scientific basis for the exploration and development of shale gas. SUMMARY

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

[0006] In order to achieve the above application purpose, the technical scheme adopted by the present application is as follows:

[0007] In the first aspect, the present application provides a method for quickly calculating the dispersion curve of the guided waves in the viscoelastic anisotropic shale reservoir, which comprises the following steps:

[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 dispersion equation of the guided waves;

[0011] Solving the dispersion equation of the guided waves to obtain the phase velocity dispersion curve and the attenuation coefficient curve of the guided waves in the viscoelastic medium, wherein the phase velocity dispersion curve and the attenuation coefficient curve of the guided waves are used to calculate the sensitivity function with respect to different parameters.

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

[0013] In some possible implementations, each element in the stiffness matrix is:

[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 zero except for the element;

[0018] The imaginary part Q ij of the element c ij in the stiffness matrix is used to represent the anisotropy of 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 symmetry of the real part of the element c ij in the stiffness matrix, and other elements are also zero.

[0023] In some possible implementations, based on the frequency domain stiffness matrix, the Thomsen-Haskell transfer matrix method is used to construct the guided wave dispersion equation, including: the plane wave solution is brought into the elastic wave equation to obtain a matrix expression of the guided wave solution function; the transfer matrix is obtained by the displacement and stress continuity of each interface; the boundary condition of zero displacement field at infinity is applied to the upper and lower half-infinite space; 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.

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

[0025]

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

[0027] In some possible implementations, the dispersion equation of the guided wave is solved by using the fast recursive method to find the root, and the phase velocity dispersion curve and the attenuation coefficient curve of the guided wave in the viscoelastic medium are obtained, that is, all values in the three directions of phase velocity, frequency and attenuation coefficient are traversed to find the point satisfying the dispersion equation equal to zero, and the specific process is as follows:

[0028] Fixing the frequency ω0, the root seeking in the three-dimensional space of the phase velocity c, the frequency ω0 and the attenuation coefficient a is converted into the root seeking in the two-dimensional plane of the phase velocity c and the attenuation coefficient a.

[0029] All values of the phase velocity c and the attenuation coefficient a in the plane are traversed, at this time the root of the dispersion equation satisfies the local minimum value, the root seeking of all modes under the frequency ω0 is completed, the root seeking of the next frequency ω0+Δω is performed until the root seeking of all frequencies is completed, the initial step length Δc and Δa are set by using the recursive method, and then the new root seeking range [c i -Δc, c i +Δc] and [a i -Δa, a i +Δa] is given according to the position of the zero point, the new step length Δc / n and Δa / n is used to continue to accurately seek the root seeking range, the number of recursions satisfies the set condition, and the above process is iterated until the required root seeking accuracy is satisfied.

[0030] After all the zero points under all frequencies are solved, the relationship between the phase velocity of the guided wave and the frequency and the relationship between the attenuation coefficient of the guided wave and the frequency in the viscoelastic medium are solved from the dispersion equation, that is, the phase velocity dispersion curve and the attenuation coefficient curve of the guided wave in the viscoelastic medium are obtained.

[0031] In some possible implementations, the sensitivity function with respect to different parameters is calculated by using the finite difference method, and the calculation formulas are as follows:

[0032]

[0033] 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, where the calculation formula of a specific parameter is:

[0034]

[0035] where, is the sensitivity function of the guided wave phase velocity to the transverse wave velocity, S ε is the sensitivity function of the guided wave phase velocity to the anisotropy parameter ε, S δ is the sensitivity function of the guided wave phase velocity to the anisotropy parameter δ, c is the guided wave phase velocity, ΔV s0 , Δε, and Δδ represent the model perturbation. is the sensitivity function of the guided wave attenuation coefficient to the transverse wave quality factor, S is the sensitivity function of the guided wave attenuation coefficient to the longitudinal wave quality factor, S and is the sensitivity function of the guided wave attenuation coefficient to the quality factor anisotropy parameter, α is the guided wave attenuation coefficient, ΔQ s0 , ΔQ p0 , Δε Q , and Δδ Q represent the model parameter perturbation.

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

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

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

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

[0040] an equation solving unit configured to solve the guided wave dispersion equation by root-finding to obtain a guided wave phase velocity dispersion curve and an attenuation coefficient curve in a viscoelastic medium, wherein the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function to different parameters.

[0041] In a third aspect, the present application provides an electronic device, comprising: at least one processor; and a memory connected with the processor in communication; wherein the memory stores instructions executable by the processor, and the instructions are executed by the processor to enable the processor to perform the method.

[0042] In a fourth aspect, the present application provides a computer readable storage medium storing one or more programs, the one or more programs comprising computer instructions for causing a computer to perform the method.

[0043] The present application has the following characteristics due to the above technical solutions:

[0044] 1. Compared with the traditional viscoelastic isotropic forward method and elastic anisotropic forward, the present application can more accurately reflect the information of the reservoir based on the viscoelastic anisotropic theory. The parameter quality factor reflecting the wave amplitude attenuation cannot be obtained based on the forward theory of the complete elastic medium, and the velocity information of the medium cannot be accurately reflected based on the forward theory of the isotropic medium. Therefore, the combination of the viscoelasticity and anisotropy of the medium can effectively reduce the error in the description of the shale reservoir, thereby improving the accuracy of the inversion of the reservoir physical property parameters based on the guided wave dispersion property, and further helping to deeply understand the mechanical properties, stress state and fracture distribution of the reservoir, so as to improve the accuracy of the reservoir evaluation and the rationality of the development decision.

[0045] 2. The present application can more efficiently and stably solve the root-seeking equation by using the fast recursive method, and even in the high-frequency high-order part, the root missing and step skipping situations will not occur. The idea of recursion is introduced in the process of root-seeking, so that the calculation efficiency of root-seeking is guaranteed.

[0046] 3. The present application can calculate the sensitivity function of the Rayleigh-type guided wave phase velocity dispersion curve to the vertical transverse wave velocity, Thomsen anisotropy parameters epsilon and delta by using the numerical solution method; the sensitivity function of the Rayleigh-type guided wave attenuation coefficient curve to the vertical longitudinal and transverse wave velocity quality factor, quality factor anisotropy parameters epsilon and delta, and the sensitivity analysis of the above parameters provides support for the parameter inversion in the viscoelastic anisotropic medium. Q Q

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

[0048] Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The accompanying drawings are included to provide a description of the preferred embodiments and are not intended to limit the scope of the application. Throughout the drawings, like reference numerals will be used to refer to like components. In the drawings:​​

[0049] Figure 1 This is a flowchart of a method for rapidly calculating the guided wave dispersion curve of a viscoelastic anisotropic shale reservoir according to an embodiment of the present invention;

[0050] Figure 2 This is a schematic diagram illustrating the fast recursive method for solving the guided wave dispersion equation according to an embodiment of the present invention.

[0051] Figure 3 Forward modeling of the dispersion curve of the phase velocity of the guided wave in the viscoelastic anisotropic medium model of this invention;

[0052] Figure 4 Forward modeling of the waveguide attenuation coefficient curve in the viscoelastic anisotropic medium model of this invention;

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

[0054] Figure 6 Sensitivity analysis of waveguide phase velocity to Thomsen anisotropy parameter ε in embodiments of the present invention;

[0055] Figure 7 Sensitivity analysis of waveguide phase velocity to Thomsen anisotropy parameter δ in embodiments of the present invention;

[0056] Figure 8 Sensitivity analysis of the waveguide attenuation coefficient to the vertical longitudinal wave velocity quality factor in this embodiment of the invention;

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

[0058] Figure 10 The waveguide attenuation coefficient in this embodiment of the invention relates to the quality factor anisotropy parameter ε. Q Sensitivity analysis;

[0059] Figure 11 The waveguide attenuation coefficient in this embodiment of the invention relates to the quality factor anisotropy parameter δ. Q Sensitivity analysis;

[0060] Figure 12 This is a structural diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0061] It is to be understood that the terminology used herein is for the purpose of describing particular example embodiments only and is not intended to be limiting. As used herein, the singular forms "a", "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. The terms "comprises", "comprising", "includes", "including" and "has" are inclusive and therefore 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 groups thereof. The method steps, processes, and operations described herein are not to be construed as necessarily requiring their performance in the particular order

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

[0063] Spatially relative terms, such as "inner", "outer", "beneath", "below", "lower", "above", "upper", and the like, can be used herein for ease of description to describe one element or feature's relationship to another element(s) or feature(s) as illustrated in the figures. The spatially relative terms are intended to encompass different orientations of the device in use or operation in addition to the orientations depicted in the figures.

[0064] Since the current research on the dispersion properties of guided waves mainly focuses on isotropic full elastic medium, the properties of guided waves cannot be accurately and comprehensively described. The method, device, equipment and medium for quickly calculating the dispersion curve of guided waves of viscoelastic anisotropic (Vertical Transversely Isotropic, VTI) shale reservoir provided by the present application, comprising: 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 frequency domain stiffness matrix; based on the frequency domain stiffness matrix, using the Thomsen-Haskell transfer matrix method to construct the guided wave dispersion equation; solving the guided wave dispersion equation by root-finding to obtain the guided wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium, wherein the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function of different parameters. Therefore, the present application derives the analytical expression of the viscoelastic VTI medium guided wave dispersion curve based on the Kelvin-Voigt viscoelastic medium model and the Thomsen-Haskell transfer matrix method, calculates the sensitivity function of the guided wave phase velocity dispersion curve relative to the vertical transverse wave velocity, Thomsen parameters ε and δ, and calculates the sensitivity function of the guided wave attenuation coefficient curve relative to the vertical transverse wave velocity quality factor, vertical longitudinal wave velocity quality factor, attenuation anisotropy parameters ε and δ, which can improve the accuracy of the reservoir physical property parameter inversion based on the guided wave dispersion properties. Q and δ Q

[0065] Exemplary embodiments of the present application will be described in greater detail below, with reference to the accompanying drawings. While exemplary embodiments of the present application are shown in the drawings, it is understood that the present application can be embodied in various forms and should not be limited by the embodiments set forth herein. Rather, these embodiments are provided so that this application will be thorough and complete, and will fully convey the scope of the application to those skilled in the art.

[0066] Embodiment one: as shown in the embodiment, the method for quickly calculating the dispersion curve of guided waves of viscoelastic anisotropic shale reservoir provided by the present application, comprising: Figure 1

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

[0068] In the present 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 parameters ε q ​​and δ Q and reservoir thickness h.

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

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

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

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

[0073] wherein σ=(σ 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 represent the displacement fields along the x, y and z directions respectively.

[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 velocity V p0 and S-wave velocity V S0 , Thomson anisotropy parameters ε, δ and γ, which are used to represent the anisotropy of velocity, as follows:

[0077]

[0078]

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

[0080] The imaginary part Q ij of the stiffness matrix element c ij can be calculated from the vertical P-wave quality factor Q p0 and S-wave quality factor Q s0 , quality factor anisotropy parameters ε Q , δQ and γ Q The anisotropy of the attenuation coefficient is calculated and expressed as follows:

[0081] Q 33 = Q p0 ;

[0082] Q 55 = Q s0 ;

[0083]

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

[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 brought into the elastic wave equation to obtain the matrix expression of the guided wave solution function; the transfer matrix is obtained from the displacement and stress continuity of each interface; the boundary condition of zero displacement field at infinity is applied to the upper and lower half-infinite space; 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] Further, the elastic wave equation is composed of the motion equation and the constitutive equation, and 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 mutual interference of P waves and SV waves, and only the wave propagating along the x-axis needs 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] The above conditions are substituted into the constitutive equation and the motion equation to obtain the following wave equation:

[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 neglected. The solution of the displacement field and stress field in the frequency domain in the above elastic wave equation is 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] Wherein, δ is the Dirac function, y = (y1, y2, y3, y4) T is a function to be solved, and the guided wave wavenumber is extended to the complex domain in the viscoelastic medium, as follows:

[0103]

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

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

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

[0107] Wherein:

[0108]

[0109] Wherein, F1, F2, F3, F4 represent the amplitude, and the elements of the matrix D(z) are calculated by the following formula:

[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] where, υ1, υ 2、 υ3, υ4 are eigenvalues, representing the vertical wave numbers of the uplink and downlink longitudinal and transverse waves, 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] where, ε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 stress field are continuous at the interface 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] Wherein, U T represents the boundary condition of the guided wave in the upper half infinite space, T1T2…T n-1 T n represents the transfer matrix of the intermediate layer, 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, the fast recursive method is used to solve and seek the root of the guided wave dispersion equation.

[0132] In the 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 the three-dimensional space (frequency, phase velocity, attenuation coefficient), and a new root-seeking method needs to be used to solve the dispersion curve of the guided wave in the viscoelastic medium.

[0133] Further, the method used by the present application is a fast recursive dispersion equation root-seeking algorithm, all values in the three directions of phase velocity, frequency and attenuation coefficient are traversed to find the point satisfying the dispersion equation equal to zero, and the efficiency and accuracy are improved through the recursive algorithm, and the specific calculation method is:

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

[0135] The second step is to traverse all values of the phase velocity c and the attenuation coefficient α in the plane, at this time the root of the dispersion equation satisfies the local minimum value, the root-seeking of all modes under the frequency ω0 is completed, and the root-seeking of the next frequency ω0+Δω is carried out until the root-seeking of all frequencies is completed. In order to improve the efficiency and accuracy of the root-seeking, the present application adopts the recursive method, the initial step length Δc and Δα can be set to a larger value, and then a new root-seeking range [c i -Δc,c i +Δc] and [α i -Δα,α i +Δα] is given according to the position of the zero point, the new step length Δc / n and Δα / n is used to continue to accurately seek the root-seeking range, the recursive number is about 10 times, and the above process is iterated until the required root-seeking accuracy is met.

[0136] When all the zero points at different frequencies are solved, the relationship between the phase velocity of the guided wave in the viscoelastic medium and the frequency and the relationship between the attenuation coefficient of the guided wave and the frequency can be solved from the dispersion equation.

[0137] S5, the sensitivity function of the phase velocity dispersion curve and the attenuation coefficient curve of the guided wave relative to different parameters is calculated based on the finite difference.

[0138] In this embodiment, the sensitivity of the phase velocity and the attenuation coefficient of the guided wave in the viscoelastic medium to different parameters can be further analyzed by means of the phase velocity dispersion curve and the attenuation coefficient curve of the guided wave in the viscoelastic anisotropic medium. The present application adopts the finite difference method to approximate the partial derivative of the phase velocity dispersion curve and the attenuation coefficient curve of the guided wave relative to different parameters, that is, the sensitivity function:

[0139]

[0140] wherein, 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] Further, after the sensitivity analysis of all the parameters affecting the phase velocity dispersion and the attenuation coefficient curve of the guided wave, the parameters that obviously affect the two curves are obtained. The parameters sensitive to the phase velocity dispersion curve of the guided wave are the vertical transverse wave velocity, the Thomsen anisotropy parameters ε and δ; the parameters sensitive to the attenuation coefficient curve of the guided wave are the vertical longitudinal wave velocity quality factor, the vertical transverse wave velocity quality factor, and the quality factor anisotropy parameters ε 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 application of the method for quickly calculating the guided wave dispersion curve of the viscoelastic anisotropic shale reservoir will be described in detail through specific embodiments.

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

[0144]

[0145] Based on the model parameters in Table 1, the Kelvin-Voigt viscoelastic medium model is used to calculate the stiffness matrix C in the frequency domain by the viscoelastic anisotropic medium parameters of the shale reservoir, the Thomsen-Haskell transfer matrix method is used to construct the guided wave dispersion equation, the fast recursive method is used to solve the guided wave dispersion equation, and the root-finding is used to obtain the phase velocity dispersion curve of the guided wave as shown in the following table:Figure 3 The attenuation coefficient curves of the guided wave are shown below. Figure 4 As shown. For Figure 3 Sensitivity analysis was performed on the fundamental order of the guided wave phase velocity dispersion curve to analyze its sensitivity to the vertical shear wave velocity, Thomsen parameters ε and δ, using the following difference formula:

[0146]

[0147]

[0148] in, It is a sensitive function of the guided wave phase velocity relative to the shear wave velocity, S ε It is a sensitive function of the waveguide phase velocity relative to the anisotropy parameter ε, S δ It is a sensitivity function of the waveguide phase velocity relative to the anisotropy parameter δ, where c is the waveguide phase velocity, f is the frequency, and ΔV is the waveguide phase velocity. s0 Δε and Δδ represent model perturbations and are used to approximate partial derivatives. The sensitivity analysis regarding the vertical shear wave velocity is obtained as follows: Figure 5 As shown, the sensitivity analysis of the Thomsen parameters ε and δ is as follows: Figure 6 and Figure 7 As shown. The waveguide phase velocity versus the transverse wave velocity structure V. s0 It is the most sensitive, exhibiting a certain sensitivity to anisotropy parameters ε and δ, and therefore can be used for the inversion of anisotropy parameters in shale reservoirs. Figure 4 Sensitivity analysis was performed on the fundamental order of the dispersion curve of the mid-guide wave attenuation coefficient. The analysis showed that the dispersion curve affected the vertical shear wave velocity quality factor, the vertical longitudinal wave velocity quality factor, and the Thomsen parameter ε describing the attenuation. Q and δ Q The sensitivity is determined using the following difference formula:

[0149]

[0150] in, It is a sensitive function of the guided wave attenuation coefficient relative to the shear wave quality factor. It is a sensitive function of the guided wave attenuation coefficient relative to the longitudinal wave quality factor. and It is a sensitivity function of the waveguide attenuation coefficient relative to the quality factor anisotropy parameter, where α is the waveguide attenuation coefficient, f is the frequency, and ΔQ is the waveguide attenuation coefficient. s0 ΔQ p0 , Δε Q and Δδ Q This represents the perturbation of the model parameters. The sensitivity analysis of the vertical P-wave velocity quality factor is obtained as follows: Figure 8 As shown, the sensitivity analysis of the vertical shear wave velocity quality factor is as follows: Figure 9The illustrated, quality factor anisotropy parameter ε Q and δ Q The sensitivity analysis as shown in Figure 10 and Figure 11 The wave attenuation coefficient is most sensitive to the quality factor parameter in the low-velocity layer, and the wave attenuation coefficient curve and the sensitivity function can provide a basis for the research and inversion of the viscoelastic properties of the shale reservoir.

[0151] Embodiment two: the above embodiment one provides a method for quickly calculating the viscoelastic anisotropic shale reservoir wave dispersion curve, and correspondingly, the present embodiment provides a device for quickly calculating the viscoelastic anisotropic shale reservoir wave dispersion curve. The device provided by the present embodiment can implement the method for quickly calculating the viscoelastic anisotropic shale reservoir wave dispersion curve of embodiment one. The device can be realized by software, hardware or a combination of software and hardware. For the convenience of description, the device is described as various units in the function. Of course, in the implementation, the functions of the units can be realized in the same or multiple software and / or hardware. For example, the device can include integrated or separate functional modules or functional units to perform the corresponding steps in the method of embodiment one. Since the device of the present embodiment is basically similar to the method embodiment, the description process of the present embodiment is relatively simple, and the related parts can be referred to the part of the description of embodiment one. The embodiment of the device for quickly calculating the viscoelastic anisotropic shale reservoir wave dispersion curve provided by the present application is only illustrative.

[0152] Specifically, the device for quickly calculating the viscoelastic anisotropic shale reservoir wave dispersion curve provided by the present embodiment comprises:

[0153] The parameter acquisition unit is configured to obtain the viscoelastic anisotropic medium parameters of the shale reservoir.

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

[0155] The equation construction unit is configured to construct the wave dispersion equation based on the frequency domain stiffness matrix by using the Thomsen-Haskell transfer matrix method.

[0156] The equation solving unit is configured to solve the wave dispersion equation to obtain the wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium, wherein the wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function of different parameters.

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

[0158] As shown in Figure 12 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 the communication between each other. The memory stores a computer program that can run on the processor, and the processor runs the computer program to execute the method of the embodiment one, which has similar principles and technical effects to the embodiment one, and will not be repeated here. Those skilled in the art can understand that Figure 12 The structure shown in the figure is only a block diagram of part of the structure related to the scheme of the present application, and does not constitute a limitation on the computing device to which the scheme of the present application is applied. The specific computing device can include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0159] In a preferred embodiment, the logical instructions in the memory described above can be realized in the form of a software function unit and sold or used as an independent product when used, which can be stored in a computer readable storage medium. Based on such understanding, the technical scheme of the present application essentially or the part that contributes to the prior art or part of the technical scheme can be embodied in the form of a software product, and the computer software product is stored in a storage medium, including a plurality of instructions for causing 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 method described in the embodiments of the present application. The foregoing storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM, Read-Only Memory), a random access memory (RAM, Random Access Memory), an optical disc and various program code storage media.

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

[0161] Embodiment four: the embodiment provides a computer readable storage medium storing one or more programs, the one or more programs including computer instructions, the computer instructions causing a computer to execute the method provided in the above embodiment one when executed by the computer.

[0162] In one preferred embodiment, the computer readable storage medium can be a tangible device that holds and stores a computer program that is executed by a processor, such as, but not limited to, an electronic storage device, a magnetic storage device, an optical storage device, a solid-state storage device, or any combination thereof. The computer readable storage medium stores the computer program instructions that cause a computer to execute the methods provided by the embodiments described above.

[0163] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing the steps in one or more flowcharts and / or 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 apparatus to function in a particular manner, such that the instructions stored in the computer readable memory produce an article of manufacture including instructions which implement the steps specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing the steps in one or more flowcharts and / or blocks

[0165] The computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart block or blocks. Figure 1 one or more flowcharts and / or blocks Figure 1 means for functionally implementing the steps in one or more flowcharts and / or blocks

[0166] Each of the embodiments in the specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other. Each of the embodiments focuses on the difference from other embodiments. In the description of the specification, the description of the terms "one preferred embodiment", "further", "specifically", "in this embodiment", and the like means 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 the specification. In the specification, the illustrative description of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any appropriate manner in any one or more embodiments or examples. Furthermore, the person skilled in the art can combine and combine the different embodiments or examples described in the specification and the features of the different embodiments or examples without contradiction.

[0167] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part 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 application.

Claims

1. A method for fast computation of viscoelastic anisotropic shale reservoir guided wave dispersion curves, characterized in that, The method comprises the following steps: obtaining shale reservoir viscoelastic anisotropic medium parameters; calculating a frequency domain stiffness matrix based on the Kelvin-Voigt viscoelastic medium model and the shale reservoir viscoelastic anisotropic medium parameters; constructing a guided wave dispersion equation based on the frequency domain stiffness matrix and the Thomsen-Haskell transfer matrix method; solving the guided wave dispersion equation by root-finding to obtain a guided wave phase velocity dispersion curve and an attenuation coefficient curve in the viscoelastic medium, wherein the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function with respect to different parameters.

2. The method for fast computation of viscoelastic anisotropic shale reservoir guided wave dispersion curves of claim 1, wherein, Shale reservoir viscoelastic anisotropic medium parameters include density , vertical P-S wave velocity and ; Thomsen anisotropy parameters , ; vertical P-S wave quality factors and ; quality factor anisotropy parameters and , and reservoir thickness .

3. The method for fast computation of viscoelastic anisotropic shale reservoir guided wave dispersion curves of claim 2, wherein, Each element in the stiffness matrix is: ； Real part of element in stiffness matrix to express the anisotropy of velocity, the calculation formula is:​ ； ； ； ； ； ； ； ； ； ； ； ； wherein all other elements are zero except for the element; imaginary parts of elements in stiffness matrix anisotropy of the attenuation coefficient, calculated as​ ； ； ； ； ； where the imaginary part of the element of the stiffness matrix is similar to the real part of the element of the stiffness matrix and the other elements are zero, , and are the quality factor anisotropy parameters.​​ 4. The method for fast calculation of viscoelastic anisotropic shale reservoir guide wave dispersion curves of claim 1, wherein, based on the frequency domain stiffness matrix, the Thomsen-Haskell transfer matrix method is used to construct the guided wave dispersion equation, including: the plane wave solution is brought into the elastic wave equation to obtain a matrix expression of the guided wave solution function; the transfer matrix is obtained by the displacement and stress continuity of each interface; the boundary condition of zero displacement field at infinity is applied to the upper and lower half-infinite space; 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.

5. The method for fast calculation of viscoelastic anisotropic shale reservoir guided wave dispersion curves of claim 4, wherein, The guided wave dispersion curve is obtained by solving the following dispersion equation: ； where ω is the angular frequency of the guided wave, , and are the angular frequency, phase velocity and attenuation coefficient of the guided wave, respectively, represents the boundary condition of the guided wave in the upper half infinite space, represents the transfer matrix of the intermediate layer, represents the boundary condition of the lower half infinite space, wherein the guided wave propagating in the layered medium needs to satisfy the boundary condition: 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 the stress field are continuous at the interface 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 fast calculation of viscoelastic anisotropic shale reservoir guided wave dispersion curves of claim 5, wherein, The guided wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium are obtained by solving the guided wave dispersion equation by root-finding using the fast recursive method, that is, all values in the three directions of phase velocity, frequency and attenuation coefficient are traversed to find the point that satisfies the dispersion equation equal to zero, and the specific process is as follows: Fixed frequency The root-finding in the three-dimensional space of phase velocity, frequency and attenuation coefficient is converted to the root-finding in the two-dimensional plane of phase velocity and attenuation coefficient. By iterating through all values ​​of the phase velocity and attenuation coefficient within the plane, the roots of the dispersion equation satisfy a local minimum, thus completing the frequency determination. Root-finding of all modes below, proceeding to the next frequency. The root search continues until all frequencies are found, using a recursive method to set the initial step size. and Then, based on the position of the zero point, a new root-finding range is given. and , with new step size and Continue to precisely search for the root range, and if the number of recursions meets the set condition, iterate the above process until the required root-searching accuracy is met; After all the zero points at different frequencies are solved, the relationship between the guided wave phase velocity and the frequency and the relationship between the guided wave attenuation coefficient and the frequency in the viscoelastic medium are solved from the dispersion equation, that is, the guided wave phase velocity dispersion curve and the attenuation coefficient curve in the viscoelastic medium are obtained.

7. The method for fast calculation of viscoelastic anisotropic shale reservoir guided wave dispersion curves of claim 5, wherein, The sensitivity function with respect to different parameters is calculated based on the finite difference method, and the calculation formulas are as follows: ； ； wherein, is a sensitivity function of the phase velocity to different parameters, is a sensitivity function of the attenuation coefficient to different parameters, is a frequency, is a model parameter, is a perturbation of the model parameter, wherein the calculation formula of a specific particular parameter is: ； ； ； ； ； ； ； where is the sensitivity function of the guided wave phase velocity to the shear wave velocity, is the sensitivity function of the guided wave phase velocity to the anisotropy parameter , is the sensitivity function of the guided wave phase velocity to the anisotropy parameter , is the guided wave phase velocity, , and denote the model perturbation; is the sensitivity function of the guided wave attenuation coefficient to the shear wave quality factor, is the sensitivity function of the guided wave attenuation coefficient to the longitudinal wave quality factor, and are the sensitivity functions of the guided wave attenuation coefficient to the quality factor anisotropy parameters, is the guided wave attenuation coefficient, , , and represent the model parameter perturbation.

8. An apparatus for fast computation of viscoelastic anisotropic shale reservoir guided wave dispersion curves, characterized by, The method comprises the following steps: a parameter acquisition unit configured to obtain shale reservoir viscoelastic anisotropic medium parameters; a stiffness matrix calculation unit configured to calculate a frequency domain stiffness matrix based on the Kelvin-Voigt viscoelastic medium model and the shale reservoir viscoelastic anisotropic medium parameters; an equation construction unit configured to construct a guided wave dispersion equation based on the frequency domain stiffness matrix and the Thomsen-Haskell transfer matrix method; an equation solving unit configured to solve the guided wave dispersion equation by root-finding to obtain a guided wave phase velocity dispersion curve and an attenuation coefficient curve in the viscoelastic medium, wherein the guided wave phase velocity dispersion curve and the attenuation coefficient curve are used to calculate the sensitivity function with respect to different parameters.

9. An electronic device, comprising: The one or more programs include computer instructions for causing a computer to execute the method according to any one of claims 1-7. The one or more programs include computer instructions for causing a computer to execute the method according to any one of claims 1-7. ​ 10. A computer-readable storage medium storing one or more programs, the one or more programs comprising instructions for: ​

Images