Inversion-assisted shale rock physical modeling method considering clay-kerogen stratification and pore complexity

Abstract

According to the inversion-assisted shale rock physical modeling method considering clay-kerogen layering and pore complexity, pores are divided into isolated micro-cracks (only related to clay and kerogen) and communicated pores (developed in various minerals), and the proportion of the communicated pores is introduced to balance the influence of the two types of pores; a layering index is introduced to describe the degree of directional arrangement of the particles. And a modeling process is established on the basis of an inclusion theory and a Brown-Korringa model. And three important auxiliary parameters, including a layering index, a communicated pore ratio and a communicated pore aspect ratio, are estimated by comparing an error between a predicted elastic parameter and actually measured well data in combination with inversion. And the established rock physical model is corrected according to the inversion result, so that the rock physical model better conforms to the actual situation, and the prediction precision of the elasticity and anisotropy parameters of the shale rock is improved. Logging data tests verify the effectiveness of the method.

Description

Technical Field

[0001] The invention belongs to the field of shale rock physical modeling, and in particular relates to an inversion-assisted shale rock physical modeling method taking clay-kerogen stratification and pore complexity into consideration. Background Art

[0002] In the early days of oil and gas exploration, conventional reservoirs (such as sandstone and mudstone reservoirs) were generally regarded as isotropic media. With the increase in global oil and gas demand, especially the diversification of energy demand 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 include shale oil, shale gas, tight oil and gas, etc. These resources are usually widely distributed, but the permeability of the reservoirs is low and it is difficult to exploit. However, the physical properties of shale reservoirs are quite different from those of conventional sandstone reservoirs, so new technical means are needed for effective exploration and development. Shale reservoirs often have significant anisotropy, especially VTI (vertically symmetric transverse isotropy) anisotropy. In VTI media, the propagation speeds of longitudinal waves (P waves) and transverse waves (S waves) in different directions are quite different, and the wave speeds vary with different directions. This anisotropic characteristic will cause more complex reflection characteristics during the propagation of seismic waves, making seismic exploration more difficult. The conventional isotropic assumption obviously cannot accurately describe the physical behavior of shale reservoirs. Since it is difficult to obtain the measured data of anisotropic parameters on the well, it is necessary to estimate the anisotropic parameters in the well through the rock physics model for subsequent research. At present, the factors considered in shale rock physics modeling are relatively simple, for example: 1. The clay-rich shale model is constructed by combining the anisotropic differential equivalent medium model (DEM) and the self-compatible model (SCA). 2. The anisotropic Backus average is used to describe the anisotropy caused by the stratification of clay and kerogen in the shale matrix. 3. The standard deviation of the orientation distribution function is used to simulate the stratification of clay-kerogen. 4. The equivalent pore aspect ratio of the rock is inverted according to the longitudinal and shear wave velocities, and the equivalent pore aspect ratio is used to predict the shear wave velocity for rock physics modeling. 5. An inversion scheme that predicts the anisotropic parameters of vertical and horizontal well logging by simultaneously estimating the stratification index and the "clay-related" pore aspect ratio. 6. A rock physics modeling approach for tight sandstone reservoirs takes into account connected and isolated pores and defines a pore connectivity parameter to characterize the overall effect of different pore types on the rock elastic modulus.

[0003] In summary, current shale modeling does not consider the combined effects of mineral particle stratification and complex pore conditions on rock elastic properties. In order to improve the accuracy of shale reservoir exploration, it is necessary to develop rock physics modeling and seismic exploration methods for VTI media. These methods need to be able to fully consider the inherent anisotropic characteristics of shale and the complex pore structure, laying the foundation for subsequent shale seismic prediction. Summary of the invention

[0004] In view of the problems existing in the prior art, the present invention provides an inversion-assisted shale rock physics modeling method taking into account clay-kerogen stratification and pore complexity, which can correct the established rock physics model with the inversion parameter results to make it more consistent with the actual situation, predict the elasticity and anisotropy parameters of shale rock, and improve the prediction accuracy.

[0005] In order to solve the above technical problems, the present invention provides the following technical solutions: an inversion-assisted shale rock physics modeling method considering clay-kerogen stratification and pore complexity, comprising the following steps:

[0006] S1. Based on the clay and kerogen contents, a clay-kerogen mixture is established and its stiffness matrix is ​​calculated;

[0007] S2, the porosity is φ iso Isolated flat microcracks are added into the clay-kerogen mixture to build clay-kerogen blocks and calculate their stiffness matrix.

[0008] S3, constructing the shale matrix by rotating and superimposing the clay-kerogen blocks according to normal distribution, and calculating its stiffness matrix;

[0009] S4. Based on the brittle mineral content, an isotropic brittle mineral mixture is established, the brittle mineral mixture is added to the shale matrix to obtain a shale skeleton, and its stiffness matrix is ​​calculated;

[0010] S5, the porosity is φ con The connected dry pores are added to the shale skeleton to obtain dry shale, and its stiffness matrix is ​​calculated;

[0011] S6. Filling the fluid or fluid mixture into the pores of the dry shale to obtain the fluid-saturated shale, and calculating its stiffness matrix;

[0012] S7, inputting the initial values ​​of the physical property parameters, calculating the velocity parameters using the stiffness matrix of the fluid-saturated shale, establishing the objective function, perturbing the auxiliary parameters to minimize the objective function value, and outputting the auxiliary parameter inversion results;

[0013] S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is modified, and the final prediction results of elasticity and anisotropy parameters are output based on the modified 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 The calculation is as follows:

[0015]

[0016]

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

[0018] The volume fraction of kerogen inclusions is obtained by normalizing the measured mineral volume fraction in the well, as follows:

[0019]

[0020] In the formula, v ker is the normalized kerogen volume fraction, also the kerogen inclusion volume fraction; is the measured kerogen volume fraction in the well, is the measured clay volume fraction in the well, is the volume fraction of the i-th brittle mineral measured in the well.

[0021] Furthermore, the aforementioned step S2 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 expressed as r iso 、r con It is expressed as follows:

[0023] r con =φ con / φ

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

[0025] Where φ is the total porosity, φ = φ con +φ iso ,φ iso is the isolated porosity, φ con is the connected porosity, r iso is the isolated pore ratio, r con is the proportion of connected pores;

[0026] S2.2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to create a clay-kerogen block containing isolated pores. Its stiffness matrix C ckb The calculation is as follows:

[0027]

[0028] Among them, C ckb is the stiffness matrix of the clay-kerogen patch, is the stiffness matrix of the nth isolated pore inclusion, G ckb is the Eshelby matrix corresponding to the clay-kerogen block, I is the identity matrix, n is the nth isolated pore inclusion, ν iso is the volume fraction of isolated pore inclusions. The volume fraction of isolated pore inclusions is calculated using porosity and isolated porosity, satisfying the following expression: 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 angles of clay-kerogen blocks 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 expectation of the rotation angle θ; LI is the standard deviation of the normal distribution, which represents the disorder of the deflection of the clay-kerogen block, that is, 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 stiffness and flexibility matrix expressions of the clay layer composed of rotated clay-kerogen blocks are estimated by using VRH average as follows:

[0037]

[0038]

[0039] Among them, C mix represents the elastic stiffness matrix of the mixture obtained by the Voigt method, Smix represents the elastic flexibility matrix of the mixture obtained by the Reuss method, and are the stiffness matrix and flexibility matrix of the clay-kerogen mixture after Bond transformation, is the azimuth;

[0040] S3.4. Calculate the VRH average, the expression is as follows:

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

[0042] Among them, C clay Represents the stiffness matrix of the clay layer composed of clay-kerogen chunks after rotation.

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

[0044] S4.1. Using the VRH average method, an isotropic brittle mineral mixture is established, and its stiffness matrix C brit Calculated by the following formula:

[0045]

[0046] in,

[0047]

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

[0049] S4.2. Using the anisotropic DEM method, the brittle mineral mixture is added to the shale matrix formed by the rotation and superposition of clay-kerogen to establish a shale skeleton, and its stiffness matrix C sk Calculated by the following formula:

[0050]

