Carbonate rock fracture prediction method and device based on rock mechanics

Abstract

The invention discloses a carbonate rock fracture prediction method and device based on rock mechanics, and relates to the technical field of oil and gas reservoir exploration, and the method is characterized in that according to drilled reservoir parameters and a lithofacies model, a function relation between lithofacies and reservoir parameters is obtained through geostatistical analysis; establishing a reservoir parameter model through phase control stochastic simulation according to the function relationship; obtaining a density model according to the reservoir parameter model, and performing integration according to the density model to obtain an overlying rock pressure model; calculating pore pressure of different stratums according to lithology; a crustal stress model is constructed based on a reservoir parameter model, an overlying rock pressure model, pore pressure and strain data based on single well crustal stress analysis, stress and strain data are imported into a fracture parameter calculation model based on the crustal stress model, and three-dimensional quantitative prediction of low-permeability reservoir fractures is achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas reservoir exploration, and more specifically, it relates to a method and device for predicting carbonate rock fractures based on rock mechanics. Background Art

[0002] In recent years, with the continuous increase in the world's demand for oil and gas energy and the continuous deepening of oil and gas exploration, global oil and gas exploration has gradually shifted from conventional to unconventional, from shallow to deep and even ultra-deep. Low-permeability and ultra-low-permeability reservoirs have received widespread attention and continuous breakthroughs. Among them, fracture prediction at different scales has become one of the current research hotspots. In low-permeability reservoirs, tectonic fractures are the key factors to improve reservoir physical properties and increase oil and gas production capacity. The accurate description and effective prediction of tectonic fractures have important practical significance for guiding the exploration and development of low-permeability oil and gas reservoirs.

[0003] Different from other regions in the world, low-permeability reservoirs are widely developed in China, and the lithologies mainly include clastic rocks and carbonate rocks, with special geological characteristics such as deep burial, strong diagenesis, low porosity and low permeability, and developed fractures. Moreover, due to frequent sea or lake level changes, there is a phenomenon of frequent interbedding of various soft and hard lithologies vertically. Under the action of multi-phase tectonic movements, the fracture development in these reservoirs is not only controlled by factors such as regional stress fields, faults and folds, but also significantly affected by mudstone interlayers and lithological interfaces, further exacerbating the heterogeneity of the fracture spatial distribution, thus bringing great difficulties to the identification and prediction of fractures. At present, good results have been achieved in describing and predicting the fracture development law from the perspective of stress fields, but there is still no effective and general method for predicting fractures in interbedded low-permeability reservoirs. There are certain deviations between the prediction results and the single-well measurement results, especially there are certain contradictions with production dynamics. The key problem is that the complex fracture mechanism of this composite rock mass is still unclear, and a reasonable strength and mechanical model has not been established to solve the problems of characterizing the "quality" (fracture generation and propagation) and "quantity" (fracture scale and quantity) of fractures.

[0004] Therefore, it is necessary to further study and explore the fracture mechanism and evolution law of interbedded rock masses based on rock mechanics experiments and fracture mechanics theories, and combine the fine geomechanics model and numerical simulation means under the constraint of rock mechanics layers to achieve three-dimensional quantitative prediction of fractures in low-permeability reservoirs. Summary of the Invention

[0005] To achieve the above object, the present invention provides the following technical solution: A method for predicting carbonate rock fractures based on rock mechanics, which includes the following steps:

[0006] S1. Establish a lithofacies model.

[0007] S2. Obtain logging data, which includes acoustic logging data and density logging data. According to the interpretation results of the logging data, obtain the reservoir parameters of the drilled wells. The reservoir parameters of the drilled wells include density, seismic wave velocity, Young's modulus, and Poisson's ratio.

[0008] S3. Based on the reservoir parameters of the drilled wells and the lithofacies model, use geostatistical analysis to obtain the functional relationship between lithofacies and reservoir parameters. Establish a reservoir parameter model through facies-controlled stochastic simulation according to this functional relationship.

[0009] S4. Obtain a density model according to the reservoir parameter model, and obtain an overlying rock pressure model by integrating according to the density model.

[0010] S5. Calculate the pore pressure of different strata by lithology.

[0011] For clastic rock strata, the calculation formula for the formation pore pressure Pp is as follows:

[0012]

[0013] In the formula, Pp is the formation pore pressure, Mpa. Sv is the overlying formation pressure, MPa, obtained through the overlying rock pressure model. Ph is the hydrostatic pressure, Mpa, which is the product of the acceleration of gravity, water density, and depth value. n is the Eaton index, dimensionless, and the Eaton coefficient of different strata is obtained based on the drilled well data. V is the seismic wave velocity and density, where Vo is obtained through the reservoir parameter model. Vn is the normal compaction value, obtained based on the normal compaction trend line of the drilled wells.

[0014] For carbonate rock strata, the calculation formula for the formation pore pressure Pp is as follows:

[0015] p p =Sv-(Ae B*DT )

[0016] In the formula, Sv is the overlying formation pressure, MPa, obtained through the overlying rock pressure model. DT is the acoustic travel time, us / m, obtained through the interpretation results of the logging data. e is the natural constant, and A and B are coefficients, which are obtained based on the carbonate rock formation data of the drilled wells.

[0017] S6. Construct a geostress model, and the formula is as follows:

[0018]

[0019]

