An inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity

Abstract

This invention proposes an inversion-assisted shale petrophysical modeling method that considers the complexity of clay-kerogen stratification and pore structure. Pores are divided into isolated microfractures (related only to clay and kerogen) and connected pores (developed in various minerals), and the proportion of connected pores is introduced to balance the influence of these two types of pores. A stratification index is introduced to describe the degree of particle orientation. A modeling process is established based on inclusion theory and the Brown-Korringa model. Three important auxiliary parameters, including the stratification index, the proportion of connected pores, and the aspect ratio of connected pores, are estimated by comparing the errors between predicted elastic parameters and actual well logging data in conjunction with the inversion. The inversion results are used to correct the established petrophysical model, making it more consistent with reality and improving the prediction accuracy of shale rock elasticity and anisotropy parameters. Well logging data testing verifies the effectiveness of this method.

Description

Technical Field

[0001] This invention belongs to the field of shale rock physical modeling, specifically relating to an inversion-assisted shale rock physical modeling method that considers clay-kerogen stratification and pore complexity. Background Technology

[0002] In the early stages of oil and gas exploration, conventional reservoirs (such as sandstone and mudstone reservoirs) were generally considered isotropic media. However, with increasing global oil and gas demand, especially the diversification of energy needs and the increasing depletion of oil resources, the focus of exploration has gradually shifted to unconventional oil and gas resources. Unconventional oil and gas mainly includes shale oil, shale gas, and tight oil and gas. These resources are generally widely distributed, but their reservoir permeability is low, making extraction difficult. However, the physical properties of shale reservoirs differ significantly from those of traditional sandstone reservoirs, thus requiring new technologies for effective exploration and development. Shale reservoirs often exhibit significant anisotropy, especially VTI (vertically axisymmetric transversely isotropic) anisotropy. In VTI media, the propagation velocities of P-waves and S-waves differ considerably in different directions, and the wave velocity varies with direction. This anisotropic characteristic causes complex reflection characteristics during seismic wave propagation, making seismic exploration more difficult. The conventional isotropic assumption clearly cannot accurately describe the physical behavior of shale reservoirs. Since well-measured data on anisotropy parameters are difficult to obtain, it is necessary to estimate these parameters using rock physics models for subsequent research. Current shale rock physics modeling considers relatively limited factors, such as: 1. Constructing clay-rich shale models using a combination of anisotropic differential equivalent medium (DEM) and self-compatible model (SCA). 2. Using anisotropic Backus averaging to describe the anisotropy caused by clay and kerogen stratification in the shale matrix. 3. Using the standard deviation of the orientation distribution function to simulate clay-kerogen stratification. 4. Inverting the equivalent porosity aspect ratio of the rock based on P and S wave velocities, and using the obtained equivalent porosity aspect ratio for rock physics modeling to predict S wave velocities. 5. An inversion scheme that simultaneously estimates the stratification index and "clay-related" porosity aspect ratio to predict anisotropy parameters in vertical and horizontal logging. 6. A rock physics modeling method for tight sandstone reservoirs, which considers connected and isolated pores and defines a pore connectivity parameter to characterize the overall influence of different pore types on the rock elastic modulus.

[0003] In summary, current shale modeling does not simultaneously consider the combined effects of mineral grain stratification and complex pore structure on the rock's elastic properties. To improve the accuracy of shale reservoir exploration, it is essential to develop rock physical modeling and seismic exploration methods specifically for VTI media. These methods need to fully consider the inherent anisotropy and complex pore structure of shale, laying the foundation for subsequent shale seismic prediction. Summary of the Invention

[0004] This invention addresses the problems existing in the prior art by providing an inversion-assisted shale rock physics modeling method that considers the complexity of clay-kerogen stratification and pores. This method can correct the established rock physics model with inversion parameter results, making it more consistent with the actual situation, predicting the elasticity and anisotropy parameters of shale rocks, and improving prediction accuracy.

[0005] To address the above technical problems, this invention provides the following technical solution: an inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity, comprising the following steps:

[0006] S1. Based on the clay and kerogen content, establish the clay-kerogen mixture and calculate its stiffness matrix;

[0007] S2, with a porosity of φ iso Isolated, flat microcracks were added to the clay-kerogen mixture to create clay-kerogen blocks, and their stiffness matrix was calculated.

[0008] S3. By rotating and stacking clay-kerogen blocks according to a normal distribution, a shale matrix is ​​constructed, and its stiffness matrix is ​​calculated.

[0009] S4. Based on the content of brittle minerals, establish an isotropic brittle mineral mixture, add the brittle mineral mixture to the shale matrix to obtain the shale skeleton, and calculate its stiffness matrix.

[0010] S5, with a porosity of φ con Connected dry pores are added to the shale skeleton to obtain dry shale, and its stiffness matrix is ​​calculated.

[0011] S6. Fill the pores of dry shale with fluid or fluid mixture to obtain fluid-saturated shale, and calculate its stiffness matrix.

