Coal reservoir ground stress prediction method and system

Abstract

The present application relates to coal reservoir stress prediction technology field, specifically to a kind of coal reservoir stress prediction method.The method includes according to the depth distribution of drilling filtrate invasion formation and cleat direction, configure array sonic logging source distance;After source distance configuration, multiple logging is collected under each source distance P-wave velocity and S-wave velocity, and formation bulk density, cleat porosity is collected;Longitudinal wave velocity ratio is calculated by source distance, and far well boundary is selected, further identify far well P-wave velocity, far well S-wave velocity, near well P-wave velocity and near well S-wave velocity;Inversion dry skeleton bulk modulus and dry skeleton shear modulus;Near well wave velocity is replaced by fluid to obtain corrected P-wave velocity, corrected S-wave velocity and original gas-water two-phase state formation bulk density;Elastic parameter is calculated and substituted into geostress calculation model, and geostress prediction result is obtained.The present application effectively reduces the near well acoustic distortion caused by drilling fluid filtrate, and improves the accuracy of coal reservoir stress prediction.

Description

Technical Field

[0001] This invention belongs to the field of geostress prediction technology, specifically relating to a method and system for predicting geostress in coal reservoirs. Background Technology

[0002] In coalbed methane exploration and development, in-situ stress prediction of coal reservoirs is a core foundation for evaluating wellbore stability, designing hydraulic fracturing schemes, and analyzing the evolution of reservoir permeability. The in-situ stress state directly affects the initiation and propagation direction of fracturing fractures, thus impacting the gas production efficiency and long-term production safety of coalbed methane wells. Currently, the widely used method in industry is to obtain the P-wave and S-wave velocities of the formation through sonic logging, combine this with the formation bulk density obtained from density logging, calculate the elastic parameters of the rock, and then substitute these parameters into an in-situ stress calculation model to obtain a continuously distributed in-situ stress profile along the wellbore.

[0003] A significant characteristic distinguishing coal reservoirs from conventional sandstone reservoirs is their highly developed cleavage system. Face cleavages and end cleavages form a complex pore-fracture network, serving as both the primary storage space for coalbed methane and a crucial seepage channel for gas production. During drilling, to maintain wellbore stability, the drilling fluid column pressure is typically designed to exceed the coal seam pore pressure. This positive pressure differential drives drilling fluid filtrate to invade the near-wellbore formation. Since cleavage permeability is much higher than that of the coal matrix, the filtrate preferentially and rapidly invades the cleavage system, displacing free gas and some adsorbed gas within the cleavages, resulting in a significant increase in the water saturation of the cleavages in the near-wellbore region. This fluid displacement alters the acoustic response of the coal rock, transforming cleavages from being saturated with highly compressible gas to being saturated with almost incompressible water, causing a shift in acoustic velocity and consequently distorting the predicted in-situ stress of the coal reservoir. Summary of the Invention

[0004] (1) Technical problems to be solved The purpose of this invention is to provide a method and system for predicting geostress in coal reservoirs, so as to solve the problem that the near-wellbore acoustic velocity is distorted due to the intrusion of drilling fluid filtrate into the coal reservoir cleavage, resulting in inaccurate geostress prediction.

[0005] (2) Technical solution To achieve the above objectives, in one aspect, the present invention provides a method for predicting geostress in coal reservoirs, the method comprising: Based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage, the source spacing of the array acoustic logging is configured; after the source spacing is configured, array acoustic logging is performed to obtain the P-wave velocity and S-wave velocity at each source spacing; density logging is performed to obtain the formation volume density; and nuclear magnetic resonance logging is performed to obtain the cleavage porosity.

[0006] The P-wave and S-wave velocity ratios are calculated per source distance, and the outer boundary of the invasion zone is selected. The P-wave and S-wave velocities corresponding to source distances exceeding the outer boundary of the invasion zone are denoted as far-well P-wave velocity and far-well S-wave velocity, while those falling within the outer boundary of the invasion zone are denoted as near-well P-wave velocity and near-well S-wave velocity.

[0007] Based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity, the bulk modulus and shear modulus of the skeleton are inverted.

[0008] It is used to obtain the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the porosity of the cleavage.

[0009] Based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, elastic parameters are calculated and substituted into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

[0010] Furthermore, the method for configuring the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage includes: Based on imaging logging, the dominant azimuth angle of the surface cleavage, surface cleavage density, and end cleavage density of the target coal reservoir are identified.

[0011] The overall invasion depth of the target coal reservoir is determined based on the radial resistivity profile obtained from array inductive resistivity logging; the overall invasion depth is allocated to the face cleavage direction and the end cleavage direction based on the face cleavage density and the end cleavage density, thus obtaining the invasion depth in the face cleavage direction and the invasion depth in the end cleavage direction.

[0012] A combination of source distances for surface cleavage azimuth is configured at the azimuth corresponding to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the surface cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the surface cleavage direction. A combination of source distances for end cleavage azimuth is configured at the azimuth perpendicular to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the end cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the end cleavage direction. The combination of source distances for surface cleavage azimuth and end cleavage azimuth is used as the configuration source distance for the azimuth array sonic logging tool.

[0013] Furthermore, the method for obtaining the cleavage porosity includes: Nuclear magnetic resonance logging was performed; the nuclear magnetic resonance porosity spectrum was divided into short relaxation components and long relaxation components according to the relaxation time distribution. The porosity corresponding to the short relaxation component was taken as the matrix micropore porosity, and the difference between the total porosity calculated from the formation volume density obtained by density logging and the matrix micropore porosity was taken as the cleavage porosity.

[0014] Furthermore, the method for calculating the P-wave and S-wave velocity ratio per source distance and selecting the boundary of the intrusion zone includes: The ratio of P-wave velocity to S-wave velocity is calculated per source distance, forming a radial sequence of P-wave and S-wave velocity ratios arranged from near to far with the acoustic detection depth. Interstitial rock interference zones in the target coal reservoir are identified based on the natural gamma logging curve. The depth points corresponding to the interstitial rock interference zones are removed from the radial sequence of P-wave and S-wave velocity ratios. The arithmetic mean of the P-wave and S-wave velocity ratios of the remaining pure coal depth points is taken per source distance to obtain the effective radial sequence of P-wave and S-wave velocity ratios.

[0015] Along the acoustic detection depth from near to far, the velocity ratio radial gradient sequence is formed by subtracting the velocity ratio of the next adjacent source from the previous source distance's P-wave velocity ratio. The velocity ratio radial gradient sequence is then smoothed by moving average to obtain a smoothed velocity ratio radial gradient sequence.

[0016] The sign of each difference value in the radial gradient sequence of the smooth velocity ratio is examined one by one along the direction from far to near the sound wave detection depth. The starting depth position where the difference value turns from non-positive to positive is determined as the outer boundary of the intrusion zone.

[0017] Furthermore, the method for identifying interbedded rocky lithological interference zones in a target coal reservoir based on natural gamma logging curves includes: Natural gamma logging curves are extracted along the depth direction of the target coal reservoir; the numerical frequency distribution of the natural gamma logging curves within the depth range of the target coal reservoir is statistically analyzed; bimodal decomposition is performed on the numerical frequency distribution, and the natural gamma value corresponding to the frequency trough between the two peaks is determined as the natural gamma boundary value between the pure coal section and the interbedded gangue section.

[0018] Based on the natural gamma ray threshold, the continuous depth intervals in the natural gamma logging curves that consistently exceed the natural gamma ray threshold are identified as candidate depth intervals for interbedded rock. All depth points whose natural gamma logging curve values ​​do not exceed the natural gamma ray threshold are designated as pure coal reference segments. The shear wave velocities at each source distance are extracted from each depth point in the pure coal reference segment, and the arithmetic mean of all pure coal depth points and all source distances is taken as the pure coal reference shear wave velocity. The arithmetic mean of the shear wave velocities at each source distance within the candidate depth intervals for interbedded rock is calculated, and the continuous depth intervals with arithmetic mean higher than the pure coal reference shear wave velocity are identified as interbedded rock lithological interference intervals.

[0019] Furthermore, the method for applying a moving average smoothing process to the velocity ratio radial gradient sequence to obtain a smoothed velocity ratio radial gradient sequence includes: For the velocity ratio radial gradient sequence, count the number of source distance spans of consecutive positive gradient segments to form a positive gradient segment source distance span distribution; from the positive gradient segment source distance span distribution, identify and separate cleavage-filled positive gradient segments that span only a single source distance and fluid displacement-type positive gradient segments that span two or more source distances.

[0020] The minimum number of odd source distances that completely cover the positive gradient segment of the cleavage filling pattern is used to determine the width of the cleavage filling interference suppression window. The velocity ratio radial gradient sequence is processed by moving average using the width of the cleavage filling interference suppression window. The arithmetic mean of the gradient values ​​in each window is calculated, and the arithmetic mean is assigned to the source distance position of the window center to replace the original gradient value, thus obtaining a smooth velocity ratio radial gradient sequence.