[0051] Among them, C skis the stiffness matrix of the shale matrix, is the stiffness matrix of the nth brittle mineral mixture inclusion, G sk is the Eshelby matrix corresponding to the shale matrix, I is the unit matrix, n is the nth brittle mineral mixture inclusion, v sk is 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 follows:

[0052]

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

[0054] Furthermore, the aforementioned step S5 is specifically as follows: using anisotropic DEM to convert the porosity to con The connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry Calculated by the following formula:

[0055]

[0056] Among them, C dry is the stiffness matrix of dry shale, is the stiffness matrix of the nth connected pore inclusion, G dry is the Eshelby matrix corresponding to dry shale, I is the unit matrix, n is the nth connected pore inclusion, v con is the percentage of interconnected pore inclusions, where the percentage of interconnected pore inclusions satisfies the following expression: con =r con =φ con / φ

[0057] Where, φ is the total porosity, φ con is the connected porosity, r con is the proportion of connected pores.

[0058] Furthermore, the aforementioned step S6 is specifically as follows: based on the BK model, the fluid or fluid mixture is filled into the pores of the dry shale to obtain the fluid-saturated shale, whose flexibility matrix S sat It is calculated by the following formula;

[0059]

[0060]

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

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

[0063] S7.1. Setting the interconnected pore aspect ratio α con , proportion of interconnected pores con The stratification index LI is used as an important auxiliary parameter to be inverted. The initial value of the auxiliary parameter is set to start the inversion and obtain α con 、r con The inversion results of LI and the inversion target parameter is α con 、r con and LI, the initial values ​​of auxiliary parameters are set as:

[0064]

[0065] Among them, m0 is the initial auxiliary parameter vector, is the initial value of the connected pore aspect ratio, is the initial value of the pore connectivity index, LI 0 is the initial value of the stratification index;

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

[0067]

[0068] Among them, m i is the auxiliary parameter value of the i-th iteration, is the prediction result when the objective function error is the smallest, is the P-wave velocity predicted by the constructed shale model, is the shear wave velocity predicted by the shale model, ρ pred (m i ) is the density predicted by the shale model, is the measured P-wave velocity in the well, is the shear wave velocity measured in the well, ρ real is the measured density 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, where: and All are calculated from the shale model flexibility matrix;

[0069] S7.3, update the auxiliary parameters in an iterative manner. When the number of iterations i = 1, the auxiliary parameters satisfy When the number of iterations i ≥ 2, the physical property parameters satisfy in, is the connected pore aspect ratio at the i-th iteration, is the pore connectivity index value at the i-th iteration, LI i is the layering index value at the i-th iteration; when the objective function value At the minimum, the corresponding The final inversion result is

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

[0071] S8.1, based on the three auxiliary parameters α con 、r con And the inversion results of LI The established shale model is modified to make it closer to the actual formation conditions. At this time, the stiffness matrix of saturated fluid shale is expressed as follows:

[0072]

[0073] Among them, m inv is the inversion result vector, To connect the pore aspect ratio inversion results, is the inversion result of the pore ratio, LI inv is the inversion result of the stratification index, G represents the shale rock physics modeling forward operator from step S1 to step S6, New stiffness matrix predicted for the corrected shale model;

[0074] S8.2. Based on the stiffness matrix of the fluid-saturated 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 physics modeling method of the present invention takes into account the stratification of clay-kerogen and introduces the stratification index LI to simulate the deflection disorder of clay-kerogen blocks. It also takes into account the pore complexity, divides the pores into isolated pores and connected pores, and introduces the pore connectivity index r conTo characterize the influence of the two pore types, the constructed model is more consistent with the actual situation of the formation.

[0077] (2) Minimize the residual between the predicted values ​​of P- and S-wave velocities and the measured data, and iteratively invert the three auxiliary parameters α con 、r con The inversion-assisted modeling strategy was applied to actual logging data, and the P-wave and S-wave velocity errors predicted by the rock physics model with fixed auxiliary parameters and auxiliary parameters corrected with depth were compared. The results showed that the modeling method that considers three auxiliary parameters at the same time has the highest prediction accuracy and the best consistency with the measured data. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 Flowchart for shale rock physics modeling to aid inversion.