[0020] In the formula, S H is the maximum horizontal principal stress, Mpa. S hσh is the minimum horizontal principal stress, in Mpa. Pp is the formation pore pressure, in MPa, obtained through S5. α is the Biot coefficient, dimensionless. Sv is the overburden pressure, in MPa, obtained through the overburden rock pressure model. v is the Poisson's ratio, dimensionless, obtained through the reservoir parameter model. E is the Young's modulus, in MPa, obtained through the reservoir parameter model. εx is the minimum tectonic stress coefficient, dimensionless. εy is the maximum tectonic stress coefficient, dimensionless. εx and εy for different formations are obtained based on the data of drilled wells.

[0021] S7. Calculate the fracture parameters, and the formula is as follows:

[0022]

[0023]

[0024]

[0025] In the formula, σ p is the rock fracture stress, in Mpa, and σ p for different formations is obtained based on the data of drilled wells. σ1, σ2, and σ3 are the maximum, intermediate, and minimum effective principal stresses, in MPa, calculated through the in-situ stress model. L1, L2, and L3 are the unit lengths along the principal stress directions of σ1, σ2, and σ3, in m, calculated based on σ1, σ2, σ3, the Young's modulus and Poisson's ratio obtained through the reservoir parameter model. μ is the Poisson's ratio, dimensionless, obtained through the reservoir parameter model. θ is the fracture rupture angle, in °, calculated through the in-situ stress model. ε3 is the minimum principal strain, dimensionless, calculated based on σ3 and the Young's modulus obtained through the reservoir parameter model. ε0 is the maximum elastic tensile strain borne by the rock, dimensionless, calculated based on the maximum horizontal principal stress obtained through the in-situ stress model and the Young's modulus obtained through the reservoir parameter model. J is the fracture surface energy, in J / m2, obtained by fitting based on the formation structure of the study area, the Young's modulus obtained through the reservoir parameter model, and the in-situ stress model. Dvf is the tectonic fracture volume density, in m2 / m3. Dlf is the tectonic fracture line density, in fractures / m. b is the effective aperture of the tectonic fracture.

[0026] The present invention is further configured such that the step S1 is:

[0027] S11. Based on the seismic fault interpretation results, using the geological stratification data of drilled wells as conditions and the marker beds interpreted by seismic as trend constraints, generate the structural surfaces of sub-layers. Then, adjust the structural surfaces according to the known geological profiles to generate a structural model of the study area that conforms to the true structural form.

[0028] S12. Based on the structural model of the study area, establish a sedimentary microfacies model of the study area according to the sedimentary microfacies interpretation results of drilled wells.

[0029] S13. Conduct lithofacies interpretation based on well logging data. On the basis of the sedimentary microfacies model in the study area, establish a lithofacies model under the control of sedimentary microfacies according to the lithofacies interpretation results.

[0030] The present invention is further configured as: The structural model of the study area established in S11 is a three-dimensional grid model, the grid direction is consistent with the line track direction of seismic data acquisition, and the lateral grid step size is consistent with the seismic data trace interval.

[0031] The present invention is further configured as: The lateral grid step size is 25m * 25m. Determine the target study formation. The vertical grid step size of the formation above the target study formation is 10m, and the vertical grid step size of the target study formation is 5m.

[0032] The present invention is further configured as: After the in-situ stress parameter model in S6 is constructed, use the rock failure criterion to simulate the wellbore shear sloughing situation, and then compare it with the wellbore condition. When the wellbore sloughing width is consistent with the wellbore diameter enlargement, perform S7.

[0033] The rock failure criterion is selected according to the in-situ stress model and the engineering mud density. The selection basis is to simulate the wellbore sloughing width under different failure criteria, compare it with the wellbore diameter enlargement, and select the rock failure criterion with the most consistent simulated wellbore sloughing width and wellbore diameter enlargement. The rock failure criteria include Mohr-Coulomb criterion, modified Lade criterion, Tresca, and Circumscribed Drucker-Prager.

[0034] The present invention is further configured as: The Mohr-Coulomb criterion is selected as the rock failure criterion.

[0035] The present invention is further configured as: For carbonate rock formations in S5, the calculation formula for the formation pore pressure Pp is as follows: p p = Sv - (224e -0.0193DT ).

[0036] The present invention is further configured as: The calculation formulas for Young's modulus and Poisson's ratio based on well logging data in S2 are as follows:

[0037]

[0038]

[0039]

[0040]

[0041] E = a + b×E dyn

[0042] v = c + d×vdyn

[0043] Among them, Δt s and Δt c are the shear wave and compressional wave slownesses respectively, obtained from the acoustic logging curve. ρ b is the density, obtained from the density logging curve. G dyn is the dynamic shear modulus, K dyn is the dynamic bulk modulus, E dyn is the dynamic Young's modulus, v dyn is the dynamic Poisson's ratio, E is the static Young's modulus, and v is the static Poisson's ratio. The Young's modulus and Poisson's ratio in S2 are E and v. a and b, c and d are obtained from the stress-strain experiment of the rock during the experiment.

[0044] The present invention is further configured as follows: When the acoustic logging data in S2 has the measured deep resistivity data, the RT curve is corrected using the Faust formula. When the measured deep resistivity of the acoustic logging data is missing, the acoustic wave slowness is reconstructed by GR normalization.

[0045] The Faust formula is as follows:

[0046]

[0047]

[0048] In the formula, V p is the compressional wave velocity, m / s. Depth is the depth, m. RT is the measured deep resistivity, Ω*m. Δt c is the compressional wave slowness, μs / m. K is the formation parameter, obtained by statistical analysis of seismic data.

[0049] The GR normalization formula is as follows:

[0050]