[0021] Furthermore, the method for inverting the bulk modulus and shear modulus of the skeleton based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity includes: Based on the ratio of porosity to cleavage porosity corresponding to immovable fluid in the short relaxation component of the nuclear magnetic resonance logging relaxation time spectrum, the initial water saturation of cleavage-bound cleavage in the target coal reservoir is determined; the cleavage-bound water porosity is calculated by the product of cleavage porosity and initial water saturation of cleavage-bound cleavage; and the cleavage gas porosity is calculated by the difference between cleavage porosity and cleavage-bound water porosity.

[0022] The equivalent bulk modulus of the cleavage-bound initial water saturation and the original formation water bulk modulus, and the ratio of cleavage porosity to the pre-set original formation gas bulk modulus, are weighted and calculated according to the Wood equation to obtain the equivalent bulk modulus of the cleavage-bound fluid.

[0023] The dry skeleton shear modulus is calculated based on the far-well shear wave velocity and formation bulk density. The calculation formula is: .

[0024] in, The density of the formation is the volumetric density. The distance is the shear wave velocity.

[0025] The saturated bulk modulus of the far-wellbore is calculated based on the far-wellbore P-wave velocity, far-wellbore S-wave velocity, and formation bulk density. The calculation formula is: .

[0026] in, The distance is the longitudinal wave velocity.

[0027] By substituting the far-well saturated bulk modulus, cleavage porosity, the pre-defined bulk modulus of the coal and petrified mineral skeleton, and the equivalent bulk modulus of the cleavage-mixed fluid into the Gassmann equation, the dry skeleton bulk modulus is obtained by back-calculation.

[0028] Furthermore, the method for obtaining corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the dry skeleton bulk modulus, dry skeleton shear modulus, and cleavage porosity includes: Based on the resistivity values ​​at each depth obtained from array inductive resistivity logging, and constrained by cleavage porosity, the cleavage water saturation at each radial position is calculated by source distance according to Archie's formula; and the original gas-water two-phase formation volume density is calculated by source distance based on the cleavage water saturation at each radial position.

[0029] Using the bulk modulus of the dry skeleton, the porosity of the cleavage, and the Gassmann equation, the saturated bulk modulus of the cleavage-restored original gas-water two-phase state is calculated for each source distance and denoted as the corrected saturated bulk modulus for each source distance.

[0030] The corrected P-wave velocity is calculated for each source distance using the corrected saturated bulk modulus, dry skeleton shear modulus, and original gas-water two-phase formation bulk density. The calculation formula is: .

[0031] in, Correct the saturated bulk modulus for each source distance; This represents the volumetric density of the formation in its original gas-water two-phase state.

[0032] The corrected shear wave velocity was calculated using the dry skeleton shear modulus and the original gas-water two-phase formation bulk density. The calculation formula is: .

[0033] Furthermore, the method for calculating the bulk density of the original gas-water two-phase formation based on the source-distance of the cleavage water saturation at each radial position includes: The original gas-water two-phase formation volume density at each source distance was calculated based on the formation volume density. The calculation methods include: .

[0034] in, For cleavage porosity, To constrain the initial water saturation of the cleavage, The original formation water density, For the water saturation of the cutting surface, The density of the drilling fluid filtrate. This represents the original formation gas density.

[0035] Based on the same inventive concept, the present invention also provides a coal reservoir in-situ stress prediction system, the system comprising: The well logging data acquisition module is used to configure the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage; after the source spacing is configured, array acoustic logging is performed to obtain the P-wave velocity and S-wave velocity at each source spacing; density logging is performed to obtain the formation volume density; and nuclear magnetic resonance logging is performed to obtain the cleavage porosity.

[0036] The intrusion zone outer boundary identification module is used to calculate the P-wave and S-wave velocity ratio per source distance and select the intrusion zone outer boundary. The P-wave velocity and S-wave velocity corresponding to the source distance exceeding the intrusion zone outer boundary are recorded as far-well P-wave velocity and far-well S-wave velocity, and those falling within the intrusion zone outer boundary are recorded as near-well P-wave velocity and near-well S-wave velocity.

[0037] The skeleton modulus inversion module is used to invert the skeleton bulk modulus and skeleton shear modulus based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity.

[0038] The fluid replacement correction module is used to replace the near-wellbore P-wave velocity and near-wellbore S-wave velocity with fluid based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the cleavage porosity to obtain the corrected P-wave velocity, the corrected S-wave velocity, and the original gas-water two-phase formation bulk density.

[0039] The geostress prediction module is used to calculate elastic parameters based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, and then substitute them into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

[0040] (3) Beneficial effects Compared with the prior art, the beneficial effects of the present invention are: 1. By configuring the source distance and distinguishing the wave velocities inside and outside the invasion zone, the bulk modulus and shear modulus of the dry skeleton are inverted using far-well information. The corrected P-wave velocity and corrected S-wave velocity are obtained by fluid substitution of the near-well wave velocity, which reduces the near-well acoustic distortion caused by drilling fluid filtrate intrusion and effectively improves the accuracy of geostress prediction results.

[0041] 2. The overall invasion depth is allocated to the face cleavage direction and the end cleavage direction by the density of face cleavage and the density of end cleavage, and the azimuth array sonic logging source distance is configured respectively. Combined with the bimodal decomposition of natural gamma logging curves, the interference interval of interbedded rock lithology is identified. The velocity ratio radial gradient sequence is processed by sliding average with the cleavage filling interference suppression window width to accurately locate the outer boundary of the invasion zone. This further reduces the differences in invasion due to cleavage anisotropy, interbedded rock lithology and local interference from cleavage filling. Attached Figure Description

[0042] Figure 1 This is a flowchart of a method for predicting geostress in a coal reservoir according to Embodiment 1 of the present invention; Figure 2 This is a schematic diagram of the module composition of a coal reservoir geostress prediction system according to Embodiment 2 of the present invention. Detailed Implementation

[0043] 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. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0044] Before providing examples, it is necessary to describe the application scenarios of this invention. This invention, through the collaborative processing of multi-source azimuth array acoustic logging and various logging data, provides a continuously distributed corrected acoustic velocity and geostress prediction profile along the wellbore, serving the design of hydraulic fracturing schemes, the evaluation of fracturing fracture initiation and propagation directions, and the wellbore stability analysis of coalbed methane wells.

[0045] Example 1: As Figure 1 As shown in the figure, this embodiment provides a method for predicting geostress in coal reservoirs, the method comprising: S1. Based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage, configure the array acoustic logging source spacing; after configuring the source spacing, perform array acoustic logging to obtain the P-wave velocity and S-wave velocity at each source spacing; perform density logging to obtain the formation volume density; and perform nuclear magnetic resonance logging to obtain the cleavage porosity.

[0046] S2. Calculate the P-wave and S-wave velocity ratios for each source distance and select the outer boundary of the invasion zone. The P-wave and S-wave velocities corresponding to the source distances exceeding the outer boundary of the invasion zone are denoted as far-well P-wave velocity and far-well S-wave velocity, while those falling within the outer boundary of the invasion zone are denoted as near-well P-wave velocity and near-well S-wave velocity.

[0047] S3. Based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity, invert the bulk modulus and shear modulus of the dry skeleton.

[0048] S4 is used to obtain the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation volume density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the porosity of the cleavage.

[0049] S5. Based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, calculate the elastic parameters and substitute them into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

[0050] For example, taking the No. 3 coal seam of a coalbed methane well as the target coal reservoir, with a depth range of 1120–1138 m, array inductive resistivity logging provides a radial resistivity profile of the No. 3 coal seam, determining the overall invasion depth to be approximately 0.52 m. Imaging logging identifies the dominant azimuth angle of face cleavage at 72° east of north, the face cleavage density at 8 cleavages / m, and the end cleavage density at 5 cleavages / m. Invasion depths are allocated according to the cleavage density ratio, with an invasion depth of approximately 0.32 m in the face cleavage direction and approximately 0.20 m in the end cleavage direction. Based on this, six source spacings are configured: 0.5 m, 1.0 m, 1.5 m, 2.0 m, 2.5 m, and 3.0 m, corresponding to acoustic detection depths of approximately 0.15 m, 0.30 m, 0.45 m, 0.60 m, 0.75 m, and 0.90 m, respectively. The correspondence between the source distance and the acoustic detection depth is determined based on the radial detection depth curve calibrated by the azimuth array acoustic logging tool under the acoustic parameters of the target formation (P-wave velocity of about 1800 m / s and S-wave velocity of about 960 m / s). The detection depth is about 30% of the source distance. In practical applications, the calibration curve provided by the tool manufacturer or the finite difference forward modeling results for the target formation parameters can be used instead.