[0012] S7. Input the initial values ​​of the physical property parameters, calculate the velocity parameters using the stiffness matrix of fluid-saturated shale, establish the objective function, perturb the auxiliary parameters to minimize the objective function value, and output the auxiliary parameter inversion results.

[0013] S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is corrected, and the final prediction results of elasticity and anisotropy parameters are output based on the corrected shale rock physics model.

[0014] Furthermore, in the aforementioned step S1, the anisotropic SCA-DEM method is used to mix clay and kerogen, and the stiffness matrix C of the clay-kerogen mixture is... ckm Calculate as follows:

[0015]

[0016]

[0017] Among them, C ckm Let C be the stiffness matrix of the clay-kerogen mixture. n Let C be the stiffness matrix of the nth phase inclusion. p Let G be the stiffness matrix of the p-th phase inclusion body. ckm To calculate the Eshelby matrix corresponding to the clay-kerogen mixture, where I is the identity matrix, n is the nth phase inclusion, and ν n ν is the volume fraction of the nth phase inclusion. ker This represents the volume fraction of kerogen inclusions.

[0018] The volume fraction of kerogen inclusions was obtained by normalizing the measured mineral volume fractions in the well, as shown in the following formula:

[0019]

[0020] In the formula, v ker This represents the normalized kerogen volume fraction, and also the kerogen inclusion volume fraction. The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

[0021] Furthermore, step S2 as described above includes the following sub-steps:

[0022] S2.1. The pore space of shale is divided into two parts: isolated pore space and connected pore space. The proportions of isolated and connected pores are respectively denoted by r. iso r con It can be expressed as follows:

[0023] r con =φ con / φ

[0024] r iso =φ is o / φ=1-r con

[0025] Where φ is the total porosity, φ = φ con +φ iso , φ iso For isolated porosity, φ con For connectivity porosity, r iso r represents the proportion of isolated pores. con The proportion of connected pores;

[0026] S2.2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to construct clay-kerogen blocks containing isolated pores, whose stiffness matrix C ckb Calculate as follows:

[0027]

[0028] Among them, C ckb Here is the stiffness matrix of the clay-kerogen chunks. Let G be the stiffness matrix of the nth isolated pore inclusion. ckb This involves calculating the Eshelby matrix corresponding to clay-kerogen fragments, where I is the identity matrix, n is the nth isolated pore inclusion, and ν... iso The volume fraction of isolated pore inclusions is calculated using porosity and isolated porosity, satisfying the following expression: v iso =r iso =φ iso / φ=1-φ con / φ.

[0029] Furthermore, the aforementioned step S3 includes the following sub-steps:

[0030] S3.1 Based on the normal distribution of the rotation angle of clay-kerogen fragments in shale, the orientation distribution function ODF is as follows:

[0031] ODF~N(mean,LI)

[0032] Where ODF represents the orientation distribution function, N represents the normal distribution, mean represents the expected value of the rotation angle θ, and LI is the standard deviation of the normal distribution, representing the disorder of the clay-kerogen block deflection, i.e., the stratification index.

[0033] S3.2 The probability density function of the preferred orientation is as follows:

[0034]

[0035] Where υ(θ,LI) represents the probability density function of the preferred orientation, and θ is the rotation angle;

[0036] S3.3 The expressions for estimating the stiffness and flexibility matrices of the clay layer composed of rotated clay-kerogen blocks using VRH averaging are as follows:

[0037]

[0038]

[0039] Among them, C mix S represents the elastic stiffness matrix of the mixture obtained by the Voigt method.mix This represents the elastic compliance matrix of the mixture obtained by the Reuss method. and These are the stiffness and compliance matrices of the clay-kerogen mixture after Bond transformation, respectively. It is the azimuth angle;

[0040] S3.4 Calculate the average VRH, as shown in the following expression:

[0041] C clay =[C mix +(S mix ) -1 ] / 2

[0042] Among them, C clay This represents the stiffness matrix of the clay layer composed of clay-kerogen blocks after rotation.

[0043] Furthermore, the aforementioned step S4 includes the following sub-steps:

[0044] S4.1. Using the VRH averaging method, establish the stiffness matrix C of an isotropic brittle mineral mixture. brit It is calculated by the following formula:

[0045]

[0046] in,

[0047]

[0048] In the formula, C brit For the stiffness matrix of a brittle mineral mixture, μ brit K is the shear modulus of a brittle mineral mixture. brit The bulk modulus of a brittle mineral mixture. The shear modulus of a brittle mineral mixture obtained by the Voigt method. The shear modulus of the brittle mineral mixture obtained by the Reuss method. The bulk modulus of a brittle mineral mixture obtained by the Voigt method. The bulk modulus of the brittle mineral mixture obtained by the Reuss method;

[0049] S4.2. Using the anisotropic DEM method, a brittle mineral mixture is added to the shale matrix formed by the rotational stacking of clay and kerogen to establish the shale skeleton, whose stiffness matrix C... sk It is calculated by the following formula:

[0050]