[0051] In the formula, Δt c is the compressional wave slowness, is the minimum value of the compressional wave slowness, is the maximum value of the compressional wave slowness. GR is the natural gamma value, GR min is the minimum value of the natural gamma, GR max is the maximum value of the natural gamma. GR min 、GR max are obtained by adjusting the scale through the overlapping display of the measured GR of the upper and lower surrounding rocks and DT. GR is obtained by actual measurement.

[0052] The present invention is further configured as follows: The value range of K is 2000 - 2500. are 40us / ft and 140us / ft respectively. GRmin 、GR max are 20 API and 300 API respectively.

[0053] The present invention is further configured that: in S2, the density logging data is corrected by using the Castagna formula or the Gardner formula on the basis of the acoustic logging data.

[0054] Castagna formula:

[0055] Gardner formula:

[0056] V p is the longitudinal wave velocity. ρ is the density. a 1 、b 1 、c 1 、d 1 、f are coefficients related to lithology, which are obtained based on seismic data statistics.

[0057] The present invention is further configured that: a 1 、b 1 、c 1 、d 1 、f for different lithologies are as follows:

[0058] Shale: a 1 =-0.0261, b 1 =0.373, c 1 =1.4 - 1.6, d 1 =1.75, f = 0.265.

[0059] Sandstone: a 1 =-0.0115, b 1 =0.261, c 1 =1.5 - 1.6, d 1 =1.66, f = 0.261.

[0060] Limestone: a 1 =-0.0296, b 1 =0.461, c 1 =0.9 - 1.1, d 1 =1.50, f = 0.225.

[0061] The present invention is further configured that: the steps for obtaining εx and εy are as follows:

[0062] Perform a mini - fracture test or a Kaiser experiment. Based on the test or experiment results, combined with the geological conditions, estimate the minimum horizontal principal stress σh and the maximum horizontal principal stress σH at discrete depth points, and calculate the effective stress ratios Kmin and Kmax at discrete depth points. The formulas are as follows:

[0063] $K_{min} = (\sigma_h - P_p) / (S_v - P_p)$

[0064] $K_{max} = (\sigma_H - P_p) / (S_v - P_p)$

[0065] Using the obtained $K_{min}$ and $K_{max}$, calculate the continuous minimum horizontal principal stress $\sigma_h$ and $\sigma_H$. The formulas are as follows:

[0066] $\sigma_h = K_{min}*(S_v - P_p)+P_p$

[0067] $\sigma_H = K_{max}*(S_v - P_p)+P_p$

[0068] Substitute the continuous maximum horizontal principal stress $\sigma_H$ and minimum horizontal principal stress $\sigma_h$ into the following formula to calculate $\epsilon_x$ and $\epsilon_y$:

[0069]

[0070]

[0071] In the formula, $\sigma$ H , Mpa. $\sigma$ h , Mpa. $P_p$ is the formation pore pressure, MPa, which is the measured value. $\alpha$ is the Biot coefficient, dimensionless. $S_v$ is the overburden pressure, MPa, which is the measured value. $v$ is the Poisson's ratio, dimensionless, which is the measured value. $E$ is the Young's modulus, MPa, which is the measured value. $\epsilon_x$ is the minimum tectonic stress coefficient, dimensionless. $\epsilon_y$ is the maximum tectonic stress coefficient, dimensionless.

[0072] The present invention is further configured as follows: The calculation formula for the formation pore pressure $P_p$ used when calculating $\epsilon_x$ and $\epsilon_y$ is as follows:

[0073] For clastic rock formations, the calculation formula for the formation pore pressure $P_p$ is as follows:

[0074]

[0075] In the formula, $S_v$ is the overburden pressure, $P_h$ is the hydrostatic pressure, $n$ is the Eaton index. $X$ is velocity, density, resistivity, where $X$ i is the measured value, $X$ Norm is the normal compaction value. $S_v$, $P_h$, $n$, $X$ i and $X$ Norm in the formula are all obtained according to the interpretation results of logging data.

[0076] For carbonate rock formations, the calculation formula for the formation pore pressure $P_p$ is as follows:

[0077] $p$ p $= S_v-(A e$ B*DT )

[0078] In the formula, Sv is the overlying formation pressure, DT is the acoustic travel time, and Sv and DT are obtained from the interpretation results of logging data. e is the natural constant, and A and B are coefficients. A and B are obtained from the statistical relationship between the effective stress σ and the acoustic travel time DT, and σ = Sv - p p , and are obtained by fitting Sv and Pp at a measured formation pressure point or a drilling mud overflow point at a certain depth. Among them, the effective stress σ, the acoustic travel time DT, and Sv and Pp at the depth of the measured formation pressure depth point or the drilling mud overflow point depth are all obtained from the interpretation results of logging data.

[0079] The present invention also provides a carbonate rock fracture prediction device based on rock mechanics, including a storage medium and a processor, and a computer program is stored on the storage medium. The processor is used to implement the above-mentioned carbonate rock fracture prediction method based on rock mechanics when executing the computer program.

[0080] In summary, the present invention has the following beneficial effects compared with the prior art: The present invention obtains the functional relationship between the lithofacies and reservoir parameters by using geostatistical analysis based on the reservoir parameters and lithofacies model of the drilled wells; establishes a reservoir parameter model through facies-controlled stochastic simulation according to this functional relationship; obtains a density model based on the reservoir parameter model, and obtains an overlying rock pressure model by integrating according to the density model; calculates the pore pressure of different formations by lithology; constructs a ground stress model based on the reservoir parameter model, the overlying rock pressure model, the pore pressure and the strain data based on single-well ground stress analysis. Based on the ground stress model, the stress and strain data are imported into the fracture parameter calculation model to realize the three-dimensional quantitative prediction of fractures in low-permeability reservoirs. Description of the Drawings