[0051] Array acoustic logging, density logging, and nuclear magnetic resonance logging were performed, yielding the following P-wave and S-wave velocities at various source distances: 0.5m source distance: P-wave velocity 2180 m / s, S-wave velocity 985 m / s; 1.0m source distance: P-wave velocity 2120 m / s, S-wave velocity 978 m / s; 1.5m source distance: P-wave velocity 2040 m / s, S-wave velocity 970 m / s; 2.0m source distance: P-wave velocity 1830 m / s, S-wave velocity 964 m / s; 2.5m source distance: P-wave velocity 1862 m / s, S-wave velocity 961 m / s; 3.0m source distance: P-wave velocity 1830 m / s, S-wave velocity 959 m / s. The formation bulk density obtained from density logging was 1.46 g / cm³. The cleavage porosity was 5.8%.

[0052] The radial sequence of P-wave and S-wave velocity ratios at six source-spacing distances is as follows: 2.21 for 0.5m source-spacing, 2.17 for 1.0m source-spacing, 2.10 for 1.5m source-spacing, 1.91 for 2.0m source-spacing, 1.94 for 2.5m source-spacing, and 1.91 for 3.0m source-spacing. The outer boundary of the intrusion zone is determined between source-spacing distances of 1.5m and 2.0m, corresponding to an acoustic detection depth of approximately 0.52m.

[0053] For source distances of 2.0m, 2.5m, and 3.0m, the detection depth exceeds 0.52m, corresponding to P-wave velocities of 1830m / s, 1862m / s, and 1830m / s, denoted as far-wellbore P-wave velocities, and corresponding to S-wave velocities of 964m / s, 961m / s, and 959m / s, denoted as far-wellbore S-wave velocities. For source distances of 0.5m, 1.0m, and 1.5m, the detection depth falls within 0.52m, corresponding to P-wave velocities of 2180m / s, 2120m / s, and 2040m / s, denoted as near-wellbore P-wave velocities, and corresponding to S-wave velocities of 985m / s, 978m / s, and 970m / s, denoted as near-wellbore S-wave velocities.

[0054] In the nuclear magnetic resonance logging relaxation time spectrum, the ratio of porosity corresponding to immobile fluid to cleavage porosity (5.8%) in the short relaxation component is approximately 0.24, and the initial water saturation of the cleavage-bound fluid is 0.24; calculated according to the Wood equation, it is 0.021 GPa. Taking the average shear wave velocity (961 m / s) and average P-wave velocity (1840 m / s) at three far-well source distances, the dry skeleton shear modulus is obtained as 1.35 GPa and the dry skeleton bulk modulus as 3.09 GPa.

[0055] Near-wellbore cleavage water saturation was calculated using Archie's formula based on source-to-source spacing from array inductive resistivity logging: 0.88 at 0.5m, 0.73 at 1.0m, and 0.55 at 1.5m. Fluid substitution was performed source-to-source spacing, and the corrected saturated bulk modulus for all three near-wellbore spacings converged to 3.14 GPa. The corrected formation bulk density for the original gas-water two-phase state was: 1.424 g / cm³ at 0.5m, 1.432 g / cm³ at 1.0m, and 1.442 g / cm³ at 1.5m. Corrected P-wave velocities were calculated source-to-source spacing as: 1862 m / s, 1857 m / s, and 1850 m / s; corrected S-wave velocities were: 974 m / s, 971 m / s, and 968 m / s.

[0056] Calculate the ground stress using a source distance of 1.5m as an example. Correct for P-wave velocity. 1850m / s, corrected shear wave velocity 968 m / s and the original gas-water two-phase formation bulk density Calculate the elastic parameters, Poisson's ratio, based on 1.442 g / cm³. = ≈0.31, Young's modulus = ≈3.54 GPa; Substituting Poisson's ratio of 0.31 and Young's modulus of 3.54 GPa into the porosity elastic horizontal in-situ stress prediction model: ; in, The minimum horizontal principal stress; This represents the maximum horizontal principal stress. These are the vertical principal stresses; This represents the effective stress coefficient of Biot. Formation pore pressure; This represents the minimum horizontal tectonic strain. This represents the maximum horizontal tectonic strain.

[0057] The vertical principal stress is estimated from the gravity of the overlying strata. =2500×9.8×1129 / 106≈27.7MPa; where, the average density of the overlying strata is... Take 2500 kg / m³, =1129m is the depth of the midpoint of coal seam No. 3.

[0058] Biot's effective stress coefficient is determined by the bulk modulus of the skeleton. Compared with the preset bulk modulus of coal and rock mineral framework Direct calculation: =1 3.09 / 4.8=0.36. The method for obtaining the pre-set bulk modulus of the coal and petrological skeleton will be described in subsequent embodiments.

[0059] Formation pore pressure =11.9 MPa. Minimum horizontal tectonic strain. =1.5×10 -4 Maximum horizontal stratigraphic strain =9.0×10 -4 All data are derived from regional geological structure analysis.

[0060] Substituting each factor into the calculation, the predicted geostress is the vertical principal stress. Approximately 27.7 MPa, maximum horizontal principal stress The minimum horizontal principal stress is approximately 18.5 MPa. It is approximately 16.5 MPa.

[0061] The method for configuring the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage includes: Based on imaging logging, the dominant azimuth angle of the surface cleavage, surface cleavage density, and end cleavage density of the target coal reservoir are identified.

[0062] The overall invasion depth of the target coal reservoir is determined based on the radial resistivity profile obtained from array inductive resistivity logging; the overall invasion depth is allocated to the face cleavage direction and the end cleavage direction based on the face cleavage density and the end cleavage density, thus obtaining the invasion depth in the face cleavage direction and the invasion depth in the end cleavage direction.

[0063] A combination of source distances for surface cleavage azimuth is configured at the azimuth corresponding to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the surface cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the surface cleavage direction. A combination of source distances for end cleavage azimuth is configured at the azimuth perpendicular to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the end cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the end cleavage direction. The combination of source distances for surface cleavage azimuth and end cleavage azimuth is used as the configuration source distance for the azimuth array sonic logging tool.

[0064] For example, the orientation of cleavage traces in the 1120–1138 m depth range was statistically analyzed using imaging logging. The peak direction of the main lobe of the rose diagram was 72° east of north, which was determined as the dominant azimuth of the surface cleavage. The number of cleavage lines per unit depth was counted along the 72° east of north azimuth and in the azimuth perpendicular to the 72° east of north azimuth (162° east of north azimuth). The surface cleavage density was 8 lines / m, and the end cleavage density was 5 lines / m.

[0065] After the cable logging tool string was lowered into the well, the array inductive resistivity logging provided real-time radial resistivity profiles for the No. 3 coal seam: 7 Ω·m for shallow probes (approximately 0.10 m depth), 11 Ω·m for medium probes (approximately 0.30 m depth), 37 Ω·m for deep probes (approximately 0.60 m depth), and 42 Ω·m for even deeper probes (approximately 0.90 m depth). The resistivity continuously increases from shallow to deep, and the gradient tends to flatten out at approximately 0.52 m, approaching the background value for deeper probes, thus determining the overall intrusion depth to be approximately 0.52 m.

[0066] The sum of the density of face cleavage (8 cleavages / m) and end cleavage (5 cleavages / m) is 13 cleavages / m. The overall penetration depth is distributed in two directions based on the ratio of cleavage densities. Penetration depth in the face cleavage direction = 0.52 × 8 / 13 ≈ 0.32 m, and penetration depth in the end cleavage direction = 0.52 × 5 / 13 = 0.20 m. The permeability of face cleavage and end cleavage is significantly anisotropic. The penetration depth of drilling fluid filtrate along the face cleavage direction is greater than that along the end cleavage direction. If a uniform overall penetration depth is used to configure the source distance in all directions, the far-end source distance in the face cleavage direction may still fall within the penetration zone, while the near-end source distance in the end cleavage direction may have already exceeded the penetration zone, leading to directional deviations in the subsequent identification of the outer boundary of the penetration zone.

[0067] At an azimuth of 72° east of north, the intrusion depth along the surface cleavage direction is 0.32m. A minimum source distance of 0.5m is selected to ensure a detection depth of approximately 0.15m falls within 0.32m, while a maximum source distance of 3.0m ensures a detection depth of approximately 0.90m exceeds 0.32m. The possible source distance combinations along the surface cleavage azimuth are 0.5m, 1.0m, 1.5m, 2.0m, 2.5m, and 3.0m. At an azimuth of 162° east of north, the intrusion depth along the end cleavage direction is 0.20m. A minimum source distance of 0.5m is selected to ensure a detection depth of approximately 0.15m falls within 0.20m, while a maximum source distance of 2.5m ensures a detection depth of approximately 0.75m significantly exceeds 0.20m. The possible source distance combinations along the end cleavage azimuth are 0.5m, 1.0m, 1.5m, 2.0m, and 2.5m. The two sets of source distances mentioned above are combined to form the source distance configuration for the azimuth array acoustic logging tool at two azimuths: 72° east of north and 162° east of north.

