A method for predicting inclined fractures and identifying oil and gas considering partial saturation attenuation mechanism

Abstract

The present invention relates to the field of oil and gas exploration and development, and specifically to a method for predicting inclined fractures and identifying oil and gas taking into account a partially saturated attenuation mechanism. The method comprises the following steps: Step 1, constructing an anisotropic stiffness matrix taking into account the attenuation mechanism of a partially saturated inclined fracture type reservoir; Step 2, constructing a frequency-variant reflection coefficient equation taking into account the attenuation mechanism of a partially saturated inclined fracture type reservoir; Step 3, conducting a five-dimensional seismic frequency-variant inversion of an inclined fracture type reservoir to achieve inclined fracture prediction and oil and gas identification. The reflection coefficient equation constructed by the present invention is more consistent with the reflection characteristics of an actual structural fracture type oil and gas reservoir, that is, the reflection characteristic characterization accuracy is higher, and the Fourier series decomposition strategy is used to separate the isotropic components and anisotropic components in the seismic data, which can effectively improve the accuracy of fracture prediction and oil and gas identification.

Description

Technical Field

[0001] The present invention relates to the field of oil and gas exploration and development, and specifically includes a method for predicting inclined fractures and identifying oil and gas by considering a partial saturation attenuation mechanism, so as to improve the accuracy of seismic prediction of fracture parameters and oil and gas distribution. Background Art

[0002] Oil and gas identification and fracture prediction are crucial components of oil and gas reservoir exploration and development. Five-dimensional seismic data contains rich azimuthal and frequency information. Inversion methods based on five-dimensional seismic data are an important approach for achieving high-precision predictions of fracture parameters and oil and gas distribution.

[0003] Currently, the prediction and fluid identification of different types of fracture groups (i.e., vertical single fracture groups, inclined single fracture groups, etc.) are primarily based on the amplitude information of five-dimensional data, using the AVAZ (Amplitude Versus Angle and Azimuth) inversion method. This method assumes that underground fluids are fully and evenly mixed and treats fractured reservoirs as equivalent elastic anisotropic media. In reality, underground reservoirs are characterized by heterogeneity and partial fluid saturation. Fluid flow in the formation can induce attenuation and dispersion of seismic waves. Ignoring the attenuation and dispersion caused by partial oil-water saturation in the fractures will reduce the accuracy of seismic prediction for fractured reservoirs.

[0004] Existing research results considering the partial saturation attenuation mechanism generally assume that the reservoir is isotropic, ignoring the anisotropic characteristics induced by directional fractures, and also reducing the accuracy of seismic prediction of fractured oil and gas reservoirs.

[0005] For example, the technical solution disclosed in the patent document with the patent name "A method, device, electronic device and storage medium for oil and gas detection" and application number: 202310363660.5 is to use the high and low frequency end information of the seismic signal to comprehensively characterize the magnitude of seismic attenuation, thereby improving the extraction accuracy of attenuation attributes and the accuracy of underground oil and gas reservoir prediction.

[0006] The technical solution disclosed in the patent document with the patent name "Method for detecting oil and gas reservoirs based on seismic attenuation intercept" and application number: 202110058904.X discloses a method for detecting oil and gas reservoirs based on seismic attenuation intercept. The attenuation intercept, in which the seismic amplitude decays with frequency, is used as a stable and reliable geophysical attribute for fluid detection and gas content prediction of the target oil and gas reservoir.

[0007] Both of the above technical solutions utilize post-stack data to achieve oil and gas monitoring by extracting the amplitude attenuation attributes in the post-stack data. They do not utilize five-dimensional seismic data and do not consider the anisotropic characteristics of the reservoir.

[0008] For example, the technical solution disclosed in the patent document with the patent name "A method for identifying the fluid saturation of carbonate reservoirs based on post-stack seismic data" and application number: 201910612941.3 discloses a method for identifying the fluid saturation of carbonate reservoirs based on post-stack seismic data.

[0009] This technical solution establishes a rock physics interpretation version and uses the attenuation properties of post-stack data to achieve oil and gas detection. It does not use five-dimensional seismic data and does not consider the anisotropic characteristics of the reservoir. Summary of the Invention

[0010] In order to solve these problems of the above-mentioned prior art, the present invention takes into account the attenuation mechanism of partially saturated reservoirs and the anisotropic characteristics of fractures, fully utilizes the azimuth and frequency information in five-dimensional seismic data, and establishes a method for inclined fracture prediction and oil and gas identification to improve the accuracy of reservoir fluid and fracture identification.

[0011] To solve the above technical problems, the present invention adopts a technical solution: a method for predicting inclined fractures and identifying oil and gas considering the partial saturation attenuation mechanism, which includes the following steps:

[0012] Step 1: construct an anisotropic stiffness matrix considering the attenuation mechanism of partially saturated inclined fractured reservoir;

[0013] Step 2: Construct a frequency-dependent reflection coefficient equation considering the attenuation mechanism of partially saturated inclined fractured reservoirs;

[0014] Step 3: Conduct five-dimensional seismic frequency-varying inversion of inclined fracture reservoirs to achieve inclined fracture prediction and oil and gas identification.

[0015] Preferably, the step 1 comprises:

[0016] For a rock with a set of vertical fractures, assuming that the fracture symmetry axis coincides with the x-axis, and considering the oscillation diffusion effect of partially saturated fluid in the fracture induced by wave motion, its complex stiffness matrix is ​​expressed as:

[0017] C H T I = [ M − M 2 m e U ˜ 3 3 l − l M m e U ˜ 3 3 l − l M m e U ˜ 3 3 0 0 0 l − l M m e U ˜ 3 3 M − l 2 m e U ˜ 3 3 l − l 2 m e U ˜ 3 3 0 0 0 l − l M m e U ˜ 3 3 l − l 2 m e U ˜ 3 3 M − l 2 m e U ˜ 3 3 0 0 0 0 0 0 m 0 0 0 0 0 0 m − m e U ˜ 1 1 0 0 0 0 0 0 m − m e U ˜ 1 1 ] (1),

[0018] (1) In the formula, Represents vertical fracture reservoir Stiffness matrix; is the compression modulus, ; and are the first Lamé constant and shear modulus for the isotropic background; is the crack density; and It is a parameter related to fracture parameters and fluid parameters;

[0019] For a reservoir with a set of inclined fractures, the stiffness matrix is ​​calculated by By rotating the coordinates, we get

[0020] (2),

[0021] (2) In the formula, The fourth-order stiffness tensor representing the inclined fracture reservoir, is the crack inclination angle, 、 、 、 is a second-order tensor The elements of i, j, k, l, s, t, p, q range from 1, 2, 3, and the second-order tensor for,

[0022] L = [ s i n g 0 − c o s g 0 1 0 c o s g 0 s i n g ] (3),

[0023] The fourth-order stiffness tensor representing the vertical fracture reservoir is Correspondingly, the fourth-order tensor element and Matrix elements The corresponding formula is as follows:

[0024] , , (4),

[0025] (4) In the formula, and is the Kronecker symbol;

[0026] Substituting (1) into (2), the stiffness matrix of the inclined fracture reservoir is obtained as follows:

