A method for predicting formation pore pressure in fractured reservoirs

Abstract

This invention discloses a method for predicting formation pore pressure in fractured reservoirs, relating to the fields of oil and gas exploration and development and geophysical logging technology. The method includes obtaining logging and geological parameters of the target formation; establishing a relationship between pore pressure and fracture aspect ratio based on the effective stress principle and fracture mechanics; calculating the theoretical P-wave velocity under different pore pressure conditions within a pressure range defined by hydrostatic pressure and normal stress on the fracture surface, using an equivalent medium model and the Gassmann equation; determining the pore pressure based on the residual between the theoretical and measured P-wave velocities; and comparing the calculated pore pressure value with the measured pressure value to establish an evaluation relationship. This invention has the advantages of reducing unnecessary and ineffective calculations, improving the stability and speed of the calculation process, and enhancing the physical rationality of the results.

Description

Technical Field

[0001] This invention relates to the fields of oil and gas exploration and development and geophysical logging technology, specifically to a method for predicting formation pore pressure in fractured reservoirs. Background Technology

[0002] In deep and ultra-deep oil and gas exploration and development projects, accurate prediction of formation pore pressure is a core prerequisite for ensuring drilling operation safety, designing reasonable mud density windows, and preventing blowouts and lost circulation accidents.

[0003] Currently, the mainstream pore pressure prediction methods in the industry (such as the Eaton method and the Bowers method) are mainly based on the "undercompactment theory" of clastic rocks. This theory holds that formation porosity decreases exponentially with the increase of overlying effective stress, and abnormally high pressure will lead to abnormally high porosity, which will be reflected in the decrease of sonic velocity on the logging curve.

[0004] However, in fractured reservoirs, acoustic velocity is controlled by the complex pore structure, not just by porosity. Especially in systems containing well-developed microfractures, the opening or closing state of fractures directly determines the formation's acoustic response. Existing prediction methods often treat the fracture aspect ratio as a fuzzy variable that varies with pressure, and use empirical attenuation functions for calculations. The limitations of this approach are: first, the range of values ​​for empirical parameters is subjective, leading to significant differences in pressure results obtained by different interpreters; second, the model fails to establish the mechanical relationship between wave velocity, fracture geometry, and effective stress.

[0005] In actual production, such predictions lacking physical support often manifest as "high pressure underreporting" or "false high pressure misreporting," severely restricting the exploration efficiency of deep and complex oil and gas reservoirs. Summary of the Invention

[0006] To address the problem that existing methods for predicting formation pore pressure in fractured reservoirs mainly rely on empirical parameters and lack physical support, often resulting in "high pressure underreporting" or "false high pressure misreporting," this invention provides a method for predicting formation pore pressure in fractured reservoirs to solve this problem.

[0007] The technical solution adopted by this invention to solve its technical problem is:

[0008] A method for predicting formation pore pressure in fractured reservoirs includes the following steps:

[0009] S1: Obtain well logging and geological parameters at the target formation depth, and calculate the hydrostatic pressure and fracture surface normal stress σ at this location.n And combined with the preset pore pressure value P set for the target formation p Calculate the effective stress σ of the normal phase acting on the crack surface. eff ;

[0010] S2: Young's modulus E of the combined pure mineral matrix s Compared with Poisson's ratio ν s Construct the relationship between the crack aspect ratio α and the theoretical crack closure pressure P. c The functional relationship;

[0011] S3: The effective stress σ of the crack surface normal phase eff With the theoretical closure pressure P of the crack c Establish an equilibrium relationship, namely the effective stress σ of the crack surface normal phase. eff Equal to the theoretical closure pressure P of the crack c And construct a preset value P for pore pressure. p The relationship between the crack aspect ratio α and the crack width ratio α;

[0012] S4: Based on the equivalent medium theoretical model and fluid substitution of the Gassmann equation, the calculation is based on the current preset pore pressure P. p The theoretical P-wave velocity V of the underlying strata p,calc ;

[0013] S5: Based on the theoretical longitudinal wave velocity V p,calc Compared with the measured longitudinal wave velocity V p,meas The difference between them is the matching target. Within the pressure range defined by hydrostatic pressure and normal stress on the crack surface, a trial-and-error method is used to match the preset pore pressure P. p The solution is performed iteratively, and the final output is the calculated pore pressure value P of the target formation. p,cal ;

[0014] S6: Verification, using the calculated pore pressure value P of the target formation. p,cal The measured pore pressure value P of the target formation p,obs Compare and establish evaluation relationships.

[0015] Preferably, the well logging and geological parameters for the target formation depth in step S1 include, but are not limited to, the mineral volume fraction f. i Total porosity φ total Formation fracture porosity φ frac pore fluid bulk modulus K fl Density ρ fl , crack direction, crack surface dip angle θ, and the angle γ between the crack direction and the direction of the maximum horizontal principal stress;

[0016] It also includes the measured P-wave velocity V at the target formation depth. p,meas Three-dimensional geostress parameters;

[0017] The three-dimensional geostress parameters include at least the vertical stress σ of the overlying strata. v Maximum horizontal principal stress σ H Minimum horizontal principal stress σ h .

[0018] Preferably, in step S1, the normal stress σn of the crack surface is calculated using tensor projection, and the formula is as follows:

[0019] ;

[0020] In the formula, σ n The normal stress along the crack surface is σ, MPa; H The maximum horizontal principal stress is σ, MPa; h The minimum horizontal principal stress is given in MPa; σ v θ is the vertical stress of the overlying strata, MPa; θ is the dip angle of the fracture surface, i.e., the angle between the fracture surface and the horizontal plane, ° (degrees); γ is the angle between the fracture direction and the direction of the maximum horizontal principal stress, ° (degrees).

[0021] Preferably, in step S1, the effective stress σ of the crack surface normal phase eff The calculation formula is:

[0022] ;

[0023] In the formula, σ eff The effective stress of the crack surface is σ (MPa); n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

[0024] Preferably, in step S2, the crack aspect ratio α and the theoretical crack closure pressure P c The functional relationship is as follows:

[0025] ;

[0026] In the formula, P c E represents the theoretical closure pressure of the crack, in MPa. s Young's modulus of pure mineral matrix, MPa, ν s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless.

[0027] Preferably, in step S3, the preset pore pressure value P p The relationship between the crack aspect ratio α and the crack aspect ratio α is as follows:

[0028] ;

[0029] In the formula, E s Young's modulus of pure mineral matrix, MPa, ν s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless; σ n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

[0030] Preferably, in step S4, the current pore pressure preset value P p The theoretical P-wave velocity V of the underlying strata p,calc The calculation process is as follows:

[0031] S41: Based on the preset pore pressure value P p And step S3 obtains the crack aspect ratio α(P) corresponding to the pore pressure. p The mineral composition, total porosity, fracture porosity, and fracture aspect ratio are input into the equivalent medium model to calculate the pore pressure preset value P based on the currently set value. p Bulk modulus K of dry rock skeleton under certain conditions dry (P p ), GPa and dry rock skeleton shear modulus μ dry (P p ), GPa;

