Carbonate reservoir pore type earthquake prediction method and device based on porous effective medium model

Abstract

The invention discloses a carbonate rock porous effective medium model established on the basis of an expanded Keys-Xu model and a Gassmann-Hill equation, a multi-porosity inversion method of a two-step method is carried out, and in addition, in the inversion process, a strategy which is different from a previous fixed pore aspect ratio and uses a variable pore aspect ratio as input is adopted. And the multi-porosity inversion precision is improved. The method not only can accurately characterize the distribution of the three pore types of the compact carbonate reservoir, but also can successfully reveal the micromechanism of oil-gas migration and aggregation, and provides a theoretical basis and a technical support for characterizing the pore structure of the ultra-deep compact reservoir by using logging and seismic data. The multi-porosity inversion strategy of the porous medium model can also be applied to quantitative prediction of pore types of other lithologic tight reservoirs, and the problem that the multi-porosity of the tight reservoir is difficult to accurately predict is solved.

Description

Technical Field

[0001] The present invention relates to the technical field of geophysical exploration, and particularly relates to a seismic inversion method and device and an electronic device for pore types of carbonate rock reservoirs based on a porous effective medium model. Background Art

[0002] Carbonate rock reservoirs account for more than 50% of the global oil and gas reserves and have great exploration potential. Different from conventional clastic rock reservoirs, carbonate rock reservoirs have the characteristic of complex pore structures due to various diagenetic processes such as cementation, compaction, dissolution, and dolomitization. The pore system of carbonate rock reservoirs usually includes cavernous pores, intragranular pores, intergranular pores, intercrystalline pores, and fractures (Zhao et al., 2013), which significantly affect the permeability, acoustic properties, and seismic interpretation of carbonate rock reservoirs (Sun et al., 2006; Weger et al., 2009; Jin et al., 2023; Guo et al., 2023; Zhang et al., 2024). Therefore, accurately characterizing the pore structure is crucial for evaluating reservoir quality, describing reservoir structure, and identifying sweet spots within carbonate rock reservoirs.

[0003] Velocity deviation logging can be used to identify fractures and rigid pores (Anselmetti and Eberli, 1999); acoustic and density logging can also be used to derive matrix porosity and total porosity through empirical relationships (Du et al., 2024). However, difficulties are faced when simultaneously evaluating multiple pore types at well locations (Wu et al., 2016; Sharifi et al., 2018). Therefore, there is an urgent need for a method that can simultaneously quantify the spatial distribution of multiple pore types within carbonate rock reservoirs.

[0004] The pore aspect ratio is defined as the ratio of the short axis to the long axis of ellipsoidal pores and is commonly used to indicate pore types (Guo et al., 2015; Fournier et al., 2018; Teillet et al., 2021). Based on the concept of classifying pore types by pore aspect ratio (rigid pores with an aspect ratio of 0.7 - 1.0 correspond to holes or modal pores; reference pores or matrix pores with an aspect ratio of 0.1 - 0.25 represent intergranular or intercrystalline pores; soft pores with an aspect ratio of 0.01 - 0.05 represent fractures), some studies have explored petrophysics-based methods, especially the effective medium theory based on inclusions, to characterize the complex pore structures of carbonate reservoirs (Xu and Payne, 2009; Zhao et al., 2013; Li, 2020; Saberi et al., 2020; Guo et al., 2021; Guo et al., 2022). In contrast, pore-type inversion methods such as the iterative method based on the DEM model and Gassmann equation (Kumar and Han, 2005), the carbonate-specific adaptation of the Xu-White model (Xun and Payne,) etc. seem to be more suitable for quantifying complex pore systems. However, these methods are limited to characterizing two pore types (reference pores / rigid pores or reference pores / fractures) at a single sampling point and are only conditional on the P-wave velocity. Therefore, it is necessary to study pore-type seismic inversion methods based on the porous hypothesis.

[0005] To solve the problem of low accuracy in pore-type inversion caused by only considering two pore types, the present invention establishes a multi-porosity model for tight carbonate reservoirs based on the extended Keys-Xu model and the Gassmann-Hill equation, and characterizes the porous types and their distributions in heterogeneous carbonate reservoirs by jointly using the P-wave and S-wave velocities of well logging and seismic data. Summary of the Invention

[0006] The purpose of the present invention is to provide a pore-type seismic inversion method, device, and electronic device based on a porous effective medium model to solve the problems raised in the above background technology.

[0007] To achieve the above purpose, the present invention adopts the following technical solutions:

[0008] In the first aspect, a pore-type seismic inversion method based on a porous effective medium model is provided, including:

[0009] Step S1: Set the initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, and reservoir rock matrix properties; Step S2: Calculate the P-wave velocities V P,wyllie 、V P,HS+ 、V P,HS- 、VP,KX ;

[0010] Step S3: Calculate the differences between the P-wave velocities calculated by the Wyllie time-average equation, the HS upper and lower boundaries, and the P-wave velocity calculated by the Keys-Xu model for a single pore type. When the differences between the P-wave velocities calculated by the Wyllie time-average equation, the HS upper and lower boundaries, and the P-wave velocity calculated by the Keys-Xu model for a single pore type all meet the cut-off condition, terminate the iteration to obtain the final reservoir pore aspect ratio. The iterative calculation formula is as follows:

[0011] |V P,KX -V P,s | < ε, s = wyllie, HS+, HS-

[0012] In the formula, V P,KX represents the P-wave velocity calculated by the Keys-Xu model for a single pore type; V P,wyllie , VP,HS+ and V P,HS- represent the P-wave velocities calculated by the Wyllie time-average equation, the HS upper boundary, and the HS lower boundary respectively; ε represents the given minimum error;

[0013] Step S4: Based on the final reservoir pore aspect ratio obtained in Step S3, and combined with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, reservoir initial reference pores, and volume fractions of the reservoir initial fractures, input them into the constructed carbonate reservoir porous effective medium model to calculate the P-wave velocity V P,cal and the S-wave velocity V S,cal ;

[0014] Step S5: Use the objective function J to calculate the errors between the P-wave velocity V P,cal and the S-wave velocity V S,cal calculated in Step S4 and their well logging measured values V P,obs , V S,obs . The expression of the objective function J is as follows:

[0015] J = |V P,obs -V P,cal | 2 +|V S,obs -V S,cal | 2

[0016] In the formula, V P,obs and V P,cal represent the well logging measured P-wave velocity and the P-wave velocity V S,obs calculated based on the carbonate reservoir porous effective medium model respectively, and V S,cal represent the well logging measured S-wave velocity and the S-wave velocity calculated based on the carbonate reservoir porous effective medium model respectively;