[0079] Figure 2 Flowchart for modeling shale rock physics considering clay-kerogen stratification and isolated / connected pores.

[0080] Figure 3 These are the measured data of the test well, where (a) is a schematic diagram of the longitudinal wave velocity, (b) is a schematic diagram of the shear wave velocity, (c) is a schematic diagram of the density, (d) is a schematic diagram of the porosity, and (e) is a schematic diagram of the mineral composition.

[0081] Figure 4 Schematic diagram of the prediction results when the three auxiliary parameters are all 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 anisotropic parameters.

[0082] Figure 5 This is a schematic diagram of the prediction results of the inversion-assisted rock physics modeling of the present invention, in which (a) is a schematic diagram of the measured and predicted elastic parameters and the auxiliary parameters of the inversion, and (b) is a schematic diagram of the predicted anisotropic parameters. DETAILED DESCRIPTION

[0083] In order to better understand the technical content of the present invention, specific embodiments are given and described as follows in conjunction with the accompanying drawings.

[0084] Various aspects of the invention are described herein with reference to the accompanying drawings, in which many illustrative embodiments are shown. The embodiments of the invention are not limited to those described in the accompanying drawings. It should be understood that the invention is implemented by 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 in the invention are not limited to any implementation. In addition, some aspects disclosed in the invention may be used alone or in any appropriate combination with other aspects disclosed in the invention.

[0085] The embodiment is described by taking a modeling test of real well logging data as an example. Figure 1 , Figure 2 As shown, the present invention discloses an inversion-assisted shale rock physics modeling method considering clay-kerogen stratification and pore complexity, comprising the following steps:

[0086] S1. Based on the clay and kerogen content, a clay-kerogen mixture is established and its stiffness matrix is ​​calculated. The details are as follows:

[0087] The anisotropic SCA-DEM method is used to mix clay and kerogen. The stiffness matrix C of the clay-kerogen mixture is ckm The calculation is as follows:

[0088]

[0089]

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

[0091] The volume fraction of kerogen inclusions is obtained by normalizing the measured mineral volume fraction in the well, as follows:

[0092]

[0093] In the formula, v ker is the normalized kerogen volume fraction, also the kerogen inclusion volume fraction; is the measured kerogen volume fraction in the well, is the measured clay volume fraction in the well, is the volume fraction of the i-th brittle mineral measured in the well.

[0094] S2, the porosity is φ iso The isolated flat microcracks are added to the clay-kerogen mixture to establish a clay-kerogen block and calculate its stiffness matrix. 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 expressed as r iso、r con It is expressed as follows:

[0096] r con =φ con / φ

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

[0098] Where φ is the total porosity, φ = φ con +φ iso ,φ iso is the isolated porosity, φ con is the connected porosity, r iso is the isolated pore ratio, r con is the proportion of connected pores.

[0099] S2.2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to create a clay-kerogen block containing isolated pores. Its stiffness matrix C ckb The calculation is as follows:

[0100]

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

[0102] S3. By rotating and superimposing clay-kerogen blocks according to normal distribution, a shale matrix is ​​constructed and its stiffness matrix is ​​calculated. The mutually coupled clay-kerogen in shale will enhance the anisotropy of shale. In order to analyze the influence of clay stratification in shale, the clay layer is regarded as a complex whole, which is composed of many clay-kerogen blocks with the same properties. The rotation angle θ of each clay-kerogen block, that is, the angle between the symmetry axis of the clay-kerogen block 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 Bond transformation on the horizontal blocks. Afterwards, the equivalent elastic properties of the superimposed clay-kerogen blocks are calculated by using VRH averaging. Specifically, it includes the following sub-steps:

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

[0104] ODF~N(mean,LI)

[0105] Among them, ODF represents the orientation distribution function, N represents the normal distribution, mean represents the expectation of the rotation angle θ, and because the main stratification direction of the stratum is horizontal, mean is usually set to zero. LI is the standard deviation of the normal distribution, which represents the disorder of the deflection of the clay-kerogen block, that is, the stratification index;

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

[0107]

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