[0032] S42: Mineral matrix bulk modulus K s and mineral matrix shear modulus μ s The elastic modulus of each mineral component can be obtained by averaging the Voigt-Reuss-Hill values:

[0033] ;

[0034] ;

[0035] In the formula, f i K represents the volume fraction of the i-th mineral, dimensionless; i and μ i Let be the bulk modulus and shear modulus of the i-th mineral, respectively, in GPa;

[0036] S43: Fluid substitution is performed using the Gassmann equation to obtain the preset pore pressure value P based on the current settings. p Saturated rock bulk modulus K under certain conditions sat (P p ) and saturated rock shear modulus μ sat (P p The calculation formula is:

[0037] ;

[0038] ;

[0039] In the formula, K sat (P p P is the preset value of pore pressure. p Saturated rock bulk modulus under certain conditions, GPa; K dry (P p P is the preset value of pore pressure. p Bulk modulus of dry rock skeleton under the given conditions, GPa; K s The equivalent bulk modulus of the mineral matrix is ​​expressed in GPa and K. fl φ is the bulk modulus of pore fluid, in GPa; total Total porosity, dimensionless; μ sat (P p P is the preset value of pore pressure. p Saturated rock shear modulus under certain conditions, GPa; μ dry (P p P is the preset value of pore pressure. p Shear modulus of dry rock skeleton under the given conditions, GPa;

[0040] S44: Calculate the density of saturated rock using the following formula:

[0041] ;

[0042] In the formula, ρ sat φ is the density of saturated rock, g / cm³. total ρ represents total porosity, dimensionless. s The equivalent density of the mineral matrix is ​​g / cm³; ρ fl f is the pore fluid density, g / cm³. i ρ is the volume fraction of the i-th mineral, dimensionless; i Let be the density of the i-th mineral, in g / cm³;

[0043] S45: Calculate the preset pore pressure P based on the saturated rock bulk modulus, shear modulus, and density. p Theoretical P-wave velocity V under certain conditions p,calc The calculation formula is:

[0044] ;

[0045] In the formula, V p,calc (P p P is the preset value of pore pressure. p Theoretical P-wave velocity under the given conditions, km / s; K sat (P p P is the preset value of pore pressure. pSaturated rock bulk modulus under certain conditions, GPa; μ sat (P p P is the preset value of pore pressure. p Saturated rock shear modulus under certain conditions, GPa; ρ sat ρ is the density of saturated rock, in g / cm³.

[0046] Preferably, in step S5, the trial algorithm is an iterative solution method based on interval interpolation, and the specific steps are as follows:

[0047] S51: Determine adaptive tolerance: Extract all measured P-wave velocities V from well logging. p,meas Set tolerance ε v= η⋅V p,meas Where η is a preset proportionality coefficient, and here η is taken as 0.01%;

[0048] S52: Determining the physical search boundary: Setting the lower bound P of the pore pressure search range. min For hydrostatic pressure, the upper limit is P. max Normal stress σ on the crack surface n ;

[0049] S53: Calculation of pore pressure at initial depth point: For the initial depth point, at the lower boundary P of the pore pressure calculation interval... min and the upper boundary P max Theoretical P-wave velocity V is calculated at this location. p,calc The theoretical longitudinal wave velocity V corresponding to the lower boundary pore pressure is obtained. p,calc (P min The theoretical longitudinal wave velocity V corresponding to the upper pore pressure. p,calc (P max ); using P min P max With V p,calc (P min V p,calc (P max The first interpolation of the corresponding pressure-wave velocity relationship yields the first prediction point P1:

[0050] ;

[0051] In the formula, P1 is the predicted pore pressure obtained from the first interpolation, in MPa; P min The lower bound of the pore pressure search range is given in MPa; P max The upper bound of the pore pressure search range is given in MPa; V p,meas For the measured P-wave velocity, km / s; V p,calc (P min P is the preset value for pore pressure. min The theoretical P-wave velocity calculated at time, km / s; Vp,calc (P max P is the preset value for pore pressure. max The theoretical P-wave velocity calculated at that time, km / s;

[0052] S54: Pore pressure subsequent interpolation update: Let P be the low-pressure side endpoint of the current interpolation interval at the current depth. L The corresponding theoretical longitudinal wave velocity is V p,calc (P L The high-pressure side endpoint is P. R The corresponding theoretical longitudinal wave velocity is V p,calc (P R Based on the pressure at both ends of the current interpolation interval and its theoretical P-wave velocity, calculate the k-th predicted pressure P. k :

[0053] ;

[0054] In the formula, P k The predicted pore pressure obtained from the k-th interpolation is expressed in MPa and V. p,meas For the measured P-wave velocity, km / s; P L The left endpoint of the current interpolation interval is the low-pressure side endpoint, in MPa; P R The right endpoint of the current interpolation interval is the high-pressure side endpoint, in MPa; V p,calc (P L P is the preset value for pore pressure. L The calculated theoretical P-wave velocity; V p,calc (P R P is the preset value for pore pressure. R The calculated theoretical longitudinal wave velocity;

[0055] S55: Pore pressure calculation at subsequent depth points: For subsequent depth points, the pore pressure result P obtained from the previous depth point is used... prev As a reference pressure at the current depth point;

[0056] First, calculate the reference pressure P under the formation parameter conditions at the current depth point. prev The corresponding theoretical longitudinal wave velocity is denoted as V. p,calc (P prev And calculate the reference residual:

[0057] r prev =V p,calc (P prev )−V p,meas ;

[0058] In the formula, r prev P represents the wave velocity residual corresponding to the reference pressure, in km / s. prevThe pore pressure result calculated at the previous depth point is in MPa; V p,calc (P prev P is the preset value of pore pressure under the current depth conditions. prev The theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity at the current depth point is in km / s;

[0059] S56: Repeat steps S54 and S55 to perform subsequent depth point pore pressure interpolation updates.

[0060] S57: Output Results: The theoretical P-wave velocity V corresponding to any predicted pressure during the calculation of the initial or subsequent depth points. p,calc Compared with the measured longitudinal wave velocity V p,meas When the residual between the two conditions meets the tolerance requirement, the predicted pressure is output as the pore pressure result at the corresponding depth point, i.e., the calculated pore pressure value P of the target formation. p,cal .

[0061] Preferably, in step S6, the evaluation relationship is established using a relative error rate, which is calculated as follows:

[0062] ;

[0063] ;

[0064] In the formula, ΔP is the difference between the predicted pore pressure value (i.e., the calculated pore pressure value) and the measured pore pressure value, in MPa; P p,cal P represents the pore pressure value calculated by this method, in MPa. p,obs δ represents the measured pore pressure value of the target formation, in MPa; δ is the relative error, in %.

[0065] Preferably, step S54 further includes:

[0066] Calculate the predicted pressure P for the kth time. k The corresponding theoretical longitudinal wave velocity V p,calc (P k And calculate the residual of the kth iteration:

[0067] r k = V p,calc (P k )-V p,meas ;

[0068] In the formula, r k V represents the wave velocity residual corresponding to the k-th predicted pressure, in km / s; p,calc (P k P is the preset value for pore pressure. kThe theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity of the target stratum is denoted as km / s.

[0069] Compared with the prior art, the present invention has the following beneficial effects:

[0070] 1. This invention uses the difference between the theoretical P-wave velocity and the measured P-wave velocity as the matching target, and employs a trial-and-error algorithm to iterate in an interpolation interval manner. This allows the subsequent predicted pressure to gradually approach the target pressure value within this interval, reducing unnecessary and invalid calculations and helping to improve the stability and speed of the calculation process as well as the physical rationality of the results.

[0071] 2. Based on the relevant parameters of the target formation depth obtained from previous geophysical exploration, the prediction method provided by this invention can effectively predict the pressure of oil and gas reservoirs in deep and complex environments. It is supported by complete relevant physical data and can effectively solve the problem of missed or false high pressure in target formations under complex environments, thereby improving the production safety factor. Attached Figure Description

[0072] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0073] Figure 1 This is a schematic diagram of the workflow of the present invention;

[0074] Figure 2 This is a schematic diagram illustrating the verification of test results for a carbonate rock section in a certain well section of Dawan as the calculation object in this invention;

[0075] Figure 3 This is a schematic diagram illustrating the verification of test results for a section of Well-1, a marine carbonate rock well in southern Sichuan, as the calculation object in this invention.

[0076] Figure 4 This is a schematic diagram illustrating the verification of test results for a well section in the deep and complex strata of the Tarim Basin, which is the object of calculation in this invention. Detailed Implementation

[0077] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0078] This invention provides a method for predicting formation pore pressure in fractured reservoirs. This method is mainly based on the relationship between fracture mechanics and rock physics to predict the formation pore pressure in fractured reservoirs.

[0079] As attached Figure 1 As shown, the present invention provides a method for predicting formation pore pressure in fractured reservoirs, comprising the following steps:

[0080] S1: Obtain well logging and geological parameters at the target formation depth, and calculate the hydrostatic pressure and fracture surface normal stress σ at this location. n And combined with the preset pore pressure value P set for the target formation p Calculate the effective stress σ of the normal phase acting on the crack surface. eff .

[0081] Specifically, well logging and geological parameters at the target formation depth, including but not limited to mineral volume fraction f. i Total porosity φ total Formation fracture porosity φ frac pore fluid bulk modulus K fl Density ρ fl The crack direction, the crack surface dip angle θ, and the angle γ between the crack direction and the direction of the maximum horizontal principal stress.

[0082] Specifically, well logging and geological parameters also include the measured P-wave velocity V at the target formation depth. p,meas Three-dimensional geostress parameters.

[0083] Specifically, the three-dimensional geostress parameters include at least the vertical stress σ of the overlying strata. v Maximum horizontal principal stress σ H Minimum horizontal principal stress σ h .

[0084] Specifically, geostress refers to the natural stress existing in the Earth's crust that has not been disturbed by engineering projects. It is also called initial stress of rock mass, absolute stress, or original rock stress. In a broader sense, it also refers to the stress within the Earth's body, including stress generated by geothermal heat, gravity, changes in the Earth's rotation speed, and other factors.

[0085] Specifically, the formula for calculating hydrostatic pressure is as follows:

[0086] ;

[0087] In the formula, P h (z) represents the hydrostatic pressure at the target depth, in MPa; ρ w The average formation water density is taken as 1.05 g / cm³; g is the acceleration due to gravity, taken as 9.81 m / s². 2 z represents the target formation depth, in meters (m).

[0088] Specifically, the normal stress σ on the crack surface n The calculation uses the tensor projection method, and its formula is:

[0089] ;

[0090] In the formula, σ n The normal stress along the crack surface is σ, MPa; H The maximum horizontal principal stress is σ, MPa; h The minimum horizontal principal stress is given in MPa; σ v θ is the vertical stress of the overlying strata, MPa; θ is the dip angle of the fracture surface, i.e., the angle between the fracture surface and the horizontal plane, ° (degrees); γ is the angle between the fracture direction and the direction of the maximum horizontal principal stress, ° (degrees).

[0091] Specifically, the effective stress σ of the crack surface normal phase eff The calculation formula is:

[0092] ;

[0093] Where: σ eff The effective stress of the crack surface is σ (MPa); n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

[0094] S2: Young's modulus E of the combined pure mineral matrix s Compared with Poisson's ratio ν s Construct the relationship between the crack aspect ratio α and the theoretical crack closure pressure P. c The functional relationship.

[0095] Specifically, the crack aspect ratio α and the theoretical crack closure pressure P c The functional relationship is as follows:

[0096] ;

[0097] In the formula, P c The theoretical closure pressure of the fracture is MPa; Es is the Young's modulus of the pure mineral matrix, MPa, ν. s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless.

[0098] Specifically, the crack aspect ratio α is used to characterize the ratio between the crack's minor axis and major axis, reflecting the crack's flattening degree; this is achieved by constructing a relationship between the crack aspect ratio α and the theoretical crack closure pressure P. c The functional relationship can link the crack geometry with the crack closure pressure, and further establish the preset value P of pore pressure through the effective stress of the crack surface normal phase. p The relationship between the crack aspect ratio α and the crack length is used for subsequent calculations.

[0099] S3: The effective stress σ of the crack surface normal phase effWith the theoretical closure pressure P of the crack c Establish an equilibrium relationship, namely the effective stress σ of the crack surface normal phase. eff Equal to the theoretical closure pressure P of the crack c And construct a preset value P for pore pressure. p The relationship between the crack aspect ratio α and the crack.

[0100] Specifically, based on the established balance relationship, that is The effective stress σ of the crack surface method in step S1 is obtained by using the above steps. eff The preset pore pressure P can be obtained from the calculation formula and the functional relationship of step S2. p The relationship between the crack aspect ratio α and the crack aspect ratio α is as follows:

[0101] ;

[0102] In the formula, E s Young's modulus of pure mineral matrix, MPa, ν s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless; σ n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

[0103] S4: Based on the equivalent medium theoretical model and fluid substitution of the Gassmann equation, the calculation is based on the current preset pore pressure P. p The theoretical P-wave velocity V of the underlying strata p,calc .

[0104] Specifically, the equivalent medium theory model is used to equate heterogeneous rocks containing mineral matrix, pores, and fractures to a macroscopically homogeneous elastic medium, and to calculate the bulk modulus Kdry(P) of the dry rock skeleton. p ) and dry rock skeleton shear modulus μdry (P p ).

[0105] Specifically, equivalent medium theory models include, but are not limited to, DEM, KTDEM, MT, and self-consistent models. When using them, you can choose one of them.

[0106] Specifically, the Gassmann equation is used to calculate the bulk modulus of saturated rock under fluid saturation conditions based on the dry rock skeleton modulus, mineral matrix modulus, pore fluid bulk modulus, and total porosity; in the Gassmann fluid substitution relation, the saturated rock shear modulus is consistent with the dry rock skeleton shear modulus.

[0107] Specifically, based on the current preset pore pressure value P p The theoretical P-wave velocity V of the underlying strata p,calc The calculation process is as follows:

[0108] S41: Based on the preset pore pressure value P p And step S3 obtains the crack aspect ratio α(P) corresponding to the pore pressure. p The mineral composition, total porosity, fracture porosity, and fracture aspect ratio are input into the equivalent medium model.

[0109] Taking the self-consistent model calculation as an example, the calculation is based on the currently set preset value of pore pressure P. p Bulk modulus K of dry rock skeleton under certain conditions dry (P p ), GPa and dry rock skeleton shear modulus μ dry (P p ), GPa.

[0110] In one implementation, a self-consistent model can be used for calculation, where the rock contains n minerals, and the volume fraction of the i-th mineral in the mineral solid phase is f. i ,and The total number of phases input to the self-consistent model is n+2, where the 1st to nth phases are mineral phases, the (n+1)th phase is a rigid porous phase, and the (n+2)th phase is a fracture porous phase. The volume fraction of each phase is:

[0111] ;

[0112] ;

[0113] ;

[0114] In the formula, φ total Total porosity; φ frac denoted as fracture porosity; j represents a specific phase among the mineral phase, rigid porous phase, and fracture porous phase.

[0115] From this we can obtain .

[0116] It should be noted that this equivalent medium model calculation includes rigid porous phase and fracture porous phase, so the proportion of each phase needs to be recalculated, i.e., the total proportion of each mineral is 1-φ. total The proportion of rigid pores is φ total -φ frac The proportion of cracks and pores is φ frac .

[0117] Therefore, x j It is the whole rock composition phase volume fraction calculated from mineral volume fraction, total porosity, and fracture porosity.

[0118] ;

[0119] ;

[0120] In the formula, K dry (P p P is the preset value of pore pressure. p Bulk modulus of dry rock skeleton under the given conditions, GPa; μ dry (P p P is the preset value of pore pressure. p The shear modulus of the dry rock skeleton under the given conditions, in GPa; n is the number of mineral types; j is the phase number of the input self-consistent equivalent medium model; n+2 represents the total number of phases in the input self-consistent equivalent medium model, where the 1st to nth phases are mineral phases, the (n+1)th phase is a rigid porous phase, and the (n+2)th phase is a fracture porous phase; x j K represents the volume fraction of the j-th constituent phase in the whole rock, dimensionless; j c and μ j c The bulk modulus and shear modulus of the j-th constituent phase are respectively, in GPa; P j and Q j These are the bulk modulus polarization factor and shear modulus polarization factor corresponding to the j-th component, respectively.

[0121] S42: Mineral matrix bulk modulus K s and mineral matrix shear modulus μ s The elastic modulus of each mineral component can be obtained by averaging the Voigt-Reuss-Hill values:

[0122] ;

[0123] ;

[0124] In the formula, f i K represents the volume fraction of the i-th mineral, dimensionless; i and μ i Let be the bulk modulus and shear modulus of the i-th mineral, respectively, in GPa.

[0125] S43: Fluid substitution is performed using the Gassmann equation to obtain the preset pore pressure value P based on the current settings. p Saturated rock bulk modulus K under certain conditions sat (P p ) and saturated rock shear modulus μ sat (P p The calculation formula is:

[0126] ;

[0127] ;

[0128] In the formula, K sat (P p P is the preset value of pore pressure. p Saturated rock bulk modulus under certain conditions, GPa; K dry (P p P is the preset value of pore pressure. p Bulk modulus of dry rock skeleton under the given conditions, GPa; K s The equivalent bulk modulus of the mineral matrix is ​​expressed in GPa and K. fl φ is the bulk modulus of pore fluid, in GPa; total Total porosity, dimensionless; μ sat (P p P is the preset value of pore pressure. p Saturated rock shear modulus under certain conditions, GPa; μ dry (P p P is the preset value of pore pressure. p Shear modulus of dry rock skeleton under the given conditions, in GPa.

[0129] S44: Calculate the density of saturated rock using the following formula:

[0130] ;

[0131] In the formula, ρ sat φ is the density of saturated rock, g / cm³. total ρ represents total porosity, dimensionless. s The equivalent density of the mineral matrix is ​​g / cm³; ρ fl f is the pore fluid density, g / cm³. i ρ is the volume fraction of the i-th mineral, dimensionless; i Let be the density of the i-th mineral, in g / cm³.

[0132] S45: Calculate the preset pore pressure P based on the saturated rock bulk modulus, shear modulus, and density. p Theoretical P-wave velocity V under certain conditions p,calc The calculation formula is:

[0133] ;

[0134] In the formula, V p,calc (P p P is the preset value of pore pressure. p Theoretical P-wave velocity under the given conditions, km / s; K sat (P p P is the preset value of pore pressure. p Saturated rock bulk modulus under certain conditions, GPa; μsat (P p ) represents the current preset pore pressure value P. p Saturated rock shear modulus under certain conditions, GPa; ρ sat ρ is the density of saturated rock, in g / cm³.

[0135] S5: Based on the theoretical longitudinal wave velocity V p,calc Compared with the measured longitudinal wave velocity V p,meas The difference between them is the matching target. Within the pressure range defined by hydrostatic pressure and normal stress on the crack surface, a trial-and-error method is used to match the preset pore pressure P. p The solution is performed iteratively, and the final output is the calculated pore pressure value P of the target formation. p,cal .

[0136] Specifically, the trial method is an iterative solution method based on interval interpolation. This process does not only take the first calculation result, but approximates it multiple times until the difference (residual) between the calculated theoretical P-wave velocity and the measured P-wave velocity is reduced and meets the preset tolerance range.

[0137] Specifically, the iterative update process for cyclically solving pore pressure adopts an iterative solution method, and the specific steps are as follows:

[0138] S51: Determine adaptive tolerance: Extract all measured P-wave velocities V from well logging. p,meas Set tolerance ε v= η⋅V p,meas , where η is a preset proportional coefficient, and here η is taken as 0.01%.

[0139] S52: Determining the physical search boundary: Setting the lower bound P of the pore pressure search range. min For hydrostatic pressure, the upper limit is P. max Normal stress σ on the crack surface n .

[0140] S53: Calculation of pore pressure at initial depth point: For the initial depth point, at the lower boundary P of the pore pressure calculation interval... min and the upper boundary P max Theoretical P-wave velocity V is calculated at this location. p,calc The theoretical longitudinal wave velocity V corresponding to the lower boundary pore pressure is obtained. p,calc (P min The theoretical longitudinal wave velocity V corresponding to the upper pore pressure. p,calc (P max ); using P min P max With V p,calc (P min V p,calc (P maxThe first interpolation of the corresponding pressure-wave velocity relationship yields the first prediction point P1:

[0141] ;

[0142] In the formula, P1 is the predicted pore pressure obtained from the first interpolation, in MPa; P min The lower bound of the pore pressure search range is given in MPa; P max The upper bound of the pore pressure search range is given in MPa; V p,meas For the measured P-wave velocity, km / s; V p,calc (P min P is the preset value for pore pressure. min The theoretical P-wave velocity calculated at time, km / s; V p,calc (P max P is the preset value for pore pressure. max The theoretical P-wave velocity, calculated at that time, is in km / s.

[0143] Calculate the theoretical longitudinal wave velocity V corresponding to the first predicted pressure P1. p,calc (P1), and calculate the first residual: r1=V p,calc (P1)-V p,meas If |r1|≤ε v Then, this method takes P1 as the pore pressure result at the current depth. If |r1|>ε v Then, this method selects the interval containing the target solution based on the residual sign for a second interpolation; when r1>0, it indicates that the theoretical P-wave velocity corresponding to the first predicted pressure P1 is greater than the measured P-wave velocity, and the target pressure is located in the direction of larger pore pressure. At this time, P1 is taken as the low-pressure side endpoint, and P... max As the high-pressure side endpoint, a new interpolation interval is formed; when r1<0, it indicates that the theoretical P-wave velocity corresponding to the first predicted pressure P1 is less than the measured P-wave velocity, and the target pressure is located in the direction of smaller pore pressure. At this time, P... min Using P1 as the low-pressure side endpoint and P1 as the high-pressure side endpoint, a new interpolation interval is formed.

[0144] S54: Initial Depth Point Pore Pressure Subsequent Interpolation Update: Let P be the low-pressure endpoint of the current interpolation interval at the initial depth. L The corresponding theoretical longitudinal wave velocity is V p,calc (P L The high-pressure side endpoint is P. R The corresponding theoretical longitudinal wave velocity is V p,calc (P R Based on the pressure at both ends of the current interpolation interval and its theoretical P-wave velocity, calculate the k-th predicted pressure P. k :

[0145] ;

[0146] In the formula, P k The predicted pore pressure obtained from the k-th interpolation is expressed in MPa and V. p,meas For the measured P-wave velocity, km / s; P L The left endpoint of the current interpolation interval is the low-pressure side endpoint, in MPa; P R The right endpoint of the current interpolation interval is the high-pressure side endpoint, in MPa; V p,calc (P L P is the preset value for pore pressure. L The theoretical P-wave velocity calculated at time; V p,calc (P R P is the preset value for pore pressure. R The theoretical longitudinal wave velocity is calculated at that time.

[0147] Specifically, calculate the k-th predicted pressure P. k The corresponding theoretical longitudinal wave velocity V p,calc (P k And calculate the residual of the kth iteration:

[0148] r k = V p,calc (P k )-V p,meas ;

[0149] In the formula, r k V represents the wave velocity residual corresponding to the k-th predicted pressure, in km / s; p,calc (P k P is the preset value for pore pressure. k The theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity of the target stratum is denoted as km / s.

[0150] If |r k |≤ε v Then this method will P k As a result of the current depth pore pressure.

[0151] If |r k ∣>ε v Then according to r k Positive and negative update interpolation interval: when r k When >0, it means P k The corresponding theoretical longitudinal wave velocity V p,calc (P k () greater than the measured longitudinal wave velocity V p,meas The target pressure is located in the direction of greater pore pressure, which will cause P k Updated to the low-pressure side endpoint, i.e., P L =Pk V p,calc (P L )= V p,calc (P k When r k When <0, it means P k The corresponding theoretical longitudinal wave velocity V p,calc (P k (less than the measured longitudinal wave velocity V) p,meas The target pressure is located in the direction of smaller pore pressure, which will cause P k Updated to the high-pressure side endpoint, i.e., P R =P k V p,calc (P R )= V p,calc (P k ).

[0152] After the update, continue in the new interpolation interval [P] L ,P R Interpolation calculations are performed within the range until the residual meets the tolerance requirements or reaches the preset maximum number of iterations.

[0153] It should be noted that both steps S53 and S54 employ interpolation, but step S53 is used to establish the initial interpolation interval for the initial depth point, i.e., with P... min and P max Using the endpoint as an example, calculate the corresponding theoretical P-wave velocity V. p,calc (P min ) and V p,calc (P max Step S54 is used to interpolate within the current interpolation interval [P] L ,P R Subsequent iterative updates will be performed within this period, based on the current endpoint pressure and the corresponding theoretical longitudinal wave velocity V. p,calc (P L V p,calc (P k Calculate the predicted pressure P for the kth time. k The interval endpoints are then updated based on the residual sign. Therefore, S53 is the initial interval construction step, and S54 is the subsequent interval update step; these two steps are not redundant.

[0154] S55: Pore pressure calculation at subsequent depth points: For subsequent depth points, the pore pressure result P obtained from the previous depth point is used... prev As a reference pressure at the current depth point.

[0155] First, calculate the reference pressure P under the formation parameter conditions at the current depth point. prev The corresponding theoretical longitudinal wave velocity is denoted as V. p,calc (P prevAnd calculate the reference residual:

[0156] r prev =V p,calc (P prev )− V p,meas ;

[0157] In the formula, r prev P represents the wave velocity residual corresponding to the reference pressure, in km / s. prev The pore pressure result calculated at the previous depth point is in MPa; V p,calc (P prev P is the preset value of pore pressure under the current depth conditions. prev The theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity at the current depth point is in km / s.

[0158] If |r prev |≤ε v Then this method will P prev As a result of the current depth pore pressure.

[0159] If |r prev ∣>ε v Then, based on the reference residual r prev The positive or negative judgment of the target pressure relative to the reference pressure P prev The direction of r, and only calculate the theoretical P-wave velocity corresponding to one boundary pressure in the corresponding direction; when r prev When the value is greater than 0, it indicates that the reference pressure P is greater than 0. prev The corresponding theoretical longitudinal wave velocity V p,calc (P prev () is greater than the measured longitudinal wave velocity V at the current depth point p,meas When the target pressure is located in the direction of greater pore pressure, only the upper limit of the pressure P is calculated. max The corresponding theoretical longitudinal wave velocity V p,calc (P max ), and P prev As the low-pressure side endpoint, P max As the high-pressure side endpoint, it forms the interpolation interval for the current subsequent points; when r prev When <0, it indicates that the reference pressure P prev The corresponding theoretical longitudinal wave velocity V p,calc (P prev (less than the measured P-wave velocity V at the current depth point) p,meas When the target pressure is located in the direction of the smaller pore pressure, only the lower limit of the pressure P is calculated. min The corresponding theoretical longitudinal wave velocity V p,calc (P min ), and P min As the low-pressure side endpoint, Pprev As the high-pressure side endpoint, it forms the interpolation interval for the current subsequent points.

[0160] S56: Repeat steps S54 and S55 to update the pore pressure at subsequent depth points via interpolation.

[0161] Specifically, steps S54 and S55 are calculated repeatedly. In each iteration, the calculated value of step S55 from the previous calculation is used as the starting interpolation value for step S54 in the next iteration, until the reference pressure at the current depth point is at P. prev Wave velocity residual r under the condition prev The absolute value is less than the tolerance ε v The iterative calculation ends when the time is right, and the next step begins.

[0162] S57: Output Results: The theoretical P-wave velocity V corresponding to any predicted pressure during the calculation of the initial or subsequent depth points. p,calc Compared with the measured longitudinal wave velocity V p,meas When the residual between the two conditions meets the tolerance requirement, the predicted pressure is output as the pore pressure result at the corresponding depth point, i.e., the calculated pore pressure value P of the target formation. p,cal .

[0163] S6: Verification, using the calculated pore pressure value P of the target formation. p,cal The measured pressure value P of the target formation p,obs Compare and establish evaluation relationships.

[0164] In step S6, the evaluation relationship is established using the relative error rate, which is calculated as follows:

[0165] ;

[0166] ;

[0167] In the formula, ΔP is the difference between the predicted pressure and the measured pressure, in MPa; P p,cal P is the pore pressure calculated by this method, in MPa; p,obs δ represents the measured pressure value of the target formation, in MPa; δ is the relative error, in %.

[0168] It should be noted that step S6 mainly participates in verification calculations in the early stages of the work. That is, when the formation pore pressure prediction method of the fractured reservoir of the present invention is used in a new area, it participates in the verification calculation of the predicted pore pressure to verify the accuracy of the formation pore pressure prediction method of the fractured reservoir provided by the present invention in this area. After the verification is successful, the prediction of pore pressure at other target formation depths in the current area can be performed using steps S1 to S5 of the present invention.

[0169] This invention provides a method for predicting formation pore pressure in fractured reservoirs. Based on the relationship between fracture mechanics and rock physics, this method establishes the correspondence between pore pressure, fracture aspect ratio, and theoretical P-wave velocity. The method uses the residual between theoretical P-wave velocity and measured P-wave velocity as the matching target and solves the predicted pore pressure at the target depth through interval interpolation iteration.

[0170] Through the above steps, the predicted pore pressure at the corresponding depth point can be obtained iteratively within the defined pore pressure search range for both the initial and subsequent depth points, constrained by the residual between the theoretical and measured P-wave velocities.

[0171] The solution method employed in this invention involves setting a preset pore pressure value to calculate the theoretical P-wave velocity, then iterating continuously until the residual between the theoretical P-wave velocity and the measured wave velocity meets the tolerance requirement, at which point the result is output. During the calculation process, the preset pore pressure value is a continuously updated pore pressure value, and the predicted pore pressure is the final output value predicted using this method, i.e., the calculated pore pressure value of the target formation.

[0172] Among them, the preset value of pore pressure P p It is the trial pressure used for theoretical P-wave velocity during the iterative solution process; when the residual between the theoretical P-wave velocity and the measured P-wave velocity corresponding to a certain candidate pore pressure meets the tolerance requirement, the candidate pore pressure is used as the predicted pore pressure output.

[0173] Specific implementation examples of the present invention:

[0174] To illustrate the applicable process of the method of the present invention and to verify the accuracy of the calculation results of the prediction method provided by the present invention, as shown in the appendix... Figure 2-4 As shown, three different depth segments were selected as examples. The formation pore pressure prediction method for fractured reservoirs described above was used to predict and calculate the pore pressure of the target formation. The pore pressure of the target formation was actually measured, and the calculated pore pressure value was compared with the actual measured pore pressure value to verify the accuracy of the formation pore pressure prediction method for fractured reservoirs provided by the present invention.

[0175] During verification, it is also necessary to actually measure the pressure of the target formation to obtain the measured pore pressure value P of the target formation. p,obs .

[0176] All three embodiments are performed according to steps S1 to S6 provided by the formation pore pressure prediction method for fractured reservoirs of the present invention. The process of steps S1 to S5 is briefly described as follows:

[0177] First, obtain the measured P-wave velocity V at the target depth point.p,meas Normal stress σ on the crack surface n Total porosity φ total φ of crack porosity frac The mineral component volume fraction and related fluid parameters are then used to establish a preset pore pressure value P based on the effective stress relationship between the normal stress on the fracture surface and the pore pressure. p The relationship between the longitudinal and transverse aspect ratio α of the crack was established; then, combining the equivalent medium model and the Gassmann equation, the theoretical longitudinal wave velocity V under different pore pressure conditions was calculated. p,calc Finally, with V p,calc With V p,meas Using the difference between them as the basis for judgment, the pore pressure is gradually updated within the preset pressure range to obtain the calculated pore pressure value P of the target formation at each depth point. p,cal .

[0178] To evaluate the prediction results, step S6 of the formation pore pressure prediction method for fractured reservoirs provided by this invention is used to compare the calculated pore pressure value of the target formation with the measured pore pressure value of the target formation.

[0179] The error calculation method is as follows:

[0180] ;

[0181] ;

[0182] In the formula, ΔP is the difference between the predicted pore pressure (i.e., the calculated pore pressure value) and the measured pore pressure value, in MPa; P p,cal P represents the pore pressure value calculated by this method, in MPa. p,obs δ represents the measured pore pressure value of the target formation, in MPa; δ is the relative error, in %. Specific Implementation Example 1:

[0184] As attached Figure 2 As shown in the figure, this specific embodiment 1 selects a carbonate rock section in a certain well section of Dawan as the calculation object.

[0185] The input parameters for this well section include the measured P-wave velocity V. p,meas Normal stress σ on the crack surface n Total porosity φ total φ frac Volume fractions of evaporite, dolomite, limestone, and mudstone, etc.

[0186] In this embodiment, measured pressure data at a depth of 4850m is available for comparative analysis.

[0187] The comparative analysis results are shown in Table 1.

[0188] Table 1. Comparison and Analysis of Measured Pressure Data and Predicted Pore Pressure at a Depth of 4850m

[0189]

[0190] As shown in Table 1, at 4850m, the pore pressure P calculated using the aforementioned method for predicting formation pore pressure in fractured reservoirs is... p,cal The measured pressure value P in the target formation is 60.39 MPa. p,obs The value is 60.20 MPa, the difference between the two is 0.19 MPa, and the relative error is approximately 0.32%.

[0191] The results show that, in this well section, the calculation results using the method provided by this invention are close to the measured pressure, and can reflect the pore pressure level near this depth point. Specific Implementation Example 2:

[0193] This embodiment selects a section of the Well-1 well in the marine carbonate rock of southern Sichuan as the calculation object. Input parameters include the measured P-wave velocity V. p,meas Normal stress σ on the crack surface n Total porosity φ total φ of crack porosity frac And mineral composition parameters.

[0194] In this embodiment, there are 3 depth points with measured pressure data. The predicted results are compared with the measured pressure as follows.

[0195] Table 2. Comparison and analysis of measured pressure data and predicted pore pressure at three depth points.

[0196]

[0197] The average absolute error of the three measured comparison points in this well section is approximately 0.32 MPa, the root mean square error is approximately 0.37 MPa, and the average relative error is approximately 0.54%. The largest absolute error occurs at 3017 m, at approximately 0.57 MPa, corresponding to a relative error of approximately 0.97%.

[0198] The comparison results show that the predicted pore pressure and the measured pressure have the same trend, and the error is within a small range. Specific Implementation Example 3:

[0200] This embodiment selects a well section in the deep, complex strata of the Tarim Basin as the calculation object. Input parameters include the measured P-wave velocity V. p,meas Normal stress σ on the crack surface n Total porosity φ total φ of crack porosityfrac And mineral composition parameters.

[0201] In this embodiment, there are 5 depth points with measured pressure data. The predicted results are compared with the measured pressure as follows.

[0202] Table 3. Comparison and analysis of measured pressure data and predicted pore pressure at 5 depth points.

[0203]

[0204] The average absolute error of the five measured comparison points in this well section is approximately 0.24 MPa, the root mean square error is approximately 0.32 MPa, and the average relative error is approximately 0.27%. The maximum absolute error is 0.57 MPa, corresponding to a relative error of approximately 0.63%.

[0205] The results show that the predicted pore pressure increases with depth, which is basically consistent with the pressure change trend reflected by the measured pressure, and the predicted results are similar to those of the measured pressure.

[0206] Comprehensive analysis:

[0207] As can be seen from the three specific embodiments, this method uses the measured longitudinal wave velocity as a constraint and performs point-by-point calculations based on the relationship between the preset value of pore pressure, the longitudinal-to-transverse ratio of the crack, and the theoretical longitudinal wave velocity.

[0208] The calculation results show that the average absolute error between the pore pressure calculated by this method and the measured pressure of the target formation is about 0.25 MPa, and the average relative error is about 0.36%. The average relative error is calculated by summing the average absolute error values ​​obtained from the three specific implementations mentioned above and then calculating the average absolute error.

[0209] In summary, these results demonstrate that the proposed method can effectively reflect the variation characteristics of pore pressure in fractured reservoirs and can be used for continuous pressure prediction at different depths. This allows for the calculation of the pore pressure value P of the target formation within the target area using steps S1 to S5 of the prediction method provided by this invention during subsequent work. p,cal The pore pressure value P of the target formation was calculated based on this. p,cal This provides theoretical guidance for subsequent production.

[0210] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A method for predicting formation pore pressure in fractured reservoirs, characterized in that, Includes the following steps: S1: Obtain well logging and geological parameters at the target formation depth, and calculate the hydrostatic pressure and fracture surface normal stress σ at this location. n And combined with the preset pore pressure value P set for the target formation p Calculate the effective stress σ of the normal phase acting on the crack surface. eff ; S2: Young's modulus E of the combined pure mineral matrix s Compared with Poisson's ratio ν s Construct the relationship between the crack aspect ratio α and the theoretical crack closure pressure P. c The functional relationship; S3: The effective stress σ of the crack surface normal phase eff With the theoretical closure pressure P of the crack c Establish an equilibrium relationship, namely the effective stress σ of the crack surface normal phase. eff Equal to the theoretical closure pressure P of the crack c And construct a preset value P for pore pressure. p The relationship between the crack aspect ratio α and the crack width ratio α; S4: Based on the equivalent medium theoretical model and fluid substitution of the Gassmann equation, the calculation is based on the current preset pore pressure P. p The theoretical P-wave velocity V of the underlying strata p,calc ; S5: Based on the theoretical longitudinal wave velocity V p,calc Compared with the measured longitudinal wave velocity V p,meas The difference between them is the matching target. Within the pressure range defined by hydrostatic pressure and normal stress on the crack surface, a trial-and-error method is used to match the preset pore pressure P. p The solution is performed iteratively, and the final output is the calculated pore pressure value P of the target formation. p,cal ; S6: Verification, using the calculated pore pressure value P of the target formation. p,cal The measured pore pressure value P of the target formation p,obs Compare and establish evaluation relationships.

2. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, The well logging and geological parameters for the target formation depth in step S1 include, but are not limited to, the mineral volume fraction f. i Total porosity φ total Formation fracture porosity φ frac pore fluid bulk modulus K fl Density ρ fl , crack direction, crack surface dip angle θ, and the angle γ between the crack direction and the direction of the maximum horizontal principal stress; It also includes the measured P-wave velocity V at the target formation depth. p,meas Three-dimensional geostress parameters; The three-dimensional geostress parameters include at least the vertical stress σ of the overlying strata. v Maximum horizontal principal stress σ H Minimum horizontal principal stress σ h .

3. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S1, the normal stress σ on the crack surface is... n The calculation uses the tensor projection method, and its formula is: ； In the formula, σ n The normal stress along the crack surface is σ, MPa; H The maximum horizontal principal stress is σ, MPa; h The minimum horizontal principal stress is given in MPa; σ v θ is the vertical stress of the overlying strata, MPa; θ is the dip angle of the fracture surface, i.e., the angle between the fracture surface and the horizontal plane, ° (degrees); γ is the angle between the fracture direction and the direction of the maximum horizontal principal stress, ° (degrees).

4. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S1, the effective stress σ of the crack surface normal phase is... eff The calculation formula is: ； In the formula, σ eff The effective stress of the crack surface is σ (MPa); n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

5. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S2, the crack aspect ratio α and the theoretical crack closure pressure P c The functional relationship is as follows: ； In the formula, P c E represents the theoretical closure pressure of the crack, in MPa. s Young's modulus of pure mineral matrix, MPa, ν s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless.

6. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S3, the preset pore pressure value P p The relationship between the crack aspect ratio α and the crack aspect ratio α is as follows: ； In the formula, E s Young's modulus of pure mineral matrix, MPa, ν s α is the Poisson's ratio of the pure mineral matrix, dimensionless; α is the aspect ratio of the fracture, dimensionless; σ n P represents the normal stress along the crack surface, in MPa. p This is the preset value for pore pressure, in MPa.

7. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S4, the current pore pressure preset value P p The theoretical P-wave velocity V of the underlying strata p,calc The calculation process is as follows: S41: Based on the preset pore pressure value P p And step S3 obtains the crack aspect ratio α(P) corresponding to the pore pressure. p The mineral composition, total porosity, fracture porosity, and fracture aspect ratio are input into the equivalent medium model to calculate the pore pressure preset value P based on the currently set value. p Bulk modulus K of dry rock skeleton under certain conditions dry (P p ), GPa and dry rock skeleton shear modulus μ dry (P p ), GPa; S42: Mineral matrix bulk modulus K s and mineral matrix shear modulus μ s The elastic modulus of each mineral component can be obtained by averaging the Voigt-Reuss-Hill values: ； ； In the formula, f i K represents the volume fraction of the i-th mineral, dimensionless; i and μ i Let be the bulk modulus and shear modulus of the i-th mineral, respectively, in GPa; S43: Fluid substitution is performed using the Gassmann equation to obtain the preset pore pressure value P based on the current settings. p Saturated rock bulk modulus K under certain conditions sat (P p ) and saturated rock shear modulus μ sat (P p The calculation formula is: ； ； In the formula, K sat (P p P is the preset value of pore pressure. p Saturated rock bulk modulus under certain conditions, GPa; K dry (P p P is the preset value of pore pressure. p Bulk modulus of dry rock skeleton under the given conditions, GPa; K s The equivalent bulk modulus of the mineral matrix is ​​expressed in GPa and K. fl φ is the bulk modulus of pore fluid, in GPa; total Total porosity, dimensionless; μ sat (P p P is the preset value of pore pressure. p Saturated rock shear modulus under certain conditions, GPa; μ dry (P p P is the preset value of pore pressure. p Shear modulus of dry rock skeleton under the given conditions, GPa; S44: Calculate the density of saturated rock using the following formula: ； In the formula, ρ sat φ is the density of saturated rock, g / cm³. total ρ represents total porosity, dimensionless. s The equivalent density of the mineral matrix is ​​g / cm³; ρ fl f is the pore fluid density, g / cm³. i ρ represents the volume fraction of the i-th mineral, dimensionless; i Let be the density of the i-th mineral, in g / cm³; S45: Calculate the preset pore pressure P based on the saturated rock bulk modulus, shear modulus, and density. p Theoretical P-wave velocity V under certain conditions p,calc The calculation formula is: ； In the formula, V p,calc (P p P is the preset value of pore pressure. p Theoretical P-wave velocity under the given conditions, km / s; K sat (P p P is the preset value of pore pressure. p Saturated rock bulk modulus under certain conditions, GPa; μ sat (P p P is the preset value of pore pressure. p Saturated rock shear modulus under certain conditions, GPa; ρ sat ρ is the density of saturated rock, in g / cm³.

8. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S5, the trial algorithm is an iterative solution method based on interval interpolation, and the specific steps are as follows: S51: Determine adaptive tolerance: Extract all measured P-wave velocities V from well logging. p,meas Set tolerance ε v= η⋅V p,meas Where η is a preset proportionality coefficient, and here η is taken as 0.01%; S52: Determining the physical search boundary: Setting the lower bound P of the pore pressure search range. min For hydrostatic pressure, the upper limit is P. max Normal stress σ on the crack surface n ; S53: Calculation of pore pressure at initial depth point: For the initial depth point, at the lower boundary P of the pore pressure calculation interval... min and the upper boundary P max Theoretical P-wave velocity V is calculated at this location. p,calc The theoretical longitudinal wave velocity V corresponding to the lower boundary pore pressure is obtained. p,calc (P min The theoretical longitudinal wave velocity V corresponding to the upper pore pressure. p,calc (P max ); using P min P max With V p,calc (P min V p,calc (P max The first interpolation of the corresponding pressure-wave velocity relationship yields the first prediction point P1: ； In the formula, P1 is the predicted pore pressure obtained from the first interpolation, in MPa; P min The lower bound of the pore pressure search range is given in MPa; P max The upper bound of the pore pressure search range is MPa; V p,meas For the measured P-wave velocity, km / s; V p,calc (P min P is the preset value for pore pressure. min The theoretical P-wave velocity calculated at time, km / s; V p,calc (P max P is the preset value for pore pressure. max The theoretical P-wave velocity calculated at that time, km / s; S54: Pore pressure subsequent interpolation update: Let P be the low-pressure side endpoint of the current interpolation interval at the current depth. L The corresponding theoretical longitudinal wave velocity is V p,calc (P L The high-pressure side endpoint is P. R The corresponding theoretical longitudinal wave velocity is V p,calc (P R Based on the pressure at both ends of the current interpolation interval and its theoretical P-wave velocity, calculate the k-th predicted pressure P. k : ； In the formula, P k The predicted pore pressure obtained from the k-th interpolation is expressed in MPa and V. p,meas For the measured P-wave velocity, km / s; P L The left endpoint of the current interpolation interval is the low-pressure side endpoint, in MPa; P R The right endpoint of the current interpolation interval is the high-pressure side endpoint, in MPa; V p,calc (P L P is the preset value for pore pressure. L The calculated theoretical P-wave velocity; V p,calc (P R P is the preset value for pore pressure. R The calculated theoretical longitudinal wave velocity; S55: Pore pressure calculation at subsequent depth points: For subsequent depth points, the pore pressure result P obtained from the previous depth point is used... prev As a reference pressure at the current depth point; First, calculate the reference pressure P under the formation parameter conditions at the current depth point. prev The corresponding theoretical longitudinal wave velocity is denoted as V. p,calc (P prev And calculate the reference residual: r prev =V p,calc (P prev )−V p,meas ; In the formula, r prev P represents the wave velocity residual corresponding to the reference pressure, in km / s. prev The pore pressure result calculated at the previous depth point is in MPa; V p,calc (P prev P is the preset value of pore pressure under the current depth conditions. prev The theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity at the current depth point is in km / s; S56: Repeat steps S54 and S55 to perform subsequent depth point pore pressure interpolation updates. S57: Output Results: The theoretical P-wave velocity V corresponding to any predicted pressure during the calculation of the initial or subsequent depth points. p,calc Compared with the measured longitudinal wave velocity V p,meas When the residual between the two conditions meets the tolerance requirement, the predicted pressure is output as the pore pressure result at the corresponding depth point, i.e., the calculated pore pressure value P of the target formation. p,cal .

9. The method for predicting formation pore pressure in fractured reservoirs according to claim 1, characterized in that, In step S6, the evaluation relationship is established using the relative error rate, which is calculated as follows: ； ； In the formula, ΔP is the difference between the predicted pore pressure value (i.e., the calculated pore pressure value) and the measured pore pressure value, in MPa; P p,cal P represents the pore pressure value calculated by this method, in MPa. p,obs δ represents the measured pore pressure value of the target formation, in MPa; δ is the relative error, in %.

10. The method for predicting formation pore pressure in fractured reservoirs according to claim 8, characterized in that, Step S54 further includes: Calculate the predicted pressure P for the kth time. k The corresponding theoretical longitudinal wave velocity V p,calc (P k And calculate the residual of the kth iteration: r k = V p,calc (P k )-V p,meas ; In the formula, r k V represents the wave velocity residual corresponding to the k-th predicted pressure, in km / s; p,calc (P k P is the preset value for pore pressure. k The theoretical P-wave velocity calculated at time, km / s; V p,meas The measured P-wave velocity of the target stratum is denoted as km / s.

Images