Shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction

Abstract

The invention provides a shale fracture attribute inversion method based on clay matrix dominant azimuth distribution correction, which utilizes an azimuth distribution function to describe the influence of directional arrangement of a clay matrix part on rock anisotropy, and meanwhile, utilizes an improved anisotropic matrix Eshelby-Cheng model to simulate the influence of a horizontal fracture. According to the method, a shale rock physical model considering directional arrangement and horizontal cracks of a clay matrix part at the same time is constructed, shale reservoir crack attributes and anisotropy parameters are estimated in combination with logging information, and in the inversion process, a simulated annealing particle swarm optimization (SA-PSO) is used for improving the solving capacity of a non-linear objective function. According to the method, the influence of anisotropy caused by directional arrangement of the clay matrix on fracture attribute inversion can be reduced, the accuracy of shale reservoir fracture attribute estimation and anisotropy parameter inversion is improved, and development of a shale gas resource exploration and development technology is assisted.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration and development technology, and more specifically, to a method for inverting shale fracture properties based on clay matrix dominance orientation distribution correction. Background Technology

[0002] Shale reservoirs are important carriers of unconventional oil and gas. Fractures are key structural elements of oil and gas reservoirs, playing a crucial role in shale oil and gas resource development as important reservoir spaces and seepage channels. Accurately quantifying fracture properties such as fracture density and aspect ratio is essential for the evaluation and development of shale gas reservoirs. Shale exhibits significant geophysical anisotropy, predominantly characterized by transverse isotropy (VTI) along the vertical axis of symmetry. This anisotropy mainly originates from two aspects: (a) inherent anisotropy: caused by the directional arrangement and platy morphology of clay minerals and kerogen in shale; (b) induced anisotropy: mainly caused by the presence of horizontally or near-horizontally developed fractures.

[0003] Existing fracture-type petrophysical models have serious shortcomings. These include non-contact approximation models, anisotropic self-compatible approximation models, the Hudson model, and the Eshelby-Cheng model, which typically assume the rock matrix is ​​isotropic during construction. However, in actual shale formations (especially in typical marine shale reservoirs such as the Sichuan Basin in China), the matrix exhibits significant inherent anisotropy due to the directional arrangement of clay minerals. Ignoring this inherent weak anisotropy leads to the erroneous attribution of the anisotropy contributed by the directional arrangement of the clay matrix to fractures during fracture attribute inversion, thus significantly underestimating the accuracy of key attribute parameters such as fracture density. Existing research has insufficient modeling of the combined effects of the directional arrangement of the clay matrix (inherent anisotropy) and horizontal fractures (induced anisotropy) on the macroscopic anisotropy of shale. Numerous problems also exist in the inversion process. Fracture attribute inversion is usually based on petrophysical models and measured data (such as well logging velocities), and the complexity of these models leads to difficulties in initialization, multi-parameter nonlinear inversion, and multiple solutions. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention provides an inversion method that can correct for the influence of the dominant orientation distribution of clay matrix and improve the accuracy of inversion of fracture property parameters in shale reservoirs.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction includes the following steps:

[0007] S1. Constructing a rock physics model:

[0008] S11. By combining the anisotropic SCA model and the anisotropic DEM model, a transversely isotropic (VTI) medium with a completely oriented vertical symmetry axis formed by the mixture of clay minerals and kerogen is simulated as the background medium.

[0009] S12. Based on the azimuth distribution function, using the standard deviation of the Gaussian distribution. The degree of directional alignment of the background medium is quantified, combined with the Legendre coefficient. and Adjust the equivalent elastic stiffness tensor of background medium one to obtain a partially oriented clay matrix model, which serves as background medium two.

[0010] S13. The brittle mineral mixture is mixed using an isotropic SCA model to simulate the brittle mineral mixture. The brittle mineral mixture and matrix pores are then gradually added to the background medium II using an anisotropic DEM model to construct an anisotropic matrix model of shale.

[0011] S14. Horizontal fractures were added to the shale anisotropic matrix model using the improved Eshelby-Cheng anisotropic matrix model, and correction terms were calculated using the equivalent anisotropic matrix. Thus, a shale skeleton model was obtained.

[0012] S15. Using Wood's formula, the pore fluid is mixed to simulate the fluid mixture. Combined with the Brown-Korringa anisotropic fluid substitution theory, the fluid is filled into the shale skeleton model to obtain a fluid-saturated shale rock physics model.

[0013] S2. Establish the inversion objective function: Using the measured P-wave velocity and S-wave velocity in the well as dual constraints, define the objective function as follows:

[0014] .

[0015] in, It is the measured longitudinal wave velocity in the well. It is the measured longitudinal wave velocity in the well. It is the longitudinal wave velocity simulated by the model. It is the shear wave velocity simulated by the model. The degree of directional alignment of clay matrix units, It is the aspect ratio of the crack. It refers to the crack density.

[0016] S3. Execution parameter inversion: The objective function is solved by using simulated annealing optimized particle swarm optimization algorithm, and the crack attribute parameters and anisotropy parameters are output.

[0017] In one embodiment of the present invention, in step S12, the function of the Gaussian distribution is defined as:

[0018] .

[0019] in, Normalization equations The solution is obtained.

[0020] In one embodiment of the present invention, in step S12, when the standard deviation of the azimuth distribution function is 0°, it can be said that the clay matrix unit is completely oriented, and its corresponding Legendre coefficient and Maximum, respectively and It exhibits the greatest degree of anisotropy; when the standard deviation of the azimuthal distribution function is 90°, it indicates that the clay matrix elements are completely randomly distributed, and its corresponding Legendre coefficient is... and The minimum values ​​are all 0, exhibiting isotropic characteristics.