[0109] S3.3. The stiffness and flexibility matrix expressions of the clay layer composed of rotated clay-kerogen blocks are estimated by using VRH average as follows:

[0110]

[0111]

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

[0113] S3.4. Calculate the VRH average, the expression is as follows:

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

[0115] Among them, C clay Represents the stiffness matrix of the clay layer composed of clay-kerogen chunks after rotation.

[0116] S4. Based on the brittle mineral content, an isotropic brittle mineral mixture is established, the brittle mineral mixture is added to the shale matrix to obtain a shale skeleton, and its stiffness matrix is ​​calculated.

[0117] S4.1. Using the VRH average method, an isotropic brittle mineral mixture is established, and its stiffness matrix C brit Calculated by the following formula:

[0118]

[0119] in,

[0120]

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

[0122] S4.2. Using the anisotropic DEM method, the brittle mineral mixture is added to the shale matrix formed by the rotation and superposition of clay-kerogen to establish a shale skeleton, and its stiffness matrix C sk Calculated by the following formula:

[0123]

[0124] Among them, C sk is the stiffness matrix of the shale matrix, is the stiffness matrix of the nth brittle mineral mixture inclusion, G skis the Eshelby matrix corresponding to the shale matrix, I is the unit matrix, n is the nth brittle mineral mixture inclusion, v sk is 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 follows:

[0125]

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

[0127] S5, the porosity is φ con The connected dry pores of φ are added to the shale skeleton to obtain dry shale and calculate its stiffness matrix. Specifically, the anisotropic DEM is used to add the porosity of φ con The connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry Calculated by the following formula:

[0128]

[0129] Among them, C dry is the stiffness matrix of dry shale, is the stiffness matrix of the nth connected pore inclusion, G dry is the Eshelby matrix corresponding to dry shale, I is the unit matrix, n is the nth connected pore inclusion, v con is the percentage of interconnected pore inclusions, where the percentage of interconnected pore inclusions satisfies the following expression: con =r con =φ con / φ, where φ is the total porosity, φ con is the connected porosity, r con is the proportion of connected pores.

[0130] S6. Fill the pores of dry shale with 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 fluid or fluid mixture to obtain fluid-saturated shale, and its flexibility matrix S sat It is calculated by the following formula;

[0131]

[0132]

[0133] Among them, S dryis the flexibility matrix of dry shale, S sat is the flexibility matrix of saturated fluid shale, S m is the equivalent elastic flexibility matrix of the constituent minerals, β fl is the compressibility of the pore fluid, β m is the compressibility of the mineral, K fl is the bulk modulus of the pore fluid, φ is the porosity, and i, j, k, l, a, b, c, and d are the coordinate axis direction indices of the flexibility matrix.

[0134] S7, input the initial values ​​of the physical property parameters, calculate the velocity parameters using the stiffness matrix of the 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, it includes the following sub-steps:

[0135] S7.1. Setting the interconnected pore aspect ratio α con , proportion of interconnected pores con The stratification index LI is used as an important auxiliary parameter to be inverted. The initial value of the auxiliary parameter is set to start the inversion and obtain α con 、r con The inversion results of LI and the inversion target parameter is α con 、r con and LI, the initial values ​​of auxiliary parameters are set as:

[0136]

[0137] Among them, m0 is the initial auxiliary parameter vector, is the initial value of the connected pore aspect ratio, is the initial value of the pore connectivity index, LI 0 is the initial value of the stratification index;

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

[0139]

[0140] Among them, m i is the auxiliary parameter value of the i-th iteration, is the prediction result when the objective function error is the smallest, is the P-wave velocity predicted by the constructed shale model, is the shear wave velocity predicted by the shale model, ρ pred (m i ) is the density predicted by the shale model, is the measured P-wave velocity in the well, is the shear wave velocity measured in the well, ρ realis the measured density 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, where: and All are calculated from the shale model flexibility matrix;

[0141] S7.2, update the auxiliary parameters in an iterative manner. When the number of iterations i = 1, the auxiliary parameters satisfy When the number of iterations i ≥ 2, the physical property parameters satisfy in, is the connected pore aspect ratio at the i-th iteration, is the pore connectivity index value at the i-th iteration, LI i is the layering index value at the i-th iteration; when the objective function value is the minimum, the corresponding The final inversion result is