[0017] Step S6: Traverse the volume fractions of the reference pores and fractures within a given range with a set step size, and terminate the iteration when the objective function J in Step S5 is less than a given error σ, so as to obtain the final volume fractions of the reference pores and fractures;

[0018] Step S7: Calculate the final volume fraction of the hard pores based on the final volume fractions of the reference pores and fractures obtained in Step S6;

[0019] Step S8: Based on the final volume fractions of the reference pores, fractures, and hard pores obtained in Step S7, and in combination with the reservoir porosity, calculate the porosities of the final reference pores, fractures, and hard pores to complete the inversion.

[0020] Further, the construction method of the carbonate reservoir porous effective medium model in Step S4 includes:

[0021] Step S41: Calculate the bulk modulus and shear modulus of the rock matrix using the Voigt-Reuss-Hill (VRH) model;

[0022] Step S42: Calculate the bulk modulus and shear modulus of the dry rock using the extended Keys-Xu model;

[0023] Step S43: Calculate the bulk modulus and shear modulus of the saturated rock using the Gassmann equation under the low-frequency assumption, and determine the bulk modulus and shear modulus of the overall saturated rock using the Gassmann-Hill equation in the carbonate reservoir.

[0024] Further, the method for calculating the bulk modulus and shear modulus of the rock matrix using the Voigt-Reuss-Hill (VRH) model in Step S41 includes:

[0025]

[0026] where f i , K i and G i represent the volume fraction, bulk modulus, and shear modulus of the i-th mineral component respectively; M represents the total number of types of mineral components contained in the rock matrix.

[0027] Further, the method for calculating the bulk modulus and shear modulus of the dry rock using the extended Keys-Xu model in Step S42 includes:

[0028] Calculate the bulk modulus and shear modulus of the dry rock using the extended Keys-Xu model applicable to hard pores, reference pores, and fractures. The formula is as follows:

[0029] K dry = K ma (1 - φ)P

[0030] G dry = G ma (1 - φ) Q

[0031] Where φ represents porosity; P and Q are geometric factors related to the pore aspect ratio:

[0032]

[0033] Where v s , v r and v c represent the volume fractions of hard pores, reference pores, and fractures, respectively; α s , α r and α c represent the pore aspect ratios of hard pores, reference pores, and fractures, respectively; T ijij (α l ) and T iijj (α l ) are functions of the pore aspect ratio.

[0034] Furthermore, the calculation method for the bulk modulus and shear modulus of the saturated rock using the Gassmann equation under the low-frequency assumption in step S43 includes:

[0035]

[0036] G sat = G dry

[0037] Where K fl represents the bulk modulus of the pore fluid.

[0038] Furthermore, the method for determining the bulk modulus and shear modulus of the overall saturated rock using the Gassmann-Hill equation in step S43 includes:

[0039]

[0040] Where the bulk modulus K sat,w of the water-saturated rock and the bulk modulus K sat,g of the gas-saturated rock are calculated from the Gassmann equation under the low-frequency assumption, and S w represents the water saturation.

[0041] Furthermore, the calculation of the P-wave velocity V P,wyllie using the Wyllie time-average equation in step S2 includes:

[0042]

[0043] where φ represents porosity; V P,fl and V P,ma respectively represent the P-wave velocities of the fluid and the rock matrix.

[0044] Furthermore, in step S2, the P-wave velocities V P,HS+ and V P，HS- are calculated using the HS upper and lower bounds. The method includes:

[0045] Step S21: Calculate the bulk modulus and shear modulus of carbonate rocks using the Hashin-Shtrikman (HS) upper and lower bounds:

[0046]

[0047] where K 1 , G 1 and f 1 respectively represent the bulk modulus, shear modulus and volume fraction of the first phase; K 2 , G 2 and f 2 represent the bulk modulus, shear modulus and volume fraction of the second phase. The calculation of the HS boundary upper and lower bounds is determined by swapping the order of the first phase and the second phase. "+" represents the upper bound and "-" represents the lower bound;

[0048] Step S22: Calculate the P-wave velocities V P，HS+ and V P，HS- of the carbonate reservoir based on the bulk modulus and shear modulus of the carbonate rocks calculated in step S21:

[0049]

[0050] Second, a seismic inversion device for carbonate reservoir pore types based on a porous effective medium model is provided, including: An initial setting unit: Set the initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, and reservoir rock matrix properties;

[0051] A first P-wave velocity calculation unit: Calculate the P-wave velocities V P，wtllie , V P，HS+ , V P，HS- , V P，KX using the Wyllie time-average equation, the HS upper and lower bounds, and the single pore type Keys-Xu model respectively;

[0052] Porosity aspect ratio calculation unit: The difference is calculated between the longitudinal wave velocities calculated by the Wyllie time-average equation, the HS upper and lower boundaries, and the longitudinal wave velocity calculated by the Keys-Xu model for a single pore type. When the differences between the longitudinal wave velocities calculated by the Wyllie time-average equation, the HS upper and lower boundaries, and the longitudinal wave velocity calculated by the Keys-Xu model for a single pore type all meet the cut-off condition, the iteration is terminated to obtain the final reservoir porosity aspect ratio. The iteration calculation formula is as follows:

[0053] |V P,KX -V P,s |<ε, s=wyllie, HS+, HS-

[0054] In the formula, V P,KX represents the longitudinal wave velocity calculated by the Keys-Xu model for a single pore type; V P,wyllie , V P,HS+ and V P,HS- respectively represent the longitudinal wave velocities calculated by the Wyllie time-average equation, the HS upper boundary, and the HS lower boundary; ε represents the given minimum error;

[0055] Second longitudinal wave velocity calculation unit: Based on the final reservoir porosity aspect ratio obtained by the porosity aspect ratio calculation unit, and combined with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, the volume fractions of the initial reference pores and initial fractures in the reservoir, they are input into the constructed carbonate reservoir porous effective medium model to calculate the longitudinal wave velocity V P，cal , the shear wave velocity V S,cal ;

[0056] Objective function unit: Use the objective function J to calculate the difference between the longitudinal wave velocity V P，cal , the shear wave velocity V S，cal calculated in the second longitudinal wave velocity calculation unit and their well logging measured values V P，obs , V S，obs . The expression of the objective function J is as follows:

[0057] J=|V P,obs -V P,cal | 2 +|V S,obs -V S,cal | 2

[0058] In the formula, V P,obs and V P,cal respectively represent the well logging measured longitudinal wave velocity and the longitudinal wave velocity V S,obs calculated based on the carbonate reservoir porous effective medium model, and V S,cal respectively represent the well logging measured longitudinal wave velocity and the shear wave velocity calculated based on the carbonate reservoir porous effective medium model;