[0081] Figure 1 is a schematic flow chart of a carbonate rock fracture prediction method based on rock mechanics;

[0082] Figure 2 is a diagram of the logging curve correction result;

[0083] Figure 3 is a structural model;

[0084] Figure 4 is a grid model;

[0085] Figure 5 is a diagram of the rock mechanics parameter calculation result;

[0086] Figure 6 is a three-dimensional maximum horizontal principal stress model;

[0087] Figure 7 is a schematic diagram of the fracture development density. Detailed Embodiments

[0088] The technical solution of the present invention will be clearly described below in conjunction with the accompanying drawings. Obviously, the described embodiments are not all embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the protection scope of the invention.

[0089] Embodiment

[0090] As Figure 1 shown, it is a schematic flow chart of a carbonate rock fracture prediction method based on rock mechanics, which includes the following steps:

[0091] S1. Establish a lithofacies model:

[0092] S11. Based on the seismic fault interpretation results, with the geological stratification data of the drilled wells as the conditions, and using the marker horizons interpreted by seismic as the trend constraints, generate the structural horizons of the small layers. Then adjust the structural horizons according to the known geological profiles. As Figure 3 shown, generate a structural model of the research area that conforms to the true structural form. As Figure 4 shown, the structural model of the research area is a three-dimensional grid model, the grid direction is consistent with the line direction of seismic data acquisition, and the lateral grid step is consistent with the seismic data trace interval. In this embodiment, the lateral grid step is 25m * 25m. Determine the target research formation. The longitudinal grid step of the formation above the target research formation is 10m, and the longitudinal grid step of the target research formation is 5m.

[0093] S12. On the basis of the structural model of the research area, establish a sedimentary microfacies model of the research area according to the sedimentary microfacies interpretation results of the drilled wells.

[0094] S13. Conduct lithofacies interpretation based on well logging data. On the basis of the sedimentary microfacies model of the research area, establish a lithofacies model controlled by sedimentary microfacies according to the lithofacies interpretation results.

[0095] S2. Obtain well logging data, and the well logging data includes acoustic logging data and density logging data. As Figure 5 shown, according to the well logging data interpretation results, obtain the reservoir parameters of the drilled wells. The reservoir parameters of the drilled wells include density, seismic wave velocity, Young's modulus, and Poisson's ratio.

[0096] The calculation formulas for Young's modulus and Poisson's ratio based on well logging data are as follows:

[0097]

[0098]

[0099]

[0100]

[0101] E = a + b×E dyn

[0102] v = c + d×v dyn

[0103] where, Δt s and Δt c are the shear wave and compressional wave slownesses respectively, obtained from the acoustic logging curve. ρ b is the density, obtained from the density logging curve. G dyn is the dynamic shear modulus, K dyn is the dynamic bulk modulus, E dyn is the dynamic Young's modulus, v dyn is the dynamic Poisson's ratio, E is the static Young's modulus, v is the static Poisson's ratio. E and v are the Young's modulus and Poisson's ratio in the reservoir parameters of the drilled wells. a and b, c and d are obtained from the stress-strain experiment of the rock during the experiment.

[0104] S3. Based on the reservoir parameters of the drilled wells and the lithofacies model, the functional relationship between lithofacies and reservoir parameters is obtained by geostatistical analysis. Based on this functional relationship, a reservoir parameter model is established by facies-controlled stochastic simulation.

[0105] S4. According to the reservoir parameter model, a density model is obtained, and based on the density model, an overburden pressure model is obtained by integration.

[0106] The integral formula for calculating the overburden pressure at the section burial depth point H is:

[0107]

[0108] In the formula, P H is the overburden pressure at point H of the section, MPa; h is the formation burial depth, m; g is the acceleration of gravity, m / s2; ρ w is the formation water density, g / cm3; H is the section burial depth, m.

[0109] S5. Calculate the pore pressure of different formations by lithology.

[0110] For clastic rock formations, the formula for calculating the formation pore pressure Pp is as follows:

[0111]

[0112] Wherein, Pp is the formation pore pressure, in Mpa; Sv is the overburden pressure, in MPa, obtained through the overburden rock pressure model; Ph is the hydrostatic pressure, in Mpa, which is the product of the acceleration of gravity, water density and depth value; n is the Eaton index, dimensionless, and the Eaton coefficient of different formations is obtained based on the analysis results of drilled wells; V is the seismic wave velocity and density, where Vo is obtained through the reservoir parameter model; Vn is the normal compaction value, obtained based on the normal compaction trend line of drilled wells.

[0113] For carbonate formations, the calculation formula for the formation pore pressure Pp is as follows:

[0114] p p = Sv - (Ae B*DT )

[0115] Wherein, Sv is the overburden pressure, in MPa, obtained through the overburden rock pressure model; DT is the acoustic travel time, in μs / m, obtained through the interpretation results of well logging data; e is the natural constant, and A and B are coefficients, which are obtained based on the analysis of carbonate formation data of drilled wells. In this embodiment, A is 224 and B is -0.0193.

[0116] S6. As Figure 6 shown, construct the in-situ stress model, and the formula is as follows:

[0117]

[0118]

[0119] Wherein, S H is the maximum horizontal principal stress, in Mpa; S h is the minimum horizontal principal stress, in Mpa; Pp is the formation pore pressure, in MPa, obtained through S5; α is the Biot coefficient, dimensionless; Sv is the overburden pressure, in MPa, obtained through the overburden rock pressure model; v is the Poisson's ratio, dimensionless, obtained through the reservoir parameter model; E is the Young's modulus, in MPa, obtained through the reservoir parameter model; εx is the minimum tectonic stress coefficient, dimensionless; εy is the maximum tectonic stress coefficient, dimensionless. The εx and εy of different formations are obtained based on the drilled well data.

