Shale reservoir pore pressure prediction method, device, equipment and medium

Abstract

The embodiment of the invention provides a shale reservoir pore pressure prediction method and device, equipment and a medium, and relates to the technical field of unconventional oil and gas shale reservoir prediction. The method comprises the following steps: collecting shale data of a target area, and constructing a rock physical model based on the shale data; removing the response of the longitudinal wave velocity to kerogen according to the rock physical model, removing the response of the longitudinal wave velocity and density to the gas content according to the rock physical model, and removing the response of the longitudinal wave velocity to anisotropy according to the rock physical model to obtain the longitudinal wave velocity after kerogen, gas content and anisotropy correction; and substituting the corrected longitudinal wave velocity into a first formula to obtain the single-well pore pressure of the target area, thereby improving the prediction precision of the single-well pore pressure.

Description

Technical Field

[0001] This invention relates to the field of unconventional oil and gas shale reservoir prediction technology, and in particular to a method, apparatus, equipment and medium for predicting pore pressure in shale reservoirs. Background Technology

[0002] Since pore pressure provides the primary driving force for fluid flow within pores, it is closely related to the recovery rate of shale reservoirs. Furthermore, pore pressure serves as a crucial parameter in drilling strategy development, providing a basis for analyses such as mud window determination and wellbore stability.

[0003] The most common method for quantitative prediction of pore pressure is through anomalies in well logging curves such as sonic velocity and resistivity. This mainly includes calculations based on empirical or semi-empirical formulas of compaction trends, and calculations using rock physics or rock mechanics models. The most representative methods include the Eaton and Bowers approaches. Eaton quantitatively characterizes pore pressure using overlying stress and the ratio of measured to normal compaction trends, based on the intrinsic relationship between effective stress and physical parameters. Bowers establishes equations for both sedimentary compaction loading and unloading processes, and predicts pore pressure based on the effective stress theorem.

[0004] However, the accuracy of pore pressure prediction based on normal compaction trends combined with empirical formulas, or by using well logging response anomalies, is affected by the normal compaction trend, while trends defined by linear or regional compaction have lower accuracy. Furthermore, due to the influence of various factors, pore pressure predictions based on well logging response anomalies related to sonic velocity, resistivity, and density contain errors in organic-rich shale formations. Therefore, pore pressure prediction based on curves influenced by multiple factors has high ambiguity, leading to low accuracy in predicted pore pressure. Summary of the Invention

[0005] To address the aforementioned technical issues, this invention provides a method, apparatus, equipment, and medium for predicting pore pressure in shale reservoirs, thereby improving the accuracy of single-well pore pressure prediction.

[0006] This invention provides a method for predicting pore pressure in shale reservoirs, the method comprising:

[0007] Collect shale data from the target area and construct a rock physical model based on the shale data;

[0008] The P-wave velocity after correction for kerogen, gas content, and anisotropy is obtained by removing the P-wave velocity response to kerogen from the rock physics model, removing the P-wave velocity and density response to gas content from the rock physics model, and removing the P-wave velocity response to anisotropy from the rock physics model.

[0009] Substituting the corrected P-wave velocity into the first formula, the single-well pore pressure in the target area is obtained;

[0010] The first formula is:

[0011] Among them, P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w Let R be the density of water, and R be the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. w.ave The average density of water in the depth range of 0 to z.

[0012] A second aspect of the present invention provides a shale reservoir pore pressure prediction device, the device comprising:

[0013] The physical model building module is used to collect shale data from the target area and build a rock physical model based on the shale data.

[0014] The velocity correction module is used to remove the response of P-wave velocity to kerogen according to the rock physics model, remove the response of P-wave velocity and density to gas content according to the rock physics model, and remove the response of P-wave velocity to anisotropy according to the rock physics model, so as to obtain the P-wave velocity after correction for kerogen, gas content and anisotropy.

[0015] The pressure prediction module is used to substitute the corrected P-wave velocity into the first formula to obtain the single-well pore pressure in the target area.

[0016] The first formula is:

[0017] Among them, P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w Let R be the density of water, and R be the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. w.ave The average density of water in the depth range of 0 to z.

[0018] A third aspect of the present invention provides an electronic device, including a processor, a memory, and a computer program stored in the memory and executable on the processor. When the computer program is executed by the processor, it implements the steps of the shale reservoir pore pressure prediction method as described in the first aspect of the present invention.

[0019] A fourth aspect of the present invention provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the steps of the shale reservoir pore pressure prediction method as described in the first aspect of the present invention.

[0020] The shale reservoir pore pressure prediction method of this invention corrects the influence of kerogen, anisotropy, and gas content on the P-wave velocity response by using a rock physics model constructed in the target area. Based on the P-wave velocity corrected for kerogen, gas content, and anisotropy, the pore pressure of a single well is predicted. This achieves the goal of the curve no longer being affected by multiple factors, and the curve reflects the pore pressure change with a simpler perturbation, thereby improving the prediction accuracy of pore pressure in the target area. This provides a basis for the optimal selection of shale reservoir engineering sites and the design of horizontal well trajectories, and also lays the foundation for the prediction of other geomechanical parameters. Attached Figure Description

[0021] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments of the present invention will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 This is a flowchart illustrating a method for predicting pore pressure in shale reservoirs according to an embodiment of the present invention;

[0023] Figure 2 This is a three-dimensional pore pressure prediction profile proposed in one embodiment of the present invention;

[0024] Figure 3 This is a flowchart of a method for determining the formation compaction coefficient and fracture density according to an embodiment of the present invention;

[0025] Figure 4 This is a single-well pore pressure prediction diagram shown in one embodiment of the present invention;

[0026] Figure 5 This is a structural block diagram of a shale reservoir pore pressure prediction device provided in an embodiment of the present invention;

[0027] Figure 6This is a schematic diagram of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0028] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0029] Reference Figure 1 , Figure 1 This is a flowchart illustrating a method for predicting pore pressure in shale reservoirs according to an embodiment of the present invention. Figure 1 As shown, the shale reservoir pore pressure prediction method of this embodiment may include the following steps:

[0030] Step S11: Collect shale data from the target area and construct a rock physical model based on the shale data.