[0068] The method for obtaining the porosity of the cleavage includes: Nuclear magnetic resonance logging was performed; the nuclear magnetic resonance porosity spectrum was divided into short relaxation components and long relaxation components according to the relaxation time distribution. The porosity corresponding to the short relaxation component was taken as the matrix micropore porosity, and the difference between the total porosity calculated from the formation volume density obtained by density logging and the matrix micropore porosity was taken as the cleavage porosity.

[0069] For example, the relaxation time spectrum obtained from nuclear magnetic resonance logging in coal seam No. 3 exhibits a bimodal distribution. The peak of the short relaxation component is located at approximately 0.8 ms, while the peak of the long relaxation component is located at approximately 120 ms. The coal matrix micropores have extremely small pore sizes, resulting in a large contact area between the fluid and the solid surface and a fast relaxation rate, corresponding to the short relaxation component. The cleavage pores have relatively large pore sizes, resulting in a slow fluid relaxation rate, corresponding to the long relaxation component. A distinct valley appears in the relaxation time spectrum within the 8–12 ms range. The relaxation time spectrum value corresponding to the valley is taken as 10 ms as the cutoff value for separating the two components. The cutoff value is determined by the position of the valley between the two peaks. The nuclear magnetic resonance porosity corresponding to the integral of the short relaxation component is 7%, denoted as the matrix micropore porosity.