[0142] S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is modified, and the final prediction results of elasticity and anisotropy parameters are output based on the modified shale rock physics model. Specifically, it includes the following sub-steps: S8.1. Based on the three auxiliary parameters α con 、r con And the inversion results of LI The established shale model is modified to make it closer to the actual formation conditions. At this time, the stiffness matrix of saturated fluid shale is expressed as follows:

[0143]

[0144] Among them, m inv is the inversion result vector, To connect the pore aspect ratio inversion results, is the inversion result of the pore ratio, LI inv is the inversion result of the stratification index, G represents the shale rock physics modeling forward operator from step S1 to step S6, New stiffness matrix predicted for the corrected shale model;

[0145] S8.2. Based on the stiffness matrix of the fluid-saturated 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 embodiment are as follows: Figure 3As shown in the figure, (a) shows the schematic diagram of longitudinal wave velocity, (b) shows the shear wave velocity, (c) shows the density, (d) shows the porosity, and (e) shows the mineral composition. The following modeling tests were carried out for comparison: (1) Using Figure 2 The shale rock physics model method, but LI, r con , α con are all set to fixed values, and the prediction results are as follows Figure 4 In the figure, (a) shows the measured and predicted elastic parameters, and (b) shows the predicted anisotropic parameters. (2) Using the inversion-assisted modeling method of the present invention, select Figure 2 The modeling method and Figure 1 The prediction process of LI and r con , α con Then the shale model is corrected with the inversion results to predict the elastic parameters and anisotropic parameters. The prediction results are as follows: Figure 5 shown. Figure 4 and Figure 5 In the figure, the black solid line represents the measured data, including P-wave velocity, S-wave velocity, and density; the blue dotted line represents the elastic parameters predicted by the model; the red solid line represents the auxiliary parameters obtained by inversion; and the blue solid line represents the 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 elastic parameter prediction error is large; the inversion assisted modeling method proposed in the present invention, such as Figure 5 In the figure, (a) shows the measured and predicted elastic parameters and the auxiliary parameters of the inversion, and (b) shows the predicted anisotropic parameters. The accuracy of the predicted elastic results is significantly improved, especially the well consistency of the P-wave velocity and S-wave velocity is significantly improved. The test shows that the modeling method proposed in the present invention can better simulate the characteristics of shale rocks and obtain more reliable anisotropic parameter simulation results.

[0147] The present invention considers the layered distribution of clay-kerogen in shale reservoirs in rock physics modeling, and uses the layering index to characterize the strength of shale layering; in order to consider the complex pore conditions of shale, the pores are divided into isolated pores and connected pores. Considering that isolated pores in shale are usually developed in non-brittle minerals, isolated microcracks are only added to non-brittle minerals such as clay-kerogen during modeling, and connected pores are added to both brittle and non-brittle minerals. The pore development conditions are simulated by considering the aspect ratio and proportion of the two pore types. In combination with inversion, three important auxiliary parameters are inverted by comparing the errors between the predicted elastic parameters and the measured well data, including the layering index, the proportion of connected pores, and the aspect ratio of connected pores. The established rock physics model is corrected with the inversion parameter results to make it more consistent with the actual situation, predict the elasticity and anisotropy parameters of shale rocks, and improve the prediction accuracy.

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

Claims

1. An inversion-assisted shale rock physics modeling method considering clay-kerogen stratification and pore complexity, characterized in that: The following steps are involved: S1. Based on the clay and kerogen contents, a clay-kerogen mixture is established and its stiffness matrix is ​​calculated; S2, the porosity is φ iso Isolated flat microcracks are added into the clay-kerogen mixture to build clay-kerogen blocks and calculate their stiffness matrix. S3, constructing the shale matrix by rotating and superimposing the clay-kerogen blocks according to normal distribution, and calculating its stiffness matrix; S4. Based on the brittle mineral content, an isotropic brittle mineral mixture is established, the brittle mineral mixture is added to the shale matrix to obtain a shale skeleton, and its stiffness matrix is ​​calculated; S5, the porosity is φ con The connected dry pores are added to the shale skeleton to obtain dry shale, and its stiffness matrix is ​​calculated; S6. Filling the fluid or fluid mixture into the pores of the dry shale to obtain the fluid-saturated shale, and calculating its stiffness matrix; S7, inputting the initial values ​​of the physical property parameters, calculating the velocity parameters using the stiffness matrix of the fluid-saturated shale, establishing the objective function, perturbing the auxiliary parameters to minimize the objective function value, and outputting the auxiliary parameter inversion results; S8. After obtaining the inversion results of the auxiliary parameters, the established shale model is modified, and the final prediction results of elasticity and anisotropy parameters are output based on the modified shale rock physics model.