[0031] In this embodiment, corresponding shale data can be collected for the target area, and a rock physics model for the target area can be constructed based on the collected shale data. The target area is the region where shale reservoir pore pressure prediction is required; for example, the target area could be a specific shale gas production area.

[0032] In one alternative example, a matrix simulation can be performed using the HS+Backus model based on shale data, and the optimal pore-throat ratio can be obtained by inverting the anisotropic SCA, followed by fluid simulation using Brown-Korringa to obtain an anisotropic rock physics model.

[0033] Step S12: Remove the response of P-wave velocity to kerogen based on the rock physics model, remove the responses of P-wave velocity and density to gas content based on the rock physics model, and remove the response of P-wave velocity to anisotropy based on the rock physics model, to obtain the P-wave velocity after correction for kerogen, gas content and anisotropy.

[0034] In this embodiment, considering the influence of factors such as kerogen, gas content, and anisotropy on the P-wave velocity, the P-wave velocity is corrected by removing the response of P-wave velocity to kerogen, the response of P-wave velocity and density to gas content, and the response of P-wave velocity to anisotropy based on the rock physics model corresponding to the target region, respectively, to obtain the P-wave velocity corrected for kerogen, gas content, and anisotropy.

[0035] Step S13: Substitute the corrected longitudinal wave velocity into the first formula to obtain the single-well pore pressure in the target area.

[0036] In this embodiment, after obtaining the P-wave velocity after correction for kerogen, gas content, and anisotropy, it can be based on Eato n The formula is obtained by substituting the P-wave velocity after correction for kerogen, gas content, and anisotropy into the first formula to solve for the single-well pore pressure in the target area.

[0037] The first formula is: P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w R is the density of water, R is the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy, and R is the elastic parameter used for pore pressure prediction. norm Here, R represents the depth trend of the elastic parameter used for pore pressure prediction. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. w.ave The average density of water in the depth range of 0 to z.

[0038] In this embodiment, a rock physics model constructed in the target area is used to correct the influence of kerogen, anisotropy, and gas content on the P-wave velocity response. That is, the influence of organic matter (such as kerogen), anisotropy, and gas content on the velocity curve response is eliminated. Based on the P-wave velocity corrected for kerogen, gas content, and anisotropy, the pore pressure of a single well is predicted. This achieves the goal that the curve is no longer affected by multiple factors, so that the velocity response anomaly can better reflect the single pore pressure. This improves the prediction accuracy of the pore pressure of a single well in the target area, providing a basis for the optimal selection of shale reservoir engineering sites and the design of horizontal well trajectories, and also laying the foundation for the prediction of other geomechanical parameters.

[0039] In conjunction with the above embodiments, in one implementation, the present invention also provides a method for predicting pore pressure in shale reservoirs. In addition to the steps described above, this method may further include steps S21 to S25:

[0040] In this embodiment, the collected shale data includes: measured P-wave velocity, S-wave velocity, and density curves; P-wave velocity, S-wave velocity, and density from well logging curves; mineral component content, mineral component modulus, and density obtained from well logging interpretation; kerogen content and fluid saturation; bulk modulus, shear modulus, and density of minerals and kerogen in the work area; fluid bulk modulus and density; formation pore pressure test results obtained from the laboratory; and P-wave velocity, S-wave velocity, and density obtained from seismic inversion.

[0041] In one optional example, the mineral components of the target area may include: clay, quartz, feldspar, calcite, dolomite, pyrite, and kerogen. Table 1 shows the bulk modulus, shear modulus, and density of each mineral component in one embodiment. Regarding the fluids, the shale reservoir in the work area contains two types of fluids: gas and water. The bulk modulus and density of the fluids are shown in Table 2, which is a table of bulk modulus and density of the fluids in one embodiment.

[0042] Table 1

[0043] Serial Number mineral Bulk modulus (GPa) Shear modulus (GPa) Density (g / cc) 1 clay 27 17 2.55 2 quartz 36.6 45 2.65 3 Feldspar 55 28 2.62 4 calcite 76.7 32 2.71 5 dolomite 94.9 45 2.87 6 Pyrite 158 149 5.02 7 kerogen 2.9 2.7 1.30

[0044] Table 2

[0045] Serial Number fluid Bulk modulus (GPa) Density (g / cc) 1 water 2.56 1 2 gas 0.038 0.15

[0046] Step S21: Difference the P-wave velocity after correction for kerogen, gas content, and anisotropy with the measured P-wave velocity, and calculate the absolute value diff.

[0047] In this embodiment, when determining the pore pressure of a single well using the first formula, a fitting coefficient is determined. Due to lateral variations in the region, 'n' may change when determining the fitting coefficient. Here, for three-dimensional pores, the fitting coefficients 'n' obtained from fitting different wells are weighted by inverse distance to obtain a planar diagram. For example... Figure 2 As shown, Figure 2 This is a three-dimensional pore pressure prediction profile proposed in one embodiment of the present invention. Figure 2 The primary target layer is the Longmaxi Formation shale gas reservoir.

[0048] In this embodiment, the absolute value diff is obtained by subtracting the obtained P-wave velocity (corrected for kerogen, gas content, and anisotropy) from the actual measured P-wave velocity. diff = |V P_correct -V P_mea |;Among them, V P_correct V is the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. P_mea This represents the measured longitudinal wave velocity.

[0049] Step S22: Using the absolute value diff as the target, fit the model parameters by using the P-wave velocity, S-wave velocity and density of the logging curve as input.

[0050] In this embodiment, after obtaining the absolute value diff, the model parameters are determined by fitting the model using the P-wave velocity, S-wave velocity, and density of the well logging curve as inputs, with this absolute value diff as the target. The fitting methods include, but are not limited to, multiple linear regression, BP neural network, RBF neural network, and random forest. In a preferred embodiment, the model with the smallest RMSE error can be selected for fitting to determine the model parameters, thus obtaining the model corresponding to those parameters.

[0051] Step S23: Using the model corresponding to the model parameters, input the P-wave velocity, S-wave velocity, and density obtained from the seismic inversion, and calculate the diff three-dimensional data volume.

[0052] In this embodiment, the P-wave velocity, S-wave velocity, and density obtained from the seismic inversion can be input into the model corresponding to the model parameters to calculate the diff three-dimensional data volume, that is, the diff three-dimensional data volume output by the model corresponding to the model parameters.