[0059] Reference pore and fracture volume fraction determination unit: Traverse the reference pore and fracture volume fractions within a given range at a set step size, and terminate the iteration when the objective function J in the objective function unit is less than a given error σ, to obtain the final reference pore and fracture volume fractions;

[0060] Hard pore volume fraction determination unit: Based on the final reference pore and fracture volume fractions obtained by the reference pore and fracture volume fraction determination unit, calculate the final hard pore volume fraction;

[0061] Three-porosity calculation unit: Based on the final reference pore, fracture, and hard pore volume fractions obtained by the hard pore volume fraction determination unit, combined with the reservoir porosity, calculate the final reference pore, fracture, and hard pore porosities to complete the inversion.

[0062] In a third aspect, an electronic device is provided, including: a memory for storing a computer program; a processor for implementing the steps of the seismic inversion method for carbonate reservoir pore types based on a porous effective medium model when executing the computer program.

[0063] In a fourth aspect, a computer-readable storage medium is provided, on which a computer program is stored, and the computer program implements the steps of the seismic inversion method for carbonate reservoir pore types based on a porous effective medium model when executed by a processor.

[0064] By adopting the above technical solutions, the accuracy of multiple porosity inversion is improved, and the distribution of the three pore types in tight carbonate reservoirs can be accurately characterized.

[0065] Compared with the prior art, the beneficial effects of the present invention are:

[0066] The method proposed by the present invention simultaneously considers the coexistence of three pore types, performs multiple porosity inversion based on the extended Keys-Xu model and Gassmann-Hill equation established carbonate rock porous effective medium model, and adopts a strategy of using variable pore aspect ratios instead of fixed pore aspect ratios as input in the inversion process, improving the accuracy of multiple porosity inversion. This method can not only accurately characterize the distribution of the three pore types in tight carbonate reservoirs, but also successfully reveal the microscopic mechanism of oil and gas migration and accumulation; it provides a theoretical basis and technical support for characterizing the pore structure of ultra-deep tight reservoirs using well logging and seismic data. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a flowchart for implementing the seismic inversion method for pore types based on a porous effective medium model according to an embodiment of the present invention;

[0068] Figure 2It is a flowchart for implementing the specific details of the two-step inversion of multiple porosities in the pore type seismic inversion method based on the porous effective medium model according to an embodiment of the present invention;

[0069] Figure 3 It is a flowchart for implementing the method of constructing a carbonate rock porous effective medium model according to an embodiment of the present invention;

[0070] Figure 4 It is a flowchart for implementing the specific details of the method of constructing a carbonate rock porous effective medium model according to an embodiment of the present invention;

[0071] Figure 5 It is a flowchart for calculating the P-wave velocity of the HS upper and lower boundaries in the pore type seismic inversion method based on the porous effective medium model according to an embodiment of the present invention;

[0072] Figure 6 It is a structural diagram of the device for pore type seismic inversion based on the porous effective medium model according to an embodiment of the present invention;

[0073] Figure 7 It is the quantitative relationship between the pore type and porosity of limestone and dolomite reservoirs and the P-wave velocity and S-wave velocity according to an embodiment of the present invention;

[0074] Figure 8 It is a result comparison diagram of the pore type seismic inversion method based on the porous effective medium model according to an embodiment of the present invention;

[0075] Figure 9 It is to predict the spatial distribution characteristics of the total porosity and multiple porosities by the pore type seismic inversion method based on the porous effective medium model according to an embodiment of the present invention.

[0076] Figure 10 It is a structural diagram of an electronic device according to an embodiment of the present invention. Detailed implementation manners

[0077] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0078] It should be noted that the terms "first", "second", etc. in the description, claims and above-mentioned drawings of the present invention are used to distinguish similar objects, and do not necessarily describe a specific order or sequence. It should be understood that the data used in this way can be interchanged under appropriate circumstances, so that the embodiments of the present invention described here can be implemented in an order other than those illustrated or described here. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device that includes a series of steps or units does not necessarily have to be limited to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0079] Figure 1-2 Figures 5 are the implementation flowcharts of the pore type seismic inversion method based on the porous effective medium model, the specific implementation flowchart of the two-step inversion of multiple porosities, and the flowchart of calculating the P-wave velocity of the HS upper and lower boundaries provided by the embodiments of the present invention. For the convenience of description, only the parts related to the embodiments of the present invention are shown and are described in detail as follows:

[0080] The pore type seismic inversion method based on the porous effective medium model includes:

[0081] Step S1: Set the initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, and reservoir rock matrix properties;

[0082] Step S2: Calculate the P-wave velocities V P,wyllie 、V P,HS+ 、V P,HS- 、V P,KX respectively by using the Wyllie time-average equation, the HS upper and lower boundaries, and the single pore type Keys-Xu model;

[0083] Step S3: Calculate the differences between the P-wave velocities calculated by the Wyllie time-average equation and the HS upper and lower boundaries and the P-wave velocity calculated by the single pore type Keys-Xu model respectively. When the differences between the P-wave velocities calculated by the Wyllie time-average equation, the HS upper and lower boundaries and the P-wave velocity calculated by the single pore type Keys-Xu model all meet the cut-off conditions, terminate the iteration to obtain the final reservoir pore aspect ratio. The iteration calculation formula is as follows:

[0084] |V P,KX -V P,s |<ε, s = wyllie, HS+, HS-

[0085] In the formula, V P,KX represents the P-wave velocity calculated by the single pore type Keys-Xu model; V P,wyllie , V P,HS+ and V P,HS-respectively represent the longitudinal wave velocities calculated using the Wiley time-averaging equation, the HS upper bound, and the HS lower bound; ε represents a given minimum error;

[0086] Step S4: Based on the final reservoir pore aspect ratio obtained in Step S3, and in combination with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, reservoir initial reference pores, and the volume fractions of the reservoir initial fractures, input them into the constructed carbonate reservoir porous effective medium model to calculate the longitudinal wave velocity V P,cal and the shear wave velocity V S,cal ;

[0087] Step S5: Use the objective function J to calculate the longitudinal wave velocity V P,cal calculated in Step S4 above, and the shear wave velocity V S,cal and their logging measurement values V P,obs , V S,obs The error between them. The expression of the objective function J is as follows:

[0088] J = |V P,obs - V P,cal | 2 + |V S,obs - V S,cal | 2

[0089] In the formula, V P,obs and V P,cal respectively represent the longitudinal wave velocity measured by logging and the longitudinal wave velocity V S,obs calculated based on the carbonate reservoir porous effective medium model, and V S,cal respectively represent the longitudinal wave velocity measured by logging and the shear wave velocity calculated based on the carbonate reservoir porous effective medium model;