[0070] The porosity of the long relaxation component of nuclear magnetic resonance (NMR) is not directly used as the cleavage porosity because NMR logging suffers from incomplete signal acquisition at extremely short relaxation times, resulting in a systematically lower total porosity calculated from NMR than from density logging. This difference is primarily reflected in the measurement of cleavage porosity. Density logging is not sensitive to the type of fluid in the pores, making its calculated total porosity more reliable. Subtracting the micropore porosity of the NMR matrix from the total porosity from the density logging allows the error from incomplete NMR signal acquisition to be eliminated from the cleavage porosity measurement. The formation bulk density obtained from density logging is 1.46 g / cm³. Combining this with the coal and petrographic skeleton density calibrated from the elemental capture spectrum logging of coal seam No. 3 (1.52 g / cm³) and the drilling fluid filtrate density (1.05 g / cm³), the calculated total porosity is: (1.52 g / cm³) = (1.52 g / cm³) / (1.05 g / cm³). 1.46) / (1.52 (1.05) = 0.06 / 0.47 ≈ 12.8%. Cleavage porosity = 12.8% 7% = 5.8%.

[0071] The method for calculating the P-wave and S-wave velocity ratio per source distance and selecting the outer boundary of the intrusion zone includes: The ratio of P-wave velocity to S-wave velocity is calculated per source distance, forming a radial sequence of P-wave and S-wave velocity ratios arranged from near to far with the acoustic detection depth. Interstitial rock interference zones in the target coal reservoir are identified based on the natural gamma logging curve. The depth points corresponding to the interstitial rock interference zones are removed from the radial sequence of P-wave and S-wave velocity ratios. The arithmetic mean of the P-wave and S-wave velocity ratios of the remaining pure coal depth points is taken per source distance to obtain the effective radial sequence of P-wave and S-wave velocity ratios.

[0072] Along the acoustic detection depth from near to far, the velocity ratio radial gradient sequence is formed by subtracting the velocity ratio of the next adjacent source from the previous source distance's P-wave velocity ratio. The velocity ratio radial gradient sequence is then smoothed by moving average to obtain a smoothed velocity ratio radial gradient sequence.

[0073] The sign of each difference value in the radial gradient sequence of the smooth velocity ratio is examined one by one along the direction from far to near the sound wave detection depth. The starting depth position where the difference value turns from non-positive to positive is determined as the outer boundary of the intrusion zone.

[0074] For example, the P-wave and S-wave velocity ratios were calculated at each depth point and source distance within the depth range of 1120–1138 m. The background readings of the natural gamma ray logging curves were 22–28 API in the 1120–1138 m range, rising to 72 API in the 1131–1133 m range, thus identifying the 1131–1133 m range as an area affected by interbedded rock lithology. This interbedded rock lithology increased the P-wave and S-wave velocity ratios at each source distance and depth point within the 1131–1133 m range to 2.35–2.45. Including representative values ​​for each source distance in the calculation would compress the difference in P-wave and S-wave velocity ratios between near-wellbore and far-wellbore areas, affecting the identification of the outer boundary of the invasion zone. After removing the depth points from 1131 to 1133 m from the radial sequence of P-wave and S-wave velocity ratios, the effective radial sequence of P-wave and S-wave velocity ratios was obtained by taking the arithmetic mean of the P-wave and S-wave velocity ratios at each depth point in the remaining pure coal section (1120–1131 m and 1133–1138 m, a total depth of 16 m). The values ​​were: 0.5 m source distance 2.21, 1.0 m source distance 2.17, 1.5 m source distance 2.10, 2.0 m source distance 1.91, 2.5 m source distance 1.94, and 3.0 m source distance 1.91.

[0075] The difference between adjacent source distance pairs is calculated from near to far for the effective P-wave and S-wave velocity ratio radial sequence, forming a velocity ratio radial gradient sequence: +0.04 between 0.5m and 1.0m, +0.07 between 1.0m and 1.5m, +0.19 between 1.5m and 2.0m, and +0.19 between 2.0m and 2.5m. 0.03, +0.03 between 2.5m and 3.0m. At a source distance of 2.5m, a localized calcite-filled cleavage is present at a depth of approximately 0.75m. This calcite infilling locally increases the P-wave velocity at the 2.5m source distance, raising the P-wave / S-wave velocity ratio to 1.94. At a source distance of 2.0m, the cleavage is not reached at a depth of 0.60m. At a source distance of 3.0m, although the cleavage is covered at a depth of 0.90m, it is homogenized and diluted within the larger probe volume; the P-wave / S-wave velocity ratio is 1.91 for both locations. The isolated peak at 2.5m exhibits a bipolar pair in the gradient sequence: the gradient near the peak (between 2.0m and 2.5m) is... The gradient at the far end of the peak (between 2.5m and 3.0m) is +0.03. After smoothing by moving average, the positive and negative values ​​of the bipolar gradient cancel each other out, and the smoothed values ​​at both positions are precisely reduced to 0. The smoothing rate ratio of the radial gradient sequence is as follows: 0 between 2.5m and 3.0m, 0 between 2.0m and 2.5m, approximately +0.08 between 1.5m and 2.0m, approximately +0.10 between 1.0m and 1.5m, and approximately +0.06 between 0.5m and 1.0m.

[0076] The signs of each difference value in the smooth velocity ratio radial gradient sequence were examined sequentially along the acoustic detection depth from far to near: 0 between 2.5m and 3.0m, 0 between 2.0m and 2.5m, +0.08 between 1.5m and 2.0m, and subsequently remained positive between 1.0m and 1.5m, and between 0.5m and 1.0m. The starting position where the difference value continuously turned from zero to positive was between the source distance of 1.5m and 2.0m, corresponding to an acoustic detection depth of approximately 0.52m, which was identified as the outer boundary of the intrusion zone.

[0077] The method for identifying interbedded rocky lithological interference zones in target coal reservoirs based on natural gamma logging curves includes: Natural gamma logging curves are extracted along the depth direction of the target coal reservoir; the numerical frequency distribution of the natural gamma logging curves within the depth range of the target coal reservoir is statistically analyzed; bimodal decomposition is performed on the numerical frequency distribution, and the natural gamma value corresponding to the frequency trough between the two peaks is determined as the natural gamma boundary value between the pure coal section and the interbedded gangue section.

[0078] Based on the natural gamma ray threshold, the continuous depth intervals in the natural gamma logging curves that consistently exceed the natural gamma ray threshold are identified as candidate depth intervals for interbedded rock. All depth points whose natural gamma logging curve values ​​do not exceed the natural gamma ray threshold are designated as pure coal reference segments. The shear wave velocities at each source distance are extracted from each depth point in the pure coal reference segment, and the arithmetic mean of all pure coal depth points and all source distances is taken as the pure coal reference shear wave velocity. The arithmetic mean of the shear wave velocities at each source distance within the candidate depth intervals for interbedded rock is calculated, and the continuous depth intervals with arithmetic mean higher than the pure coal reference shear wave velocity are identified as interbedded rock lithological interference intervals.

[0079] For example, in the No. 3 coal seam at a depth of 1120–1138 m, the natural gamma logging curve readings are stably distributed between 22–28 API in the 1120–1131 m and 1133–1138 m ranges, rising to approximately 72 API in a continuous interval of 1131–1133 m. Statistical analysis of the frequency distribution of natural gamma readings within the 1120–1138 m depth range shows a bimodal distribution: the low-value peak is concentrated between 22–28 API, and the high-value peak is concentrated between 68–76 API, with a distinct frequency trough at approximately 45 API between the low and high values. Coal has extremely weak natural radioactivity, resulting in low natural gamma readings; the content of radioactive elements such as potassium, uranium, and thorium in the mudstone interbeds is much higher than in coal, with readings concentrated in the high-value range. The two lithologies form their own independent numerical clusters in terms of frequency distribution, and the frequency troughs objectively mark the lithological boundary between the pure coal section and the interbedded rock section. The 45API value, corresponding to the frequency trough between the low and high peaks, is taken as the natural gamma boundary between the pure coal section and the interbedded gangue section.

[0080] The continuous depth range where the natural gamma reading consistently exceeds 45 API is 1131–1133 m, and this range is identified as the candidate depth range for interbedded coal. The 1120–1131 m and 1133–1138 m (a total depth range of 16 m) depths where the reading does not exceed 45 API are used as pure coal reference ranges. The shear wave velocity at each depth point is extracted for all source distances, and the arithmetic mean of all pure coal depth points and all source distances is approximately 962 m / s, which is denoted as the pure coal reference shear wave velocity.

[0081] Within the candidate depth range of 1131–1133 m for interbedded rock, the arithmetic mean of the shear wave velocity at each depth and source distance is approximately 1068 m / s, significantly higher than the reference shear wave velocity of 962 m / s for pure coal. The shear wave velocity is not sensitive to changes in pore fluid; saturation of drilling fluid filtrate by near-wellbore cleavage does not cause an increase in shear wave velocity. The elastic modulus of the argillaceous interbedded rock mineral skeleton (mainly composed of clay minerals) is higher than that of the coal rock skeleton, and the increased shear stiffness of the skeleton results in a higher shear wave velocity than in the pure coal section. The higher shear wave velocity at 1131–1133 m compared to the reference shear wave velocity for pure coal is attributed to changes in skeleton lithology; therefore, 1131–1133 m is identified as the interbedded rock lithological interference range.

[0082] The method for applying a moving average smoothing process to the velocity ratio radial gradient sequence to obtain a smoothed velocity ratio radial gradient sequence includes: For the velocity ratio radial gradient sequence, count the number of source distance spans of consecutive positive gradient segments to form a positive gradient segment source distance span distribution; from the positive gradient segment source distance span distribution, identify and separate cleavage-filled positive gradient segments that span only a single source distance and fluid displacement-type positive gradient segments that span two or more source distances.

[0083] The minimum number of odd source distances that completely cover the positive gradient segment of the cleavage filling pattern is used to determine the width of the cleavage filling interference suppression window. The velocity ratio radial gradient sequence is processed by moving average using the width of the cleavage filling interference suppression window. The arithmetic mean of the gradient values ​​in each window is calculated, and the arithmetic mean is assigned to the source distance position of the window center to replace the original gradient value, thus obtaining a smooth velocity ratio radial gradient sequence.

[0084] For example, the gradient values ​​of the five adjacent source-distance pairs in the velocity-to-radial gradient sequence are as follows: +0.04 between 0.5m and 1.0m, +0.07 between 1.0m and 1.5m, +0.19 between 1.5m and 2.0m, and +0.19 between 2.0m and 2.5m. 0.03, between 2.5m and 3.0m +0.03.

[0085] In the statistical velocity ratio radial gradient sequence, the continuous positive gradient segment consists of three consecutive positive gradient values ​​between 0.5m and 1.0m and between 1.5m and 2.0m, forming a continuous positive gradient segment with a span of 3; the gradient value between 2.0m and 2.5m is... At 0.03, the continuous positive value segment is interrupted; the gradient value between 2.5m and 3.0m +0.03 constitutes an isolated positive gradient segment with a span of 1. The source distance span distribution of the positive gradient segments is as follows: one segment with a span of 3 and another segment with a span of 1.

[0086] The continuous positive gradient segment with a span of 3 is a fluid displacement type positive gradient segment. It is formed because drilling fluid filtrate intrudes into near-wellbore cleavages, causing the fluid within the cleavages to change from a gas-saturated state to a water-saturated state. This leads to a continuous increase in the P-wave velocity ratio at multiple source distances within the intrusion zone, thus forming a fluid displacement type positive gradient segment. The isolated positive gradient segment with a span of 1 is a cleavage-filled type positive gradient segment. It is formed because a locally calcite-filled cleavage is present at a detection depth of approximately 0.75m at a 2.5m source distance. The calcite filling causes a local increase in P-wave velocity at this source distance, forming an isolated peak with a apex at 2.5m in the radial sequence of P-wave and S-wave velocity ratios. At a 2.0m source distance (detection depth 0.60m, not reaching the cleavage-filled location) and a 3.0m source distance (detection depth 0.90m, the cleavage-filled location is diluted in a larger detection volume), the P-wave and S-wave velocity ratios are both the background value of 1.91. Isolated peaks are represented in the gradient sequence as adjacent bipolar pairs: the near-end gradient (between 2.0m and 2.5m) is... The gradient at the far end (between 2.5m and 3.0m) is +0.03. The positive value spans a single source distance, which is the positive gradient segment of the identified cleavage filling type.

[0087] The positive gradient segment of the cleavage filling pattern spans a single source distance. The minimum odd number required to completely cover the span of a single source distance is 3. Therefore, the width of the cleavage filling interference suppression window is determined to be 3 source distances.

[0088] A moving average of the velocity ratio radial gradient sequence was performed with a window width of three source distances. At the near-end boundary (between 0.5m and 1.0m), extending only one position further outward, the mean of two points was (0.04 + 0.07) / 2 ≈ +0.06; at the inner location between 1.0m and 1.5m, the mean of three points was (0.04 + 0.07 + 0.19) / 3 = +0.10; and at the inner location between 1.5m and 2.0m, the mean of three points was (0.07 + 0.19 + (…). 0.03)) / 3≈+0.08; Gradient value between 2.0m and 2.5m ( 0.03) was identified as a bipolar gradient pair of cleavage-filler type before smoothing (between 2.0m and 2.5m: The near-end negative value member (between 0.03, 2.5m and 3.0m: +0.03) is considered. If a full 3-point window extending proximally is applied, the adjacent positions on the intrusion side (between 1.5m and 2.0m: +0.19) will be included in the mean calculation, causing the bipolar pairs to not cancel each other out precisely. To ensure that the positive and negative values ​​within the identified cleavage-filled bipolar gradient pairs cancel each other out precisely in the smoothing calculation, only the 2-point mean of the bipolar pair's near-end negative value member and its only adjacent position is taken. 0.03) + (+0.03) / 2 = 0; At the far boundary position (between 2.5m and 3.0m), there are no more distant neighboring points. Take the average of the two points between itself and the only available neighbor: (+0.03) + ( 0.03) / 2=0.

[0089] The smoothed gradient values ​​at each location are: approximately +0.06 between 0.5m and 1.0m, approximately +0.10 between 1.0m and 1.5m, approximately +0.08 between 1.5m and 2.0m, 0 between 2.0m and 2.5m, and 0 between 2.5m and 3.0m, resulting in a smoothed radial gradient sequence. The bipolar gradient is aligned to +0.03 and... The 0.03 values ​​cancel each other out precisely in the calculation of their respective window means, and the positive gradient contribution of the cleavage filling model drops to 0, thus not constituting a sustained positive starting signal in subsequent scans from far to near.

[0090] The method for inverting the bulk modulus and shear modulus of the skeleton based on far-wellbore P-wave velocity, far-wellbore S-wave velocity, formation bulk density, and cleavage porosity includes: Based on the ratio of porosity to cleavage porosity corresponding to immovable fluid in the short relaxation component of the nuclear magnetic resonance logging relaxation time spectrum, the initial water saturation of cleavage-bound cleavage in the target coal reservoir is determined; the cleavage-bound water porosity is calculated by the product of cleavage porosity and initial water saturation of cleavage-bound cleavage; and the cleavage gas porosity is calculated by the difference between cleavage porosity and cleavage-bound water porosity.

[0091] The equivalent bulk modulus of the cleavage-bound initial water saturation and the original formation water bulk modulus, and the ratio of cleavage porosity to the pre-set original formation gas bulk modulus, are weighted and calculated according to the Wood equation to obtain the equivalent bulk modulus of the cleavage-bound fluid.

[0092] The dry skeleton shear modulus is calculated based on the far-well shear wave velocity and formation bulk density. The calculation formula is: .

[0093] in, The density of the formation is the volumetric density. The distance is the shear wave velocity.

[0094] The saturated bulk modulus of the far-wellbore is calculated based on the far-well P-wave velocity, far-well S-wave velocity, and formation bulk density. The formula for calculating the saturated bulk modulus of the far-wellbore is as follows: .

[0095] in, The distance is the longitudinal wave velocity.

[0096] By substituting the far-well saturated bulk modulus, cleavage porosity, the pre-defined bulk modulus of the coal and petrified mineral skeleton, and the equivalent bulk modulus of the cleavage-mixed fluid into the Gassmann equation, the dry skeleton bulk modulus is obtained by back-calculation.

[0097] For example, the outer boundary of the intrusion zone is located at a detection depth of approximately 0.52 m. The acoustic detection depths at the three source distances of 2.0 m, 2.5 m, and 3.0 m all exceed 0.52 m, corresponding to P-wave velocities of 1830 m / s, 1862 m / s, and 1830 m / s and S-wave velocities of 964 m / s, 961 m / s, and 959 m / s, denoted as the far-well P-wave velocity and far-well S-wave velocity.

[0098] The far-well area was not invaded by drilling fluid filtrate, but the cleavage of the No. 3 coal seam was not purely gas-saturated in its original state. After the formation of the cleavage network, it was in long-term contact with formation water. There was original formation water bound by capillary forces at the micropore throats of the cleavage walls and at the cleavage intersections. After natural gas was injected, the displacement pressure difference in the cleavage was less than the capillary force, and this part of the formation water could not be driven out by the gas, constituting the original bound water of the cleavage. The elastic compressibility of the formation water participated in the overall response of the cleavage pore fluid. If the pre-set original formation gas bulk modulus is directly used as the cleavage pore fluid modulus and substituted into the Gassmann equation, the contribution of this part of the formation water (whose bulk modulus is much higher than that of gas) will be completely ignored, the equivalent bulk modulus of the cleavage pore fluid will be underestimated, resulting in a lower bulk modulus of the dry skeleton obtained by inversion, and a directional deviation will occur when performing fluid replacement on the near-well water saturation rate.

[0099] In the nuclear magnetic resonance logging relaxation time spectrum, the porosity corresponding to the immobile fluid in the short relaxation component is 1.39%, which is approximately 0.24 compared to the cleavage porosity of 5.8%. This yields the initial water saturation of the cleavage-bound fluid. =0.24.

[0100] The porosity of water bound by cleavage is 5.8% × 0.24 = 1.39%; the porosity of gas cleavage is 5.8%. 1.39% = 4.41%; the ratio of cleavage porosity to cleavage porosity = 4.41% / 5.8% ≈ 0.76.

[0101] The original formation temperature of coal seam No. 3 was approximately 60°C, and the formation pressure was approximately 11.9 MPa. The corresponding original formation water bulk modulus under these conditions was... =2.25 GPa, preset original formation gas bulk modulus. Weighted according to Wood's equation: ; in, For cleavage porosity; The equivalent bulk modulus of the cleavage-mixed fluid; For cleavage porosity, To predetermine the initial formation gas bulk modulus, it was calculated using the Batzle-Wang equation based on formation temperature, formation pressure, and coalbed methane composition. The methane content of coalbed methane in coal seam No. 3 is approximately 97%, and under formation temperature of 60°C and formation pressure of 11.9 MPa, the calculated value is 0.016 GPa.

[0102] Solving The equivalent bulk modulus of the cleavage pore fluid in the far-well region is 0.021 GPa.

[0103] It should be noted that the formula for calculating the shear modulus of the skeleton is... In this process, the substituted formation volume density should, strictly speaking, be consistent with the fluid state at the radial detection location of the shear wave velocity used. The shear wave velocity used for inversion is taken from the far-wellhead source distance, reflecting the original gas-water two-phase state of the formation outside the invasion zone. The corresponding substituted formation volume density should also be the density under this original state, not the near-wellbore invasion state density reflected by density logging (density logging detection depth is approximately 0.1m, only reflecting the near-wellbore invasion area). For ease of demonstration and calculation, this embodiment uses the formation volume density obtained from density logging, 1.46 g / cm³, as the substituted density. In practical applications, the original formation volume density can be obtained by determining the coal matrix density and cleavage porosity through core experiments, combined with the original gas-water saturation reconstructed from the interpretation of deep lateral resistivity logging curves, to eliminate the fluid state mismatch error between the invasion state density and the original state velocity.

[0104] The average shear wave velocity at the source distance of the three far-wellbore samples is approximately (964 + 961 + 959) / 3 ≈ 961 m / s. The formation bulk density is 1.46 g / cm³. Substituting these values ​​into the formula for calculating the dry skeleton shear modulus, we obtain the dry skeleton shear modulus. It is approximately 1.35 GPa.

[0105] The average P-wave velocity at the source distance of the three far wells is approximately (1830 + 1862 + 1830) / 3 ≈ 1840 m / s. Substituting this into the formula for calculating the saturated bulk modulus of the far well, we obtain the saturated bulk modulus of the far well. The value is 3.14 GPa.

[0106] The wellbore saturated bulk modulus is 3.14 GPa, the cleavage porosity is 0.058, and the pre-set bulk modulus of the coal and petrological framework is... and equivalent bulk modulus of cleavage-mixed fluid Substituting 0.021 GPa into the Gassmann equation, the bulk modulus of the skeleton is calculated inversely: ; The pre-defined bulk modulus of the coal and petrological framework was obtained by conducting ultrasonic experiments on core samples from coal seam No. 3 in a dry state at 60°C and a confining pressure of 70 MPa (far exceeding the cleavage closure pressure, as both cleavage and microfractures were closed). The measured longitudinal wave velocity was 2390 m / s, the transverse wave velocity was 1390 m / s, and the core framework density was 1.52 g / cm³. The bulk modulus of the coal and petrological framework was then calculated. =1520×(2390² (4 / 3 × 1390²) ≈ 4.8 GPa. Solving for the bulk modulus of the skeleton yields... =3.09 GPa.

[0107] The method for obtaining corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the cleavage porosity includes: Based on the resistivity values ​​at each depth obtained from array inductive resistivity logging, and constrained by cleavage porosity, the cleavage water saturation at each radial position is calculated by source distance according to Archie's formula; and the original gas-water two-phase formation volume density is calculated by source distance based on the cleavage water saturation at each radial position.

[0108] Using the bulk modulus of the dry skeleton, the porosity of the cleavage, and the Gassmann equation, the saturated bulk modulus of the cleavage-restored original gas-water two-phase state is calculated for each source distance and denoted as the corrected saturated bulk modulus for each source distance.

[0109] The corrected P-wave velocity is calculated for each source distance using the corrected saturated bulk modulus, dry skeleton shear modulus, and original gas-water two-phase formation bulk density. The calculation formula is: .

[0110] in, Correct the saturated bulk modulus for each source distance; This represents the volumetric density of the formation in its original gas-water two-phase state.

[0111] The corrected shear wave velocity was calculated using the dry skeleton shear modulus and the original gas-water two-phase formation bulk density. The calculation formula is: .

[0112] For example, the drilling fluid continuously permeates into the cleavage near the wellbore under a positive pressure differential of approximately 1.0 MPa. The displacement pressure differential is greatest near the wellbore, and the filtrate has displaced most of the free gas in the cleavage. As the radial distance increases, the displacement pressure differential decreases, and the amount of filtrate penetrating gradually decreases. The water saturation of the cleavage forms a radial gradient distribution that continuously decreases from the wellbore outward in the near-wellbore region. The acoustic detection depths at the three near-wellbore source distances of 0.5 m, 1.0 m, and 1.5 m are approximately 0.15 m, 0.30 m, and 0.45 m, respectively, falling within different invasion depth ranges. If a fixed water saturation level is used to represent all near-wellbore locations, the actual fluid state at each source distance cannot be accurately reflected, resulting in errors in the calculation of the fluid modulus at each source distance in different directions.

[0113] Resistivity at detection depths corresponding to source spacings of 0.5m, 1.0m, and 1.5m using array inductive resistivity logging. The resistivity values ​​are 7.4 Ω·m, 10.7 Ω·m, and 19.0 Ω·m, respectively. Coal cleavage uses fracture surfaces as the main conduction channels, and the conductive path is relatively regular. Using the Archie lithology coefficient a=1, cementation index m=1.5, and saturation index n=2, the formation water resistivity of coal seam No. 3 at approximately 60°C is... =0.08Ω·m, cleavage porosity =0.058, where =0.058 1.5 ≈0.0140. The cleavage water saturation is calculated by back-calculating the source-by-source distance using Archie's formula: ; The cleavage water saturation was obtained at a source distance of 0.5m. Approximately 0.88; water saturation at a source distance of 1.0m. Approximately 0.73; water saturation at a source distance of 1.5m. Approximately 0.55. A radial distribution sequence of cleavage water saturation is formed in the near-wellbore area: 0.88 at 0.5m, 0.73 at 1.0m, and 0.55 at 1.5m.

[0114] The water saturation of each source distance will be used for subsequent correction calculations of the original gas-water two-phase formation volume density at each source distance; the corrected saturated bulk modulus required for fluid replacement is directly determined by the dry skeleton modulus and the original gas-water two-phase mixed fluid modulus at each source distance.

[0115] The bulk modulus of the dry skeleton is 3.09 GPa, the cleavage porosity is 0.058, and the equivalent bulk modulus of the cleavage-mixed fluid is... Substituting 0.021 GPa into the Gassmann equation, the corrected saturated bulk modulus for each source distance in the original gas-water two-phase state under cleavage recovery is calculated on a source-distance-by-source basis: .

[0116] Solve for the corrected saturated bulk modulus of each source distance Both are 3.14 GPa, consistent with the saturated bulk modulus of the far well obtained directly from the acoustic velocity of the far well.

[0117] The corrected saturated bulk modulus of 3.14 GPa and the dry skeleton shear modulus of 1.35 GPa were compared with the original gas-water two-phase formation bulk density at a source distance of 0.5 m. Approximately 1.424 g / cm³; bulk density of the original gas-water two-phase formation at a source distance of 1.0 m. Approximately 1.432 g / cm³; bulk density of the original gas-water two-phase formation at a source distance of 1.5 m. Approximately 1.442 g / cm³, source-by-source distance is used to correct the P-wave velocity. The calculation formula yields the corrected P-wave velocity at a source distance of 0.5m. Approximately 1862 m / s; corrected P-wave velocity at a source distance of 1.0 m. Approximately 1857 m / s; corrected P-wave velocity at a source distance of 1.5 m. It is approximately 1850 m / s.

[0118] Substituting the shear modulus of the dry skeleton (1.35 GPa) and the original gas-water two-phase formation bulk density at each source distance into the corrected shear wave velocity... The calculation formula yields the corrected shear wave velocity at a source distance of 0.5m. Approximately 974 m / s; corrected shear wave velocity at a source distance of 1.0 m. Approximately 971 m / s; corrected shear wave velocity at a source distance of 1.5 m. It is approximately 968 m / s.

[0119] The method for calculating the bulk density of the original gas-water two-phase formation based on the source-distance source-by-source cleavage water saturation at each radial position includes: The original gas-water two-phase formation volume density at each source distance was calculated based on the formation volume density. The calculation methods include: ; in, For cleavage porosity, To constrain the initial water saturation of the cleavage, The original formation water density, To determine the water saturation of the cutting surface, The density of the drilling fluid filtrate. This represents the original formation gas density.

[0120] For example, the formation bulk density obtained from density logging is 1.46 g / cm³. The cleavage pores are filled with a mixture of drilling fluid filtrate and residual gas, and the degree of filtrate filling varies at each source distance. The sonic velocity has been restored to the original gas-water two-phase state through fluid replacement. If the density still retains the contribution of the invading filtrate, it will not match the actual invading state and needs to be corrected for each source distance according to the actual invading state. For ease of example calculation and demonstration, this embodiment uses the formation bulk density obtained from density logging of 1.46 g / cm³ as the standard for all source distances. Substitute the values ​​into the calculation; in practical applications, the actual measured formation volume density at the radial position corresponding to each source distance should be substituted into the calculation one source distance at a time to obtain a more accurate original gas-water two-phase formation volume density.

[0121] The original formation gas density was determined by the coalbed methane composition obtained from well tests and the formation temperature and pressure conditions. The No. 3 coal seam has a methane content of approximately 97%. According to the Peng-Robinson equation of state, under formation temperature of 60°C and formation pressure... Methane compressibility factor at 11.9 MPa The original formation gas density is approximately 0.90. ;in, The open temperature is 333K; The universal gas constant is 8.314 J / (mol·K); The molar mass of coalbed methane is given. The methane content of coalbed methane in coal seam No. 3 is approximately 97%, with a methane molar mass of 0.016 kg / mol. The content of other components (N2, CO2, etc.) is extremely low, and their influence on the molar mass of the mixed gas is negligible. Therefore, a pure methane molar mass of 0.016 kg / mol is used. The original formation gas density is then obtained. Approximately 0.076 g / cm³ 3 .

[0122] Original formation water density of coal seam No. 3 1.05 g / cm³, compared to the density of drilling fluid filtrate =1.05g / cm³ is equal, the formula simplifies to: ; The water saturation at each source distance is: 0.88 at 0.5m, 0.73 at 1.0m, and 0.55 at 1.5m. =0.24, =0.058, substituting the source distance into the equation, we obtain the gas-saturated formation bulk density at a source distance of 0.5m. Approximately 1.424 g / cm³; gas-saturated formation bulk density at a source distance of 1.0 m. Approximately 1.432 g / cm³; gas-saturated formation bulk density at a source distance of 1.5 m. It is approximately 1.442 g / cm³.

[0123] Example 2: Based on the same inventive concept, such as Figure 2 As shown in the figure, this embodiment also provides a coal reservoir in-situ stress prediction system, the system comprising: The well logging data acquisition module is used to configure the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage; after the source spacing is configured, array acoustic logging is performed to obtain the P-wave velocity and S-wave velocity at each source spacing; density logging is performed to obtain the formation volume density; and nuclear magnetic resonance logging is performed to obtain the cleavage porosity.

[0124] The intrusion zone outer boundary identification module is used to calculate the P-wave and S-wave velocity ratio per source distance and select the intrusion zone outer boundary. The P-wave velocity and S-wave velocity corresponding to the source distance exceeding the intrusion zone outer boundary are recorded as far-well P-wave velocity and far-well S-wave velocity, and those falling within the intrusion zone outer boundary are recorded as near-well P-wave velocity and near-well S-wave velocity.

[0125] The skeleton modulus inversion module is used to invert the skeleton bulk modulus and skeleton shear modulus based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity.

[0126] The fluid replacement correction module is used to replace the near-wellbore P-wave velocity and near-wellbore S-wave velocity with fluid based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the cleavage porosity to obtain the corrected P-wave velocity, the corrected S-wave velocity, and the original gas-water two-phase formation bulk density.

[0127] The geostress prediction module is used to calculate elastic parameters based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, and then substitute them into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

[0128] It should be noted that the specific ways in which each module performs operations in the system described in the above embodiments have been described in detail in the embodiments related to the method, and will not be elaborated here.

[0129] Finally, it should be noted that although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for predicting geostress in coal reservoirs, characterized in that, The method includes: Based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage, the source spacing of the array acoustic logging is configured; after the source spacing is configured, array acoustic logging is performed to obtain the P-wave velocity and S-wave velocity at each source spacing; density logging is performed to obtain the formation volume density; nuclear magnetic resonance logging is performed to obtain the cleavage porosity. The P-wave and S-wave velocity ratios are calculated per source distance, and the outer boundary of the invasion zone is selected. The P-wave and S-wave velocities corresponding to source distances exceeding the outer boundary of the invasion zone are denoted as far-well P-wave velocity and far-well S-wave velocity, and those falling within the outer boundary of the invasion zone are denoted as near-well P-wave velocity and near-well S-wave velocity. Based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density, and cleavage porosity, the bulk modulus and shear modulus of the skeleton are inverted. It is used to obtain the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the porosity of the cleavage. Based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, elastic parameters are calculated and substituted into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

2. The method for predicting in-situ stress in coal reservoirs according to claim 1, characterized in that, The method for configuring the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage includes: Based on imaging logging, the dominant azimuth angle of the surface cleavage, surface cleavage density, and end cleavage density of the target coal reservoir are identified. The overall invasion depth of the target coal reservoir is determined based on the radial resistivity profile obtained from array inductive resistivity logging; the overall invasion depth is allocated to the face cleavage direction and the end cleavage direction based on the face cleavage density and the end cleavage density, thus obtaining the invasion depth in the face cleavage direction and the invasion depth in the end cleavage direction. A combination of source distances for surface cleavage azimuth is configured at the azimuth corresponding to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the surface cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the surface cleavage direction. A combination of source distances for end cleavage azimuth is configured at the azimuth perpendicular to the dominant azimuth of the surface cleavage. The detection depth of the minimum source distance of the combination falls within the intrusion depth of the end cleavage direction, while the detection depth of the maximum source distance exceeds the intrusion depth of the end cleavage direction. The combination of source distances for surface cleavage azimuth and end cleavage azimuth is used as the configuration source distance for the azimuth array sonic logging tool.

3. The method for predicting in-situ stress in coal reservoirs according to claim 1, characterized in that, The method for obtaining the porosity of the cleavage includes: Nuclear magnetic resonance logging was performed; the nuclear magnetic resonance porosity spectrum was divided into short relaxation components and long relaxation components according to the relaxation time distribution. The porosity corresponding to the short relaxation component was taken as the matrix micropore porosity, and the difference between the total porosity calculated from the formation volume density obtained by density logging and the matrix micropore porosity was taken as the cleavage porosity.

4. The method for predicting in-situ stress in coal reservoirs according to claim 1, characterized in that, The method for calculating the P-wave and S-wave velocity ratio per source distance and selecting the outer boundary of the intrusion zone includes: The ratio of P-wave velocity to S-wave velocity is calculated per source distance to form a radial sequence of P-wave and S-wave velocity ratios arranged from near to far with the acoustic detection depth. Interstitial rock interference zones in the target coal reservoir are identified based on the natural gamma logging curve. The depth points corresponding to the interstitial rock interference zones are removed from the radial sequence of P-wave and S-wave velocity ratios. The arithmetic mean of the P-wave and S-wave velocity ratios of the remaining pure coal depth points is taken per source distance to obtain the effective radial sequence of P-wave and S-wave velocity ratios. Along the direction of acoustic detection depth from near to far, the P-wave and S-wave velocity ratios of the previous source distance are subtracted from the P-wave and S-wave velocity ratios of the next adjacent source distance to form a velocity ratio radial gradient sequence; the velocity ratio radial gradient sequence is then smoothed by moving average to obtain a smooth velocity ratio radial gradient sequence. The sign of each difference value in the radial gradient sequence of the smooth velocity ratio is examined one by one along the direction from far to near the sound wave detection depth. The starting depth position where the difference value turns from non-positive to positive is determined as the outer boundary of the intrusion zone.

5. The method for predicting in-situ stress in coal reservoirs according to claim 4, characterized in that, The method for identifying interbedded rocky lithological interference zones in target coal reservoirs based on natural gamma logging curves includes: Natural gamma logging curves are extracted along the depth direction of the target coal reservoir; the numerical frequency distribution of the natural gamma logging curves within the depth range of the target coal reservoir is statistically analyzed; bimodal decomposition is performed on the numerical frequency distribution, and the natural gamma value corresponding to the frequency valley between the two peaks is determined as the natural gamma boundary value between the pure coal section and the interbedded gangue section. Based on the natural gamma ray threshold, the continuous depth intervals in the natural gamma logging curves that consistently exceed the natural gamma ray threshold are identified as candidate depth intervals for interbedded rock. All depth points whose natural gamma logging curve values ​​do not exceed the natural gamma ray threshold are designated as pure coal reference segments. The shear wave velocities at each source distance are extracted from each depth point in the pure coal reference segment, and the arithmetic mean of all pure coal depth points and all source distances is taken as the pure coal reference shear wave velocity. The arithmetic mean of the shear wave velocities at each source distance within the candidate depth intervals for interbedded rock is calculated, and the continuous depth intervals with arithmetic mean higher than the pure coal reference shear wave velocity are identified as interbedded rock lithological interference intervals.

6. The method for predicting in-situ stress in coal reservoirs according to claim 4, characterized in that, The method for applying a moving average smoothing process to the velocity ratio radial gradient sequence to obtain a smoothed velocity ratio radial gradient sequence includes: For the velocity ratio radial gradient sequence, count the number of source distance spans of consecutive positive gradient segments to form a positive gradient segment source distance span distribution; from the positive gradient segment source distance span distribution, identify and separate cleavage-filled positive gradient segments that span only a single source distance and fluid displacement positive gradient segments that span two or more source distances. The minimum number of odd source distances that completely cover the positive gradient segment of the cleavage filling pattern is used to determine the width of the cleavage filling interference suppression window. The velocity ratio radial gradient sequence is processed by moving average using the width of the cleavage filling interference suppression window. The arithmetic mean of the gradient values ​​in each window is calculated, and the arithmetic mean is assigned to the source distance position of the window center to replace the original gradient value, thus obtaining a smooth velocity ratio radial gradient sequence.

7. The method for predicting in-situ stress in coal reservoirs according to claim 1, characterized in that, The method for inverting the bulk modulus and shear modulus of the skeleton based on far-wellbore P-wave velocity, far-wellbore S-wave velocity, formation bulk density, and cleavage porosity includes: Based on the ratio of porosity to cleavage porosity corresponding to immobile fluid in the short relaxation component of the nuclear magnetic resonance logging relaxation time spectrum, the initial water saturation of cleavage-bound cleavage in the target coal reservoir is determined; the cleavage-bound water porosity is calculated by the product of cleavage porosity and initial water saturation of cleavage-bound cleavage; and the cleavage gas porosity is calculated by the difference between cleavage porosity and cleavage-bound water porosity. The equivalent bulk modulus of the cleavage-bound initial water saturation and the original formation water bulk modulus, and the ratio of cleavage porosity to cleavage porosity and the pre-set original formation gas bulk modulus are used to calculate the equivalent bulk modulus of the cleavage-bound fluid using the Wood equation. The dry skeleton shear modulus is calculated based on the far-well shear wave velocity and formation bulk density. The calculation formula is: ； in, The density of the formation is the volumetric density. The velocity of the shear wave at the far well; The saturated bulk modulus of the far-wellbore is calculated based on the far-wellbore P-wave velocity, far-wellbore S-wave velocity, and formation bulk density. The calculation formula is: ； in, For the longitudinal wave velocity at the far end of the well; By substituting the far-well saturated bulk modulus, cleavage porosity, the pre-defined bulk modulus of the coal and petrified mineral skeleton, and the equivalent bulk modulus of the cleavage-mixed fluid into the Gassmann equation, the dry skeleton bulk modulus is obtained by back-calculation.

8. The method for predicting in-situ stress in coal reservoirs according to claim 1, characterized in that, The method for obtaining corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density by fluid substitution of near-wellbore P-wave velocity and near-wellbore S-wave velocity based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the cleavage porosity includes: Based on the resistivity values ​​of array inductive resistivity logging at each detection depth, and constrained by cleavage porosity, the cleavage water saturation at each radial position is calculated back by source distance according to Archie's formula; and the original gas-water two-phase formation volume density is calculated by source distance based on the cleavage water saturation at each radial position. Using the bulk modulus of the dry skeleton, the porosity of the cleavage, and the Gassmann equation, the saturated bulk modulus of the cleavage-restored original gas-water two-phase state is calculated for each source distance and denoted as the corrected saturated bulk modulus for each source distance. The corrected P-wave velocity is calculated for each source distance using the corrected saturated bulk modulus, dry skeleton shear modulus, and original gas-water two-phase formation bulk density. The calculation formula is: ； in, Correct the saturated bulk modulus for each source distance; The volumetric density of the formation in its original gas-water two-phase state; The corrected shear wave velocity was calculated using the dry skeleton shear modulus and the original gas-water two-phase formation bulk density. The calculation formula is: 。 9. The method for predicting in-situ stress in coal reservoirs according to claim 8, characterized in that, The method for calculating the bulk density of the original gas-water two-phase formation based on the source-distance source-by-source cleavage water saturation at each radial position includes: The original gas-water two-phase formation volume density at each source distance was calculated based on the formation volume density. The calculation methods include: ； in, For cleavage porosity, To constrain the initial water saturation of the cleavage, The original formation water density, To determine the water saturation of the cutting surface, The density of the drilling fluid filtrate. This represents the original formation gas density.

10. A coal reservoir in-situ stress prediction system, used for performing the method described in any one of claims 1 to 9, characterized in that, The system includes: The well logging data acquisition module is used to configure the array acoustic logging source spacing based on the depth distribution of drilling fluid intrusion into the formation and the direction of cleavage; after the source spacing is configured, array acoustic logging is performed to obtain the P-wave velocity and S-wave velocity at each source spacing; density logging is performed to obtain the formation volume density; and nuclear magnetic resonance logging is performed to obtain the cleavage porosity. The intrusion zone outer boundary identification module is used to calculate the P-wave and S-wave velocity ratio per source distance and select the intrusion zone outer boundary; the P-wave velocity and S-wave velocity corresponding to the source distance exceeding the intrusion zone outer boundary are recorded as far-well P-wave velocity and far-well S-wave velocity, and those falling within the intrusion zone outer boundary are recorded as near-well P-wave velocity and near-well S-wave velocity. The skeleton modulus inversion module is used to invert the skeleton bulk modulus and skeleton shear modulus based on the far-well P-wave velocity, far-well S-wave velocity, formation bulk density and cleavage porosity. The fluid replacement correction module is used to replace the near-wellbore P-wave velocity and near-wellbore S-wave velocity with fluid based on the bulk modulus of the dry skeleton, the shear modulus of the dry skeleton, and the porosity of the cleavage to obtain the corrected P-wave velocity, the corrected S-wave velocity, and the original gas-water two-phase formation bulk density. The geostress prediction module is used to calculate elastic parameters based on the corrected P-wave velocity, corrected S-wave velocity, and original gas-water two-phase formation bulk density, and then substitute them into the geostress calculation model to obtain the geostress prediction results of the target coal reservoir.

Images