[0027] C T T I = [ C 1 1 C 1 2 C 1 3 0 C 1 5 0 C 1 2 C 2 2 C 2 3 0 C 2 5 0 C 1 3 C 2 3 C 3 3 0 C 3 5 0 0 0 0 C 4 4 0 C 4 6 C 1 5 C 2 5 C 3 5 0 C 5 5 0 0 0 0 C 4 6 0 C 6 6 ] (5),

[0028] In formula (5),

[0029] C 1 1 = [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 4 g + [ 4 e ( M U ˜ 3 3 − m U ˜ 1 1 ) ] c o s 2 g + M ( 1 − M m e U ˜ 3 3 ) ,

[0030] ,

[0031] C 1 3 = [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 4 g + [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 2 g + l ( 1 − M m e U ˜ 3 3 ) ,

[0032] C 1 5 = { [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 2 g + 2 e ( M U ˜ 3 3 − m U ˜ 1 1 ) } s i n g c o s g , ,

[0033] C 2 3 = ( − 2 l e U ˜ 3 3 ) c o s 2 g + l [ 1 − l m e U ˜ 3 3 ] , ,

[0034] C 3 3 = [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 4 g − [ 4 e ( l U ˜ 3 3 + m U ˜ 1 1 ) ] c o s 2 g + ( M − l 2 m e U ˜ 3 3 ) ,

[0035] C 3 5 = { [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 2 g + 2 e ( l U ˜ 3 3 + m U ˜ 1 1 ) } s i n g c o s g ,

[0036] , ,

[0037] C 5 5 = [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 4 g − [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 2 g + m ( 1 − e U ˜ 1 1 ) ,

[0038] ,

[0039] in,

[0040] U ˜ 1 1 = 1 6 ( l + 2 m ) 3 ( 3 l + 4 m ) [ 1 + Oh ˜ 1 1 ( oh ) ] , (6a),

[0041] , (6b),

[0042] In formulas (6a) and (6b),

[0043] Oh ˜ 1 1 ( oh ) = 4 α f π i oh m [ q l or l + ( 1 − q l ) or g ] l + 2 m 3 l + 4 m ,

[0044] ,

[0045] Oh ˜ 3 3 2 = oh ( l + 2 m ) π m ( l + m ) α f 3 ( 1 K l − 1 K g ) 2 ( q l K l + 1 − q l K g ) − 2 × [ or l F l ( q l ) + or g F g ( 1 − q l ) ] , (7),

[0046] In formula (7), is the crack aspect ratio; is the angular frequency, , is the frequency; and are the liquid viscosity and gas viscosity respectively; and are the saturations of gas and liquid in the fracture under equilibrium state, ; and denote the bulk modulus of liquid and gas, respectively; F l ( q ) ≈ 0 . 0 5 3 ( 1 − q ) [ 1 + c o s π ( 1 − q ) ] , F g ( q ) ≈ 0 . 0 5 8 ( 1 − q ) [ 1 + c o s π ( 1 − q ) ] ;

[0047] Will and Simplified to:

[0048] , , (8),

[0049] Discard the formula (6b) , perform Taylor expansion on it and retain the first-order terms to obtain:

[0050] (9),

[0051] In formula (9), ; ;

[0052] neglect Imaginary part,

[0053] (10),

[0054] In the earthquake frequency band, ignore The dispersion and attenuation of , Equation (6a) can be approximated as:

[0055] (11),

[0056] Substituting equations (10) and (11) into equation (5), we can obtain the approximate formula for the complex stiffness coefficient of the inclined fracture reservoir considering the partial saturation attenuation effect:

[0057] C 1 1 ≈ [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 ( − m d T + M g d N − oh 2 M g ϒ ) ] c o s 2 g + M ( 1 − d N + oh 2 ϒ ) ,

[0058] C 1 2 ≈ [ 2 l g ( d N − oh 2 ϒ ) ] c o s 2 g + l ( 1 − d N + oh 2 ϒ ) ,

[0059] C 1 3 ≈ [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + l ( 1 − d N + oh 2 ϒ ) ,

[0060] C 1 5 ≈ { [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + 2 ( − m d T + M g d N − oh 2 M g ϒ ) } s i n g c o s g ,

[0061] , C 2 3 ≈ [ − 2 l g ( d N − oh 2 ϒ ) ] c o s 2 g + l [ 1 − x d N + oh 2 x ϒ ] ,

[0062] C 2 5 ≈ [ 2 l g ( d N − oh 2 ϒ ) ] s i n g c o s g ,

[0063] C 3 3 ≈ [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g − [ 4 ( m d T + l g d N − oh 2 l g ϒ ) ] c o s 2 g + ( M − l x d N + oh 2 l x ϒ ) ,

[0064] C 3 5 ≈ { [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + 2 ( m d T + l g d N − oh 2 l g ϒ ) } s i n g c o s g ,

[0065] , ,

[0066] C 5 5 ≈ [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + m ( 1 − d T ) ,

[0067] ,(12),

[0068] In formula (12), ; is the crack normal weakness parameter, is the fracture tangential weakness parameter, and , ; ,Will Defined as the indicator factor of oil-bearing fractured reservoir.

[0069] Preferably, the step 2 comprises:

[0070] For the case where two inclined fractured reservoirs are separated by a horizontal interface, under the assumption that the medium properties on both sides of the interface are weakly different, by discarding the term containing the product of the fracture weakness parameter and the perturbation, the perturbation of the complex stiffness coefficient is expressed as:

[0071] D C 1 1 = 4 m [ D ( d T c o s 4 g ) − g D ( d N c o s 4 g ) + oh 2 g D ( ϒ c o s 4 g ) ] + 4 [ M g D ( d N c o s 2 g ) − m D ( d T c o s 2 g ) − oh 2 M g D ( ϒ c o s 2 g ) ] + D M − M D d N + oh 2 M D ϒ ,

[0072] D C 1 2 = 2 l g [ D ( d N c o s 2 g ) − oh 2 D ( ϒ c o s 2 g ) ] + D l − l D d N + oh 2 l D ϒ ,

[0073] D C 1 3 = 4 m [ g D ( d N c o s 4 g ) − D ( d T c o s 4 g ) − oh 2 g D ( ϒ c o s 4 g ) + D ( d T c o s 2 g ) − g D ( d N c o s 2 g ) + oh 2 g D ( ϒ c o s 2 g ) ] + D l − l D d N + oh 2 l D ϒ ,

[0074] D C 1 5 = 4 m [ D ( d T s i n g c o s 3 g ) − g D ( d N s i n g c o s 3 g ) + oh 2 g D ( ϒ s i n g c o s 3 g ) ] + 2 [ M g D ( d N s i n g c o s g ) − m D ( d T s i n g c o s g ) − oh 2 M g D ( ϒ s i n g c o s g ) ] ,

[0075] ,

[0076] D C 2 3 = − 2 l g [ D ( d N c o s 2 g ) − oh 2 D ( ϒ c o s 2 g ) ] + D l − l x D d N + oh 2 l x D ϒ ,

[0077] D C 2 5 = 2 l g [ D ( d N s i n g c o s g ) − oh 2 D ( ϒ s i n g c o s g ) ] ,

[0078] D C 3 3 = 4 m [ D ( d T c o s 4 g ) − g D ( d N c o s 4 g ) + oh 2 g D ( ϒ c o s 4 g ) ] − 4 [ m D ( d T c o s 2 g ) + l g D ( d N c o s 2 g ) − oh 2 l g D ( ϒ c o s 2 g ) ] + D M − l x D d N + oh 2 l x D ϒ ,

[0079] D C 3 5 = 4 m [ g D ( d N s i n g c o s 3 g ) − D ( d T s i n g c o s 3 g ) − oh 2 g D ( ϒ s i n g c o s 3 g ) ] + 2 [ m D ( d T s i n g c o s g ) + l g D ( d N s i n g c o s g ) − oh 2 l g D ( ϒ s i n g c o s g ) ] ,

[0080] , ,

[0081] D C 5 5 = 4 m [ g D ( d N c o s 4 g ) − D ( d T c o s 4 g ) − oh 2 g D ( ϒ c o s 4 g ) + D ( d T c o s 2 g ) − g D ( d N c o s 2 g ) + oh 2 g D ( ϒ c o s 2 g ) ] + D m − m D d T ,

[0082] , (13),

[0083] In formula (13), △ means the difference in the parameters of the upper and lower layers of the reflection interface;

[0084] Combining Born's phase-stabilization method with inverse scattering theory, the seismic reflection coefficient of any anisotropic medium can be written as:

[0085] (14),

[0086] In formula (14), is the mass density; the horizontal line on the variable means the average value of the upper and lower medium parameters; is the incident angle of the seismic wave; ;and,

[0087] ,

[0088] ,

[0089] ,

[0090] ,

[0091] ,

[0092] ,

[0093] ,

[0094] ,

[0095] ,

[0096] ,

[0097] ,

[0098] ,

[0099] (15),

[0100] In formula (15), is the average value of the background P-wave velocity without cracks above and below the reflection interface; is the observation azimuth;

[0101] Substituting Equations (13) and (15) into Equation (14), we obtain the frequency-dependent seismic reflection coefficient equation for partially saturated inclined fracture reservoirs:

[0102] (16),

[0103] In formula (16):

[0104] ,

[0105] , , ,

[0106] ,

[0107] ,

[0108] ,

[0109] ,

[0110] C 2 ( i , f ) = g [ − 2 s i n 2 i t a n 2 i c o s 4 f + ( t a n 2 i − 4 s i n 2 i ) c o s 2 f + 1 − 2 c o s 2 i ] ,

[0111] ,

[0112] The reflection coefficient in equation (16) Expand to , is the crack azimuth.

[0113] Preferably, the step 3 comprises:

[0114] Rewrite the frequency-dependent reflection coefficient equation (16) of the inclined fracture reservoir into the Fourier series form:

[0115] (17),

[0116] In formula (17), is the 2h-order Fourier coefficient, and its expression is:

[0117] (18),

[0118] (19),

[0119] (20),

[0120] In formula (18):

[0121] (twenty one),

[0122] The weight coefficient is:

[0123] , ,

[0124] ,

[0125] , ,

[0126] ,

[0127] , ,

[0128] ,

[0129] , ,

[0130] , ,

[0131] , ,

[0132] , , ,

[0133] Rewrite Equation (17) as the sum of sine and cosine:

[0134] (twenty two),

[0135] In formula (22), , , directly calculated using discrete Fourier transform:

[0136] (twenty three),

[0137] Crack azimuth Under the constraints of fracture azimuth logging interpretation results, and It is estimated that:

[0138] (twenty four),

[0139] The second-order Fourier coefficients and the fourth-order Fourier coefficients are calculated by equation (25):

[0140] r 2 h ( i , oh ) = s i g n [ u 2 h ( i , oh ) c o s 2 f s y m ] [ u 2 h ( i , oh ) ] 2 + [ v 2 h ( i , oh ) ] 2 (25),

[0141] In formula (25), h = 1 or 2, To obtain the sign function, the reflection coefficients in equations (23) to (25) are Replaced with the corresponding five-dimensional seismic data, the result calculated by formula (25) is the 2h-order component of the five-dimensional seismic data;

[0142] The five-dimensional seismic data with N incident angles and M sampling points are transformed by continuous wavelet transform to obtain the two frequency components of the data. The seismic data with different frequency components are transformed by discrete Fourier transform to obtain the second-order components of the seismic data. The forward equation of the second-order components of the five-dimensional seismic data with different frequency components is obtained by convolving Equation (19) with the seismic wavelet:

[0143] (26),

[0144] In formula (26),

[0145] d 2 = [ d 2 ( i 1 , oh 1 ) ⋯ d 2 ( i N , oh 1 ) d 2 ( i 1 , oh 2 ) ⋯ d 2 ( i N , oh 2 ) ] T ,

[0146] m 2 = [ R ϒ 4 R ϒ 2 R d N 4 R d N 2 R d T 4 R d T 2 ] T ,

[0147] W = [ w ( i 1 , oh 1 ) 0 0 0 0 0 0 ⋱ 0 ⋮ ⋮ ⋮ ⋮ 0 w ( i N , oh 1 ) 0 ⋮ ⋮ ⋮ ⋮ 0 w ( i 1 , oh 2 ) 0 ⋮ ⋮ ⋮ ⋮ 0 ⋱ 0 0 0 0 0 0 w ( i N , oh 2 ) ] ,

[0148] A 2 = [ a 2 4 ( i 1 , oh 1 ) a 2 2 ( i 1 , oh 1 ) b 2 4 ( i 1 ) b 2 2 ( i 1 ) c 2 4 ( i 1 ) c 2 2 ( i 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ a 2 4 ( i N , oh 1 ) a 2 2 ( i N , oh 1 ) b 2 4 ( i N ) b 2 2 ( i N ) c 2 4 ( i N ) c 2 2 ( i N ) a 2 4 ( i 1 , oh 2 ) a 2 2 ( i 1 , oh 2 ) b 2 4 ( i 1 ) b 2 2 ( i 1 ) c 2 4 ( i 1 ) c 2 2 ( i 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ a 2 4 ( i N , oh 2 ) a 2 2 ( i N , oh 2 ) b 2 4 ( i N ) b 2 2 ( i N ) c 2 4 ( i N ) c 2 2 ( i N ) ] ,

[0149] In the above formula, The incident angle is , the angular frequency is The second-order component of the five-dimensional seismic data at time ; The angle of incidence is , the angular frequency is The azimuthally averaged wavelet matrix at ; and,

[0150] R ϒ 4 = [ D ( ϒ s i n 4 g ) | t = t 1 , D ( ϒ s i n 4 g ) | t = t 2 , ⋯ , D ( ϒ s i n 4 g ) | t = t M ] T R ϒ 2 = [ D ( ϒ s i n 2 g ) | t = t 1 , D ( ϒ s i n 2 g ) | t = t 2 , ⋯ , D ( ϒ s i n 2 g ) | t = t M ] T R d N 4 = [ D ( d N s i n 4 g ) | t = t 1 , D ( d N s i n 4 g ) | t = t 2 , ⋯ , D ( d N s i n 4 g ) | t = t M ] T R d N 2 = [ D ( d N s i n 2 g ) | t = t 1 , D ( d N s i n 2 g ) | t = t 2 , ⋯ , D ( d N s i n 2 g ) | t = t M ] T R d T 4 = [ D ( d T s i n 4 g ) | t = t 1 , D ( d T s i n 4 g ) | t = t 2 , ⋯ , D ( d T s i n 4 g ) | t = t M ] T R d T 2 = [ D ( d T s i n 2 g ) | t = t 1 , D ( d T s i n 2 g ) | t = t 2 , ⋯ , D ( d T s i n 2 g ) | t = t M ] T

[0151] (27),

[0152] In formula (27), Representative establishment The diagonal matrix of ; Representative Sampling moments;

[0153] Formula (26) can be expressed as:

[0154] (28),

[0155] In formula (28), , , represents the decorrelation matrix, which is estimated from the well logging interpretation results of the parameters to be inverted;

[0156] Under the Bayesian theory framework, the posterior probability density function of the parameter to be inverted is and the prior probability density function and likelihood function Positively correlated,

[0157] (29),

[0158] Assuming that the noise of the second-order component of the five-dimensional seismic data obeys a Gaussian distribution, the likelihood function is:

[0159] P ( d 2 | m ′ 2 ) = 1 2 π s d 2 e x p { − [ d 2 − G ′ 2 m ′ 2 ] T [ d 2 − G ′ 2 m ′ 2 ] 2 s d 2 2 } (30),

[0160] In formula (30), is the noise variance. Assuming that the inverted parameter obeys the Cauchy distribution, the prior probability density function is:

[0161] P ( m ' 2 ) = 1 ( π s m ' 2 2 ) 6 M ∏ i = 1 6 M 1 1 + [ m ' 2 ] 2 / s m ' 2 2 (31),

[0162] In formula (31), is the variance of the parameter to be inverted after decorrelation, estimated using the well logging interpretation results;

[0163] Substitute Equations (30) and (31) into Equation (29), and take the maximum posterior probability solution of Equation (29) as the goal, while introducing the low-frequency model constraint, and the final objective function is:

[0164] (32),

[0165] In formula (32), the low-frequency model constraint is:

[0166] (33),

[0167] In formula (33), 、 、 、 、 and is the regularization coefficient; 、 、 、 、 and In order 、 、 、 、 and Low-frequency model of 、 、 、 、 and is the integral matrix;

[0168] Using the reweighted iterative least squares algorithm to solve equation (32) we get ,but , using the Dow integral method to obtain Parameters to be inverted at sampling time:

[0169] (34),

[0170] In formula (34), is the starting time of each channel, The oil-bearing fracture reservoir indicator factor, fracture normal weakness parameter, fracture tangential weakness parameter and fracture dip angle at each sampling moment are calculated by the following formula:

[0171] ϒ ( t i i ) = [ ϒ s i n 2 g ( t i i ) ] 2 ϒ s i n 4 g ( t i i ) ， d N ( t i i ) = [ d N s i n 2 g ( t i i ) ] 2 d N s i n 4 g ( t i i ) ,

[0172] d T ( t i i ) = [ d T s i n 2 g ( t i i ) ] 2 d T s i n 4 g ( t i i ) , g ( t i i ) = a r c s i n ( d N 4 ( t i i ) d N 2 ( t i i ) ) ,(35),

[0173] Using the normal weakness parameter The high value area indicates the development area of ​​inclined fractures. The high-value areas are used to delineate the distribution range of oil-bearing fracture reservoirs.

[0174] The beneficial technical effects brought about by the present invention are:

[0175] The present invention fully considers the heterogeneity of oil and gas distribution, namely the partial saturation characteristics of the fluid, and the anisotropic characteristics induced by inclined fractures. The reflection coefficient equation constructed by the present invention is more consistent with the reflection characteristics of actual structural fracture-type oil and gas reservoirs, that is, the reflection characteristic characterization is more accurate. The Fourier series decomposition strategy is used to separate the isotropic and anisotropic components in the seismic data, which can effectively improve the accuracy of fracture prediction and oil and gas identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0176] Figure 1 This is a flow chart of an embodiment of the present invention.

[0177] Figure 2 Indicator factor for oil-bearing fractured reservoirs Schematic diagram of inversion results.

[0178] Figure 3 Schematic diagram of the inversion results of the fracture normal weakness parameters.

[0179] Figure 4 Schematic diagram of the inversion results of the fracture tangential weakness parameters.

[0180] Figure 5 Statistical histograms of the fracture dip angle around the wellbore, where (a) is the statistical histogram of the logging interpretation results of the fracture dip angle around the wellbore, and (b) is the statistical histogram of the inversion results of the fracture dip angle around the wellbore. DETAILED DESCRIPTION

[0181] The detailed description and technical contents of the present invention are described below with reference to the accompanying drawings. However, the drawings are only provided for reference and explanation and are not intended to limit the present invention.

[0182] like Figure 1 As shown, an embodiment of the present invention includes the following steps:

[0183] Step 1: construct an anisotropic stiffness matrix considering the attenuation mechanism of partially saturated inclined fractured reservoir;

[0184] Step 2: Construct a frequency-dependent reflection coefficient equation considering the attenuation mechanism of partially saturated inclined fractured reservoirs;

[0185] Step 3: Conduct five-dimensional seismic frequency-varying inversion of inclined fracture reservoirs to achieve inclined fracture prediction and oil and gas identification.

[0186] Step 1: The specific steps of constructing the anisotropic stiffness matrix considering the attenuation mechanism of partially saturated inclined fractured reservoir are as follows.

[0187] For a rock with a set of vertical fractures (assuming the fracture symmetry axis coincides with the x-axis), considering the oscillation diffusion effect of partially saturated fluid in the fracture induced by wave motion, its complex stiffness matrix can be expressed as:

[0188] C H T I = [ M − M 2 m e U ˜ 3 3 l − l M m e U ˜ 3 3 l − l M m e U ˜ 3 3 0 0 0 l − l M m e U ˜ 3 3 M − l 2 m e U ˜ 3 3 l − l 2 m e U ˜ 3 3 0 0 0 l − l M m e U ˜ 3 3 l − l 2 m e U ˜ 3 3 M − l 2 m e U ˜ 3 3 0 0 0 0 0 0 m 0 0 0 0 0 0 m − m e U ˜ 1 1 0 0 0 0 0 0 m − m e U ˜ 1 1 ] (1),

[0189] (1) In the formula, Represents vertical fracture reservoir Stiffness matrix; is the compression modulus, ; and are the first Lamé constant and shear modulus for the isotropic background; is the crack density; and It is a parameter related to fracture parameters (fracture aspect ratio, etc.) and fluid parameters (fluid viscosity and bulk modulus, etc.);

[0190] For a reservoir with a set of inclined fractures, the stiffness matrix can be obtained by By rotating the coordinates, we get

[0191] (2),

[0192] (2) In the formula, The fourth-order stiffness tensor representing the inclined fracture reservoir, is the crack inclination angle, 、 、 、 is a second-order tensor The value ranges of i, j, k, l, s, t, p, and q are 1, 2, and 3, for example, is the sth row and ith column element of the tensor L, and the second-order tensor for,

[0193] L = [ s i n g 0 − c o s g 0 1 0 c o s g 0 s i n g ] (3),

[0194] The fourth-order stiffness tensor representing the vertical fracture reservoir is Correspondingly, it needs to be explained that the fourth-order tensor element in the present invention is and Matrix elements The corresponding formula is as follows:

[0195] , , (4),

[0196] (4) In the formula, and is the Kronecker symbol;

[0197] Substituting equation (1) into equation (2), the stiffness matrix of the inclined fracture reservoir can be obtained as follows:

[0198] C T T I = [ C 1 1 C 1 2 C 1 3 0 C 1 5 0 C 1 2 C 2 2 C 2 3 0 C 2 5 0 C 1 3 C 2 3 C 3 3 0 C 3 5 0 0 0 0 C 4 4 0 C 4 6 C 1 5 C 2 5 C 3 5 0 C 5 5 0 0 0 0 C 4 6 0 C 6 6 ] (5),

[0199] In formula (5),

[0200] C 1 1 = [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 4 g + [ 4 e ( M U ˜ 3 3 − m U ˜ 1 1 ) ] c o s 2 g + M ( 1 − M m e U ˜ 3 3 ) ,

[0201] ,

[0202] C 1 3 = [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 4 g + [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 2 g + l ( 1 − M m e U ˜ 3 3 ) ,

[0203] C 1 5 = { [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 2 g + 2 e ( M U ˜ 3 3 − m U ˜ 1 1 ) } s i n g c o s g , ,

[0204] C 2 3 = ( − 2 l e U ˜ 3 3 ) c o s 2 g + l [ 1 − l m e U ˜ 3 3 ] , ,

[0205] C 3 3 = [ 4 m e ( U ˜ 1 1 − U ˜ 3 3 ) ] c o s 4 g − [ 4 e ( l U ˜ 3 3 + m U ˜ 1 1 ) ] c o s 2 g + ( M − l 2 m e U ˜ 3 3 ) ,

[0206] C 3 5 = { [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 2 g + 2 e ( l U ˜ 3 3 + m U ˜ 1 1 ) } s i n g c o s g ,

[0207] , ,

[0208] C 5 5 = [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 4 g − [ 4 m e ( U ˜ 3 3 − U ˜ 1 1 ) ] c o s 2 g + m ( 1 − e U ˜ 1 1 ) ,

[0209] ,

[0210] in,

[0211] U ˜ 1 1 = 1 6 ( l + 2 m ) 3 ( 3 l + 4 m ) [ 1 + Oh ˜ 1 1 ( oh ) ] ,(6a),

[0212] , (6b),

[0213] In formulas (6a) and (6b),

[0214] Oh ˜ 1 1 ( oh ) = 4 α f π i oh m [ q l or l + ( 1 − q l ) or g ] l + 2 m 3 l + 4 m ,

[0215] ,

[0216] Oh ˜ 3 3 2 = oh ( l + 2 m ) π m ( l + m ) α f 3 ( 1 K l − 1 K g ) 2 ( q l K l + 1 − q l K g ) − 2 × [ or l F l ( q l ) + or g F g ( 1 − q l ) ] , (7),

[0217] In formula (7), is the crack aspect ratio; is the angular frequency, , is the frequency; and are the liquid viscosity and gas viscosity respectively; and are the saturations of gas and liquid in the fracture under equilibrium state, ; and denote the bulk modulus of liquid and gas, respectively; F l ( q ) ≈ 0 . 0 5 3 ( 1 − q ) [ 1 + c o s π ( 1 − q ) ] , F g ( q ) ≈ 0 . 0 5 8 ( 1 − q ) [ 1 + c o s π ( 1 − q ) ] ; If the compression resistance of the gas is ignored, then and , and assuming that the gas viscosity is very small, then and Can be simplified to:

[0218] , , (8),

[0219] Therefore, discard the , and Taylor expansion of it and retaining the first-order terms yields:

[0220] (9),

[0221] In formula (9), ; ;

[0222] because ,so The imaginary part is significantly lower than the real part, so it can be ignored The imaginary part, that is,

[0223] (10),

[0224] There is a theory that proves that within the earthquake frequency band, The dispersion and attenuation of can be ignored, so Equation (6a) can be approximated as:

[0225] (11),

[0226] Substituting Equations (10) and (11) into Equation (5), we can obtain the approximate formula for the complex stiffness coefficient of the inclined fracture reservoir considering the partial saturation attenuation effect:

[0227] C 1 1 ≈ [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 ( − m d T + M g d N − oh 2 M g ϒ ) ] c o s 2 g + M ( 1 − d N + oh 2 ϒ ) ,

[0228] C 1 2 ≈ [ 2 l g ( d N − oh 2 ϒ ) ] c o s 2 g + l ( 1 − d N + oh 2 ϒ ) ,

[0229] C 1 3 ≈ [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + l ( 1 − d N + oh 2 ϒ ) ,

[0230] C 1 5 ≈ { [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + 2 ( − m d T + M g d N − oh 2 M g ϒ ) } s i n g c o s g ,

[0231] , C 2 3 ≈ [ − 2 l g ( d N − oh 2 ϒ ) ] c o s 2 g + l [ 1 − x d N + oh 2 x ϒ ] ,

[0232] C 2 5 ≈ [ 2 l g ( d N − oh 2 ϒ ) ] s i n g c o s g ,

[0233] C 3 3 ≈ [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g − [ 4 ( m d T + l g d N − oh 2 l g ϒ ) ] c o s 2 g + ( M − l x d N + oh 2 l x ϒ ) ,

[0234] C 3 5 ≈ { [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + 2 ( m d T + l g d N − oh 2 l g ϒ ) } s i n g c o s g ,

[0235] , ,

[0236] C 5 5 ≈ [ − 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 4 g + [ 4 m ( d T − g d N + oh 2 g ϒ ) ] c o s 2 g + m ( 1 − d T ) ,

[0237] ,(12),

[0238] In formula (12), ; is the crack normal weakness parameter, is the fracture tangential weakness parameter, and , , both are linear increasing functions of crack density and can well indicate the spatial distribution of cracks; , It is a monotonically increasing function of liquid viscosity, liquid saturation and fracture density. That is, the greater the fracture density, the higher the oil saturation. The larger the Defined as the indicator factor of oil-bearing fractured reservoir.

[0239] Step 2: The specific steps for constructing the frequency-dependent reflection coefficient equation considering the attenuation mechanism of partially saturated inclined fractured reservoirs are as follows.

[0240] For the case where two inclined fractured reservoirs are separated by a horizontal interface, under the assumption that the medium properties on both sides of the interface are weakly different, by discarding the term containing the product of the fracture weakness parameter and the perturbation, the perturbation of the complex stiffness coefficient can be expressed as:

[0241] D C 1 1 = 4 m [ D ( d T c o s 4 g ) − g D ( d N c o s 4 g ) + oh 2 g D ( ϒ c o s 4 g ) ] + 4 [ M g D ( d N c o s 2 g ) − m D ( d T c o s 2 g ) − oh 2 M g D ( ϒ c o s 2 g ) ] + D M − M D d N + oh 2 M D ϒ ,

[0242] D C 1 2 = 2 l g [ D ( d N c o s 2 g ) − oh 2 D ( ϒ c o s 2 g ) ] + D l − l D d N + oh 2 l D ϒ ,

[0243] D C 1 3 = 4 m [ g D ( d N c o s 4 g ) − D ( d T c o s 4 g ) − oh 2 g D ( ϒ c o s 4 g ) + D ( d T c o s 2 g ) − g D ( d N c o s 2 g ) + oh 2 g D ( ϒ c o s 2 g ) ] + D l − l D d N + oh 2 l D ϒ ,

[0244] D C 1 5 = 4 m [ D ( d T s i n g c o s 3 g ) − g D ( d N s i n g c o s 3 g ) + oh 2 g D ( ϒ s i n g c o s 3 g ) ] + 2 [ M g D ( d N s i n g c o s g ) − m D ( d T s i n g c o s g ) − oh 2 M g D ( ϒ s i n g c o s g ) ] ,

[0245] ,

[0246] D C 2 3 = − 2 l g [ D ( d N c o s 2 g ) − oh 2 D ( ϒ c o s 2 g ) ] + D l − l x D d N + oh 2 l x D ϒ ,

[0247] D C 2 5 = 2 l g [ D ( d N s i n g c o s g ) − oh 2 D ( ϒ s i n g c o s g ) ] ,

[0248] D C 3 3 = 4 m [ D ( d T c o s 4 g ) − g D ( d N c o s 4 g ) + oh 2 g D ( ϒ c o s 4 g ) ] − 4 [ m D ( d T c o s 2 g ) + l g D ( d N c o s 2 g ) − oh 2 l g D ( ϒ c o s 2 g ) ] + D M − l x D d N + oh 2 l x D ϒ ,

[0249] D C 3 5 = 4 m [ g D ( d N s i n g c o s 3 g ) − D ( d T s i n g c o s 3 g ) − oh 2 g D ( ϒ s i n g c o s 3 g ) ] + 2 [ m D ( d T s i n g c o s g ) + l g D ( d N s i n g c o s g ) − oh 2 l g D ( ϒ s i n g c o s g ) ] ,

[0250] , ,

[0251] D C 5 5 = 4 m [ g D ( d N c o s 4 g ) − D ( d T c o s 4 g ) − oh 2 g D ( ϒ c o s 4 g ) + D ( d T c o s 2 g ) − g D ( d N c o s 2 g ) + oh 2 g D ( ϒ c o s 2 g ) ] + D m − m D d T ,

[0252] , (13),

[0253] In formula (13), △ means the difference in the parameters of the upper and lower layers of the reflection interface;

[0254] Combining Born's phase-stabilization method with inverse scattering theory, the seismic reflection coefficient of any anisotropic medium can be written as:

[0255] (14),

[0256] In formula (14), is the mass density; the horizontal line on the variable means the average value of the upper and lower medium parameters; is the incident angle of the seismic wave; ;and,

[0257] ,

[0258] ,

[0259] ,

[0260] ,

[0261] ,

[0262] ,

[0263] ,

[0264] ,

[0265] ,

[0266] ,

[0267] ,

[0268] ,

[0269] (15),

[0270] In formula (15), is the average value of the background P-wave velocity without cracks above and below the reflection interface; is the observation azimuth;

[0271] Substituting Equations (13) and (15) into Equation (14), we can obtain the frequency-dependent seismic reflection coefficient equation for partially saturated inclined fracture reservoirs:

[0272] (16),

[0273] In formula (16):

[0274] ,

[0275] , , ,

[0276] ,

[0277] ,

[0278] ,

[0279] ,

[0280] C 2 ( i , f ) = g [ − 2 s i n 2 i t a n 2 i c o s 4 f + ( t a n 2 i − 4 s i n 2 i ) c o s 2 f + 1 − 2 c o s 2 i ] ,

[0281] .

[0282] In order to be applicable to inclined fracture reservoirs with arbitrary fracture dip angles, the reflection coefficient in Eq. (16) is replaced by Expand to , is the crack azimuth.

[0283] Step 3: Carry out five-dimensional seismic frequency-varying inversion of inclined fracture reservoirs. The specific steps for achieving inclined fracture prediction and fluid identification are as follows.

[0284] Rewrite the frequency-dependent reflection coefficient equation (16) of the inclined fracture reservoir into the Fourier series form:

[0285] (17),

[0286] In formula (17), is the 2h-order Fourier coefficient, and its expression is:

[0287] (18),

[0288] (19),

[0289] (20),

[0290] In formula (18):

[0291] (twenty one),

[0292] The weight coefficient is:

[0293] , ,

[0294] ,

[0295] , ,

[0296] ,

[0297] , ,

[0298] ,

[0299] , ,

[0300] , ,

[0301] , ,

[0302] , , ,

[0303] To calculate the Fourier coefficients, rewrite Equation (17) in the form of a sum of sines and cosines:

[0304] (twenty two),

[0305] In formula (22), , , which can be directly calculated using discrete Fourier transform:

[0306] (twenty three),

[0307] Crack azimuth Under the constraints of the fracture azimuth logging interpretation results, and It is estimated that:

[0308] (twenty four),

[0309] The second-order Fourier coefficients and the fourth-order Fourier coefficients can be calculated by formula (25):

[0310] r 2 h ( i , oh ) = s i g n [ u 2 h ( i , oh ) c o s 2 f s y m ] [ u 2 h ( i , oh ) ] 2 + [ v 2 h ( i , oh ) ] 2 (25),

[0311] In formula (25), h = 1 or 2, To obtain the sign function, the reflection coefficients in equations (23) to (25) are Replaced with the corresponding five-dimensional seismic data, the result calculated by formula (25) is the 2h-order component of the five-dimensional seismic data;

[0312] Continuous wavelet transform is performed on the five-dimensional seismic data with N incident angles and M sampling points to obtain the two frequency components of the data. Then, discrete Fourier transform is performed on the seismic data with different frequency components to obtain the second-order components of the seismic data. Finally, the forward equation of the second-order components of the five-dimensional seismic data with different frequency components is obtained by convolving Equation (19) with the seismic wavelet:

[0313] (26),

[0314] In formula (26),

[0315] d 2 = [ d 2 ( i 1 , oh 1 ) ⋯ d 2 ( i N , oh 1 ) d 2 ( i 1 , oh 2 ) ⋯ d 2 ( i N , oh 2 ) ] T ,

[0316] m 2 = [ R ϒ 4 R ϒ 2 R d N 4 R d N 2 R d T 4 R d T 2 ] T ,

[0317] W = [ w ( i 1 , oh 1 ) 0 0 0 0 0 0 ⋱ 0 ⋮ ⋮ ⋮ ⋮ 0 w ( i N , oh 1 ) 0 ⋮ ⋮ ⋮ ⋮ 0 w ( i 1 , oh 2 ) 0 ⋮ ⋮ ⋮ ⋮ 0 ⋱ 0 0 0 0 0 0 w ( i N , oh 2 ) ] ,

[0318] A 2 = [ a 2 4 ( i 1 , oh 1 ) a 2 2 ( i 1 , oh 1 ) b 2 4 ( i 1 ) b 2 2 ( i 1 ) c 2 4 ( i 1 ) c 2 2 ( i 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ a 2 4 ( i N , oh 1 ) a 2 2 ( i N , oh 1 ) b 2 4 ( i N ) b 2 2 ( i N ) c 2 4 ( i N ) c 2 2 ( i N ) a 2 4 ( i 1 , oh 2 ) a 2 2 ( i 1 , oh 2 ) b 2 4 ( i 1 ) b 2 2 ( i 1 ) c 2 4 ( i 1 ) c 2 2 ( i 1 ) ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ a 2 4 ( i N , oh 2 ) a 2 2 ( i N , oh 2 ) b 2 4 ( i N ) b 2 2 ( i N ) c 2 4 ( i N ) c 2 2 ( i N ) ] ,

[0319] In the above formula, The incident angle is , the angular frequency is The second-order component of the five-dimensional seismic data at time ; The angle of incidence is , the angular frequency is The azimuthally averaged wavelet matrix at ; and,

[0320] R ϒ 4 = [ D ( ϒ s i n 4 g ) | t = t 1 , D ( ϒ s i n 4 g ) | t = t 2 , ⋯ , D ( ϒ s i n 4 g ) | t = t M ] T R ϒ 2 = [ D ( ϒ s i n 2 g ) | t = t 1 , D ( ϒ s i n 2 g ) | t = t 2 , ⋯ , D ( ϒ s i n 2 g ) | t = t M ] T R d N 4 = [ D ( d N s i n 4 g ) | t = t 1 , D ( d N s i n 4 g ) | t = t 2 , ⋯ , D ( d N s i n 4 g ) | t = t M ] T R d N 2 = [ D ( d N s i n 2 g ) | t = t 1 , D ( d N s i n 2 g ) | t = t 2 , ⋯ , D ( d N s i n 2 g ) | t = t M ] T R d T 4 = [ D ( d T s i n 4 g ) | t = t 1 , D ( d T s i n 4 g ) | t = t 2 , ⋯ , D ( d T s i n 4 g ) | t = t M ] T R d T 2 = [ D ( d T s i n 2 g ) | t = t 1 , D ( d T s i n 2 g ) | t = t 2 , ⋯ , D ( d T s i n 2 g ) | t = t M ] T

[0321] (27),

[0322] In formula (27), Representative establishment The diagonal matrix of ; Representative Sampling moments;

[0323] In order to overcome the instability of the inversion results caused by the correlation of the inverted parameters, Equation (26) is expressed as:

[0324] (28),

[0325] In formula (28), , ; Represents the decorrelation matrix, which can be estimated from the well logging interpretation results of the parameters to be inverted and is a widely used method in the industry;

[0326] Under the Bayesian theory framework, the posterior probability density function of the parameter to be inverted is and the prior probability density function and likelihood function is positively correlated, that is

[0327] (29),

[0328] Assuming that the noise of the second-order component of the five-dimensional seismic data obeys a Gaussian distribution, the likelihood function is:

[0329] P ( d 2 | m ′ 2 ) = 1 2 π s d 2 e x p { − [ d 2 − G ′ 2 m ′ 2 ] T [ d 2 − G ′ 2 m ′ 2 ] 2 s d 2 2 } (30),

[0330] In formula (30), is the noise variance. Assuming that the inverted parameter obeys the Cauchy distribution, the prior probability density function is:

[0331] P ( m ' 2 ) = 1 ( π s m ' 2 2 ) 6 M ∏ i = 1 6 M 1 1 + [ m ' 2 ] 2 / s m ' 2 2 (31),

[0332] In formula (31), is the variance of the parameter to be inverted after decorrelation, which can be estimated using the well logging interpretation results;

[0333] Substitute equations (30) and (31) into equation (29), and take the maximum posterior probability solution of equation (29) as the goal, while introducing the low-frequency model constraint, then the final objective function is:

[0334] (32),

[0335] In formula (32), the low-frequency model constraint is:

[0336] (33),

[0337] In formula (33), 、 、 、 、 and is the regularization coefficient; 、 、 、 、 and In order 、 、 、 、 and Low-frequency model of 、 、 、 、 and is the integral matrix;

[0338] Using the reweighted iterative least squares algorithm to solve equation (32) we get ,but , and then use the Dow integral method to get the Parameters to be inverted at sampling time:

[0339] (34),

[0340] In formula (34), is the starting time of each track; on this basis, The oil-bearing fracture reservoir indicator factor, fracture normal weakness parameter, fracture tangential weakness parameter and fracture dip angle at each sampling moment can be calculated by the following formula:

[0341] ϒ ( t i i ) = [ ϒ s i n 2 g ( t i i ) ] 2 ϒ s i n 4 g ( t i i ) ， d N ( t i i ) = [ d N s i n 2 g ( t i i ) ] 2 d N s i n 4 g ( t i i ) ,

[0342] d T ( t i i ) = [ d T s i n 2 g ( t i i ) ] 2 d T s i n 4 g ( t i i ) , g ( t i i ) = a r c s i n ( d N 4 ( t i i ) d N 2 ( t i i ) ) ,(35),

[0343] Finally, using the normal weakness parameter The high value area indicates the development area of ​​inclined fractures; using the indicator factor The high-value areas are used to delineate the distribution range of oil-bearing fracture reservoirs.

[0344] Testing of the method of the present invention:

[0345] The method proposed in the present invention is applied to the prediction of a tight fractured oil reservoir in southwestern China to verify the rationality and feasibility of the method. Figure 2-Figure 4 The predicted results for the oil-bearing fracture reservoir indicator factor, fracture normal weakness parameter, and fracture tangential weakness parameter are shown. The figure shows that the predicted values ​​for these factors, fracture normal weakness parameter, and fracture tangential weakness parameter are all high in the range indicated by the white arrows, indicating high oil saturation and fracture development within this range. Horizontal well drilling logs show that the formation indicated by the white arrows represents a well-developed fractured reservoir, validating the predictions. Figure 5 Figures (a) and (b) show the well logging interpretation results and the inversion results statistical histograms of the fracture dip angle around the wellbore, respectively. By comparing the data of the two histograms, it can be found that the fracture dip inversion results are reasonable.

[0346] The above embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Any equivalent substitution or modification made by those skilled in the art based on the present invention is within the protection scope of the present invention.

Claims

1. A method for predicting inclined fractures and identifying oil and gas considering the partial saturation attenuation mechanism, characterized in that: It includes the following steps: Step 1, construct an anisotropic stiffness matrix considering the attenuation mechanism of partially saturated inclined fractured reservoir; Step 2, construct a frequency-dependent reflection coefficient equation considering the attenuation mechanism of partially saturated inclined fractured reservoir; Step 3: Conduct five-dimensional seismic frequency-variable inversion of inclined fracture reservoirs to achieve inclined fracture prediction and oil and gas identification; In step 2, the frequency-dependent reflection coefficient equation is: （16）， In formula (16): ， ， ， ， ， ， ， ， ， ， in, is the mass density; is the compression modulus, ; and are the first Lamé constant and shear modulus of the isotropic background; the horizontal line on the variable means the average value of the upper and lower medium parameters; △ means the difference between the upper and lower medium parameters of the reflection interface; , is the ratio of the shear modulus to the compression modulus of the isotropic background; ; is the incident angle of the seismic wave; is the observation azimuth; is the angular frequency, , is the frequency; It is the indicator factor of oil-bearing fracture reservoir; is the crack inclination; is the crack normal weakness parameter, is the fracture tangential weakness parameter; In formula (16), Expand to , is the crack azimuth.

2. The inclined fracture prediction and oil and gas identification method considering the partial saturation attenuation mechanism as claimed in claim 1 is characterized in that: The step 1 comprises: For a rock with a set of vertical fractures, assuming that the fracture symmetry axis coincides with the x-axis, and considering the oscillation diffusion effect of partially saturated fluid in the fracture induced by wave motion, its complex stiffness matrix is ​​expressed as: （1）， (1) In the formula, Represents a vertical fracture reservoir Stiffness matrix; is the compression modulus, ; and are the first Lamé constant and shear modulus for the isotropic background; is the crack density; and It is a parameter related to fracture parameters and fluid parameters; For a reservoir with a set of inclined fractures, the stiffness matrix is ​​calculated by By rotating the coordinates, we get （2）， (2) In the formula, The fourth-order stiffness tensor representing the inclined fractured reservoir, is the crack inclination angle, , , , is a second-order tensor The elements of i, j, k, l, s, t, p, q range from 1, 2, 3, and are second-order tensors for, （3）， The fourth-order stiffness tensor representing the vertical fracture reservoir is Correspondingly, the fourth-order tensor element and Matrix Elements The corresponding formula is as follows: ， ，（4）， (4) In the formula, and is the Kronecker symbol; Substituting equation (1) into equation (2), the stiffness matrix of the inclined fracture reservoir is obtained as follows: （5）， In formula (5), ， ， ， ， ， ， ， ， ， ， ， ， ， in, , （6a）， ， （6b）， In formula (6a) and (6b), ， ， ，（7）， In formula (7), is the crack aspect ratio; is the angular frequency, , is the frequency; and are the liquid viscosity and gas viscosity respectively; and are the saturations of gas and liquid in the fracture under equilibrium state, ; and denote the bulk modulus of liquid and gas, respectively; , ; Will and Simplified to: ， ，（8）， Discard the , perform Taylor expansion on it and retain the first-order terms to obtain: （9）， In formula (9), ; ; neglect Imaginary part, （10）， In the earthquake frequency band, ignore The dispersion and attenuation of , equation (6a) can be approximated as: （11）， Substituting equations (10) and (11) into equation (5), we can obtain the approximate formula of the complex stiffness coefficient of the inclined fracture reservoir considering the partial saturation attenuation effect: , , , , ， , , , , ， ， , ,（12）， In formula (12), ; is the crack normal weakness parameter, is the fracture tangential weakness parameter, and , ; ,Will Defined as the indicator factor of oil-bearing fractured reservoir.

3. The inclined fracture prediction and oil and gas identification method considering the partial saturation attenuation mechanism as claimed in claim 2 is characterized in that: The step 2 comprises: For the case where two inclined fractured reservoirs are separated by a horizontal interface, under the assumption that the medium properties on both sides of the interface are weakly different, by discarding the term containing the product of the fracture weakness parameter and the disturbance amount, the disturbance amount of the complex stiffness coefficient is expressed as: ， ， ， ， ， ， ， ， ， ， ， ， ，（13）， In formula (13), △ means the difference in the parameters of the upper and lower layers of the reflection interface; Combining Born's phase stability method with inverse scattering theory, the seismic reflection coefficient of any anisotropic medium can be written as: （14）， In formula (14), is the mass density; the horizontal line on the variable means the average value of the upper and lower medium parameters; is the incident angle of the seismic wave; ;and, ， ， ， ， ， ， ， ， ， ， ， ， ，（15）， In formula (15), is the average value of the crack-free background P-wave velocity above and below the reflection interface; is the observation azimuth; Substituting equation (13) and equation (15) into equation (14), the frequency-dependent reflection coefficient equation (16) is obtained.

4. The inclined fracture prediction and oil and gas identification method considering the partial saturation attenuation mechanism as claimed in claim 3 is characterized in that: The step 3 comprises: Rewrite the frequency-dependent reflection coefficient equation (16) into Fourier series form: （17）， In formula (17), is the 2h-order Fourier coefficient, and its expression is: （18）， （19）， （20）， In formula (18): （21）， The weight coefficient is: , ， ， , ， ， , ， ， , ， , ， , ， , , ， Rewrite equation (17) as the sum of sine and cosine: （22）， In formula (22), , , directly calculated using discrete Fourier transform: （23）， Crack azimuth Under the constraints of fracture azimuth logging interpretation results, and It is estimated that: （24）， In formula (24), , , The second-order Fourier coefficients and the fourth-order Fourier coefficients are calculated by equation (25): （25）， In formula (25), h = 1 or 2, To obtain the sign function, the reflection coefficient in equation (23) to equation (25) is Replaced with the corresponding five-dimensional seismic data, the result calculated by equation (25) is the 2h-order component of the five-dimensional seismic data; The five-dimensional seismic data with N incident angles and B sampling points are transformed by continuous wavelet transform to obtain the two frequency components of the data. The seismic data with different frequency components are transformed by discrete Fourier transform to obtain the second-order components of the seismic data. The forward equations of the second-order components of seismic data with different frequency components are established as follows: （26）， In formula (26), ， ， ， ， In the above formula, The incident angle is , the angular frequency is The second-order component of the five-dimensional seismic data at ; The incident angle is , the angular frequency is The azimuthally averaged wavelet matrix at ; and, ， ，（27）， In formula (27), Representative establishment The diagonal matrix of ; Representative sampling moments; Formula (26) can be expressed as: （28）， In formula (28), , , represents the decorrelation matrix, which is estimated from the well logging interpretation results of the parameters to be inverted; Under the Bayesian theory framework, the posterior probability density function of the inverted parameter is With the prior probability density function and the likelihood function Positively correlated, （29）， Assuming that the noise of the second-order component of the five-dimensional seismic data obeys a Gaussian distribution, the likelihood function is: （30）， In formula (30), is the noise variance. Assuming that the inverted parameters obey the Cauchy distribution, the prior probability density function is: （31）， In formula (31), is the variance of the parameter to be inverted after decorrelation, which is estimated using the well logging interpretation results; Substitute equations (30) and (31) into equation (29), and take the maximum a posteriori probability solution of equation (29) as the goal, while introducing low-frequency model constraints, the final objective function is: （32）， In formula (32), the low-frequency model constraint is: （33）， In formula (33), , , , , and is the regularization coefficient; , , , , and In order , , , , and Low frequency model of , , , , and is the integral matrix; Using the reweighted iterative least squares algorithm to solve equation (32) we get ,but , using the Dow integral method, we get Parameters to be inverted at sampling time: （34）, In formula (34), is the start time of each channel, The oil-bearing fracture reservoir indicator factor, fracture normal weakness parameter, fracture tangential weakness parameter and fracture dip angle at each sampling time are calculated by the following formula: , ,（35）, Using the normal weakness parameter The high value area indicates the development area of ​​inclined fractures. The high-value areas are used to delineate the distribution range of oil-bearing fracture reservoirs.

Images