[0051] Among them, C skHere is the stiffness matrix of the shale matrix. Let G be the stiffness matrix of the nth brittle mineral inclusion. sk This calculates the Eshelby matrix corresponding to the shale matrix, where I is the identity matrix, n is the nth brittle mineral inclusion, and v sk This represents the percentage content of inclusions in brittle mineral mixtures; the percentage content of inclusions in brittle mineral mixtures is obtained by normalizing the volume fraction of brittle minerals, as shown in the following formula:

[0052]

[0053] In the formula, The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

[0054] Furthermore, step S5 described above specifically involves: using an anisotropic DEM to determine the porosity φ. con Connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry It is calculated by the following formula:

[0055]

[0056] Among them, C dry Here is the stiffness matrix of dry shale. Let G be the stiffness matrix of the nth connected pore inclusion. dry This calculates the Eshelby matrix corresponding to dry shale, where I is the identity matrix, n is the nth connected pore inclusion, and v con The percentage content of connected pore inclusions, where the percentage content of connected pore inclusions satisfies the following expression: v con =r con =φ con / φ

[0057] Where φ is the total porosity, φ con For connectivity porosity, r con The percentage of connected pores.

[0058] Furthermore, step S6 described above specifically involves: filling the pores of dry shale with fluid or a fluid mixture based on the BK model to obtain fluid-saturated shale, whose compliance matrix S sat It is calculated by the following formula;

[0059]

[0060]

[0061] Among them, S dry S is the compliance matrix of dry shale. sat S is the compliance matrix of saturated fluid shale. m β is the equivalent elastic compliance matrix of the constituent minerals. fl For the compressibility of pore fluids, β m For the compressibility of minerals, K fl Let φ be the bulk modulus of the porous fluid, and φ be the porosity, where i, j, k, l, a, b, c, and d are the coordinate axis indices of the compliance matrix.

[0062] Furthermore, the aforementioned step S7 includes the following sub-steps:

[0063] S7.1, Set the aspect ratio α of the connecting pores con The proportion of interconnected pores r con The layering index LI is used as an important auxiliary parameter to be inverted. Initial values ​​for the auxiliary parameter are set to begin the inversion process, and α is obtained. con r con The inversion results of LI, with the inversion objective parameter being α. con r con And LI, the initial values ​​of the auxiliary parameters are set as follows:

[0064]

[0065] Where m0 is the initial auxiliary parameter vector, The initial value of the aspect ratio of the connected pores. LI is the initial value of the pore connectivity index. 0 This is the initial value of the stratification index;

[0066] S7.2. Set the inversion objective function as follows:

[0067]

[0068] Where, m i Let be the auxiliary parameter value for the i-th iteration. This is the prediction result when the objective function error is minimized. The predicted P-wave velocity for the constructed shale model. ρ represents the predicted shear wave velocity from the constructed shale model. pred (m i The density is predicted by the established shale model. The measured P-wave velocity in the well. ρ is the transverse wave velocity measured in the well. real This represents the density measured in the well. σ1 is the weight value of the P-wave velocity data item, σ2 is the weight value of the S-wave velocity data item, and σ3 is the weight value of the density data item. and All were calculated from the compliance matrix of the shale model;

[0069] S7.3. Update the auxiliary parameters using an iterative method. When the iteration number i = 1, the auxiliary parameters satisfy... When the iteration number i ≥ 2, the physical property parameters satisfy in, The aspect ratio of the connected pores at the i-th iteration is given. Let LI be the pore connectivity index value at the i-th iteration. i The stratification index value at the i-th iteration; when the objective function value When it is at its minimum, the corresponding This is the final inversion result.

[0070] Furthermore, the aforementioned step S8 includes the following sub-steps:

[0071] S8.1, Based on obtaining three auxiliary parameters α con r con Inversion results with LI The established shale model is modified to better reflect real formation conditions. The stiffness matrix of the saturated fluid shale is then expressed as follows:

[0072]

[0073] Where, m inv The inversion result vector, To connect the aspect ratio inversion results of the pores, For the inversion results of the proportion of interconnected pores, LI inv For the stratification index inversion result, G represents the forward modeling operator for shale petrophysical modeling in steps S1 to S6. To correct the new stiffness matrix predicted by the shale model;

[0074] S8.2. Based on the stiffness matrix of saturated fluid shale, the final elastic and anisotropic parameters are obtained, including the P-wave velocity, S-wave velocity, density, and three Thomsen anisotropic parameters predicted by the modified shale rock physics model.

[0075] Compared with the prior art, the beneficial technical effects of the present invention using the above technical solution are as follows:

[0076] (1) The shale rock physical modeling method of the present invention considers the stratification of clay-kerogen and introduces a stratification index LI to simulate the degree of deflection disorder of clay-kerogen fragments. It also considers pore complexity, dividing pores into isolated pores and connected pores, and introduces a pore connectivity index r. conThis is used to characterize the influence of the two pore types, making the constructed model more consistent with the actual formation conditions.

[0077] (2) Minimize the residuals between the predicted values ​​of P-wave and S-wave velocities and the measured logging data, and iteratively invert the three auxiliary parameters α. con r con The initial rock physics model was then corrected using LI (Liquidity Inversion). The inversion-assisted modeling strategy was applied to actual well logging data, and the errors in P-wave and S-wave velocities predicted by rock physics models with fixed auxiliary parameters and those with auxiliary parameters corrected with depth were compared. The results show that the modeling method considering three auxiliary parameters simultaneously yields the highest prediction accuracy and the best consistency with the measured data. Attached Figure Description

[0078] Figure 1 Flowchart for inversion-aided shale rock physical modeling.

[0079] Figure 2 Flowchart for modeling a shale rock physics model that takes into account clay-kerogen stratification and isolated / connected pores.

[0080] Figure 3 The diagram shows the measured data of the test well. In the diagram, (a) is a schematic diagram of P-wave velocity, (b) is a schematic diagram of S-wave velocity, (c) is a schematic diagram of density, (d) is a schematic diagram of porosity, and (e) is a schematic diagram of mineral composition.

[0081] Figure 4 The figure shows the prediction results when all three auxiliary parameters are constant. In the figure, (a) is a schematic diagram of the measured and predicted elastic parameters, and (b) is a schematic diagram of the predicted anisotropy parameters.

[0082] Figure 5 The figure shows the prediction results of the inversion-assisted rock physics modeling of the present invention. In the figure, (a) is a schematic diagram of the measured and predicted elastic parameters and the inverted auxiliary parameters, and (b) is a schematic diagram of the predicted anisotropy parameters. Detailed Implementation

[0083] To better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0084] In this invention, various aspects of the invention are described with reference to the accompanying drawings, in which numerous illustrative embodiments are shown. Embodiments of the invention are not limited to those depicted in the drawings. It should be understood that the invention is implemented through any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed herein are not limited to any particular implementation. Furthermore, some aspects of the invention disclosed may be used alone or in any suitable combination with other aspects of the invention disclosed.