[0053] Step S24: Subtract the P-wave velocity obtained from the seismic inversion from the diff three-dimensional data volume to obtain the corrected three-dimensional P-wave velocity volume V. P_3D_correct ;

[0054] In this embodiment, the P-wave velocity obtained from seismic inversion (such as the P-wave velocity volume obtained from seismic inversion) can be subtracted from the diff three-dimensional data volume output by the model to obtain the corrected three-dimensional P-wave velocity volume V. P_3D_correct .

[0055] Step S25: Transfer the corrected three-dimensional longitudinal wave velocity volume V P_3D_correct The P-wave velocity R after correction for kerogen, gas content, and anisotropy will be the result obtained by fitting the measured P-wave velocity based on well data from the work area. norm Substituting into the first formula, we obtain P in the first formula. P At this time, P P The three-dimensional pore pressure of the target region.

[0056] In this embodiment, after obtaining the corrected three-dimensional longitudinal wave velocity volume V P_3D_correct Then, the corrected three-dimensional longitudinal wave velocity volume V P_3D_correctR in the first formula refers to the corrected three-dimensional P-wave velocity volume as the P-wave velocity R in the first formula after correction for kerogen, gas content, and anisotropy; and R in the first formula refers to the result obtained by fitting the measured P-wave velocity based on the well data of the work area. norm Substitute the values ​​into the first formula, calculate the first formula, and obtain P in the first formula. P At this time, P P The three-dimensional pore pressure of the target region is the three-dimensional pore pressure body of the target region.

[0057] In this embodiment, based on single-well pore pressure prediction, a three-dimensional pore pressure prediction method is provided. By correcting the velocity volume obtained from seismic inversion, the influence of other non-pressure factors on the velocity volume is removed. That is, the correction amount is calculated by subtracting the logging curves before and after correction, and the model is used to predict the change. After prediction, the velocity obtained from seismic inversion is adjusted to obtain the corrected velocity, thereby completing the prediction of three-dimensional pore pressure and improving the accuracy of three-dimensional pore pressure prediction.

[0058] In addition, in another optional embodiment, when defining trends during three-dimensional density extension, different trends can be defined for different lithologies based on the seismic inversion results, as well as different three-dimensional extension approaches.

[0059] In conjunction with the above embodiments, in one implementation, the present invention also provides a method for predicting pore pressure in shale reservoirs. In addition to the steps described above, this method may further include step S31:

[0060] Step S31: Based on the density three-dimensional volume extended to TVD=0, calculate the overlying stress according to the second formula.

[0061] In this embodiment, the overburden stress in the first formula can be obtained by solving the second formula based on a density three-dimensional volume extended to TVD=0, or a density curve extended to TVD=0. The second formula in this embodiment is: ρb is the density of the rock, g is the acceleration due to gravity, z is the depth, and ρ b.ave The average density of the rock in the depth range of 0 to z.

[0062] In one implementation, based on the above embodiments, since the measured density curve only contains data for a certain stratum, while calculating the overburden stress requires density data from TVD=0 (true vertical depth) to the predicted depth, density curve compensation is necessary. The work area contains three lithofacies: shale, sandstone, and limestone. The intersection of logging impedance and P-wave / S-wave velocity ratios is performed, with the colors representing the two data volumes of P-wave impedance and P-wave / S-wave velocity ratios obtained from pre-stack inversion within the work area. Based on the pre-stack inversion results, the intersection analysis determines the polygons for the three lithofacies (shale, sandstone, and limestone), and each polygon is replaced with the density-depth trend of the corresponding lithofacies. Using the XGBoost network, the relative density is predicted based on preferred sensitive seismic attributes, and the relative density and depth trends are merged in the frequency domain to obtain the density curve compensated to TVD=0.

[0063] Therefore, this embodiment proposes a density curve extension method (steps ① to ④). That is, the present invention also provides a method for predicting pore pressure in shale reservoirs. In this method, in addition to the steps described above, steps ① to ④ may also be included:

[0064] Step 1: Using the density data of the non-overpressured strata, fit a density variation trend with vertical depth according to the exponential relationship. The third formula is as follows:

[0065] X(z)=X Top +(X Base -X Top )×k×e cz .

[0066] In this embodiment, density data from the non-overpressured strata are first used to perform a rock physical depth trend analysis. A density variation trend with vertical depth is then fitted using a third formula based on an exponential relationship. In the third formula, c is the compaction coefficient (in 1 / m), determined based on the geological conditions and actual data of the work area; z is the burial depth (in meters); k is a correction coefficient, determined based on actual data; and X refers to elastic parameters (i.e., any elastic coefficient), including at least density and velocity. Top To fit the elastic parameter value of the trend peak, X Base Let X(z) be the elastic parameter value at the bottom of the fitted trend, and let X(z) be the trend of the fitted elastic parameter.

[0067] In one embodiment, R is substituted into the first formula when calculating the three-dimensional pore pressure. norm The result is obtained by fitting the measured P-wave velocity based on the well data of the work area in the manner described in step ① above.

[0068] Step 2: Based on the fourth formula, perform trend decomposition on the density curve to obtain the relative density component rel(X);

[0069] The fourth formula is: rel(X) = ln(X) - trend(lnX);

[0070] Take the logarithm of the density curve to obtain lnX, where X is the density. Then, use the third formula to calculate trend(lnX).

[0071] In this embodiment, the density curve needs to be decomposed based on the fourth formula (rel(X)=ln(X)-trend(lnX)) to obtain the relative density component rel(X). Then, the logarithm of the density curve is taken to obtain lnX (X refers to density here). The trend(lnX) is calculated using the method in step ① (i.e., the third formula), and then the relative density component rel(X) is obtained according to the fourth formula.

[0072] It is important to note that the work area may contain various lithologies, requiring density trend analysis for each lithology. Specifically: for well logging data, if lithology information is available, fit trends based on the density data of different lithologies. For seismic data, based on cross-plot analysis of seismic inversion results, circle polygons representing different lithologies in the cross-plot of different data volumes according to their distribution patterns. Each polygon can then be replaced with the density-depth trend of the corresponding lithology.