[0021] In one embodiment of the present invention, in step S13, the brittle mineral includes one or more of quartz, calcite and pyrite.

[0022] In one embodiment of the present invention, in step S14, the equivalent stiffness matrix expression of the improved anisotropic matrix Eshelby-Cheng model is:

[0023] .

[0024] in, It is the rigid matrix of the anisotropic matrix. These are the improved corrections. It refers to crack porosity.

[0025] In one embodiment of the present invention, the improved modification item The required matrix modulus is given by the following formula:

[0026] .

[0027] in, , , , , It is the density of the matrix. and It is the velocity of an equivalent isotropic matrix.

[0028] In one embodiment of the present invention, the relationship between fracture density and fracture longitudinal and transverse directions in the rock physics model is as follows:

[0029] .

[0030] in, It is the crack density. It is the aspect ratio of the crack. It refers to crack porosity.

[0031] In rock physics models, reservoir pore space consists of fractures and spherical pores, with the following relationship:

[0032] .

[0033] in, It is the total porosity of the reservoir. It refers to the porosity of spherical pores. It refers to crack porosity.

[0034] In one embodiment of the present invention, step S3 of the annealing-optimized particle swarm optimization algorithm includes the following steps:

[0035] 31) Initialize the particles.

[0036] 32) Calculate particle fitness.

[0037] 33) Update the temporary position of the particle.

[0038] 34) Determine if the particle movement rules are met. If yes, update the particle position and velocity and proceed to step 35. If no, restore the particle's original position and velocity and then proceed to step 35.

[0039] 35) Cooling operation.

[0040] 36) Determine if the termination condition is met. If yes, output the optimal particle position; otherwise, return to step 32.

[0041] The beneficial effects of this invention are as follows:

[0042] (1) Accurately separate inherent anisotropy from induced anisotropy: By introducing the orientation distribution function (ODF), the dominant distribution of clay minerals and kerogen (the degree of directional arrangement of clay matrix units) is clearly quantified. This study contributes to the inherent anisotropy of shale, avoids its influence being incorrectly attributed to fractures, and solves the problem of conventional methods underestimating fracture property parameters such as fracture density.

[0043] (2) Accurate simulation of the impact of fractures on anisotropy: The innovative improved Eshelby-Cheng model, by using the equivalent isotropic modulus derived from the anisotropic matrix VTI property to calculate the fracture correction term, can more accurately describe the contribution of horizontal fractures to the macroscopic anisotropy of shale matrix with inherent anisotropy.

[0044] (3) The model has clear physical meaning and strong applicability: It combines the advantages of anisotropic SCA model (connectivity problem), anisotropic DEM model (mixed sequence adaptation), orientation distribution function (directional randomness), and improved Eshelby-Cheng (crack response of weak anisotropic matrix) to construct a rock physics model that is more consistent with the actual microstructure of shale (oriented clay matrix + horizontal cracks + spherical pores + mineral mixing).