[0085] The example implementation uses a modeling test based on real well logging data as an example for illustration. Figure 1 , Figure 2 As shown, this invention discloses an inversion-assisted shale petrophysical modeling method that considers clay-kerogen stratification and pore complexity, comprising the following steps:

[0086] S1. Based on the clay and kerogen content, establish the clay-kerogen mixture structure and calculate its stiffness matrix. Details are as follows:

[0087] The stiffness matrix C of the clay-kerogen mixture was obtained by mixing clay and kerogen using the anisotropic SCA-DEM method. ckm Calculate as follows:

[0088]

[0089]

[0090] Among them, C ckm Let C be the stiffness matrix of the clay-kerogen mixture. n Let C be the stiffness matrix of the nth phase inclusion. p Let G be the stiffness matrix of the p-th phase inclusion body. ckm To calculate the Eshelby matrix corresponding to the clay-kerogen mixture, where I is the identity matrix, n is the nth phase inclusion, and ν n ν is the volume fraction of the nth phase inclusion. ker This represents the volume fraction of kerogen inclusions.

[0091] The volume fraction of kerogen inclusions was obtained by normalizing the measured mineral volume fractions in the well, as shown in the following formula:

[0092]

[0093] In the formula, v ker This represents the normalized kerogen volume fraction, and also the kerogen inclusion volume fraction. The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

[0094] S2, with a porosity of φ iso Isolated, flattened microcracks are added to a clay-kerogen mixture to create clay-kerogen blocks, and their stiffness matrix is ​​calculated. The specific steps include the following:

[0095] S2.1. The pore space of shale is divided into two parts: isolated pore space and connected pore space. The proportions of isolated and connected pores are respectively denoted by r. isor con It can be expressed as follows:

[0096] r con =φ con / φ

[0097] r iso =φ is o / φ=1-r con

[0098] Where φ is the total porosity, φ = φ con +φ iso , φ iso For isolated porosity, φ con For connectivity porosity, r iso r represents the proportion of isolated pores. con The percentage of connected pores.

[0099] S2.2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to construct clay-kerogen blocks containing isolated pores, whose stiffness matrix C ckb Calculate as follows:

[0100]

[0101] Among them, C ckb Here is the stiffness matrix of the clay-kerogen chunks. Let G be the stiffness matrix of the nth isolated pore inclusion. ckb This involves calculating the Eshelby matrix corresponding to clay-kerogen fragments, where I is the identity matrix, n is the nth isolated pore inclusion, and ν... iso The volume fraction of isolated pore inclusions is calculated using porosity and isolated porosity, satisfying the following expression: v iso =r iso =φ iso / φ=1-φ con / φ.