[0073] Step 3: Perform seismic transformation on the seismic data. The seismic transformation includes at least phase shift, attribute extraction, and seismic inversion. At the well point, the relative density curve is used as the output of the XGBoost model, and the seismic transformation is used as the input of the XGBoost model. Based on the importance of XGBoost model features, the seismic transformation features are selected.

[0074] Using the seismic transformation features as input and the relative density curve as output, the model is trained based on a machine learning model. The optimal parameters of the model are obtained using grid search, resulting in a trained network.

[0075] In this embodiment, seismic data needs to undergo seismic transformations such as phase shifting, attribute extraction, and seismic inversion. At the well points, the relative density curve is used as the output of the XGBoost model, and the seismic transformation is used as the input. Sensitive features are selected based on the feature importance of the XGBoost model. Feature importance can be calculated by assigning a weight to each feature's contribution to the result. Features with a cumulative contribution of 90% are selected as preferred features, thus obtaining the seismic transformation features. Then, the selected seismic transformation features are used as input, and the relative density curve is used as output. The model is trained using machine learning models such as XGBoost, and the optimal parameters are obtained using grid search, resulting in a trained network. Step ③ here can be understood as label creation and network training.

[0076] Step 4: Based on the trained network described in Step 3, the seismic transformation features are used as output to obtain the relative density three-dimensional volume;

[0077] The relative density three-dimensional volume is combined with the relative density components of different lithologies obtained in step ② to obtain the density three-dimensional volume from TVD=0 to the vertical depth to be predicted: X=e (rel(x)+trend(lnX)) .

[0078] In this embodiment, based on the network trained in step ③, seismic transformation features can be used as output to obtain a relative density three-dimensional volume. This relative density three-dimensional volume is then combined with the density trends of different lithologies (relative components of density for different lithologies) obtained in step ② to obtain a density three-dimensional volume from TVD = 0 to the vertical depth to be predicted: X = e (rel(x)+trend(lnX)) Therefore, the overlying stress is calculated based on this density three-dimensional volume.

[0079] In an alternative implementation, a rock physics model (i.e., a shale rock physics matrix model) can be constructed through the following steps:

[0080] Step (1): If the kerogen is in an over-mature stage, its development tends to be isotropic. Therefore, it can be mixed with quartz, carbonate minerals, etc. during the construction of the background matrix. Here, quartz, calcite, dolomite, pyrite, and kerogen can be mixed using the Hashin-Shritkman average boundary model.

[0081] K HS+ =Λ(μ max ), K HS- =Λ(μ min (1)

[0082] μ HS+ =Γ(ζ(K) max μ max )), μ HS- =Γ(ζ(K) min μ mix (2)

[0083]

[0084] Where K represents the bulk modulus in GPa; μ represents the shear modulus in GPa; and <.> represents the average based on the volume fraction of each mineral.

[0085] If kerogen is in its immature or mature stage, scanning electron microscopy reveals it to be a load-bearing mineral, predominantly in an elongated form. Increased kerogen content enhances matrix stratification; therefore, kerogen at this stage is added using the Backus average method. The Backus average calculation formula is as follows:

[0086]

[0087] C 11 =λ+2G (11)

[0088] C 33 =λ+2G (12)

[0089] C 13 =λ (13)

[0090] C 44 =G (14)

[0091] C 66 =G (15)

[0092]

[0093] Among them, C 11 C 13 C 33 C 44 C 66 Let be five independent elastic constants, representing the corresponding position coefficients of the elastic stiffness matrix. λ is the Lamé coefficient, K is the bulk modulus, and G is the shear modulus.

[0094] Step (2): Shale clay undergoes directional alignment under compaction. As compaction intensity increases, the degree of directional alignment also increases. This directional alignment of the clay increases the sedimentary stratification of the formation, leading to increased VTI anisotropy. With increased sedimentary stratification and VTI anisotropy, the P-wave velocity and S-wave velocity along the direction perpendicular to the formation decrease. Here, the clay compaction coefficient is used to characterize the formation compaction coefficient CF:

[0095] V p_clay_vertical (CF)=V P_clay_iso ×(1-CF) (17)

[0096] V s_clay_vertical (CF)=V s_clay_iso ×(1-CF) (18)

[0097] The bulk modulus and shear modulus of clay become:

[0098]

[0099] Among them, the Backus average formula (substituting (19) and (20) into formula (6)-(16)) can be used to add clay to the mixture obtained by modeling in step (1) above, and obtain the matrix stiffness coefficient matrix C without pores and cracks. 33m C 44m Based on C 33m C44m The matrix bulk modulus (K) was calculated. m ), shear modulus (G) m ) and density (ρ m ):in,

[0100]

[0101] G m =C 44m

[0102]

[0103] Step (3): Calculate the fluid properties in the shale reservoir: Using the type, saturation, and bulk modulus of the pore fluids from the collected shale data, calculate the mixed fluid properties based on the Wood model. The mixed fluid properties include the bulk modulus (K) of the mixed fluid. fluid ) and density (ρ fluid ).in, f i represents the fluid saturation level for each fluid.

[0104] Step (4): Adding cracks to the matrix of step (2): Adding cracks with pores of a certain shape and a certain density to the matrix to simulate the elastic modulus of a matrix containing pores and cracks:

[0105]

[0106] G d =C 44d

[0107]

[0108] ρ d =ρ m

[0109] Where, d c C is the crack density. 33d C 44d K represents the stiffness coefficient matrix parameters of dry rock. m G is the bulk modulus of the matrix. m For the matrix shear modulus, φ t K represents the total porosity. d G represents the bulk modulus of dry rock. d ρ is the shear modulus of dry rock. d This is the density of dry rock.

[0110] Step (5): Perform fluid displacement and calculate the P-wave and S-wave velocities of the rock: After modeling the dry rock, use the Brown-Korringa model to perform fluid displacement, filling the pores and fractures in this step with the fluid calculated in step (3), and calculate the bulk modulus (K) of the fluid-containing shale reservoir. sat ), shear modulus (G) sat ) and density (ρ b Based on the bulk modulus, shear modulus, and density, the P-wave and S-wave velocities (Vp) of shale reservoirs can be calculated. Psat V Ssat ).

[0111] Based on steps (1) to (5), a rock physics model (i.e., a shale rock physics matrix model) is constructed.