[0120] After the in-situ stress parameter model is constructed, use the rock failure criterion to simulate the wellbore shear sloughing situation, and then compare it with the wellbore condition. If the wellbore sloughing width is consistent with the wellbore diameter enlargement, it means that the in-situ stress parameter model fits the actual situation. If the wellbore sloughing width is inconsistent with the wellbore diameter enlargement, then return to any step in S1 - S5 until the wellbore sloughing width is consistent with the wellbore diameter enlargement.

[0121] The rock fracture criterion is selected based on the in-situ stress model and the engineering mud density. The selection basis is to simulate the borehole breakout width under different fracture criteria and compare it with the borehole diameter enlargement, and select the rock fracture criterion with the most consistent simulated borehole breakout width and borehole diameter enlargement. The rock fracture criteria include the Mohr-Coulomb criterion, the modified Lade criterion, Tresca, and the Circumscribed Drucker-Prager. In this embodiment, the Mohr-Coulomb criterion is selected as the rock fracture criterion.

[0122] S7. As Figure 7 shown, calculate the fracture parameters, and the formula is as follows:

[0123]

[0124]

[0125]

[0126] In the formula, σ p is the rock fracture stress, in Mpa, and σ p of different formations is obtained based on the analysis of the drilled well data. σ1, σ2, and σ3 are the maximum, intermediate, and minimum effective principal stresses, respectively, in MPa. L1, L2, and L3 are the unit lengths along the principal stress directions σ1, σ2, and σ3, in m, and are calculated based on σ1, σ2, σ3, and the Young's modulus and Poisson's ratio obtained through the reservoir parameter model. μ is the Poisson's ratio, dimensionless, and is obtained through the reservoir parameter model. θ is the fracture rupture angle, in degrees, and is calculated based on the in-situ stress model. ε3 is the minimum principal strain, dimensionless, and is calculated based on σ3 and the Young's modulus obtained through the reservoir parameter model. ε0 is the maximum elastic tensile strain borne by the rock, dimensionless, and is calculated based on the maximum horizontal principal stress obtained through the in-situ stress model and the Young's modulus obtained through the reservoir parameter model. J is the fracture surface energy, in J / m2, and is obtained by fitting based on the formation structure of the study area, the Young's modulus obtained through the reservoir parameter model, and the in-situ stress model. Dvf is the tectonic fracture volume density, in m2 / m3. Dlf is the tectonic fracture line density, in fractures / m. b is the effective aperture of the tectonic fracture.

[0127] Specifically, σ1, σ2, and σ3 are used to draw Mohr's circles in the stress space according to the in-situ stress model. The diameter of Mohr's circle is the difference between the maximum principal stress and the minimum principal stress. According to the effective stress principle, the maximum effective principal stress is half of the diameter of Mohr's circle plus the normal stress value at the origin, and the minimum effective principal stress is half of the diameter of Mohr's circle minus the normal stress value at the origin. The intermediate effective principal stress is equal to the result of dividing the sum of the maximum and minimum effective principal stresses by two through the Mohr-Coulomb criterion.

[0128] Specifically, L1, L2, and L3 are calculated based on the shear strains corresponding to σ1, σ2, and σ3 through L = (σ1 - σ3) / (2*G*γ), where L is the element length; γ is the shear strain, obtained based on the strain measurement data of the rock or soil; L is the element length, and G is the shear modulus, obtained by calculating G = E / (2*(1 + ν)) from the Young's modulus E and Poisson's ratio ν obtained from the reservoir parameter model.

[0129] Specifically, θ = arctan((σ1 - σ3) / 2c), where c is the cohesion of the rock or soil, obtained through experiments.

[0130] Specifically, as Figure 2 shown, in this embodiment, based on the effective medium theory, the upper and lower limit values of the composite rock modulus and velocity are simulated or estimated according to different mineral compositions and pore fluid properties. Based on the rock physics boundary theory, the logging data of the work area is diagnostically analyzed through the crossplot method according to the rock physics diagnostic model, and the problematic logging data is corrected. This makes the corrected logging data more in line with the rock physics variation law. The correction steps are as follows:

[0131] When the acoustic logging data has measured deep resistivity data, the RT curve is corrected using the Faust formula. When the measured deep resistivity of the acoustic logging data is missing, the acoustic travel time is reconstructed using GR normalization.

[0132] The Faust formula is as follows:

[0133]

[0134]

[0135] In the formula, V p is the longitudinal wave velocity, m / s. Depth is the depth, m. RT is the measured deep resistivity, Ω*m. Δt c is the longitudinal wave travel time, μs / m. K is the formation parameter, obtained by statistical analysis of seismic data. In this embodiment, the value range of K is 2000 - 2500.

[0136] The GR normalization formula is as follows:

[0137]

[0138] In the formula, Δt c is the longitudinal wave travel time, is the minimum longitudinal wave travel time, is the maximum longitudinal wave travel time. GR is the natural gamma value, GR min is the minimum natural gamma value, GR max is the maximum natural gamma value. GRmax Obtained by adjusting the scale through the overlapping display of the measured GR and DT of the upper and lower surrounding rocks. GR is obtained through actual measurement. In this embodiment are 40 us / ft and 140 us / ft respectively. GR min 、GR max are 20 API and 300 API respectively.