2. An inversion-assisted shale rock physics 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 The calculation is as follows: Among them, C ckm is the stiffness matrix of the clay-kerogen mixture, C n is the stiffness matrix of the n-th phase inclusion, C p is the stiffness matrix of the p-th phase inclusion, G ckm To calculate the Eshelby matrix corresponding to the clay-kerogen mixture, I is the identity matrix, n is the nth phase inclusion, ν n is the volume fraction of the nth phase inclusion, ν ker is the volume fraction of kerogen inclusions; The volume fraction of kerogen inclusions is obtained by normalizing the measured mineral volume fraction in the well, as follows: In the formula, v ker is the normalized kerogen volume fraction, also the kerogen inclusion volume fraction; is the measured kerogen volume fraction in the well, is the measured clay volume fraction in the well, is the volume fraction of the i-th brittle mineral measured in the well.

3. The inversion-assisted shale rock physics 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 expressed as r iso 、r con It is expressed as follows: r con =φ con / f r iso =φ iso / φ=1-r con Where φ is the total porosity, φ = φ con +φ iso ,φ iso is the isolated porosity, φ con is the connected porosity, r iso is the isolated pore ratio, r con is the proportion of connected pores; S2.

2. Using the anisotropic DEM method, isolated pores are added to the clay-kerogen to create a clay-kerogen block containing isolated pores. Its stiffness matrix C ckb The calculation is as follows: Among them, C ckb is the stiffness matrix of the clay-kerogen patch, is the stiffness matrix of the nth isolated pore inclusion, G ckb is the Eshelby matrix corresponding to the clay-kerogen block, I is the identity matrix, n is the nth isolated pore inclusion, ν iso is the volume fraction of isolated pore inclusions. The volume fraction of isolated pore inclusions is calculated using porosity and isolated porosity, satisfying the following expression: iso =r iso =φ iso / φ=1-φ con / φ.

4. The inversion-assisted shale rock physics 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 angles of clay-kerogen blocks 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 expectation of the rotation angle θ; LI is the standard deviation of the normal distribution, which represents the disorder of the deflection of the clay-kerogen block, that is, 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 stiffness and flexibility matrix expressions of the clay layer composed of rotated clay-kerogen blocks are estimated by using VRH average as follows: Among them, C mix represents the elastic stiffness matrix of the mixture obtained by the Voigt method, S mix represents the elastic flexibility matrix of the mixture obtained by the Reuss method, and are the stiffness matrix and flexibility matrix of the clay-kerogen mixture after Bond transformation, is the azimuth; S3.

4. Calculate the VRH average, the expression is as follows: C clay =[C mix +(S mix ) -1 ]] / 2 Among them, C clay Represents the stiffness matrix of the clay layer composed of clay-kerogen chunks after rotation.

5. The inversion-assisted shale rock physics 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 average method, an isotropic brittle mineral mixture is established, and its stiffness matrix C brit Calculated by the following formula: in, In the formula, C brit is the stiffness matrix of the brittle mineral mixture, μ brit is the shear modulus of the brittle mineral mixture, K brit is the bulk modulus of the brittle mineral mixture, is the shear modulus of brittle mineral mixture obtained by Voigt method, is the shear modulus of the brittle mineral mixture obtained by the Reuss method, is the bulk modulus of brittle mineral mixture obtained by Voigt method, is the bulk modulus of the brittle mineral mixture obtained by the Reuss method; S4.