[0112] In conjunction with the above embodiments, the present invention also provides a method for predicting pore pressure in shale reservoirs. In this method, the step S12 above, "removing the response of P-wave velocity to kerogen according to the rock physics model," specifically includes steps S51 to S52:

[0113] Step S51: If the kerogen is in the over-mature stage, mix the quartz, calcite, dolomite, pyrite and kerogen using the Hashin-Shritkman average boundary model.

[0114] In this embodiment, for over-mature shale, kerogen is incorporated into the matrix using a Hashin-Shritkman boundary model. Specifically, when the kerogen is in an over-mature stage, quartz, calcite, dolomite, pyrite, and kerogen are mixed using a Hashin-Shritkman average boundary model. The above model is then used to simulate the constructed rock physics model to obtain P-wave velocity, S-wave velocity, and density.

[0115] Step S52: If the kerogen is in an immature or mature stage, the kerogen is added to the other non-clay mineral mixture in a Backus-like manner.

[0116] In this embodiment, for immature to mature shale, kerogen is added to the mixture of other non-clay minerals in a Backus-average manner. Specifically, the kerogen may be added to the mixture of other non-clay minerals in a Backus-average manner when it is in an immature or mature stage, thereby removing the response of P-wave velocity and density to gas content based on step S52 and / or step S51.

[0117] In one example, the target area is over-mature shale, and minerals 2-6 in Table 1 can be mixed using the Hashin-Shritkman boundary model, while clay mineral 1 can be mixed with the above non-clay mixture using Backus averaging.

[0118] In an alternative implementation, the influence of kerogen can be removed when constructing the rock physics model based on steps (1) to (5) above, and kerogen is no longer added during matrix modeling. The volume fractions of the various mineral components used in the calculations are adjusted as follows:

[0119]

[0120] Among them, f i For each original mineral component, f i_new The adjusted mineral composition was achieved by mixing non-clay minerals such as quartz, calcite, dolomite, and pyrite using the Hashin-Shritkman average boundary model. The response of P-wave velocity and density to gas content was removed.

[0121] In conjunction with the above embodiments, the present invention also provides a method for predicting pore pressure in shale reservoirs. In this method, the response of P-wave velocity to anisotropy is removed based on the rock physics model. This can be achieved by eliminating the influence of compaction parameters on clay velocity. Then, the VTI matrix modulus is calculated, and the fracture modulus parameters are replaced with matrix modulus and density using the Hudson-Cheng model based on fracture density to calculate the velocity. Specifically, this can be achieved through the following steps:

[0122] Considering that the factors that cause anisotropy include the depositional stratification described in step (2) and the crack density described in steps (4) and (5), the anisotropy is considered to be caused by factors including the depositional stratification described in step (2) and the crack density described in steps (4) and (5).

[0123] For the anisotropy caused by sedimentary stratification, step (2) is characterized by the compaction coefficient CF. Here, CF = 0, and the collected clay bulk modulus, clay modulus and density are used to add clay into the matrix using the Backus average formula.

[0124] To eliminate the anisotropy caused by crack development, the formula for calculating the elastic modulus of the matrix containing pores and cracks described in step (4) can be transformed into:

[0125]

[0126] This achieves the removal of the anisotropic response of the longitudinal wave velocity.

[0127] In conjunction with the above embodiments, the present invention also provides a method for predicting pore pressure in shale reservoirs. In this method, the step S12 above, "removing the response of P-wave velocity and density to gas content according to the rock physical model," may specifically include step S61, and in addition to the above steps, may also include steps S62 to S63:

[0128] Step S61: Use the Brown-Korringa model to perform fluid replacement, change the bulk modulus Kf of the fluid, and replace the fluid in the original formation with the bulk modulus of water.

[0129] In this embodiment, since the shale reservoir contains oil or gas, fluid replacement is performed using the Brown-Korringa model in step (5) to change the bulk modulus K of the fluid. f The bulk modulus of water is replaced by fluids (oil or gas) in the original formation, thereby eliminating the response of P-wave velocity and density to gas content.

[0130] Step S62: Calculate the bulk modulus, shear modulus, and density of the fluid-containing shale reservoir.

[0131] In this embodiment, after dry rock modeling, the Brown-Korringa model can be used to perform fluid replacement, filling the pores and fractures with the fluid calculated in the aforementioned step (3), and calculating the bulk modulus, shear modulus and density of the fluid-containing shale reservoir.

[0132] Step S63: Calculate the P-wave velocity and S-wave velocity of the shale reservoir based on the bulk modulus, shear modulus, and density of the fluid-containing shale reservoir.

[0133] In this embodiment, the longitudinal wave velocity and transverse wave velocity of the shale reservoir can be further obtained based on the calculated bulk modulus, shear modulus, and density of the fluid-containing shale reservoir.

[0134] In conjunction with the above embodiments, the present invention also provides a method for predicting pore pressure in shale reservoirs. In this method, step S62 may specifically include step S71, and step S63 may specifically include step S72.

[0135] Step S71: Calculate the bulk modulus K of the fluid shale reservoir using the fifth formula. sat Shear modulus G sat and density ρ b .

[0136] The fifth formula in this embodiment includes:

[0137]

[0138] G sat=G d ;

[0139] ρ b =(1-φ t )ρ d +ρ fluid φ t ;

[0140] Step S72: Calculate the P-wave velocity V of the shale reservoir using the sixth formula. Psat and transverse wave velocity V Ssat ;

[0141] The sixth formula in this embodiment includes:

[0142]

[0143] Where B is the equivalent pore pressure coefficient or effective stress coefficient, with a value between 0 and 1, and in shale reservoirs, the most common value is between 0.75 and 0.9; K d φ is the bulk modulus of dry rock. t G represents the total porosity. d K represents the shear modulus of dry rock. m ρ is the bulk modulus of the matrix. d ρ is the density of dry rock. fluid K is the density of the mixed fluid. f K is the bulk modulus of the fluid. S This represents the bulk modulus of the solid matrix.

[0144] Furthermore, in one embodiment, a method for predicting pore pressure in shale reservoirs is provided, which further includes determining the formation compaction coefficient and fracture density (i.e., obtaining the optimal formation compaction parameters and fracture density by comparing the measured P-wave velocity and the predicted P-wave velocity):