[0139] The density logging data is corrected using the Castagna formula or the Gardner formula based on the acoustic logging data.

[0140] Castagna formula:

[0141] Gardner formula:

[0142] V p is the longitudinal wave velocity. ρ is the density. a 1 、b 1 、c 1 、d 1 、f are coefficients related to lithology, obtained by statistical analysis of seismic data. In this embodiment, a 1 、b 1 、c 1 、d 1 、f for different lithologies are as follows:

[0143] Shale: a 1 =-0.0261, b 1 =0.373, c 1 =1.4 - 1.6, d 1 =1.75, f = 0.265.

[0144] Sandstone: a 1 =-0.0115, b 1 =0.261, c 1 =1.5 - 1.6, d 1 =1.66, f = 0.261.

[0145] Limestone: a 1 =-0.0296, b 1 =0.461, c 1 =0.9 - 1.1, d 1 =1.50, f = 0.225.

[0146] Specifically, the steps for obtaining εx and εy used in S6 are as follows:

[0147] Conduct small-scale fracturing tests or Kaiser experiments. Based on the test or experiment results and combined with geological conditions, through empirical formulas such as the Mohr-Coulomb criterion and -Brown criterion, estimate the minimum horizontal principal stress σh and the maximum horizontal principal stress σH at discrete depth points, and calculate the effective stress ratios Kmin and Kmax at discrete depth points. The formulas are as follows:

[0148] Kmin = (σh - Pp) / (Sv - Pp)

[0149] Kmax = (σH - Pp) / (Sv - Pp)

[0150] Using the obtained Kmin and Kmax, calculate the continuous minimum horizontal principal stress σh and σH. The formulas are as follows:

[0151] σh = Kmin * (Sv - Pp) + Pp

[0152] σH = Kmax * (Sv - Pp) + Pp

[0153] Substitute the continuous maximum horizontal principal stress σH and minimum horizontal principal stress σh into the following formula to calculate εx and εy:

[0154]

[0155]

[0156] In the formula, σ H , Mpa. σ n , Mpa. Pp is the formation pore pressure, MPa, which is the measured value. α is the Biot coefficient, dimensionless. Sv is the overburden pressure, MPa, which is the measured value. v is the Poisson's ratio, dimensionless, which is the measured value. E is the Young's modulus, MPa, which is the measured value. εx is the minimum tectonic stress coefficient, dimensionless. εy is the maximum tectonic stress coefficient, dimensionless.

[0157] Specifically, the calculation formula for the formation pore pressure Pp used when calculating εx and εy is as follows:

[0158] For clastic rock formations, the calculation formula for the formation pore pressure Pp is as follows:

[0159]

[0160] In the formula, Sv is the overburden pressure, Ph is the hydrostatic pressure, n is the Eaton index. X is velocity, density, resistivity, where X i is the measured value, X Norm is the normally compacted value. Sv, Ph, n, X i and X NormAll are obtained according to the interpretation results of logging data.

[0161] For carbonate formations, the calculation formula for formation pore pressure Pp is as follows:

[0162] p p = Sv - (A e B*DT )

[0163] In the formula, Sv is the overburden pressure, DT is the acoustic travel time. Sv and DT are obtained according to the interpretation results of logging data. e is the natural constant, and A and B are coefficients. A and B are obtained according to the statistical relationship between the effective stress σ and the acoustic travel time DT, and σ = Sv - p p , and are obtained by fitting Sv and Pp at the measured formation pressure point or the drilling mud overflow point at a certain depth. Among them, the effective stress σ, the acoustic travel time DT, the Sv and Pp at the measured formation pressure depth point or the drilling mud overflow point depth are all obtained according to the interpretation results of logging data. The A here is the same as the A used in S5, and the B here is the same as the B used in S5.

[0164] This embodiment also provides a carbonate fracture prediction device based on rock mechanics, including a storage medium and a processor. A computer program is stored on the storage medium. The processor is configured to implement the above-mentioned carbonate fracture prediction method based on rock mechanics when executing the computer program.

[0165] This embodiment obtains the functional relationship between the lithofacies and reservoir parameters by using geostatistical analysis based on the reservoir parameters and lithofacies model of the drilled wells; establishes a reservoir parameter model through facies-controlled stochastic simulation based on this functional relationship; obtains a density model according to the reservoir parameter model, and obtains an overlying rock pressure model by integrating based on the density model; calculates the pore pressure of different formations by lithology; constructs a in-situ stress model based on the reservoir parameter model, the overlying rock pressure model, the pore pressure and the strain data based on single-well in-situ stress analysis. Based on the in-situ stress model, the stress and strain data are imported into the fracture parameter calculation model, realizing the three-dimensional quantitative prediction of fractures in low-permeability reservoirs.