[0102] S3. By rotating and stacking clay-kerogen blocks according to a normal distribution, a shale matrix is ​​constructed, and its stiffness matrix is ​​calculated. The coupled clay-kerogen in shale enhances its anisotropy. To analyze the influence of clay stratification in shale, the clay layer is considered as a complex whole composed of many clay-kerogen blocks with identical properties. The rotation angle θ of each clay-kerogen block (the angle between the block's axis of symmetry and the vertical direction) is changed to simulate the layered distribution of clay in actual shale formations. The elastic properties of clay-kerogen blocks with different rotation angles can be obtained by performing a Bond transformation on the horizontal blocks. Then, VRH averaging is used to calculate the equivalent elastic properties of the stacked clay-kerogen blocks. Specifically, the following sub-steps are included:

[0103] S3.1 Based on the normal distribution of the rotation angle of clay-kerogen fragments in shale, the orientation distribution function ODF is as follows:

[0104] ODF~N(mean,LI)

[0105] Where ODF represents the orientation distribution function, N represents the normal distribution, and mean represents the expected value of the rotation angle θ. Since the main stratification direction of the strata is horizontal, mean is usually set to zero. LI is the standard deviation of the normal distribution, representing the degree of disorder in the deflection of clay-kerogen parcels, i.e., the stratification index;

[0106] S3.2 The probability density function of the preferred orientation is as follows:

[0107]

[0108] Where υ(θ,LI) represents the probability density function of the preferred orientation, and θ is the rotation angle.

[0109] S3.3 The expressions for estimating the stiffness and flexibility matrices of the clay layer composed of rotated clay-kerogen blocks using VRH averaging are as follows:

[0110]

[0111]

[0112] Among them, C mix S represents the elastic stiffness matrix of the mixture obtained by the Voigt method. mix This represents the elastic compliance matrix of the mixture obtained by the Reuss method. and These are the stiffness and compliance matrices of the clay-kerogen mixture after Bond transformation, respectively. This is the azimuth angle.

[0113] S3.4 Calculate the average VRH, as shown in the following expression:

[0114] C clay =[C mix +(S mix ) -1 ]] / 2

[0115] Among them, C clay This represents the stiffness matrix of the clay layer composed of clay-kerogen blocks after rotation.

[0116] S4. Based on the content of brittle minerals, establish an isotropic brittle mineral mixture, add the brittle mineral mixture to the shale matrix to obtain the shale skeleton, and calculate its stiffness matrix.

[0117] S4.1. Using the VRH averaging method, establish the stiffness matrix C of an isotropic brittle mineral mixture. brit It is calculated by the following formula:

[0118]

[0119] in,

[0120]

[0121] In the formula, C brit For the stiffness matrix of a brittle mineral mixture, μ brit K is the shear modulus of a brittle mineral mixture. brit The bulk modulus of a brittle mineral mixture. The shear modulus of a brittle mineral mixture obtained by the Voigt method. The shear modulus of the brittle mineral mixture obtained by the Reuss method. The bulk modulus of a brittle mineral mixture obtained by the Voigt method. The bulk modulus of the brittle mineral mixture obtained by the Reuss method;

[0122] S4.2. Using the anisotropic DEM method, a brittle mineral mixture is added to the shale matrix formed by the rotational stacking of clay and kerogen to establish the shale skeleton, whose stiffness matrix C... sk It is calculated by the following formula:

[0123]

[0124] Among them, C sk Here is the stiffness matrix of the shale matrix. Let G be the stiffness matrix of the nth brittle mineral inclusion. skThis calculates the Eshelby matrix corresponding to the shale matrix, where I is the identity matrix, n is the nth brittle mineral inclusion, and v sk This represents the percentage content of inclusions in brittle mineral mixtures; the percentage content of inclusions in brittle mineral mixtures is obtained by normalizing the volume fraction of brittle minerals, as shown in the following formula:

[0125]

[0126] In the formula, The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

[0127] S5, with a porosity of φ con Connected dry pores are added to the shale skeleton to obtain dry shale, and its stiffness matrix is ​​calculated. Specifically, an anisotropic DEM is used to define the stiffness matrix of shale with porosity φ. con Connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry It is calculated by the following formula:

[0128]

[0129] Among them, C dry Here is the stiffness matrix of dry shale. Let G be the stiffness matrix of the nth connected pore inclusion. dry This calculates the Eshelby matrix corresponding to dry shale, where I is the identity matrix, n is the nth connected pore inclusion, and v con The percentage content of connected pore inclusions, where the percentage content of connected pore inclusions satisfies the following expression: v con =r con =φ con / φ, where φ is the total porosity, φ con For connectivity porosity, r con The percentage of connected pores.

[0130] S6. Fill the pores of dry shale with a fluid or fluid mixture to obtain fluid-saturated shale, and calculate its stiffness matrix. Specifically: Based on the BK model, fill the pores of dry shale with a fluid or fluid mixture to obtain fluid-saturated shale, and calculate its compliance matrix S. sat It is calculated by the following formula;

[0131]

[0132]

[0133] Among them, S dryS is the compliance matrix of dry shale. sat S is the compliance matrix of saturated fluid shale. m β is the equivalent elastic compliance matrix of the constituent minerals. fl For the compressibility of pore fluids, β m For the compressibility of minerals, K fl Let φ be the bulk modulus of the porous fluid, and φ be the porosity, where i, j, k, l, a, b, c, and d are the coordinate axis indices of the compliance matrix.

[0134] S7. Input the initial values ​​of the physical property parameters, calculate the velocity parameters using the stiffness matrix of fluid-saturated shale, establish the objective function, perturb the auxiliary parameters to minimize the objective function value, and output the auxiliary parameter inversion results. Specifically, this includes the following sub-steps:

[0135] S7.1, Set the aspect ratio α of the connecting pores con The proportion of interconnected pores r con The layering index LI is used as an important auxiliary parameter to be inverted. Initial values ​​for the auxiliary parameter are set to begin the inversion process, and α is obtained. con r con The inversion results of LI, with the inversion objective parameter being α. con r con And LI, the initial values ​​of the auxiliary parameters are set as follows:

[0136]

[0137] Where m0 is the initial auxiliary parameter vector, The initial value of the aspect ratio of the connected pores. LI is the initial value of the pore connectivity index. 0 This is the initial value of the stratification index;

[0138] S7.2. Set the inversion objective function as follows:

[0139]

[0140] Where, m i Let be the auxiliary parameter value for the i-th iteration. This is the prediction result when the objective function error is minimized. The predicted P-wave velocity for the constructed shale model. ρ represents the predicted shear wave velocity from the constructed shale model. pred (m i The density is predicted by the established shale model. The measured P-wave velocity in the well. ρ is the transverse wave velocity measured in the well. realThis represents the density measured in the well. σ1 is the weight value of the P-wave velocity data item, σ2 is the weight value of the S-wave velocity data item, and σ3 is the weight value of the density data item. and All were calculated from the compliance matrix of the shale model;

[0141] S7.2. Update the auxiliary parameters using an iterative method. When the iteration number i = 1, the auxiliary parameters satisfy... When the iteration number i ≥ 2, the physical property parameters satisfy in, The aspect ratio of the connected pores at the i-th iteration is given. Let LI be the pore connectivity index value at the i-th iteration. i This represents the stratification index value at the i-th iteration; when the objective function value is minimized, the corresponding... This is the final inversion result.

[0142] S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is corrected, and the final predicted results of elasticity and anisotropy parameters are output based on the corrected shale petrophysical model. Specifically, this includes the following sub-steps: S8.1. Based on the obtained three auxiliary parameters α... con r con Inversion results with LI The established shale model is modified to better reflect real formation conditions. The stiffness matrix of the saturated fluid shale is then expressed as follows:

[0143]

[0144] Where, m inv The inversion result vector, To connect the aspect ratio inversion results of the pores, For the inversion results of the proportion of interconnected pores, LI inv For the stratification index inversion result, G represents the forward modeling operator for shale petrophysical modeling in steps S1 to S6. To correct the new stiffness matrix predicted by the shale model;

[0145] S8.2. Based on the stiffness matrix of saturated fluid shale, the final elastic and anisotropic parameters are obtained, including the P-wave velocity, S-wave velocity, density, and three Thomsen anisotropic parameters predicted by the modified shale rock physics model.

[0146] The well data used for testing in the embodiments are as follows: Figure 3As shown in the figure, (a) shows the longitudinal wave velocity, (b) shows the transverse wave velocity, (c) shows the density, (d) shows the porosity, and (e) shows the mineral composition. The following modeling tests were conducted for comparison: (1) Using... Figure 2 The shale rock physical modeling method, but LI, r con α con All are set to fixed values, and the prediction results are as follows: Figure 4 As shown in the figure, (a) shows the measured and predicted elasticity parameters, and (b) shows the predicted anisotropy parameters. (2) Using the inversion-assisted modeling method of the present invention, the following parameters are selected: Figure 2 Modeling methods and Figure 1 The prediction process, while inverting LI and r con α con Then, the inversion results are used to correct the shale model before predicting elastic parameters and anisotropy parameters. The prediction results are as follows: Figure 5 As shown. Figure 4 and Figure 5 In the diagram, the black solid line represents measured data, including P-wave velocity, S-wave velocity, and density; the blue dashed line represents elastic parameters predicted by the model; the red solid line represents auxiliary parameters obtained through inversion; and the blue solid line represents anisotropic parameters predicted by the model, including ε, δ, and γ. It can be seen that when the three parameters are set to constant values, such as... Figure 4 The figure shows that the prediction error of the elastic parameters is relatively large; the inversion-assisted modeling method proposed in this invention, such as... Figure 5 In the figure, (a) shows the measured and predicted elastic parameters and the inverted auxiliary parameters, and (b) shows the predicted anisotropy parameters. The accuracy of the predicted elastic results is significantly improved, especially the well consistency of P-wave velocity and S-wave velocity. The test results demonstrate that the modeling method proposed in this invention can better simulate the characteristics of shale rocks and obtain relatively reliable anisotropy parameter simulation results.

[0147] This invention considers the stratification distribution of clay-kerogen in shale reservoirs during rock physics modeling, using a stratification index to characterize the strength of shale stratification. To account for the complex porosity of shale, pores are divided into isolated pores and connected pores. Given that isolated pores in shale typically develop in non-brittle minerals, isolated microcracks are only added to non-brittle minerals such as clay-kerogen during modeling, while connected pores are added to both brittle and non-brittle minerals. The porosity development is simulated by considering the aspect ratio and proportion of the two pore types. Furthermore, by comparing the errors between predicted elastic parameters and actual well data, three important auxiliary parameters—stratification index, connected pore proportion, and connected pore aspect ratio—are retrieved through inversion. The retrieved parameter results are used to correct the established rock physics model, making it more consistent with reality, predicting the elasticity and anisotropy parameters of shale rocks, and improving prediction accuracy.

[0148] While the present invention has been described above with reference to preferred embodiments, it is not intended to limit the invention. Those skilled in the art can make various modifications and refinements without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention shall be determined by the claims.

Claims

1. A method for inverting and assisting shale petrophysical modeling that considers clay-kerogen stratification and pore complexity, characterized in that, Includes the following steps: S1. Based on the clay and kerogen content, establish the clay-kerogen mixture and calculate its stiffness matrix; S2, with a porosity of φ iso Isolated, flat microcracks were added to the clay-kerogen mixture to create clay-kerogen blocks, and their stiffness matrix was calculated. S3. By rotating and stacking clay-kerogen blocks according to a normal distribution, a shale matrix is ​​constructed, and its stiffness matrix is ​​calculated. S4. Based on the content of brittle minerals, establish an isotropic brittle mineral mixture, add the brittle mineral mixture to the shale matrix to obtain the shale skeleton, and calculate its stiffness matrix. S5, with a porosity of φ con Connected dry pores are added to the shale skeleton to obtain dry shale, and its stiffness matrix is ​​calculated. S6. Fill the pores of dry shale with fluid or fluid mixture to obtain fluid-saturated shale, and calculate its stiffness matrix. S7. Input the initial values ​​of the physical property parameters, calculate the velocity parameters using the stiffness matrix of fluid-saturated shale, establish the objective function, perturb the auxiliary parameters to minimize the objective function value, and output the auxiliary parameter inversion results. S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is corrected, and the final prediction results of elasticity and anisotropy parameters are output based on the corrected shale rock physics model.

2. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, In step S1, clay and kerogen are mixed using the anisotropic SCA-DEM method. The stiffness matrix C of the clay-kerogen mixture is... ckm Calculate as follows: Among them, C ckm Let C be the stiffness matrix of the clay-kerogen mixture. n Let C be the stiffness matrix of the nth phase inclusion. p Let G be the stiffness matrix of the p-th phase inclusion body. ckm To calculate the Eshelby matrix corresponding to the clay-kerogen mixture, where I is the identity matrix, n is the nth phase inclusion, and ν n ν is the volume fraction of the nth phase inclusion. ker This represents the volume fraction of kerogen inclusions. The volume fraction of kerogen inclusions was obtained by normalizing the measured mineral volume fractions in the well, as shown in the following formula: In the formula, v ker This represents the normalized kerogen volume fraction, and also the kerogen inclusion volume fraction. The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

3. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S2 includes the following sub-steps: S2.

1. The pore space of shale is divided into two parts: isolated pore space and connected pore space. The proportions of isolated and connected pores are respectively denoted by r. iso r con It can be expressed as follows: r con =φ con / f r iso =φ iso / φ=1-r con Where φ is the total porosity, φ = φ con +φ iso , φ iso For isolated porosity, φ con For connectivity porosity, r iso r represents the proportion of isolated pores. con The proportion of connected pores; S2.

2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to construct clay-kerogen blocks containing isolated pores, whose stiffness matrix C ckb Calculate as follows: Among them, C ckb Here is the stiffness matrix of the clay-kerogen chunks. Let G be the stiffness matrix of the nth isolated pore inclusion. ckb This involves calculating the Eshelby matrix corresponding to clay-kerogen fragments, where I is the identity matrix, n is the nth isolated pore inclusion, and ν... iso The volume fraction of isolated pore inclusions is calculated using porosity and isolated porosity, satisfying the following expression: v iso =r iso =φ iso / φ=1-φ con / φ.

4. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S3 includes the following sub-steps: S3.1 Based on the normal distribution of the rotation angle of clay-kerogen fragments in shale, the orientation distribution function ODF is as follows: ODF~N(mean,LI) Where ODF represents the orientation distribution function, N represents the normal distribution, mean represents the expected value of the rotation angle θ, and LI is the standard deviation of the normal distribution, representing the disorder of the clay-kerogen block deflection, i.e., the stratification index. S3.2 The probability density function of the preferred orientation is as follows: Where υ(θ,LI) represents the probability density function of the preferred orientation, and θ is the rotation angle; S3.3 The expressions for estimating the stiffness and flexibility matrices of the clay layer composed of rotated clay-kerogen blocks using VRH averaging are as follows: Among them, C mix S represents the elastic stiffness matrix of the mixture obtained by the Voigt method. mix This represents the elastic compliance matrix of the mixture obtained by the Reuss method. and These are the stiffness and compliance matrices of the clay-kerogen mixture after Bond transformation, respectively. It is the azimuth angle; S3.4 Calculate the average VRH, as shown in the following expression: C clay =[C mix +(S mix ) -1 ]] / 2 Among them, C clay This represents the stiffness matrix of the clay layer composed of clay-kerogen blocks after rotation.

5. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S4 includes the following sub-steps: S4.

1. Using the VRH averaging method, establish the stiffness matrix C of an isotropic brittle mineral mixture. brit It is calculated by the following formula: in, In the formula, C brit For the stiffness matrix of a brittle mineral mixture, μ brit K is the shear modulus of a brittle mineral mixture. brit The bulk modulus of a brittle mineral mixture. The shear modulus of a brittle mineral mixture obtained by the Voigt method. The shear modulus of the brittle mineral mixture obtained by the Reuss method. The bulk modulus of a brittle mineral mixture obtained by the Voigt method. The bulk modulus of the brittle mineral mixture obtained by the Reuss method; S4.

2. Using the anisotropic DEM method, a brittle mineral mixture is added to the shale matrix formed by the rotational stacking of clay and kerogen to establish the shale skeleton, whose stiffness matrix C... sk It is calculated by the following formula: Among them, C sk Here is the stiffness matrix of the shale matrix. Let G be the stiffness matrix of the nth brittle mineral inclusion. sk This calculates the Eshelby matrix corresponding to the shale matrix, where I is the identity matrix, n is the nth brittle mineral inclusion, and v sk This represents the percentage content of inclusions in brittle mineral mixtures; the percentage content of inclusions in brittle mineral mixtures is obtained by normalizing the volume fraction of brittle minerals, as shown in the following formula: In the formula, The measured volume fraction of kerogen in the well. The measured clay volume fraction in the well. Let be the volume fraction of the i-th brittle mineral measured in the well.

6. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S5 specifically involves: using an anisotropic DEM to determine the porosity φ con Connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry It is calculated by the following formula: Among them, C dry The stiffness matrix of dry shale. Let G be the stiffness matrix of the nth connected pore inclusion. dry This calculates the Eshelby matrix corresponding to dry shale, where I is the identity matrix, n is the nth connected pore inclusion, and v con The percentage content of connected pore inclusions, where the percentage content of connected pore inclusions satisfies the following expression: v con =r con =φ con / φ where φ is the total porosity, φ con For connectivity porosity, r con The percentage of connected pores.

7. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S6 specifically involves: filling the pores of dry shale with fluid or a fluid mixture based on the BK model to obtain fluid-saturated shale, whose compliance matrix S sat It is calculated by the following formula; Among them, S dry S is the compliance matrix of dry shale. sat S is the compliance matrix of saturated fluid shale. m β is the equivalent elastic compliance matrix of the constituent minerals. fl For the compressibility of pore fluids, β m For the compressibility of minerals, K fl Let φ be the bulk modulus of the porous fluid, and φ be the porosity, where i, j, k, l, a, b, c, and d are the coordinate axis indices of the compliance matrix.

8. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S7 includes the following sub-steps: S7.1, Set the aspect ratio α of the connecting pores con The proportion of interconnected pores r con The layering index LI is used as an important auxiliary parameter to be inverted. Initial values ​​for the auxiliary parameter are set to begin the inversion process, and α is obtained. con r con The inversion results of LI, The inversion objective parameter is α con r con And LI, the initial values ​​of the auxiliary parameters are set as follows: Where m0 is the initial auxiliary parameter vector, The initial value of the aspect ratio of the connected pores. LI is the initial value of the pore connectivity index. 0 This is the initial value of the stratification index; S7.

2. Set the inversion objective function as follows: Where, m i Let be the auxiliary parameter value for the i-th iteration. This is the prediction result when the objective function error is minimized. The predicted P-wave velocity for the constructed shale model. ρ represents the predicted shear wave velocity from the constructed shale model. pred (m i The density is predicted by the established shale model. The measured P-wave velocity in the well. ρ is the transverse wave velocity measured in the well. real Here, σ1 represents the density measured in the well, σ2 represents the weight value of the P-wave velocity data item, σ3 represents the weight value of the S-wave velocity data item, and σ4 represents the weight value of the density data item. and All were calculated from the compliance matrix of the shale model; S7.

3. Update the auxiliary parameters using an iterative method. When the iteration number i = 1, the auxiliary parameters satisfy... When the iteration number i ≥ 2, the physical property parameters satisfy in, The aspect ratio of the connected pores at the i-th iteration is given. Let LI be the pore connectivity index value at the i-th iteration. i The stratification index value at the i-th iteration; when the objective function value When it is at its minimum, the corresponding This is the final inversion result.

9. The inversion-assisted shale petrophysical modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that, Step S8 includes the following sub-steps: S8.1, Based on obtaining three auxiliary parameters α con r con Inversion results with LI The established shale model is modified to better reflect real formation conditions. The stiffness matrix of the saturated fluid shale is then expressed as follows: Where, m inv The inversion result vector, To connect the aspect ratio inversion results of the pores, For the inversion results of the proportion of interconnected pores, LI inv For the stratification index inversion result, G represents the forward modeling operator for shale petrophysical modeling in steps S1 to S6. To correct the new stiffness matrix predicted by the shale model; S8.

2. Based on the stiffness matrix of saturated fluid shale, the final elastic and anisotropic parameters are obtained, including the P-wave velocity, S-wave velocity, density, and three Thomsen anisotropic parameters predicted by the modified shale rock physics model.

Images