[0145] The formation compaction coefficient CF described in step (2) is disturbed at intervals of 0.01 between 0 and 1, and the fracture density d c The rock physics model was perturbed at intervals of 0.01 between 0 and 0.3 to calculate the P-wave velocity (Vp) of the shale reservoir. P_pre ) and transverse wave velocity (V S_pre ).

[0146] Construct an objective equation that minimizes the sum of the squared difference between the predicted P-wave velocity and the measured P-wave velocity, and the squared difference between the predicted S-wave velocity and the measured S-wave velocity.

[0147] f(CF,d c )=min((V P_pre -V P_mea ) 2 +(VS_pre -V S_mea ) 2 );

[0148] This allows us to obtain the clay compaction coefficient (CF) curve and fracture density curve of the shale reservoir. For example... Figure 3 As shown, Figure 3 This is a flowchart illustrating a method for determining the formation compaction coefficient and fracture density according to an embodiment of the present invention. Figure 3 In this process, the formation compaction coefficient and fracture density are first perturbed and then introduced into the rock physics model. Then, the P-wave and S-wave velocities are predicted. Based on the objective function and the measured P-wave and S-wave velocities, it is determined whether the optimal combination of compaction coefficient and fracture density can be output.

[0149] In conjunction with the above embodiments, in one embodiment, such as Figure 4 As shown, Figure 4 This is a single-well pore pressure prediction diagram illustrated in one embodiment of the present invention. Figure 4 In the image, from left to right, are the overburden stress, uncorrected (solid line), and corrected (dashed line) P-wave velocity, S-wave velocity, density, and single-well pore pressure predicted based on measured and corrected P-wave velocities. The black points represent formation pore pressure points measured experimentally. Figure 4 It can be seen that the pore pressure obtained after correction calculation has a higher degree of agreement with the measured value. This embodiment was tested in a shale reservoir development area. By comparing with the measured pore pressure, it can be seen that pore pressure prediction based on the corrected curve (in this embodiment) can effectively improve the accuracy of pore pressure prediction, avoiding the problem in related technologies where the input curve is affected by multiple factors rather than a single pore pressure, resulting in low accuracy of predicted pore pressure.

[0150] It should be noted that, for the sake of simplicity, the method embodiments are all described as a series of actions. However, those skilled in the art should understand that the embodiments of the present invention are not limited to the described order of actions, because according to the embodiments of the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to the embodiments of the present invention.

[0151] Based on the same inventive concept, one embodiment of the present invention provides a shale reservoir pore pressure prediction device 500. (Reference) Figure 5 , Figure 5 This is a structural block diagram of a shale reservoir pore pressure prediction device provided in an embodiment of the present invention. Figure 5 As shown, the device 500 includes:

[0152] The physical model construction module 501 is used to collect shale data from the target area and construct a rock physical model based on the shale data.

[0153] The velocity correction module 502 is used to remove the response of P-wave velocity to kerogen according to the rock physics model, remove the response of P-wave velocity and density to gas content according to the rock physics model, and remove the response of P-wave velocity to anisotropy according to the rock physics model, so as to obtain the P-wave velocity after correction for kerogen, gas content and anisotropy.

[0154] Pressure prediction module 503 is used to substitute the corrected longitudinal wave velocity into the first formula to obtain the single-well pore pressure in the target area.

[0155] The first formula is:

[0156] Among them, P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w Let R be the density of water, and R be the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. wave The average density of water in the depth range of 0 to z.

[0157] Optionally, the shale data includes at least one of the following: measured P-wave velocity, S-wave velocity, and density curves; P-wave velocity, S-wave velocity, and density from well logging curves; mineral composition content, kerogen content, and fluid saturation obtained from well logging interpretation; bulk modulus, shear modulus, and density of minerals and kerogen in the work area; bulk modulus and density of fluids; formation pore pressure test results obtained from the laboratory; and P-wave velocity, S-wave velocity, and density obtained from seismic inversion.

[0158] The device 500 further includes:

[0159] The absolute value determination module is used to calculate the absolute value diff by subtracting the P-wave velocity after correction for kerogen, gas content and anisotropy from the measured P-wave velocity.

[0160] diff = |V P_correct -V P_mea |;V P_correct V is the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. P_mea The measured longitudinal wave velocity;

[0161] The model parameter determination module is used to determine the model parameters by fitting the P-wave velocity, S-wave velocity and density of the logging curve as inputs, with the absolute value diff as the target.

[0162] The three-dimensional data determination module is used to calculate the diff three-dimensional data volume by using the model corresponding to the model parameters, inputting the P-wave velocity, S-wave velocity, and density obtained by the seismic inversion;

[0163] The three-dimensional velocity correction module is used to subtract the P-wave velocity obtained from the seismic inversion from the diff three-dimensional data volume to obtain the corrected three-dimensional P-wave velocity volume V. P_3D_correct ;

[0164] The three-dimensional pore pressure determination module is used to determine the corrected three-dimensional longitudinal wave velocity volume V. P_3D_correct The P-wave velocity R after correction for kerogen, gas content, and anisotropy will be the result obtained by fitting the measured P-wave velocity based on the well data from the work area. norm Substituting into the first formula, we obtain P in the first formula. P At this time, P P The three-dimensional pore pressure of the target region.

[0165] Optionally, the device 500 further includes:

[0166] A stress determination module is used to calculate the overlying stress based on a density three-dimensional volume extended to TVD=0, according to a second formula.

[0167] The second formula is:

[0168] Where, ρ b Let ρ be the density of the rock, g be the acceleration due to gravity, z be the depth, and ρ be the velocity of the rock. b.ave The average density of the rock in the depth range of 0 to z.

[0169] Optionally, the device 500 further includes:

[0170] The fitting module is used to fit a density trend with vertical depth using density data from non-overpressured strata, following an exponential relationship. The third formula is as follows:

[0171] X(z)=X Top +(X Base -X Top )×k×e cz ;

[0172] In the third formula, c is the compaction coefficient in units of 1 / m, z is the burial depth in units of meters, k is the correction coefficient, and X refers to elastic parameters, including at least density and velocity. TopTo fit the elastic parameter value of the trend peak, X Base Let X(z) be the elastic parameter value at the bottom of the fitted trend, and let X(z) be the trend of the fitted elastic parameter.