[0090] Step S6: Terminate the iteration when the objective function J in Step S5 is less than the given error σ by traversing the volume fractions of the reference pores and fractures within a given range in a set step size, and obtain the final volume fractions of the reference pores and fractures;

[0091] Step S7: Based on the final volume fractions of the reference pores and fractures obtained in Step S6, calculate the final volume fraction of the hard pores;

[0092] Step S8: Based on the final volume fractions of the reference pores, fractures, and hard pores obtained in Step S7, and in combination with the reservoir porosity, calculate the porosities of the final reference pores, fractures, and hard pores to complete the inversion.

[0093] In particular, based on the two-step porous type seismic inversion method process:

[0094] (1) Estimation of the aspect ratios of the three pore types

[0095] Set the initial pore aspect ratio and its search range, and calculate the longitudinal wave velocity of carbonate rocks containing rigid pores, reference pores and fractures through the Wyllie time-average equation:

[0096]

[0097] where φ represents porosity; V P,fl and V P,ma represent the longitudinal wave velocities of the fluid and the rock matrix, respectively.

[0098] Use the Hashin-Shtrikman (HS) bounds to calculate the elastic modulus of carbonate rocks containing hard pores, reference pores and fractures:

[0099]

[0100] where K 1 , G 1 and f 1 represent the bulk modulus, shear modulus and volume fraction of the first phase, respectively; K 2 , G 2 and f 2 represent the bulk modulus, shear modulus and volume fraction of the second phase. The calculation of the upper and lower bounds of the HS bounds is determined by swapping the order of the first phase and the second phase.

[0101] Use the HS bounds to calculate the longitudinal wave velocity of carbonate rock reservoirs containing rigid pores, reference pores and fractures:

[0102]

[0103] By assigning the initial pore aspect ratio, apply the Keys-Xu model for a single pore type to calculate the longitudinal wave velocity of carbonate rocks, and take the difference between the longitudinal wave velocities calculated by the Wyllie time-average equation and the HS bounds to update and iterate the pore aspect ratio. When the error between the longitudinal wave velocities calculated by the Wyllie time-average equation and the HS upper and lower bounds and the longitudinal and transverse wave velocities calculated by the Keys-Xu model for a single pore type satisfies the following condition, terminate the iteration to obtain the pore aspect ratio.

[0104] |V P,KX -V P,s |<ε, s = wyllie, HS+, HS-

[0105] where V P,KX represents the longitudinal wave velocity calculated by the Keys-Xu model for a single pore type; V P,wyllie , V P,HS+ and V P,HS- represent the longitudinal wave velocities calculated by the Wyllie time-average equation, the HS upper bound and the HS lower bound, respectively; ε represents the given minimum error;

[0106] Specifically, V P,KX represents the P-wave velocity calculated by the Keys-Xu model for a single pore type. The Keys-Xu model applicable to a single pore type can calculate the bulk modulus and shear modulus of dry rock:

[0107] K dry = K ma (1 - φ) P

[0108] G dry = G ma (1 - φ) Q

[0109] where φ represents porosity; P 0 and Q 0 are geometric factors related to the pore aspect ratio:

[0110]

[0111] where α represents the pore aspect ratio; T ijij (α) and T iijj (α) are functions of the pore aspect ratio.

[0112] Based on this, the bulk modulus and shear modulus of saturated rock are calculated using the Gassmann equation under the low-frequency assumption:

[0113]

[0114] G sat = G dry

[0115] where K fl represents the bulk modulus of the pore fluid.

[0116] The method for determining the bulk modulus and shear modulus of the overall saturated rock using the Gassmann-Hill equation includes:

[0117]