2. Using the anisotropic DEM method, the brittle mineral mixture is added to the shale matrix formed by the rotation and superposition of clay-kerogen to establish a shale skeleton, and its stiffness matrix C sk Calculated by the following formula: Among them, C sk is the stiffness matrix of the shale matrix, is the stiffness matrix of the nth brittle mineral mixture inclusion, G sk is the Eshelby matrix corresponding to the shale matrix, I is the unit matrix, n is the nth brittle mineral mixture inclusion, v sk is 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 follows: In the formula, is the measured kerogen volume fraction in the well, is the measured clay volume fraction in the well, is the volume fraction of the i-th brittle mineral measured in the well.

6. The inversion-assisted shale rock physics modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that: Step S5 is specifically as follows: using anisotropic DEM to convert the porosity to con The connected dry pores are added to the shale skeleton to obtain the shale dry skeleton, whose stiffness matrix C dry Calculated by the following formula: Among them, C dry is the stiffness matrix of dry shale, is the stiffness matrix of the nth connected pore inclusion, G dry is the Eshelby matrix corresponding to dry shale, I is the unit matrix, n is the nth connected pore inclusion, v con is the percentage of interconnected pore inclusions, where the percentage of interconnected pore inclusions satisfies the following expression: con =r con =φ con / φwhere φ is the total porosity, φ con is the connected porosity, r con is the proportion of connected pores.

7. An inversion-assisted shale rock physics modeling method considering clay-kerogen stratification and pore complexity according to claim 1, characterized in that: Step S6 is as follows: based on the BK model, the fluid or fluid mixture is filled into the pores of the dry shale to obtain the fluid-saturated shale, whose flexibility matrix S sat It is calculated by the following formula; Among them, S dry is the flexibility matrix of dry shale, S sat is the flexibility matrix of saturated fluid shale, S m is the equivalent elastic flexibility matrix of the constituent minerals, β fl is the compressibility of the pore fluid, β m is the compressibility of the mineral, K fl is the bulk modulus of the pore fluid, φ is the porosity, and i, j, k, l, a, b, c, and d are the coordinate axis direction indices of the flexibility matrix.

8. An inversion-assisted shale rock physics 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. Setting the interconnected pore aspect ratio α con , proportion of interconnected pores con The stratification index LI is used as an important auxiliary parameter to be inverted. The initial value of the auxiliary parameter is set to start the inversion and obtain α con 、r con And the inversion results of LI, The inversion target parameter is α con 、r con and LI, the initial values ​​of auxiliary parameters are set as: Among them, m0 is the initial auxiliary parameter vector, is the initial value of the connected pore aspect ratio, is the initial value of the pore connectivity index, LI 0 is the initial value of the stratification index; S7.

2. Set the inversion objective function to: Among them, m i is the auxiliary parameter value of the i-th iteration, is the prediction result when the objective function error is the smallest, is the P-wave velocity predicted by the constructed shale model, is the shear wave velocity predicted by the shale model, ρ pred (m i ) is the density predicted by the shale model, is the measured P-wave velocity in the well, is the shear wave velocity measured in the well, ρ real is the measured density 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, where: and All are calculated from the shale model flexibility matrix; S7.3, update the auxiliary parameters in an iterative manner. When the number of iterations i = 1, the auxiliary parameters satisfy When the number of iterations i ≥ 2, the physical property parameters satisfy in, is the connected pore aspect ratio at the i-th iteration, is the pore connectivity index value at the i-th iteration, LI i is the layering index value at the i-th iteration; when the objective function value At the minimum, the corresponding The final inversion result is 9. An inversion-assisted shale rock physics 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 the three auxiliary parameters α con 、r con And the inversion results of LI The established shale model is modified to make it closer to the actual formation conditions. At this time, the stiffness matrix of saturated fluid shale is expressed as follows: Among them, m inv is the inversion result vector, To connect the pore aspect ratio inversion results, is the inversion result of the pore ratio, LI inv is the inversion result of the stratification index, G represents the shale rock physics modeling forward operator from step S1 to step S6, New stiffness matrix predicted for the corrected shale model; S8.

2. Based on the stiffness matrix of the fluid-saturated 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