[0173] The component determination module is used to perform trend decomposition processing on the density curve based on the fourth formula to obtain the relative component rel(X) of the density;

[0174] The fourth formula is: rel(X) = ln (X -trend(lnX);

[0175] The logarithmic processing module is used to take the logarithm of the density curve to obtain lnX, where X is the density, and to calculate trend(lnX) using the third formula.

[0176] The feature transformation module is used to perform seismic transformation on seismic data. The seismic transformation includes at least phase shift, attribute extraction and seismic inversion. At the well point, the relative density curve is used as the output of the XGBoost model, and the seismic transformation is used as the input of the XGBoost model. Based on the feature importance of the XGBoost model, the seismic transformation features are selected.

[0177] The network training module is used to take the seismic transformation features as input and the relative density curve as output, train the network based on the machine learning model, and obtain the optimal parameters of the model by grid search to obtain the trained network.

[0178] The density determination module is used to obtain a relative density three-dimensional volume based on the trained network in the network training module and using the seismic transformation features as output.

[0179] The density three-dimensional volume determination module is used to merge the relative density three-dimensional volume with the relative components of densities of different lithologies obtained by the component determination module, to obtain the density three-dimensional volume from TVD=0 to the vertical depth to be predicted: X=e (rel(X)+trend( l nx)) .

[0180] Optionally, the speed correction module 502 includes:

[0181] The first correction module is used to mix quartz, calcite, dolomite, pyrite and kerogen using the Hashin-Shritkman average boundary model if the kerogen is in an over-mature stage.

[0182] The second correction module is used to add kerogen to the mixture of other non-clay minerals in a Backus-average manner if the kerogen is in an immature or mature stage.

[0183] Optionally, the speed correction module 502 includes:

[0184] The bulk modulus replacement module is used for fluid displacement using the Brown-Korringa model to change the bulk modulus K of the fluid. f The bulk modulus of water is used to replace the fluid in the original formation.

[0185] The device 500 further includes:

[0186] The first determining module is used to calculate the bulk modulus, shear modulus, and density of fluid-containing shale reservoirs;

[0187] The second determining module is used to determine the longitudinal wave velocity and transverse wave velocity of the shale reservoir based on the bulk modulus, shear modulus, and density of the fluid-containing shale reservoir.

[0188] Optionally, the first determining module includes:

[0189] The first determining submodule is used to calculate the bulk modulus K of the fluid shale reservoir using the fifth formula. sat Shear modulus G sat and density ρ b ;

[0190] The fifth formula includes:

[0191]

[0192]

[0193] G sat =G d ;

[0194] ρ b =(1-φ t )ρ d +ρ fluid φ t ;

[0195] The second determining module includes:

[0196] The second determining submodule is used to calculate the P-wave velocity V of the shale reservoir using the sixth formula. Psat and transverse wave velocity V Ssat ;

[0197] The sixth formula includes:

[0198]

[0199] Where B is the equivalent pore pressure coefficient or effective stress coefficient, with a value between 0 and 1; K d φ is the bulk modulus of dry rock. t G represents the total porosity.d K represents the shear modulus of dry rock. m ρ is the bulk modulus of the matrix. d ρ is the density of dry rock. fluid K is the density of the mixed fluid. f K is the bulk modulus of the fluid. S This represents the bulk modulus of the solid matrix.

[0200] Based on the same inventive concept, another embodiment of the present invention provides an electronic device 600, such as... Figure 6 As shown. Figure 6 This is a schematic diagram of an electronic device according to an embodiment of the present invention. The electronic device includes a processor 601, a memory 602, and a computer program stored in the memory 602 and executable on the processor 601. When the computer program is executed by the processor, it implements the steps in the shale reservoir pore pressure prediction method described in any of the above embodiments of the present invention.

[0201] Based on the same inventive concept, another embodiment of the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps in the shale reservoir pore pressure prediction method as described in any of the above embodiments of the present invention.

[0202] As the device embodiment is basically similar to the method embodiment, the description is relatively simple, and relevant parts can be found in the description of the method embodiment.

[0203] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0204] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, apparatus, or computer program products. Therefore, embodiments of the present invention can take the form of entirely hardware embodiments, entirely software embodiments, or embodiments combining software and hardware aspects. Furthermore, embodiments of the present invention can take the form of computer program products implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0205] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, terminal devices (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing terminal device to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing terminal device, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0206] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing terminal device to operate in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0207] These computer program instructions can also be loaded onto a computer or other programmable data processing terminal equipment, causing a series of operational steps to be performed on the computer or other programmable terminal equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable terminal equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0208] Although preferred embodiments of the present invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of the embodiments of the present invention.

[0209] Finally, it should be noted that in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or terminal device that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or terminal device. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or terminal device that includes said element.

[0210] The present invention provides a detailed description of a method, apparatus, equipment, and medium for predicting pore pressure in shale reservoirs. Specific examples have been used to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. At the same time, those skilled in the art will recognize that there will be changes in the specific implementation methods and application scope based on the ideas of the present invention. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for predicting pore pressure in shale reservoirs, characterized in that, The method includes: Collect shale data from the target area and construct a rock physical model based on the shale data; The P-wave velocity after correction for kerogen, gas content, and anisotropy is obtained by removing the P-wave velocity response to kerogen from the rock physics model, removing the P-wave velocity and density response to gas content from the rock physics model, and removing the P-wave velocity response to anisotropy from the rock physics model. Substituting the corrected P-wave velocity into the first formula, the single-well pore pressure in the target area is obtained; The first formula is: Among them, P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w Let R be the density of water, and R be the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. wave The average density of water in the depth range of 0 to z.