[0118] where the bulk modulus K sat,w of the water-saturated rock and the bulk modulus K sat,g of the gas-saturated rock are calculated from the Gassmann equation under the low-frequency assumption, and S w represents the water saturation. (Specifically, according to the above formula, the bulk modulus K sat,w of the water-saturated rock can be obtained (i.e., when K fl = K w in the Gassmann equation under the above low-frequency assumption). According to the above formula, the bulk modulus K sat,g(i.e., calculated when K in the Gassmann equation under the above low-frequency assumption fl = K g ).

[0119] Thus, the P-wave velocity calculated based on the Keys-Xu model with a single pore type is as follows:

[0120]

[0121] (2) Inversion of the percentages of three pore types

[0122] Adopt the following modeling process of the triple-porosity effective medium for carbonate rocks, combined with the pore aspect ratio, total porosity, rock matrix, and fluid properties calculated above as inputs to calculate the P-wave velocity and S-wave velocity of the rock. Use the objective function J to calculate the errors between the P-wave velocity and S-wave velocity calculated by the triple-porosity effective medium model for carbonate rocks and their measured values. If this error is greater than or equal to the given error σ, then search for the reference pore and fracture volume fractions that make the error less than the given error σ by traversing the reference pore and fracture volume fractions within a given range with a certain step size (grid search method). The objective function is as follows:

[0123] J = |V P,obs - V P,cal | 2 + |V S,obs - V S,cal | 2

[0124] In the formula, V P,obs and V P,cal respectively represent the P-wave velocity measured by logging and the P-wave velocity calculated based on the TPEM model;

[0125] V S,obs and V S,cal respectively represent the S-wave velocity measured by logging and the S-wave velocity calculated based on the TPEM model.

[0126] Specifically, first, set the initial pore aspect ratio (the initial pore aspect ratio, 0.5 for hard pores, 0.1 for reference pores, and 0.001 for fractures), porosity (subject to the actual situation, about 2% - 12% in the case of this article), fluid properties (subject to the actual situation, including the elastic modulus, density, and volume fraction of water and gas in the present invention), and rock matrix properties (subject to the actual situation, including the elastic modulus, density, and volume fraction of dolomite, calcite, quartz, and clay in the present invention). Using the Wyllie time-average equation, the HS boundary, and the single-pore-type Keys-Xu model to calculate the P-wave velocity, find the difference between the P-wave velocity calculated by the single-pore-type Keys-Xu model and the P-wave velocity calculated by the Wyllie time-average equation and the HS boundary. When the difference is greater than 0.1 km / s, traverse the given range of pore aspect ratio with a step size of 0.001 to find the pore aspect ratio that can make the error less than 0.1 km / s. Based on the fluid properties (including fluid elastic modulus, density, and volume fraction) and rock matrix properties (including rock matrix elastic modulus, density, and volume fraction) subject to the actual situation, the calculated pore aspect ratio, total porosity (about 2% - 12% in the present invention), and the given initial reference pore and fracture volume fractions (reference pore volume fraction is 1, fracture is 0), calculate the error according to the objective function. When the error is greater than 0.1 km / s, traverse the combination of reference pore and fracture volume fractions within the range of 0 - 1 with a step size of 0.01 to make the error meet the requirements, and subtract the volume fractions of the reference pore and fracture from 1 to obtain the volume fraction of solution pores; multiply them by the total porosity respectively to obtain three porosities. The mineral components and fluid volume fractions involved in this case are all obtained from well logging curves, and their elastic modulus and density are shown in Table 1.

[0127] Table 1 Elastic Modulus and Density of Each Mineral Component

[0128]

[0129]

[0130] Figure 3-4 The flowchart and specific implementation flowchart of the method for constructing a rock physical model of a hydrate reservoir with dual occurrence forms according to the embodiments of the present invention are shown. For ease of description, only the parts related to the embodiments of the present invention are shown and are described in detail as follows:

[0131] The method for constructing a porous effective medium model of a carbonate reservoir includes:

[0132] Step S41: Calculate the bulk modulus and shear modulus of the rock matrix using the Voigt-Reuss-Hill (VRH) model;

[0133] Step S42: Calculate the bulk modulus and shear modulus of the dry rock using the extended Keys-Xu model;

[0134] Step S43: Calculate the bulk modulus and shear modulus of the saturated rock using the Gassmann equation under the low-frequency assumption, and determine the bulk modulus and shear modulus of the overall saturated rock using the Gassmann-Hill equation in carbonate reservoirs.

[0135] Specifically, the rock physics model based on the porous assumption first uses VRH average mixing of dolomite, calcite, quartz, and clay to form the rock matrix; adds hard pores, reference pores, and fractures to it using the extended Keys-Xu model to form the dry rock skeleton, and finally uses the Gassmann-Hill equation to fill the pores with fluid to form the saturated rock.

[0136] Specifically, the rock physics modeling process for porous media in tight carbonate reservoirs:

[0137] (1) Calculation of the bulk modulus and shear modulus of the rock matrix

[0138] Carbonate reservoirs are mainly composed of calcite and dolomite, with a small amount of low quartz, anhydrite, and clay minerals. Given the bulk modulus of each mineral component, the shear modulus can be calculated using the Voigt-Reuss-Hill (VRH) average (Hill, 1952) to obtain the bulk modulus and shear modulus of the rock mechanism:

[0139]

[0140] In the formula, f i , K i and G i respectively represent the volume fraction, bulk modulus, and shear modulus of the i-th mineral component; M represents the total number of mineral components included in the rock matrix.

[0141] (2) Calculation of the bulk modulus and shear modulus of the dry rock

[0142] For carbonate reservoirs, the extended Keys-Xu model applicable to hard pores, reference pores, and fractures can calculate the bulk modulus and shear modulus of the dry rock:

[0143] K dry = K ma (1 - φ) P

[0144] G dry = G ma (1 - φ) Q

[0145] In the formula, φ represents the porosity; P and Q are geometric factors related to the pore aspect ratio:

[0146]

[0147] wherein, v s , v r and v c respectively represent the volume fractions of rigid pores, reference pores and cracks; α s , α r and α c respectively represent the pore aspect ratios of rigid pores (hard pores), reference pores and cracks; T ijij (α l ) and T iijj (α l ) are functions of the pore aspect ratio (Berryman, 1980), and their expressions are as follows:

[0148] T iijj (α) = 3F 1 / F 2

[0149]

[0150]

[0151] A = G j / G ma -1

[0152]

[0153] R = (1 - 2γ m ) / 2(1 - γ m )

[0154]

[0155] wherein, F 1～9 , θ and are all functions of the pore aspect ratio; A, B, R and γ m are all functions of the elastic modulus; K ma and G ma respectively represent the bulk modulus and shear modulus of the rock matrix; K j and G j respectively represent the bulk modulus and shear modulus of the inclusions; α represents the pore aspect ratio.

[0156] (3) Calculation of the bulk modulus and shear modulus of saturated rock

[0157] To determine the bulk modulus and shear modulus of fluid-saturated rock, the Gassmann equation under the low-frequency assumption is used to calculate the bulk modulus and shear modulus of saturated rock:

[0158]

[0159] G sat = G dry

[0160] Wherein, K fl represents the bulk modulus of the pore fluid.

[0161] According to the above formula, the bulk modulus K of the water-saturated rock can be obtained sat,w (that is, when K in the Gassmann equation under the above low-frequency assumption fl = K w is calculated).

[0162] According to the above formula, the bulk modulus K of the gas-saturated rock can be obtained sat,g (that is, when K in the Gassmann equation under the above low-frequency assumption fl = K g is calculated).

[0163] In a tight carbonate reservoir, oil, gas and water are not evenly distributed in the pore space. The Gassmann-Hill equation can determine the bulk modulus and shear modulus of the overall saturated rock:

[0164]

[0165] Wherein, the bulk modulus (K sat,w ) of the water-saturated rock and the bulk modulus (K sat,g ) of the gas-saturated rock are calculated by the Gassmann equation under the low-frequency assumption, and S w represents the water saturation.

[0166] Then, the longitudinal wave velocity calculated by the extended Keys-Xu model applicable to hard pores, reference pores and fractures can be obtained as:

[0167]

[0168] The shear wave velocity calculated by the extended Keys-Xu model applicable to hard pores, reference pores and fractures is:

[0169]

[0170] Figure 6 is the structural diagram of the pore type seismic inversion device based on the porous effective medium model provided by the embodiment of the present invention; for the convenience of description, only the parts related to the embodiment of the present invention are shown and are described in detail as follows:

[0171] The pore type seismic inversion device based on the porous effective medium model includes:

[0172] Initial setting unit: Set the initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, and reservoir rock matrix properties;

[0173] First P-wave velocity calculation unit: Calculate the P-wave velocities V P，wyllie , V P，HS+ , V P，HS- , V P，KX respectively by using the Wyllie time-average equation, the HS upper and lower bounds, and the single pore type Keys-Xu model;

[0174] Pore aspect ratio calculation unit: Subtract the P-wave velocities calculated by the Wyllie time-average equation and the HS upper and lower bounds from the P-wave velocity calculated by the single pore type Keys-Xu model respectively. When the differences between the P-wave velocities calculated by the Wyllie time-average equation, the HS upper and lower bounds and the P-wave velocity calculated by the single pore type Keys-Xu model all meet the cut-off condition, terminate the iteration to obtain the final reservoir pore aspect ratio. The iteration calculation formula is as follows:

[0175] |V P,KX - V P，i | < ε, i = wyllie, HS+, HS-

[0176] In the formula, V P,KX represents the P-wave velocity calculated by the single pore type Keys-Xu model; V P,wyllie , V P,HS+ and V P,HS- represent the P-wave velocities calculated by the Wyllie time-average equation, the HS upper bound and the HS lower bound respectively; ε represents the given minimum error;

[0177] Second P-wave velocity calculation unit: Based on the final reservoir pore aspect ratio obtained by the pore aspect ratio calculation unit, and combined with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, the volume fractions of the reservoir initial reference pores and reservoir initial fractures, input them into the constructed carbonate reservoir porous effective medium model to calculate the P-wave velocity V P,cal , shear wave velocity V S,cal ;

[0178] Objective function unit: Use the objective function J to calculate the error between the P-wave velocity V P,cal , shear wave velocity V S,cal calculated in the second P-wave velocity calculation unit and their well logging measured values V P,obs , V S,obs . The expression of the objective function J is as follows:

[0179] J = |V P,obs - V P,cal | 2 + |V S,obs - VS,cal | 2

[0180] Wherein, V P,obs and V P,cal respectively represent the longitudinal wave velocity measured by logging and the longitudinal wave velocity V calculated based on the porous effective medium model of carbonate reservoirs S,obs and V S,cal respectively represent the longitudinal wave velocity measured by logging and the shear wave velocity calculated based on the porous effective medium model of carbonate reservoirs;

[0181] Reference pore and fracture volume fraction determination unit: Traverse the reference pore and fracture volume fractions within a given range in a set step size, and terminate the iteration when the objective function J in the objective function unit is less than a given error σ, to obtain the final reference pore and fracture volume fractions;

[0182] Hard pore volume fraction determination unit: Calculate the final hard pore volume fraction based on the final reference pore and fracture volume fractions obtained by the reference pore and fracture volume fraction determination unit;

[0183] Three-porosity calculation unit: Based on the final reference pore, fracture, and hard pore volume fractions obtained by the hard pore volume fraction determination unit, and combined with the reservoir porosity, calculate the final reference pore, fracture, and hard pore porosities to complete the inversion.

[0184] Figure 7 is the quantitative relationship between the pore types and porosities of limestone and dolomite reservoirs and the longitudinal wave velocity and shear wave velocity provided by the embodiments of the present invention; for ease of description, only the parts related to the embodiments of the present invention are shown and are described in detail as follows:

[0185] The quantitative relationship between pore types and porosities and longitudinal wave velocity and shear wave velocity measured through experimental tests is as Figure 7 to verify the feasibility of the established model. Both the longitudinal wave and shear wave velocities decrease with the increase of porosity. However, the sensitivity of wave velocity to porosity varies with pore types. In the pore system dominated by hard pores and reference pores, the velocity decreases linearly with porosity but increases with the increase of the proportion of hard pores. In contrast, in the pore system containing fractures and reference pores, the velocity decreases sharply with the increase of porosity and fracture fraction. These findings indicate that fractures have a more significant impact on elastic properties compared to hard pores. Generally speaking, template analysis shows that the main pore types include reference pores, and fractures contribute significantly to the pores.

[0186] Figure 8 is the result comparison diagram of the pore type seismic inversion method based on the porous effective medium model provided by the embodiments of the present invention; for ease of description, only the parts related to the embodiments of the present invention are shown and are described in detail as follows:

[0187] The target layer for carbonate pore structure prediction research is taken from 2600 to 2755 meters of Well Mitan 1 in the Ordos Basin in the northwest of China. Based on the above modeling process, multiple porosity inversion is carried out with the given initial pore aspect ratios (the aspect ratio of rigid pores is 1.0, the aspect ratio of reference pores is 0.1, and the aspect ratio of fractures is 0.01).

[0188] Figure 8 Results of multiple porosity inversion of Well Mitan 1 based on the TPEM model. (a) P-wave velocity; (b) S-wave velocity; (c) Pore aspect ratio; (d) Rigid pore porosity; (e) Reference pore porosity; (f) Fracture porosity. In (a) and (b), the black solid lines represent the wave velocities measured by logging; the gray solid lines represent the wave velocities inverted with variable aspect ratios; the gray dashed lines represent the wave velocities inverted with fixed aspect ratios. In (c), the dark gray solid line represents the aspect ratio of rigid pores; the gray solid line represents the aspect ratio of reference pores; the light gray solid line represents the aspect ratio of fractures. In (d), (e) and (f), the black solid lines represent the corresponding porosities interpreted by logging; the gray solid lines represent the corresponding porosities inverted with variable aspect ratios; the gray dashed lines represent the corresponding porosities inverted with fixed aspect ratios.

[0189] Figure 8 The results of predicting multiple porosities based on the TPEM model with fixed pore aspect ratios and variable pore aspect ratios are shown. Whether the specified pore aspect ratio or the variable pore aspect ratio is used, the P-wave velocity and S-wave velocity predicted by this method are very close to the measured values of logging. The rigid pore porosity inverted by the TPEM with fixed pore aspect ratios is slightly higher than the result interpreted by logging. The error between the reference pore porosity and the logging interpretation value is small, and the fracture porosity cannot match the high-value part of the logging interpretation result; while the predicted results of the reference pores, rigid pores and fracture porosities based on the TPEM model with variable pore aspect ratios can basically match the logging interpretation results. This shows that adopting the strategy of variable pore aspect ratios in the process of multiple porosity prediction can effectively improve the prediction accuracy.

[0190] Figure 9 To predict the spatial distribution characteristics of total porosity and multiple porosities by the pore type seismic inversion method based on the porous effective medium model according to the embodiments of the present invention. For the sake of description, only the parts related to the embodiments of the present invention are shown and are described in detail as follows:

[0191] Figure 9 Spatial distribution characteristics of total porosity and multiple porosities predicted based on the porous effective medium model (TPEM model). (a) Total pore space distribution; (b) Rigid pore space distribution; (c) Reference pore space distribution; (d) Fracture space distribution.

[0192] Figure 9The distribution characteristics of the total porosity and multiple porosities inverted based on the triple-porosity effective medium model (TPEM model) in three-dimensional space. The results show that the total porosity is relatively high around Well Mitan 1, and reference pores and fractures are developed. Fracture development is also observed in faults and uplift zones, which is consistent with the geological conditions and logging measurement results in the Gaojiabao area of the Ordos Basin, proving that the pore type inversion method based on the TPEM model is effective for identifying favorable reservoirs.

[0193] Figure 10 The structural diagram of an electronic device provided according to an embodiment of the present invention is as Figure 10 shown. The device includes: a memory 21 for storing computer programs; a processor 22 for implementing the steps of the pore type seismic inversion method based on the triple-porosity effective medium model when executing the computer programs. Among them, the processor 22 may include one or more processing cores, such as a 4-core processor, an 8-core processor, etc. The processor 22 may be implemented in at least one hardware form of a digital signal processor (DSP), a field-programmable gate array (FPGA), or a programmable logic array (PLA). The processor 22 may also include a main processor and a coprocessor. The main processor is a processor for processing data in the wake state, also known as the central processing unit (CPU); the coprocessor is a low-power processor for processing data in the standby state. In some embodiments, the processor 22 may be integrated with a graphics processing unit (GPU), and the GPU is responsible for rendering and drawing the content to be displayed on the display screen. In some embodiments, the processor 22 may further include an artificial intelligence (AI) processor for processing computational operations related to machine learning.

[0194] The memory 21 may include one or more computer-readable storage media, which may be non-transitory. The memory 21 may also include high-speed random access memory and non-volatile memory, such as one or more magnetic disk storage devices and flash storage devices. In this embodiment, the memory 21 is at least used to store the following computer program 211. After the computer program is loaded and executed by the processor 22, it can implement the relevant steps of the pore type seismic inversion method based on the porous effective medium model disclosed in any of the foregoing embodiments. In addition, the resources stored in the memory 21 may also include an operating system 212 and data 213, etc., and the storage method may be temporary storage or permanent storage. Among them, the operating system 212 may include Windows, Unix, Linux, etc. The data 213 may include, but is not limited to, the data involved in the pore type seismic inversion method based on the porous effective medium model, etc.

[0195] In some embodiments, the electronic device may further include a display screen 23, an input / output interface 24, a communication interface 25, a power supply 26, and a communication bus 27.

[0196] Those skilled in the art can understand that Figure 10 the structure shown in

[0197] does not constitute a limitation on the electronic device, and may include more or fewer components than shown in the figure.

[0198] For the introduction of an electronic device provided by the present invention, please refer to the foregoing method embodiments. The present invention will not be elaborated herein again, and it has the same beneficial effects as the steps of the pore type seismic inversion method based on the porous effective medium model described above.

[0199] Furthermore, the present invention also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by the processor 22, it implements the steps of the pore type seismic inversion method based on the porous effective medium model as described above.

[0200] It can be understood that if the methods in the above embodiments are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solutions of the present invention, in essence, or the parts that contribute to the prior art, or all or part of the technical solutions, can be embodied in the form of a software product. This computer software product is stored in a storage medium and executes all or part of the steps of the methods of the various embodiments of the present invention. The aforementioned storage media include: various media such as USB flash drives, mobile hard disks, read-only memories (ROM), random access memories (RAM), magnetic disks, or optical discs that can store program codes.

[0201] For the introduction of a computer-readable storage medium provided by the present invention, please refer to the above method embodiments. The present invention will not repeat it here. It has the same beneficial effects as the steps of the pore type seismic inversion method based on the porous effective medium model.

[0202] The above has provided a detailed introduction to a pore type seismic inversion method, device, equipment, and medium based on a porous effective medium model of the present invention. The various embodiments in the specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. The same or similar parts among the various embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple. The relevant parts can be referred to the description in the method part. It should be noted that for those of ordinary skill in the art in this technical field, without departing from the principle of the present invention, several improvements and modifications can still be made to the present invention, and these improvements and modifications also fall within the protection scope of the present invention.

[0203] It should also be noted that in this specification, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise", or any other variant thereof are intended to cover non-exclusive inclusion, so that a process, method, article, or device including a series of elements not only includes those elements but also includes other elements not explicitly listed, or further includes elements inherent to such process, method, article, or device. Without further limitation, an element defined by the statement "including one..." does not exclude the existence of another identical element in the process, method, article, or device including the element.

Claims

1. A carbonate rock pore type seismic inversion method based on a porous effective medium model, characterized in that: include: Step S1: setting the initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, and reservoir rock matrix properties; Step S2: Calculate the P-wave velocity V using the Willy time-averaged equation, HS upper and lower boundaries, and the Keys-Xu model for a single pore type P，wyllie 、V P，HS+ 、V P，HS- 、V P，KX ; Step S3: The P-wave velocity calculated by the Willy time average equation, the HS upper and lower boundaries and the P-wave velocity calculated by the Keys-Xu model of a single pore type are respectively subtracted. When the difference between the P-wave velocity calculated by the Willy time average equation, the HS upper and lower boundaries and the P-wave velocity calculated by the Keys-Xu model of a single pore type meets the cutoff condition, the iteration is terminated to obtain the final reservoir pore aspect ratio. The iterative calculation formula is as follows: |V P,KX -V P,s |<ε,s=wyllie,HS+,HS- Where V P,KX represents the P-wave velocity calculated by the Keys-Xu model for a single pore type; V P,wyllie , V P,HS+ and V P,HS- denote the P-wave velocities calculated using the Wiley time-averaged equation, the HS upper bound, and the HS lower bound, respectively; ε represents the given minimum error; Step S4: Based on the final reservoir pore aspect ratio obtained in step S3, combined with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, reservoir initial reference pores and reservoir initial fracture volume fractions, it is input into the constructed carbonate reservoir porous effective medium model to calculate the P-wave velocity V P，cal , shear wave velocity V S,cal ; Step S5: Calculate the longitudinal wave velocity V calculated in step S4 using the objective function J P,cal , shear wave velocity V S,cal Its logging measurement value V P，obs 、V S，ovs The error between them, the expression of the objective function J is as follows: J=|V P,obs -V P,cal | 2 +|V S,obs -V S,cal | 2 Where V P,obs and V P,cal are the P-wave velocity measured by well logging and the P-wave velocity V calculated based on the porous effective medium model of carbonate reservoirs. S,obs and V S,cal They represent the P-wave velocity measured by well logging and the S-wave velocity calculated based on the porous effective medium model of carbonate reservoirs; Step S6: traversing the reference hole and crack volume fractions within a given range with a set step length so that the iteration is terminated when the objective function J in step S5 is less than a given error σ, and the final reference hole and crack volume fractions are obtained; Step S7: Based on the volume fractions of the final reference holes and cracks obtained in step S6, the volume fraction of the final hard holes is calculated; Step S8: Based on the volume fractions of the final reference holes, fractures, and hard holes obtained in step S7 and combined with the reservoir porosity, the porosity of the final reference holes, fractures, and hard holes is calculated to complete the inversion.

2. The carbonate rock pore type seismic inversion method based on the porous effective medium model according to claim 1 is characterized by: The method for constructing the porous effective medium model of the carbonate reservoir in step S4 includes: Step S41: using the Voigt-Reuss-Hill (VRH) model to calculate the rock matrix bulk modulus and shear modulus; Step S42: using the extended Keys-Xu model to calculate the bulk modulus and shear modulus of dry rock; Step S43: The bulk modulus and shear modulus of saturated rock are calculated using the Gassmann equation under the low-frequency assumption, and the bulk modulus and shear modulus of the overall saturated rock are determined using the Gassmann-Hill equation in the carbonate reservoir.

3. The carbonate rock pore type seismic inversion method based on the porous effective medium model according to claim 2 is characterized by: The method for calculating the rock matrix bulk modulus and shear modulus using the Voigt-Reuss-Hill (VRH) model in step S41 includes: In the formula, f i , K i and G i They represent the volume fraction, bulk modulus and shear modulus of the i-th mineral component respectively; M represents the total number of mineral components contained in the rock matrix.

4. The carbonate rock pore type seismic inversion method based on the porous effective medium model according to claim 3 is characterized by: The method for calculating the bulk modulus and shear modulus of dry rock using the extended Keys-Xu model in step S42 includes: The extended Keys-Xu model for hard holes, reference holes, and cracks calculates the bulk modulus and shear modulus of dry rock as follows: K dry =K ma (1-φ) P G dry =G ma (1-φ) Q Where φ represents the porosity; P and Q are geometric factors related to the pore aspect ratio: In the formula, v s , v r and v c Represent the volume fractions of hard pores, reference pores and cracks respectively; α s , α r and α c represent the pore aspect ratios of hard pores, reference pores and cracks, respectively; T ijij (α l ) and T iijj (α l ) is a function of the pore aspect ratio.

5. The carbonate rock pore type seismic inversion method based on the porous effective medium model according to claim 4 is characterized by: The method for calculating the bulk modulus and shear modulus of saturated rock using the Gassmann equation under the low-frequency assumption in step S43 includes: G sat =G dry In the formula, K fl represents the bulk modulus of the pore fluid.

6. The carbonate rock pore type seismic inversion method based on the porous effective medium model according to claim 5 is characterized by: The method for determining the bulk modulus and shear modulus of the overall saturated rock using the Gassmann-Hill equation in step S43 includes: Where, the bulk modulus K of water-saturated rock is sat,w and the bulk modulus K of gas-saturated rock sat,g Calculated by the Gassmann equation under the low frequency assumption, S w Indicates water saturation.

7. The carbonate rock pore type seismic inversion method based on porous effective medium model according to claim 1, characterized in that: In step S2, the longitudinal wave velocity V is calculated using the Willy time-averaged equation P，wyllie include: Where φ represents porosity; V P,fl and V P,ma represent the longitudinal wave velocities of fluid and rock matrix, respectively.

8. The carbonate rock pore type seismic inversion method based on porous effective medium model according to claim 1, characterized in that: In step S2, the longitudinal wave velocity V is calculated using the upper and lower boundaries of HS. P，HS+ 、V P，HS- Methods include: Step S21: Calculate the bulk modulus and shear modulus of carbonate rock using the Hashin-Shtrikman (HS) upper and lower boundaries: Where K1, G1 and f1 are the bulk modulus, shear modulus and volume fraction of the first phase respectively; K2, G2 and f2 are the bulk modulus, shear modulus and volume fraction of the second phase. The calculation of the upper and lower bounds of the HS boundary is determined by swapping the order of the first and second phases. "+" represents the upper bound and "-" represents the lower bound. Step S22: Calculate the carbonate reservoir P-wave velocity V based on the bulk modulus and shear modulus of the carbonate rock calculated in step S21 P，HS+ 、V P，HS- :

9. A carbonate rock pore type seismic inversion device based on a porous effective medium model, characterized in that: include: Initial setting unit: setting initial reservoir pore aspect ratio, reservoir porosity, reservoir fluid properties, reservoir rock matrix properties; The first P-wave velocity calculation unit: P-wave velocity V is calculated using the Willy time-averaged equation, HS upper and lower boundaries, and the Keys-Xu model for a single pore type. P，wyllie 、V P，HS+ 、V P，HS- 、V P，KX ; Pore ​​aspect ratio calculation unit: The P-wave velocity calculated by the Willy time average equation, the HS upper and lower boundaries, and the P-wave velocity calculated by the Keys-Xu model of a single pore type are respectively subtracted. When the difference between the P-wave velocity calculated by the Willy time average equation, the HS upper and lower boundaries, and the P-wave velocity calculated by the Keys-Xu model of a single pore type meets the cutoff condition, the iteration is terminated to obtain the final reservoir pore aspect ratio. The iterative calculation formula is as follows: |V P,KX -V P,s |<ε,s=wyllie,HS+,HS- Where V P,KX represents the P-wave velocity calculated by the Keys-Xu model for a single pore type; V P,wyllie , V P,HS+ and V P,HS- denote the P-wave velocities calculated using the Wiley time-averaged equation, the HS upper bound, and the HS lower bound, respectively; ε represents the given minimum error; Second P-wave velocity calculation unit: Based on the pore aspect ratio calculation unit, the final reservoir pore aspect ratio is obtained, and combined with the reservoir fluid properties, reservoir rock matrix properties, reservoir porosity, reservoir initial reference pores and reservoir initial fracture volume fractions, it is input into the constructed carbonate reservoir porous effective medium model to calculate the P-wave velocity V P，cal , shear wave velocity V S,cal ; Objective function unit: Calculate the longitudinal wave velocity V calculated in the second longitudinal wave velocity calculation unit using the objective function J P，cal , shear wave velocity V S，cal Its logging measurement value V P,obs 、V S,obs The error between them, the expression of the objective function J is as follows: J=|V P,obs -V P,cal | 2 +|V S,obs -V S,cal | 2 Where V P,obs and V P,cal are the P-wave velocity measured by well logging and the P-wave velocity V calculated based on the porous effective medium model of carbonate reservoirs. S,obs and V S,cal They represent the P-wave velocity measured by well logging and the S-wave velocity calculated based on the porous effective medium model of carbonate reservoirs; A reference hole and crack volume fraction determination unit: a method of traversing the reference hole and crack volume fractions within a given range with a set step length so that the iteration is terminated when the objective function J in the objective function unit is less than a given error σ, and the final reference hole and crack volume fractions are obtained; A hard hole volume fraction determination unit: based on the final reference hole and crack volume fractions obtained by the reference hole and crack volume fraction determination unit, the final hard hole volume fraction is calculated; Three porosity calculation units: Based on the volume fraction of the hard pores, the final reference pores, fractures, and hard pores obtained by the unit are determined, and combined with the reservoir porosity, the porosity of the final reference pores, fractures, and hard pores is calculated to complete the inversion.

10. An electronic device, characterized in that: include: Memory for storing computer programs; A processor is used to implement the steps of the carbonate rock pore type seismic inversion method based on a porous effective medium model as described in any one of claims 1 to 8 when executing the computer program.

Images