[0045] (4) Significantly improves inversion accuracy: When applied to actual data (such as the Longmaxi Formation shale of Well A in the Jiaoshiba area), the fracture attributes inverted by this invention ( , ), the degree of directional alignment σ of the clay matrix units and the anisotropy parameters ( The results are stable, have clear physical meaning, and correspond well with actual geological features (such as layered structure, organic matter distribution, and mineral content variations).

[0046] (5) Supporting exploration and development: Providing more accurate fracture properties of shale reservoirs ( , ) and anisotropy parameters ( This provides more reliable geophysical interpretation and engineering parameter support for key aspects of shale reservoirs, such as seismic anisotropy imaging, reservoir heterogeneity description, fracturing design optimization, "sweet spot" prediction, and reserve assessment, thus promoting the development of efficient shale gas exploration and development technologies. Attached Figure Description

[0047] Figure 1 Scanning electron microscope image of typical microstructure of Longmaxi Formation shale reservoir;

[0048] Where a represents a crack in the parallel bedding, b represents interlocking clay minerals and kerogen, c represents intergranular pores of clay minerals, and d represents intergranular pores between mineral particles.

[0049] Figure 2 This is a schematic diagram of a shale clay matrix.

[0050] Figure 3 The azimuth distribution characteristics of clay matrix corresponding to Gaussian distribution functions with different standard deviations; among them, in (a) =0°, in (b) =30°, in (c) =90°.

[0051] Figure 4 This is a process for rock physics modeling of fractured shale based on the dominant orientation distribution correction of clay matrix.

[0052] Figure 5 The variation of velocity and anisotropic parameters with the degree of directional alignment of clay;

[0053] Among them, (a) shows the relationship between the degree of directional alignment of the clay matrix and the variation of the longitudinal wave velocity (Vp) and transverse wave velocity (vs) of the model, and (b) shows the relationship between the degree of directional alignment of the clay matrix and the variation of the Thomsen anisotropy parameter.

[0054] Figure 6 The values ​​represent the changes in velocity and anisotropy parameters with the crack aspect ratio; where (a) represents the relationship between the longitudinal wave velocity (Vp) and transverse wave velocity (vs) in the crack aspect ratio model, and (b) represents the relationship between the crack aspect ratio and the Thomsen anisotropy parameters.

[0055] Figure 7 The values ​​represent the changes in velocity and anisotropy parameters with crack density; where (a) shows the relationship between crack density and the model's P-wave velocity (Vp) and S-wave velocity (vs), and (b) shows the relationship between crack density and Thomsen's anisotropy parameters.

[0056] Figure 8 The values ​​represent the changes in velocity and anisotropy parameters with spherical porosity; where (a) represents the relationship between the porosity of the spherical pore and the model longitudinal wave velocity (Vp) and transverse wave velocity (vs), and (b) represents the relationship between the porosity of the spherical pore and the Homsen anisotropy parameter.

[0057] Figure 9 This is a flowchart of the inversion of fracture attribute parameters in shale reservoirs based on the dominant orientation distribution of clay matrix.

[0058] Figure 10 A cross-plot of trial calculation results for the inversion model of reservoir fracture attribute parameters;

[0059] In the diagram, the black dashed line represents the true value of the parameter, the blue circle represents the reservoir fracture attribute inversion result based on the clay matrix dominant orientation distribution correction, and the red circle represents the reservoir fracture attribute inversion result using the conventional Eshelby-Cheng model.

[0060] Figure 11 This is a logging image of the Longmaxi Formation shale from Well A in the Jiaoshiba area.

[0061] Figure 12 These are the actual and simulated values ​​of the longitudinal and transverse wave velocities in the well.

[0062] Figure 13 The results are the inversion results of the directional arrangement of the clay matrix and the reservoir fracture property parameters.

[0063] Figure 14 Statistical histogram of parameter inversion results.

[0064] Figure 15This is the result of reservoir anisotropy parameter inversion. Detailed Implementation

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

[0066] 1. Theory and Methodology

[0067] 1.1 Elastic parameters of anisotropic media

[0068] Anisotropy is the most significant geophysical response characteristic of shale reservoirs. Compared with sandstone and carbonate rocks, shale generally exhibits stronger velocity anisotropy (>15%), predominantly characterized by transverse isotropy (VTI) along its vertical axis of symmetry. Its elastic properties can be represented by a stiffness matrix with five independent components, as shown in the following expression:

[0069] (1).

[0070] in It can also be represented using the Thomsen anisotropy parameter:

[0071] (2).

[0072] in, It is the vertical longitudinal wave velocity. It is the vertical transverse wave velocity. These are three dimensionless anisotropic parameters. Well logging methods can obtain accurate P-wave velocity and S-wave velocity, but cannot directly measure anisotropic parameters. Therefore, they can be indirectly calculated by establishing the relationship between the internal properties of the rock and the P-wave velocity, S-wave velocity, and anisotropic parameters through a rock physics model.

[0073] 1.2 Anisotropic Self-Compatible Approximation (SCA) Model

[0074] The self-compatible anisotropic approximation (SCA) model can simulate the stiffness matrix of a perfectly oriented medium. This model treats the inclusion body as an ellipsoid and can simulate the state of perfectly oriented mineral grains, as expressed below:

[0075] (3).

[0076] in, It is the equivalent medium stiffness tensor obtained using the SCA method. This corresponds to the volume content of the contained substances. This is the stiffness tensor of the containing object, where n and p refer to the containing object, and I is the unit tensor matrix. It is an elastic tensor that controls the shape of the contents.

[0077] 1.3 Anisotropic Differential Equivalent Medium (DEM) Model

[0078] The anisotropic differential equivalent medium (DEM) model can sequentially add other rock components as inclusions into the anisotropic medium, as shown in the following expression:

[0079] .

[0080] in, It is the equivalent medium stiffness tensor obtained using the DEM method. Let v be the stiffness tensor of the i-th inclusion, v be the volume fraction of the inclusion, and I be the unit tensor matrix. It is an elastic tensor that controls the shape of the contents.

[0081] 1.4 Azimuth Distribution Function

[0082] The orientation distribution function (ODF) can simulate the effect of mineral grain orientation on rock anisotropy. The equivalent stiffness matrix elements based on Voigt's iso-strain averaging are:

[0083] (5).

[0084] The anisotropic intermediate and the isotropic intermediate are represented as follows:

[0085] (6).

[0086] in, Represents the elements of the equivalent stiffness matrix. and There are two Legendre coefficients. Represents the stiffness matrix elements of a completely oriented medium.

[0087] For a transversely isotropic (VTI) medium with a perpendicular axis of symmetry, its Legendre coefficient can be expressed as:

[0088] (7).

[0089] in, , , ; In the Euler Cartesian coordinate system Shaft and inclusions Interaxial angle.

[0090] The Gaussian distribution can describe the directional arrangement of minerals and is used to calculate the equivalent stiffness matrix. Its expression is as follows:

[0091] (8).

[0092] in, Normalization equations The solution is obtained.

[0093] 1.5 Eshelby-Cheng Model

[0094] The Eshelby-Cheng model proposes a novel expansion based on the Pade approximation, avoiding the problem that the second-order expansion of the Hudson model is not a single convergent sequence and that using first-order corrections yields better results than inappropriately using second-order corrections. It treats the fractures in the rock as idealized ellipsoidal inclusions arranged in a perfectly oriented manner, with the fractures isolated and fluids unable to flow between them. The equivalent stiffness matrix expression of the Eshelby-Cheng model is:

[0095] (9).

[0096] in, It is the porosity of the crack. It is the stiffness matrix of isotropic rocks. This is a correction for the effect of cracks on isotropic rocks, and its expression is:

[0097] (10).

[0098] The expression for its intermediate parameter is as follows:

[0099] (11).

[0100] in, It's the aspect ratio of the crack. and It is the modulus of isotropic rocks.

[0101] 1.6 Brown-Korringa Theory

[0102] Gassmann theory is an important tool for studying fluid substitution problems, applicable to isotropic rocks or equivalent media. Brown and Korringa extended Gassmann theory, establishing the Brown-Korringa theory applicable to fluid substitution in anisotropic media. This theory is also applicable to low-frequency conditions, i.e., the anisotropic form of Gassmann theory, expressed as follows:

[0103] (12).

[0104] in, It is the effective elastic flexibility tensor of the dry rock skeleton. It is the effective elastic flexibility tensor of pore fluid-saturated rock. It is the effective elastic flexibility tensor of the rock matrix. It is the compressibility coefficient of pore fluid. It is the compressibility coefficient of the rock matrix. It is porosity. It is the bulk modulus of pore fluids. It is the bulk modulus of rocks and minerals.

[0105] 2. Rock physics modeling and model trial calculations

[0106] 2.1 Petrophysical Modeling of Fractured Shale Based on Clay Matrix Dominant Orientation Distribution Correction

[0107] Shale oil and gas reservoirs are strata composed of dark, organic-rich mudstone and shale, along with other rocks interbedded in thin layers. They consist of diagenetic minerals (clay and non-clay minerals), multi-scale and multi-type pore spaces, organic kerogen, and multi-phase pore-filling materials. Scanning electron microscopy (SEM) is an important tool for studying the microscopic porosity and structural characteristics of rocks. Typical SEM images of the microstructure of Longmaxi Formation shale reservoir cores are shown below. Figure 1 As shown. Figure 1 In the diagram, 'a' represents a crack parallel to the bedding planes. These cracks are mostly horizontal or nearly horizontal in shape, and brittle minerals are embedded in the clay matrix as dispersed particles. Figure 1 In the diagram, 'b' represents interlocking clay minerals and kerogen. The clay minerals are mostly layered, while the kerogen is embedded within them, exhibiting an overall directional distribution pattern. Figure 1 In the diagram, 'c' represents the intergranular pores of clay minerals. Due to the high compressibility of clay minerals, and influenced by sedimentation and diagenetic compaction, these pores are flattened, fracture-like structures. Figure 1 In this context, d represents the intergranular pores between mineral grains. These pores are formed by incomplete cementation of mineral grains or by later diagenetic processes. They are mostly spherical or near-spherical in shape, but some are also elliptical or irregular.

[0108] Shale reservoirs exhibit significant seismic anisotropy. Clay minerals are one of the main factors causing this inherent anisotropy. As the primary load-bearing matrix, plastic clay minerals reduce the mechanical strength of shale rocks. Due to sedimentary processes and subsequent diagenetic compaction, the mixture of clay minerals and kerogen is typically oriented, while numerous horizontally distributed fractures parallel to bedding planes are also present. These factors are the main causes of anisotropy in shale reservoirs. Describing the oriented arrangement of clay minerals and kerogen, as well as the horizontally distributed fractures, using appropriate rock physics models is both a key focus and a challenge in shale rock physics modeling.

[0109] Clay particles are the most basic units constituting the solid clay matrix of shale reservoirs. The spatial distribution of multiple clay particle units follows a dominant azimuth distribution function, with specific statistical characteristics in their azimuth and aspect ratio. They exhibit a certain macroscopic trend in depth, depending on the degree of compaction of the shale reservoir, resulting in inherent anisotropic characteristics in the elastic properties of shale. Kerogen has a relatively low elastic modulus, close to that of fluids, and is embedded in clay minerals, exhibiting a similar directional arrangement and distribution characteristics to clay particles. Therefore, it needs to be considered when constructing a laterally isotropic clay matrix. A schematic diagram is shown below. Figure 2 As shown.

[0110] The plate-like clay matrix unit composed of clay minerals and kerogen can be equivalently represented as a transversely isotropic (VTI) medium with a vertical axis of symmetry. Kerogen can be considered as an ellipsoidal inclusion parallel to the distribution of clay minerals. For any porosity, the clay minerals and kerogen in the plate-like dual-phase system are interconnected. The equivalent elastic stiffness tensor of the transversely isotropic clay matrix is ​​the basis for simulating the elastic properties of the solid matrix in shale reservoirs. Due to the significant difference in elastic parameters between kerogen and clay minerals, to avoid the connectivity defects present in the anisotropic SCA model, the equivalent elastic parameters of the plate-like clay matrix unit must be simulated in conjunction with the anisotropic DEM model. This ensures the interconnection of clay minerals and kerogen, avoids the asymmetry of the elastic modulus of the mixture caused by the order of addition, and considers the directional arrangement and distribution characteristics of clay minerals and kerogen. Influenced by sedimentation and subsequent diagenetic compaction, the orientation of platy clay matrix units varies, but they are generally arranged along a dominant direction. This results in multiple platy clay matrix units with different orientations constituting the clay matrix in shale reservoirs, giving shale inherent anisotropic elastic properties. The azimuth distribution function quantifies the degree of directional alignment of platy clay matrix units, reflecting the degree of shale compaction; this geological significance is reflected in the coefficients of the azimuth distribution function. Compared to Fisher or Gamma distributions, the standard deviation and mean of the Gaussian distribution are more intuitive. Figure 3This diagram illustrates the azimuth distribution function simulated for Gaussian distributions with different standard deviations and the corresponding clay particle distributions. When the standard deviation of the azimuth distribution function is within 0 ± 1°, it indicates that the clay matrix elements are completely oriented, and the corresponding Legendre coefficient... and Maximum, respectively and This exhibits the greatest degree of anisotropy; when the standard deviation of the azimuthal distribution function is less than 90 ± 1°, it can be considered that the clay matrix elements are completely randomly distributed, and the corresponding Legendre coefficient is... and The minimum values ​​are all 0, exhibiting isotropic characteristics.

[0111] Brittle minerals such as quartz, calcite, and pyrite are important physical properties for improving the fracturing resilience of shale reservoirs. They are typically distributed dispersedly within the clay matrix. Compared to clay minerals and kerogen, they have higher elastic moduli and are more brittle. To avoid differences in modeling results due to varying order of addition of brittle minerals, an isotropic SCA model can be used to simulate the equivalent elastic modulus of the brittle mineral mixture. Considering that the pores in shale are mainly non-connected spherical pores, when modeling anisotropic shale matrices, an anisotropic DEM model can be used to progressively add brittle mineral mixtures and spherical pores into the clay matrix to obtain an anisotropic shale matrix model.

[0112] In common rock physics models containing horizontal fractures, the matrix is ​​typically isotropic. However, due to the inherent anisotropy of shale, the shale matrix exhibits weak anisotropy. Rock physics experiments have demonstrated the applicability of fracture-type rock physics models like the Eshelby-Cheng model to weakly anisotropic matrices. By adding horizontal fractures to the anisotropic shale matrix using the Eshelby-Cheng model with anisotropic matrices, an anisotropic shale skeleton model is constructed. The equivalent stiffness matrix expression of this model is:

[0113] (12).

[0114] in, It is an anisotropic matrix. These are the improved corrections. This refers to crack porosity. In the Eshelby-Cheng model, the correction term is calculated from the matrix modulus and crack properties. For weakly anisotropic matrices, the matrix modulus required for the correction term is given by the following formula:

[0115] (13).

[0116] in, , , , , It is the density of the matrix. and It is the velocity of an equivalent isotropic matrix.

[0117] The pore fluids in shale are mainly gas and water. Since shale has significant transverse isotropic characteristics and in order to be applicable to low-frequency conditions, the equivalent elastic modulus of the fluid mixture can be calculated first using the Wood formula. Then, the Brown-Korringa model can be used to replace the pores with fluid, thus obtaining a fluid-saturated shale rock physics model.

[0118] In summary, the equivalent elastic modulus of horizontally fractured shale, corrected for the dominant orientation distribution of the clay matrix, can be simulated. Its longitudinal wave velocity, transverse wave velocity, and anisotropic parameters can be expressed by the above equivalent elastic modulus:

[0119] (14).

[0120] The rock physics model of fractured shale proposed in this invention, based on the correction of the dominant orientation distribution of clay matrix, is as follows: Figure 4 As shown. The specific modeling process is as follows: (1) The anisotropic SCA-DEM model is used to simulate the completely oriented clay and kerogen mixture (clay matrix unit) to obtain the completely oriented clay matrix model, i.e., the anisotropic background medium one.

[0121] (2) The background medium one is adjusted by using the orientation distribution function (ODF) to obtain the clay matrix model with partial orientation of clay and kerogen mixture (clay matrix unit), i.e., the anisotropic background medium two.

[0122] (3) Use the isotropic SCA model to mix brittle minerals (including quartz, calcite and pyrite, etc.) to simulate brittle mineral mixtures. Then use the anisotropic DEM model to gradually add brittle mineral mixtures and matrix pores (spherical pores) to the background medium 2 to obtain the shale anisotropic matrix model.

[0123] (4) Horizontal fractures were added to the anisotropic matrix using the improved Eshelby-Cheng model to obtain the shale dry skeleton model.

[0124] (5) The pore fluid was mixed using the Wood formula to simulate the fluid mixture. The fluid was replaced under anisotropic conditions using the Brown-Korringa model to obtain a fluid-saturated shale rock physics model.

[0125] 2.2 Relationship between elastic properties and reservoir fracture properties

[0126] In the above modeling process, the influence of the partial orientation of the clay matrix on anisotropy was simulated using the azimuth distribution function, and the influence of horizontal fractures was simulated using the improved Eshelby-Cheng model. This model can accurately construct the relationship between reservoir fracture properties and elastic parameters, and combined with the actual P-wave velocity and S-wave velocity in the well, it can estimate the fracture properties of shale reservoirs.

[0127] The relationship between crack density and crack width / height is as follows:

[0128] (15).

[0129] in, It is the crack density. It is the aspect ratio of the crack. It refers to crack porosity.

[0130] In the model, the reservoir pore space is composed of fractures and spherical pores, and their relationship is as follows:

[0131] (16).

[0132] in, It is the total porosity of the reservoir. It refers to the porosity of spherical pores. It refers to crack porosity.

[0133] To analyze the relationship between the model's elastic properties and reservoir fracture attributes, the standard deviation of the angle between the major axis of the clay matrix element and the formation can be used. The degree of directional alignment of clay matrix units is described. A standard deviation approaching 0° indicates that the clay matrix units are completely oriented, while a standard deviation approaching 90° indicates that the overall clay matrix exhibits isotropic characteristics. Reservoir fracture properties can be described by fracture density and fracture aspect ratio. This invention simulates the effects of clay matrix directional alignment, fracture aspect ratio, fracture density, and spherical porosity on the model's P-wave velocity, S-wave velocity, and Thomsen anisotropy parameter, respectively. The results are as follows: Figure 5-8 As shown, the crack density can be obtained by converting the above parameters, and the Thomsen anisotropy parameter... These characterize the strength of anisotropy in P-wave velocity, the strength of anisotropy in S-wave velocity, and the second derivative of the P-wave phase velocity under perpendicular incidence. When the crack aspect ratio is 0.3 and the spherical porosity is 0.06, [the following is a separate, unrelated sentence:] ... Figure 5 It can be seen that the longitudinal wave velocity and transverse wave velocity are negatively correlated with the degree of directional alignment of the clay matrix, and the anisotropy of longitudinal wave velocity and transverse wave velocity are also negatively correlated, that is, they increase with standard deviation. Increase, longitudinal wave velocity and transverse wave velocity It also increases the anisotropy parameter. The decrease indicates that the change in standard deviation is related to the change in standard deviation. As the matrix increases, the model transforms into an isotropic medium, resulting in increased P and S wave velocities and decreased anisotropy parameters. The second derivative of the P wave phase velocity... As the degree of directional alignment increases, it first increases and then decreases. When the standard deviation is 30° and the spherical porosity is 0.06, it is determined by... Figure 6 It can be seen that the P-wave velocity and S-wave velocity are positively correlated with the fracture aspect ratio, while the P-wave velocity anisotropy and S-wave velocity anisotropy are negatively correlated, that is, they increase with the fracture aspect ratio. The increase of longitudinal wave velocity and transverse wave velocity It also increases the anisotropy parameter. As the crack morphology decreases, the change in P-wave velocity is more pronounced than that in S-wave velocity, indicating that as the crack morphology in the model changes towards a spherical shape, both P-wave and S-wave velocities increase, the anisotropy of the model decreases, and the effect on P-wave velocity is more significant. When the standard deviation is 30° and the crack aspect ratio is 0.3, [the following text appears to be incomplete and requires further context: "from..."] Figure 7-8 It can be seen that the P-wave velocity and S-wave velocity are negatively correlated with the fracture density, while the P-wave velocity anisotropy and S-wave velocity anisotropy are positively correlated, that is, they increase with fracture density. The increase of longitudinal wave velocity and transverse wave velocity It also reduces the anisotropy parameter. As the porosity increases, both longitudinal and transverse wave velocities show a positive correlation with spherical porosity, while longitudinal wave velocity anisotropy and transverse wave velocity anisotropy show a negative correlation. That is, with increasing spherical porosity... The increase of longitudinal wave velocity and transverse wave velocity It also increases the anisotropy parameter. As the velocity decreases, the change in P-wave velocity is more pronounced than that in S-wave velocity, indicating that variations in pore size within the model affect both the model's velocity and anisotropy. The influence of spherical pores is greater than that of fractures, and their impact on P-wave velocity anisotropy is greater than that on S-wave velocity anisotropy. This demonstrates that reservoir fracture properties and the partial directional arrangement of the clay matrix significantly influence the elastic properties of shale, and both must be considered in rock physics modeling.

[0134] 2.3 Inversion of Reservoir Fracture Attribute Parameters

[0135] Quantifying fracture properties in shale reservoirs requires inversion estimation using rock physics models and measured data. While complex rock physics models can more accurately simulate the impact of actual shale reservoir fracture properties on elastic properties, they suffer from difficulties in initialization and ambiguity in multi-parameter inversion. Compared to seismic exploration methods, well logging provides more accurate and abundant prior information such as rock physics parameters, rock velocities, and core data, allowing for thorough initialization of the rock physics model and improving the accuracy of fracture property inversion. Incorporating constraint information can reduce the impact of ambiguity during the inversion process. By using wellbore P-wave and S-wave velocities as constraints, the simulated P-wave and S-wave velocities can be matched to suppress the influence of ambiguity. This invention utilizes the rich rock physics parameter information provided by well logging methods, combined with a rock physics model, to establish the relationship between shale microstructure and well logging curves, and then inverts the fracture property parameters of the rock based on P-wave and S-wave velocities. To account for the influence of clay matrix orientation on rock elastic properties, the degree of clay matrix orientation is also used as an inversion parameter. The objective function for model parameter inversion, combining well logging data, is constructed as follows:

[0136] (17).

[0137] in It is the measured longitudinal wave velocity in the well. It is the actual measured shear wave velocity in the well. It is the longitudinal wave velocity simulated by the model. It is the shear wave velocity simulated by the model. It refers to the degree of directional arrangement of clay matrix units. It is the aspect ratio of the crack. It refers to the crack density.

[0138] To simultaneously solve for the three unknown parameters in the objective function, a global optimization algorithm with strong capabilities for solving nonlinear, multi-parameter problems and high computational efficiency is required. This invention employs the Simulated Annealing Particle Swarm Optimization (SA-PSO) algorithm, which combines the advantages of simulated annealing and particle swarm optimization. This algorithm can quickly and effectively find the optimal solution and has a certain ability to escape local optima.

[0139] The process of inverting fracture attribute parameters in shale reservoirs based on clay matrix dominant orientation distribution correction is as follows: Figure 9 As shown, the specific process is as follows:

[0140] 1. Input prior information such as well logging data and core data.

[0141] 2. Construct a rock physics model of fractured shale based on the dominant orientation distribution correction of clay matrix and initialize the model using prior information.

[0142] 2.1 By comprehensively utilizing the anisotropic SCA-DEM model to simulate the completely oriented clay and kerogen mixture (clay matrix unit), a completely oriented clay matrix model is obtained.

[0143] 2.2 The orientation distribution function (ODF) is used to adjust the fully oriented clay matrix to obtain a clay matrix model in which the clay and kerogen mixture (clay matrix unit) is partially oriented.

[0144] 2.3 Using an isotropic SCA model, brittle minerals (including quartz, calcite, and pyrite) are mixed to simulate a brittle mineral mixture. Then, using an anisotropic DEM model, the brittle mineral mixture and matrix pores (spherical pores) are gradually added to the background medium 2 to obtain the shale anisotropic matrix model.

[0145] 2.4 The improved Eshelby-Cheng model was used to add horizontal fractures to the anisotropic matrix to obtain the shale dry skeleton model.

[0146] 2.5 The Wood formula was used to mix the pore fluids to simulate the fluid mixture. Then, the Brown-Korringa model was used to perform fluid replacement under anisotropic conditions to obtain a fluid-saturated shale rock physics model.

[0147] 3. Using the P-wave velocity and S-wave velocity in the well as constraints, the simulated P-wave velocity and S-wave velocity are matched with them to suppress the influence of multiple solutions.

[0148] 4. A model parameter inversion objective function combining well logging data was constructed: .

[0149] 5. Introduce the simulated annealing optimized particle swarm optimization (SA-PSO) algorithm to improve the ability to solve nonlinear and multi-parameter problems and the computational efficiency of the inversion process.

[0150] Specifically, the annealing optimization particle swarm optimization algorithm includes the following processes:

[0151] 5.1 Initialize particles.

[0152] 5.2 Calculate particle fitness.

[0153] 5.3 Update the temporary positions of particles.

[0154] 5.4 Determine if the particle movement rules are met. If yes, update the particle position and velocity, and proceed to process 5.5. If no, restore the particle's original position and velocity, and then proceed to process 5.5.

[0155] 5.5 Cooling operation.

[0156] 5.6 Determine if the termination condition is met. If yes, output the optimal particle position; otherwise, return to process 5.2.

[0157] 6. Output reservoir fracture attribute parameters and anisotropy parameters.

[0158] The inversion results of reservoir fracture attribute parameters are as follows Figure 10 As shown in the figure, the black dashed line represents the true value of the parameter, the blue circle represents the reservoir fracture attribute inversion result based on the clay matrix dominant orientation distribution correction, and the red circle represents the reservoir fracture attribute inversion result using the conventional Eshelby-Cheng model. It can be clearly seen from the figure that the reservoir fracture attribute inversion result based on the clay matrix dominant orientation distribution correction has a higher degree of fit with the true value of the parameter, which is far better than the inversion result of the conventional method, thus verifying the effectiveness and accuracy of the method of the present invention.

[0159] Example 1

[0160] The Sichuan Basin, a major favorable area for marine shale gas exploration in China, is located at the transition zone between the Tethys-Himalayan tectonic domain and the Circum-Pacific tectonic domain. The Jiaoshiba Longmaxi Formation shale gas field is situated in the transitional zone between the northeast-trending and near-north-south trending structures in southeastern Sichuan Basin. Two global marine transgressions occurred at the end of the Ordovician and the beginning of the Silurian, forming the organic-rich Longmaxi Formation shale. This invention selects the shale at the bottom of Well A in the Jiaoshiba Longmaxi Formation shale gas field in the Sichuan Basin, China, as the research object. This stratum is a gas-producing stratum, characterized by stable shale distribution, high organic carbon content, high organic matter maturity, and suitable brittleness.

[0161] Well logging data of Longmaxi Formation shale in Well A in Jiaoshiba area as follows: Figure 11 As shown, it provides a deep, natural gamma spectrum ( ), P-wave velocity ( ), transverse wave velocity ( ),density( ), porosity ( ), water saturation ( Information on the mineral composition of the shale matrix minerals and pore fluids is presented. Well logging data shows a positive correlation between the P-wave and S-wave velocity curves at the studied well location, and a negative correlation between the velocity curves and the natural gamma-ray spectrum curve. Core analysis indicates that the Longmaxi Formation rocks in Well A are primarily composed of clay and quartz, with minor amounts of calcite and pyrite, and also contain organic kerogen. The pore fluids are mainly water and gas. Table 1 presents the elastic modulus and density of the shale matrix minerals and pore fluids.

[0162] Table 1 Elastic modulus and density of shale matrix minerals and pore fluids

[0163]

[0164] The rock physics model of shale reservoir based on clay matrix dominance distribution correction, proposed in this invention, is initialized using well logging curves, mineral composition information, and the elastic modulus and density of matrix minerals and pore fluids as prior information. The simulated P-wave and S-wave velocities are matched to the model using actual wellbore P-wave and S-wave velocities as constraints. The inversion process is optimized using simulated annealing optimized particle swarm optimization (SA-PSO) algorithm to solve the objective function, inverting the clay matrix orientation, fracture aspect ratio, and fracture density. Simultaneously, the Thomsen anisotropy parameters are estimated. The inversion results are as follows: Figure 11-13 As shown.

[0165] Figure 12 The figure shows the matching between the P-wave and S-wave velocities measured in actual well logging and those obtained from the rock physics model. The black curve represents the actual values, and the red curve represents the simulated values. The simulated and actual curves show a high degree of fit and small global error, verifying that the rock physics model constructed in this invention is suitable for the well logging data. However, because both P-wave and S-wave velocities are used for constraints, local errors exist in the results. In actual production, due to limitations in well location conditions, logging technology, and acquisition costs, S-wave velocity information in well logging data is often scarce. The model of this invention can be used to predict S-wave velocity in well logging data lacking S-wave velocity information, providing more accurate S-wave velocity data for subsequent reservoir exploration and development.

[0166] Figure 13 The results represent the inversion of the directional arrangement of clay matrix and reservoir fracture attribute parameters in shale reservoirs, where the reservoir fracture attributes are described by fracture aspect ratio and fracture density. Figure 14 Statistical histograms of the inversion results of the directional arrangement of clay matrix and reservoir fracture attribute parameters. Figure 15 The results show the inversion of reservoir anisotropy parameters. The results indicate that the directional arrangement of the clay matrix in this section exhibits periodic jumps, the fracture aspect ratio is generally low (around 0.1), and the fracture density is around 0.02. The Thomsen anisotropy parameters show an overall positive correlation, with ε and γ showing a strong correlation, while γ and δ show a negative correlation in some areas. Combined with well logging data and geological background analysis, this section exhibits significant logging response characteristics of high natural gamma, high porosity, and low density. The shale reservoir has a distinct layered structure and is generally transversely isotropic, but this layered structure causes periodic jumps in the inversion results of the directional arrangement of the clay matrix. Plastic clay minerals, as the main load-bearing matrix, reduce the mechanical strength of the shale rock. As the clay matrix tends towards isotropic characteristics, the morphology of horizontal fractures tends to be spherical, and the fracture density decreases accordingly. With increasing depth, the organic matter kerogen content increases, and the rock quartz mineral content also increases, resulting in a decreasing fracture density and an increasing fracture aspect ratio.

[0167] This invention focuses on accurately quantifying fracture attribute parameters to achieve effective fracture characteristic description. Starting from the geological and geophysical characteristics of shale reservoirs, it comprehensively analyzes the influence of directional and horizontally distributed fractures in the clay matrix on shale reservoir anisotropy. A reasonable anisotropic petrophysical model is constructed using rock physics theory. Combined with well logging data, a method for inverting fracture attributes and anisotropic parameters in shale reservoirs based on correction for the dominant azimuth distribution of the clay matrix is ​​proposed. This method uses the azimuth distribution function (ODF) to describe the influence of partially directional clay matrix arrangement on rock anisotropy. Simultaneously, a modified Eshelby-Cheng model is used to simulate the influence of horizontal fractures, avoiding attributing the influence of directional clay matrix arrangement on anisotropy to fractures, thus improving the accuracy of shale reservoir fracture attribute description. Model trials have verified the effectiveness and accuracy of this method, and its application to marine shale reservoirs in the Sichuan Basin of China has yielded ideal results. This proposed method effectively improves the accuracy of fracture attribute estimation and anisotropy parameter inversion in shale reservoirs, providing richer and more accurate supporting information for shale reservoir seismic imaging and reservoir characterization, and promoting the exploration and development of shale gas resources. Shale has a complex microstructure, with many factors affecting its anisotropy. Although the directional arrangement of the clay matrix is ​​a common and significant factor, further research is needed to establish a more accurate relationship between fracture attributes and anisotropy and to improve the accuracy of fracture characterization. This requires the development of shale petrophysical models that simulate more factors and corresponding inversion methods.

[0168] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for inverting shale fracture properties based on clay matrix dominant orientation distribution correction, characterized in that, Includes the following steps: S1. Constructing a rock physics model: S11. By comprehensively utilizing the anisotropic SCA model and the anisotropic DEM model, a transversely isotropic medium with a completely oriented vertical symmetry axis formed by the mixture of clay minerals and kerogen is simulated as background medium one. S12. Based on the azimuth distribution function, using the standard deviation of the Gaussian distribution. The degree of directional alignment of the background medium is quantified, combined with the Legendre coefficient. and Adjust the equivalent elastic stiffness tensor of background medium one to obtain a partially oriented clay matrix model, which serves as background medium two. S13. The brittle mineral mixture is mixed using an isotropic SCA model to simulate the brittle mineral mixture. The brittle mineral mixture and matrix porosity are gradually added to the background medium II using an anisotropic DEM model to construct an anisotropic matrix model of shale. S14. Horizontal fractures were added to the shale anisotropic matrix model using the improved Eshelby-Cheng anisotropic matrix model, and correction terms were calculated using the equivalent anisotropic matrix. A shale skeleton model was obtained; S15. Using Wood's formula to mix pore fluids, simulate fluid mixtures, and combine the Brown-Korringa anisotropic fluid substitution theory to fill the shale skeleton model with fluid to obtain a fluid-saturated shale rock physics model. S2. Establish the inversion objective function: Using the measured P-wave velocity and S-wave velocity in the well as dual constraints, define the objective function as follows: ； in, It is the measured longitudinal wave velocity in the well. It is the actual measured shear wave velocity in the well. It is the longitudinal wave velocity simulated by the model. It is the shear wave velocity simulated by the model. It's an angle. Standard deviation of the distribution It is the aspect ratio of the crack. It is the crack density; S3. Execution parameter inversion: The objective function is solved by using simulated annealing optimized particle swarm optimization algorithm, and the crack attribute parameters and anisotropy parameters are output.

2. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, In step S12, the Gaussian distribution function is defined as: ； in, Normalization equations The solution is obtained.

3. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, In step S12, when the standard deviation of the azimuth distribution function is 0°, it can be said that the clay matrix element is completely oriented, and its corresponding Legendre coefficient is... and Maximum, respectively and This indicates the highest degree of anisotropy; when the standard deviation of the azimuthal distribution function is 90°, it represents a completely random distribution of clay matrix elements, and the corresponding Legendre coefficient... and The minimum values ​​are all 0, exhibiting isotropic characteristics.

4. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, In step S13, the brittle mineral includes any one or more of quartz, calcite, and pyrite.

5. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, In step S14, the equivalent stiffness matrix expression of the improved anisotropic matrix Eshelby-Cheng model is as follows: ； in, It is the rigid matrix of the anisotropic matrix. These are the improved corrections. It refers to crack porosity.

6. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 4, characterized in that, Improved corrections The required matrix modulus is given by the following formula: ； in, , , , , It is the density of the matrix. and It is the velocity of an equivalent isotropic matrix.

7. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, The relationship between fracture density and fracture cross-section in the rock physics model is as follows: ； in, It is the crack density. It is the aspect ratio of the crack. It refers to crack porosity; In rock physics models, reservoir pore space consists of fractures and spherical pores, with the following relationship: ； in, It is the total porosity of the reservoir. It refers to the porosity of spherical pores. It refers to crack porosity.

8. The shale fracture attribute inversion method based on clay matrix dominant orientation distribution correction as described in claim 1, characterized in that, In step S3, the annealing optimization particle swarm algorithm includes the following steps: 31) Initialize the particles; 32) Calculate particle fitness; 33) Update the temporary positions of the particles; 34) Determine if the particle movement rules are met. If yes, update the particle position and velocity, and proceed to step 35); otherwise, restore the particle's original position and velocity, and then proceed to step 35). 35) Cooling operation; 36) Determine if the termination condition is met. If yes, output the optimal particle position; otherwise, return to step 32.

Images