2. The method according to claim 1, characterized in that, The shale data includes: measured P-wave velocity, S-wave velocity, and density curves; P-wave velocity, S-wave velocity, and density from well logging curves; mineral composition content, kerogen content, and fluid saturation obtained from well logging interpretation; bulk modulus, shear modulus, and density of minerals and kerogen in the work area; bulk modulus and density of fluids; formation pore pressure test results obtained from the laboratory; and P-wave velocity, S-wave velocity, and density obtained from seismic inversion. The method further includes: The absolute value diff is calculated by subtracting the P-wave velocity after correction for kerogen, gas content, and anisotropy from the measured P-wave velocity. diff = |V P_correct -V P_mea |;V P_correct V is the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. P_mea The measured longitudinal wave velocity; Using the absolute value diff as the target, the P-wave velocity, S-wave velocity, and density of the well logging curve are used as inputs for fitting to determine the model parameters; Using the model corresponding to the model parameters, input the P-wave velocity, S-wave velocity, and density obtained from the earthquake inversion, and calculate the diff three-dimensional data volume; Subtracting the P-wave velocity obtained from the seismic inversion from the diff three-dimensional data volume yields the corrected three-dimensional P-wave velocity volume V. P_3D_correct ; The corrected three-dimensional longitudinal wave velocity body V P_3D_correct The P-wave velocity R after correction for kerogen, gas content, and anisotropy will be the result obtained by fitting the measured P-wave velocity based on the well data from the work area. norm Substituting into the first formula, we obtain P in the first formula. P At this time, P P The three-dimensional pore pressure of the target region.

3. The method according to claim 1 or 2, characterized in that, The method further includes: Based on a density three-dimensional volume extended to TVD=0, the overlying stress is calculated according to the second formula; The second formula is: Where, ρ b Let ρ be the density of the rock, g be the acceleration due to gravity, z be the depth, and ρ be the velocity of the rock. b.ave The average density of the rock in the depth range of 0 to z.

4. The method according to claim 3, characterized in that, The method further includes: Step ①: Using the density data of the non-overpressured strata, fit a density variation trend with vertical depth according to the exponential relationship. The third formula is as follows: X(z)=X Top +(X Base -X Top )×k×e cz ; In the third formula, c is the compaction coefficient in units of 1 / m, z is the burial depth in units of meters, k is the correction coefficient, and X refers to elastic parameters, including at least density and velocity. T op is the elasticity parameter value for fitting the peak of the trend, X Base Let X(z) be the elastic parameter value at the bottom of the fitted trend, and let X(z) be the trend of the fitted elastic parameter. Step 2: Based on the fourth formula, perform trend decomposition on the density curve to obtain the relative density component rel(X); The fourth formula is: rel(X) = ln(X) - trend(lnX); Take the logarithm of the density curve to obtain lnX, where X is the density. Calculate trend(lnX) using the third formula. Step 3: Perform seismic transformation on the seismic data. The seismic transformation includes at least phase shift, attribute extraction, and seismic inversion. At the well point, the relative density curve is used as the output of the XGBoost model, and the seismic transformation is used as the input of the XGBoost model. Based on the importance of XGBoost model features, the seismic transformation features are selected. Using the seismic transformation features as input and the relative density curve as output, the model is trained based on a machine learning model. The optimal parameters of the model are obtained by grid search, resulting in a trained network. Step 4: Based on the trained network described in Step 3, the seismic transformation features are used as output to obtain the relative density three-dimensional volume; The relative density three-dimensional volume is combined with the relative density components of different lithologies obtained in step ② to obtain the density three-dimensional volume from TVD=0 to the vertical depth to be predicted: X=e (rel(X)+trend(lnX)) 。 5. The method according to claim 1, characterized in that, The removal of the P-wave velocity response to kerogen based on the rock physics model includes: If the kerogen is in an over-mature stage, quartz, calcite, dolomite, pyrite and kerogen are mixed using the Hashin-Shritkman average boundary model. If the kerogen is in an immature or mature stage, it is added to the other non-clay mineral mixture in a Backus-like manner.

6. The method according to claim 1, characterized in that, The step of removing the response of P-wave velocity and density to gas content based on the rock physics model includes: The Brown-Korringa model was used to perform fluid replacement, changing the bulk modulus Kf of the fluid to replace the bulk modulus of water in the original formation. The method further includes: Calculate the bulk modulus, shear modulus, and density of fluid-containing shale reservoirs; Based on the bulk modulus, shear modulus, and density of the fluid-containing shale reservoir, the P-wave velocity and S-wave velocity of the shale reservoir are determined.

7. The method according to claim 6, characterized in that, The calculation of the bulk modulus, shear modulus, and density of fluid-containing shale reservoirs includes: The bulk modulus K of the fluid shale reservoir is calculated using the fifth formula. sat Shear modulus G sat and density ρ b ; The fifth formula includes: G sat =G d ; r b =(1-φ t )r d +r fluid f t ; The step of determining the P-wave velocity and S-wave velocity of the shale reservoir based on its bulk modulus, shear modulus, and density includes: The P-wave velocity V of the shale reservoir is calculated using the sixth formula. Psat and transverse wave velocity V Ssat ; The sixth formula includes: Where B is the equivalent pore pressure coefficient or effective stress coefficient, with a value between 0 and 1; K d φ is the bulk modulus of dry rock. t G represents the total porosity. d K represents the shear modulus of dry rock. m ρ is the bulk modulus of the matrix. d ρ is the density of dry rock. fluid K is the density of the mixed fluid. f K is the bulk modulus of the fluid. S This represents the bulk modulus of the solid matrix.

8. A device for predicting pore pressure in shale reservoirs, characterized in that, The device includes: The physical model building module is used to collect shale data from the target area and build a rock physical model based on the shale data. The velocity correction module is used to remove the response of P-wave velocity to kerogen according to the rock physics model, remove the response of P-wave velocity and density to gas content according to the rock physics model, and remove the response of P-wave velocity to anisotropy according to the rock physics model, so as to obtain the P-wave velocity after correction for kerogen, gas content and anisotropy. The pressure prediction module is used to substitute the corrected P-wave velocity into the first formula to obtain the single-well pore pressure in the target area. The first formula is: Among them, P Pnorm For hydrostatic pressure, P P For the pore pressure of a single well, σ v For overburden stress, ρ w Let R be the density of water, and R be the longitudinal wave velocity after correction for kerogen, gas content, and anisotropy. norm The longitudinal wave velocity represents the normal compaction depth trend, where n is the fitting coefficient, g is the gravitational acceleration, z is the depth, and ρ is the velocity. w.ave The average density of water in the depth range of 0 to z.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the computer program is executed by the processor, it implements the shale reservoir pore pressure prediction method as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program thereon, characterized in that, When the computer program is executed by the processor, it implements the shale reservoir pore pressure prediction method as described in any one of claims 1 to 7.

Images