[0166] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention can have various changes and modifications. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A carbonate rock fracture prediction method based on rock mechanics, characterized in that: It includes the following steps: S1. Establish a lithofacies model; S2. Obtain logging data, which includes acoustic logging data and density logging data; according to the logging data interpretation results, obtain the reservoir parameters of the drilled wells, and the reservoir parameters of the drilled wells include density, seismic wave velocity, Young's modulus, and Poisson's ratio; S3. Based on the reservoir parameters of the drilled wells and the lithofacies model, use geostatistical analysis to obtain the functional relationship between lithofacies and reservoir parameters; based on this functional relationship, establish a reservoir parameter model through facies-controlled stochastic simulation; S4. Obtain a density model according to the reservoir parameter model, and integrate according to the density model to obtain an overlying rock pressure model; S5. Calculate the pore pressure of different strata by lithology; For clastic rock strata, the calculation formula of formation pore pressure Pp is as follows: In the formula, Pp is the formation pore pressure, Mpa; Sv is the overlying formation pressure, MPa, obtained through the overlying rock pressure model; Ph is the hydrostatic pressure, Mpa, which is the product of gravitational acceleration, water density, and depth value; n is the Eaton index, dimensionless, and the Eaton coefficient of different strata is obtained according to the drilled well data; V is the seismic wave velocity and density, where Vo is obtained through the reservoir parameter model; Vn is the normal compaction value, obtained according to the normal compaction trend line of the drilled wells; For carbonate rock strata, the calculation formula of formation pore pressure Pp is as follows: p p = Sv - (Ae B*DT ) In the formula, Sv is the overlying formation pressure, MPa, obtained through the overlying rock pressure model; DT is the acoustic travel time, us / m, obtained through the logging data interpretation results; e is the natural constant, and A and B are coefficients, and A and B are obtained according to the carbonate rock formation data of the drilled wells; S6. Construct a geostress model, and the formula is as follows: Where S H is the maximum horizontal principal stress, Mpa; S h is the minimum horizontal principal stress, Mpa; Pp is the formation pore pressure, MPa, obtained through S5; α is the Biot coefficient, dimensionless; Sv is the overburden pressure, MPa, obtained through the overburden rock pressure model; v is the Poisson's ratio, dimensionless, obtained through the reservoir parameter model; E is the Young's modulus, MPa, obtained through the reservoir parameter model; εx is the minimum tectonic stress coefficient, dimensionless; εy is the maximum tectonic stress coefficient, dimensionless; εx and εy of different strata are obtained according to the drilled well data; S7. Calculate the fracture parameters, and the formula is as follows: where, σ p is the rock fracture stress, in Mpa, and σ of different formations p is obtained based on the drilled well data; σ1, σ2, and σ3 are the maximum, intermediate, and minimum effective principal stresses, in MPa, obtained by calculating according to the in-situ stress model; L1, L2, and L3 are the unit lengths along the principal stress directions σ1, σ2, and σ3, in m, obtained by calculating based on σ1, σ2, σ3, as well as the Young's modulus and Poisson's ratio obtained through the reservoir parameter model; μ is the Poisson's ratio, dimensionless, obtained through the reservoir parameter model; θ is the fracture rupture angle, in °, obtained by calculating according to the in-situ stress model; ε3 is the minimum principal strain, dimensionless, obtained by calculating based on σ3 and the Young's modulus obtained through the reservoir parameter model; ε0 is the maximum elastic tensile strain borne by the rock, dimensionless, obtained by calculating based on the maximum horizontal principal stress obtained through the in-situ stress model and the Young's modulus obtained through the reservoir parameter model; J is the fracture surface energy, in J / m2, obtained by fitting based on the formation structure of the study area, the Young's modulus obtained through the reservoir parameter model, and the in-situ stress model; Dvf is the tectonic fracture volume density, in m2 / m3; Dlf is the tectonic fracture line density, in fractures / m; b is the effective aperture of the tectonic fracture.

2. A carbonate rock fracture prediction method based on rock mechanics according to claim 1, characterized in that: The S1 step is: S11. Based on the seismic fault interpretation results, taking the geological stratification data of the drilled wells as the condition, using the marker bed of the seismic interpretation as the trend constraint, generate the structural plane of the small layer; then adjust the structural plane according to the known geological section to generate a research area structural model that conforms to the true structural form; S12. On the basis of the research area structural model, establish a research area sedimentary microfacies model according to the sedimentary microfacies interpretation results of the drilled wells; S13. Conduct lithofacies interpretation based on logging data, and on the basis of the research area sedimentary microfacies model, establish a lithofacies model controlled by sedimentary microfacies according to the lithofacies interpretation results.

3. A carbonate rock fracture prediction method based on rock mechanics according to claim 2, characterized in that: The research area structural model established in S11 is a three-dimensional grid model, the grid direction is consistent with the line track direction of seismic data acquisition, and the lateral grid step length is consistent with the seismic data trace interval.

4. A carbonate rock fracture prediction method based on rock mechanics according to claim 3, characterized in that: The lateral grid step is 25m * 25m; the target research formation is determined, and the vertical grid step of the formation above the target research formation is 10m, while the vertical grid step of the target research formation is 5m.

5. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: After the in-situ stress parameter model in S6 is constructed, the borehole shear sloughing situation is simulated using the rock failure criterion, and then compared with the borehole condition. When the borehole sloughing width is consistent with the borehole diameter enlargement, S7 is carried out; The rock failure criterion is selected according to the in-situ stress model and the engineering mud density. The selection basis is to simulate the borehole sloughing width under different failure criteria and compare it with the borehole diameter enlargement, and select the rock failure criterion with the most consistent simulated borehole sloughing width and borehole diameter enlargement; the rock failure criteria include Mohr-Coulomb criterion, modified Lade criterion, Tresca, and Circumscribed Drucker-Prager.

6. A method for predicting carbonate rock fractures based on rock mechanics according to claim 5, characterized in that: The Mohr-Coulomb criterion is selected as the rock failure criterion.

7. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: In S5, for carbonate rock formations, the calculation formula for formation pore pressure Pp is as follows: p p = Sv - (224e -0.0193DT ).

8. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: The calculation formulas for Young's modulus and Poisson's ratio in S2 based on well logging data are as follows: E = a + b×E dyn V = c + d×v dyn where, Δt s and Δt c are the shear wave and compressional wave slownesses respectively, obtained from the acoustic logging curve; ρ b is the density, obtained from the density logging curve; G dyn is the dynamic shear modulus, K dyn is the dynamic bulk modulus, E dyn is the dynamic Young's modulus, V dyn is the dynamic Poisson's ratio, E is the static Young's modulus, v is the static Poisson's ratio; the Young's modulus and Poisson's ratio in S2 are E and v; a and b, c and d are obtained from the stress-strain experiment of the rock during the experiment.

9. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: When the acoustic well logging data in S2 has measured deep resistivity data, the RT curve is corrected using the Faust formula; when the measured deep resistivity of the acoustic well logging data is missing, the acoustic travel time is reconstructed using GR normalization; The Faust formula is as follows: where, V p is the longitudinal wave velocity, m / s; Depth is the depth, m; RT is the measured deep resistivity, Ω*m; Δt c is the longitudinal wave travel time difference, μs / m; K is the formation parameter, obtained by statistical analysis of seismic data; The GR normalization formula is as follows: where Δt c is the compressional wave travel time difference, is the minimum value of the compressional wave travel time difference, is the maximum value of the compressional wave travel time difference; GR is the natural gamma value, GR min is the minimum value of the natural gamma, GR max is the maximum value of the natural gamma; GR max is obtained by adjusting the scale through the overlay display of the measured GR of the upper and lower surrounding rocks and DT; GR is obtained through actual measurement.

10. A method for predicting carbonate rock fractures based on rock mechanics according to claim 9, characterized in that: The value range of K is 2000 - 2500; They are 40 us / ft and 140 us / ft respectively; GR min , GR max They are 20 API and 300 API respectively.

11. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: The density well logging data in S2 is corrected using the Castagna formula or the Gardner formula based on the acoustic well logging data; Castagna formula: Gardner formula: V p is the longitudinal wave velocity; ρ is the density; a 1 , b 1 , c 1 , d 1 , and f are coefficients related to lithology, which are obtained based on seismic data statistics.

12. A method for predicting carbonate rock fractures based on rock mechanics according to claim 11, characterized in that: a of different lithologies 1 、b 1 、c 1 、d 1 、f are as follows: Shale: a 1 = -0.0261, b 1 = 0.373, c 1 = 1.4 - 1.6, d 1 = 1.75, f = 0.265; Sandstone: a 1 = -0.0115, b 1 = 0.261, c 1 = 1.5 - 1.6, d 1 = 1.66, f = 0.261; Limestone: a 1 = -0.0296, b 1 = 0.461, c 1 = 0.9 - 1.1, d 1 = 1.50, f = 0.

225.

13. A method for predicting carbonate rock fractures based on rock mechanics according to any one of claims 1 - 4, characterized in that: The steps for obtaining εx and εy are as follows: A small-scale fracturing test or Kaiser experiment is carried out. Based on the test or experiment results, combined with the geological conditions, the minimum horizontal principal stress σh and the maximum horizontal principal stress σH at discrete depth points are estimated, and the effective stress ratios Kmin and Kmax at discrete depth points are calculated. The formulas are as follows: Kmin = (σh - Pp) / (Sv - Pp) Kmax = (σH - Pp) / (Sv - Pp) Using the obtained Kmin and Kmax, calculate the continuous minimum horizontal principal stress σh and σH. The formulas are as follows: σh = Kmin * (Sv - Pp) + Pp σH = Kmax * (Sv - Pp) + Pp Substitute the continuous maximum horizontal principal stress σH and minimum horizontal principal stress σh into the following formula to calculate εx and εy: where σ H , Mpa; σ h , Mpa; Pp is the formation pore pressure, MPa, which is the measured value; a is the Biot coefficient, dimensionless; Sv is the overburden pressure, MPa, which is the measured value; v is the Poisson's ratio, dimensionless, which is the measured value; E is the Young's modulus, MPa, which is the measured value; εx is the minimum tectonic stress coefficient, dimensionless; εy is the maximum tectonic stress coefficient, dimensionless.

14. A carbonate rock fracture prediction method based on rock mechanics according to claim 13, characterized in that: The calculation formula for the formation pore pressure Pp used when calculating εx and εy is as follows: For clastic rock formations, the calculation formula for the formation pore pressure Pp is as follows: Wherein, Sv is the overlying formation pressure, Ph is the hydrostatic pressure, and n is the Eaton exponent; X is velocity, density, resistivity, where X i is the measured value, and X Nonu is the normally compacted value; Sv, Ph, n, X i and X Norm in the formula are all obtained based on the interpretation results of logging data; For carbonate rock formations, the calculation formula for the formation pore pressure Pp is as follows: p p = Sv - (A e B*DT ) Wherein, Sv is the overlying formation pressure, DT is the acoustic travel time, and Sv and DT are obtained according to the interpretation results of logging data; e is the natural constant, A and B are coefficients, and A and B are obtained according to the statistical relationship between the effective stress σ and the acoustic travel time DT and σ = Sv - p p , and are obtained by fitting Sv and Pp at a measured formation pressure point or a drilling mud overflow point at a certain depth; wherein the effective stress σ, the acoustic travel time DT, and Sv and Pp at the depth of the measured formation pressure depth point or the drilling mud overflow point are all obtained according to the interpretation results of logging data.

15. A carbonate rock fracture prediction device based on rock mechanics, characterized in that: It includes a storage medium and a processor. A computer program is stored on the storage medium; the processor is used to implement the carbonate rock fracture prediction method based on rock mechanics according to any one of claims 1-14 when executing the computer program.

Images