A method for characterizing the migration zonation of overburden coalbed methane

Abstract

This invention relates to the field of mine safety and resource development technology, and particularly to a method for characterizing the migration and zonation of overburden coalbed methane, comprising an acquisition module, a simulation module, a processing module, an inversion module, a modeling module, and an output module. This invention combines physical similarity simulation with discrete element numerical simulation to achieve multi-source data synergy, significantly improving fracture identification accuracy and obtaining a permeability field with continuous spatial distribution characteristics. This data can be used for coalbed methane migration path analysis, enrichment zone identification, and extraction scheme optimization.

Description

Technical Field

[0001] This invention relates to the field of mine safety and resource development technology, and in particular to a method for characterizing the migration and zonation of overburden coalbed methane. Background Technology

[0002] After coal mining, the overlying strata in the goaf gradually form three typical zones during the long-term evolution of mining operations and post-mining periods: caving zones, fracture zones, and flexural subsidence zones. In caving zones, the rock blocks are broken and piled up, with extremely high fracture connectivity, allowing coalbed methane to flow primarily through rapid fracture conduction. In fracture zones, fractures are well-developed and exhibit inter-layer connectivity, with coalbed methane also migrating mainly along the fracture network. In flexural subsidence zones, although macroscopic fractures are not well-developed, the pressure relief caused by mining and long-term displacement adjustments create areas of increased permeability, where coalbed methane mainly flows through pores and diffuses. Due to the fundamental differences in the migration mechanisms of these three overlying zones, accurately describing the differentiated migration patterns of coalbed methane in each zone is crucial for identifying coalbed methane enrichment areas and optimizing extraction strategies.

[0003] However, existing coalbed methane migration models struggle to accurately reflect the differences in migration mechanisms across the three zones of the overburden. Specifically, many existing models simplify the entire overburden by treating it as a homogeneous or quasi-homogeneous medium. This approach ignores the significant differences in highly connected fracture networks in the caving zone, the interlayer connectivity of fracture zones, and the pore flow mechanisms in the tortuous subsidence zone. Consequently, these models fail to accurately depict the differentiated migration behavior of coalbed methane in each zone. Regarding fracture information extraction in the caving and fracture zones, existing methods often rely on simple binary statistics, which are insufficient for the demands of fracture image extraction, connectivity identification, and conductivity characterization under complex backgrounds. Their accuracy is particularly low under conditions of illumination disturbance and fine fracture noise interference. In constructing the permeability field in the tortuous subsidence zone, existing methods often rely on empirical relationships based on single stresses, neglecting the influence of displacement gradients and long-term evolution on pore flow capacity, making it difficult to obtain a continuous and reasonable permeability distribution.

[0004] The inadequacy of this zoning mechanism, the low accuracy of fracture identification, and the discontinuity of permeability field inversion severely restrict the accuracy of coalbed methane migration path analysis and enrichment zone identification. Summary of the Invention

[0005] The purpose of this invention is to provide a method for characterizing the zonation of gas migration in overburden coalbed methane, in order to solve the technical problems of insufficient expression of zonation mechanism, low accuracy of fracture identification, and discontinuity in the inversion of parameters of bending subsidence zone in existing overburden coalbed methane migration models.

[0006] The present invention provides a method for characterizing the migration and zonation of overburden coalbed methane. The method is executed by a system, which includes at least a data connection acquisition module, a simulation module, a processing module, an inversion module, a modeling module, and an output module. The acquisition of basic parameters of coal seam and overlying rock includes at least geometric parameters and physical and mechanical parameters, and at the same time determines the coal seam mining sequence; The simulation module constructs a physical similarity model based on the principle of similarity to simulate the overburden evolution process after coal seam mining and cessation of mining, and collects overburden fracture images at different mining and cessation stages; The processing module identifies and denoises the acquired rock fracture images, extracts fracture skeleton density, average aperture, and local connectivity index, and constructs a fracture image flow function. The inversion module is based on the discrete element numerical model to invert the continuous permeability field of the bent subsidence zone, and constructs a permeability function including stress relief term, displacement gradient term and displacement rate term to obtain a permeability field with continuous spatial distribution characteristics. The modeling module couples the fracture-dominated seepage zone of the collapse zone and fracture zone with the pore seepage zone of the tortuous subsidence zone through a smooth transition function to establish a zonal transport model. The output module outputs the coalbed methane migration path and enrichment zone based on the zonal migration model.

[0007] In some embodiments, the processing module includes multi-directional enhancement processing, which uses a directional filter bank to convolve the original crack image, and the expression for the directional response function is: ( , )= ( , ) ( , ) in, ( , () represents the grayscale value of the original image. ( , ) is the direction The enhancement filter kernel is used; the maximum envelope is taken for the response in all directions to obtain the enhanced image: ( , )=max( ).

[0008] In some embodiments, an adaptive binary segmentation is achieved by constructing a local contrast-grayscale joint threshold function based on the enhanced image. The expression for the threshold is: ( , )= + +

[0009] in, This represents the average gray level of a local window. For local standard deviation, For gradient terms, and This is an adjustment coefficient; when the pixel value is below the threshold, it is determined to be a crack region.

[0010] In some embodiments, a morphological-topological joint denoising method is used to screen candidate crack regions based on the binary segmentation results, and a crack candidate discrimination index is constructed, including thinness length: = 2 /

[0011] in, For the target skeleton length, The target projected area; simultaneously, the local connectivity length is calculated. and direction consistency Only retain those that meet the requirements. > 0、 > 0、 > The area of ​​0 is used to remove short-branch noise, isolated noise points, and non-cracked patches, among which S 0、 L 0、 Q 0 is the preset threshold.

[0012] In some embodiments, the preserved region is reconstructed using a skeletonization method, and parameters such as fracture skeleton density, average aperture, maximum connected cluster ratio, and orientation consistency are extracted. A fracture image flow guidance function is then constructed, the expression of which is: ( , )=1+ 1 ( , )+ 2 ( , )+ 3 ( , ) in, The density of the fracture skeleton, For average opening, The local connectivity index. 1. 2. 3 is the correction factor; subscript f For cracks, F f The physical meaning of is a dimensionless number used to quantitatively characterize the amplification effect of fracture networks on local permeability, and its value is greater than or equal to 1.

[0013] In some embodiments, the inversion module includes the establishment of a discrete element numerical model, which is consistent with the physical similarity model, and obtains continuous stress field and displacement field data of the bending subsidence zone at different stages through numerical simulation.

[0014] In some embodiments, the expression for the permeability function is:

[0015] in, The original stratum permeability, For average stress, For the original stress, For displacement field, For displacement gradient, For displacement rate, and These are displacement and velocity characteristic quantities, respectively. , , These are parameters to be determined.

[0016] In some embodiments, for the collapse zone and fracture zone, based on the original formation permeability k0, the local permeability is corrected using the fracture image conductivity function to obtain the permeability of the fracture-dominated seepage zone, expressed as: ( , )= 0 ( , ); For the tortuous subsidence zone, the continuous permeability field obtained by inversion is used. ( , , ) is used as the permeability of the pore seepage zone.

[0017] In some embodiments, a smooth transition function is constructed at the interface between the two zones, and the expression of the smooth transition function is:

[0018] in, The distance to the zonal boundary, in the fracture zone region. d Negative values ​​are found in the curved subsidence zone region. d Positive value; at the transition boundary d=d 0, d 0 represents the transition control distance. d 0 is determined based on the unit grid size and band thickness, and is usually taken as >3 grid unit sizes; This is a smoothing coefficient that controls the steepness of the transition curve; the larger the value, the more pronounced the transition. Based on a smooth transition function, the permeability fields of the caving zone, fracture zone, and tortuous subsidence zone are continuously spliced ​​together to obtain a unified permeability field, the expression of which is: ( , , )=Φ ( , )+(1 Φ) ( , , ); By incorporating a unified permeability field into the coalbed methane migration control equation, a zonal migration model that can reflect the combined effects of fracture conduction and pore diffusion is established.

[0019] Compared with the prior art, the present invention has the following beneficial effects: 1. This invention organically combines physical similarity simulation and discrete element numerical simulation, giving full play to the advantages of physical simulation in crack visualization and numerical simulation in continuous field acquisition, thereby realizing the utilization of multi-source data; 2. In the crack identification process, a multi-directional enhancement, adaptive binary segmentation and morphological topology joint denoising method is adopted, which significantly improves the automatic identification accuracy of crack images in collapse zones and crack zones, and overcomes the problem of insufficient identification accuracy of traditional methods in complex backgrounds; 3. In the permeability field construction stage, a stress-displacement joint inversion model is introduced to comprehensively consider the influence of multiple factors such as stress relief, displacement gradient, and displacement rate, so as to overcome the problem that a single stress empirical formula is difficult to accurately describe the permeability enhancement law of the bending subsidence zone and obtain a permeability field with continuous spatial distribution characteristics. 4. In the model coupling stage, by constructing the fracture image conduction function, continuous permeability field and zonal smooth transition function, the unified coupling between the fracture conduction zone and the pore diffusion zone is realized, and a coalbed methane zonal migration model that can truly reflect the differentiated migration mechanism of the three overburden zones is established. 5. This system can be used for coalbed methane migration path analysis, enrichment zone identification, and extraction scheme optimization, and has promising applications in the development of coalbed methane resources in abandoned mines. Attached Figure Description

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

[0021] Figure 1 This is a schematic diagram of the implementation architecture of the overburden coalbed methane migration zonation characterization method of the present invention. Detailed Implementation

[0022] The following will be based on embodiments of the present invention. Figure 1 The technical solutions in the embodiments of the present invention will be clearly and completely described. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0023] After coal mining, the overlying strata in the goaf gradually form three typical zones during the long-term evolution of mining operations and post-mining cessation: caving zones, fracture zones, and flexural subsidence zones. In caving zones, the rock blocks are broken and piled up with extremely high fracture connectivity, allowing coalbed methane (CBM) to flow primarily through rapid fracture conduction. In fracture zones, fractures are well-developed and exhibit inter-layer connectivity, with CBM primarily migrating along the fracture network. In flexural subsidence zones, although macroscopic fractures are not well-developed, they create permeability-enhancing areas under the effects of mining-induced pressure relief and long-term displacement adjustments, with CBM primarily flowing through pores via seepage and diffusion. Existing CBM migration models often simplify the entire overlying strata by treating them as homogeneous or quasi-homogeneous media, neglecting the differences in migration mechanisms among the three zones. This results in models failing to accurately depict the differentiated migration behavior of CBM in each zone. In extracting fracture information from collapse zones and fracture zones, existing methods mostly rely on simple binary statistics, which cannot meet the high-precision extraction requirements of fracture images under complex backgrounds, especially with low recognition accuracy under conditions of illumination disturbance and fine fracture noise interference. Regarding the construction of permeability fields in bent subsidence zones, existing methods often use empirical formulas based on a single stress for estimation, neglecting the influence of displacement gradients and long-term evolution on pore seepage capacity, making it difficult to obtain continuous and reasonable permeability distribution results.

[0024] Example This embodiment provides a system for the zonal migration of overburden coalbed methane, considering zonal migration patterns. It includes at least a data acquisition module, a simulation module, a processing module, an inversion module, a modeling module, and an output module. Specifically: In the acquisition module, firstly, geometric parameters such as coal seam thickness, spacing, dip angle, lithology, and overlying strata structure of the study area are obtained through geological structural statistics. Coal seam thickness directly affects the spatial extent of the goaf, and thus the height of the caving zone and the degree of fracture zone development; coal seam spacing determines the mutual influence between coal seams under multi-coal seam mining conditions; coal seam dip angle affects the asymmetric development characteristics of overlying strata zoning; lithological parameters determine the fracture pattern and fragmentation characteristics of the overlying strata under mining operations; and the overlying strata structure controls the vertical expansion and penetration characteristics of fractures.

[0025] Secondly, physical and mechanical parameters such as coal and rock layer density, elastic modulus, strength, and initial permeability were obtained through field data collection and laboratory measurements. These parameters serve as the foundational inputs for subsequent physical similarity simulations and numerical simulations. Coal and rock layer density affects the self-weight stress distribution of the overburden, while elastic modulus and strength determine the deformation and fracture behavior of the rock layers under mining operations. Initial permeability is an important indicator for assessing coalbed methane production capacity and migration ability. Simultaneously, based on the coal seam mining sequence and overburden zoning characteristics, the spatial boundaries and evolution patterns of each zoning zone were determined, providing geometric constraints for subsequent simulations.

[0026] In the simulation module, a physically similar model is constructed based on the principle of similarity, maintaining geometric, kinematic, and dynamic similarity to the prototype geological conditions. The selection of the similarity ratio needs to comprehensively consider factors such as model size, loading conditions, and measurement accuracy. Typically, a geometric similarity ratio of 1:50 to 1:100 is chosen, and a density similarity ratio of approximately 1:1.5. The time similarity ratio is not an independent parameter but is derived from the geometric similarity ratio based on the dynamic similarity criterion under a gravitational field. Its expression is:

[0027] In the formula, For time similarity ratio, This represents the geometric similarity ratio. For example, a geometric similarity ratio of 1:100 means that one day of evolution in the model corresponds to 10 days of long-term evolution in the prototype. This relationship ensures the similarity between the gravity-driven overburden movement, caving, and compaction processes in the model and the prototype.

[0028] During the model laying process, similar materials are prepared according to the design ratio, laid layer by layer and compacted to ensure that the mechanical properties of each rock layer maintain a similar relationship with the prototype. Displacement sensors, stress sensors and acoustic emission probes can be pre-embedded in the model to monitor the deformation, stress distribution and fracture characteristics of the overlying rock during mining.

[0029] The coal seam mining simulation employs a step-by-step excavation method, allowing sufficient time for the overlying strata to fully respond after each excavation stage to simulate the long-term evolution process after actual mining. High-definition cameras or industrial CCDs are used to acquire images of overlying strata fractures during different mining and shutdown phases. The image acquisition system is equipped with multi-angle light sources to reduce the impact of shadows, and a calibration board is used for spatial calibration to obtain the true spatial location information of the fractures.

[0030] The processing module performs four steps on the acquired crack images: multi-directional enhancement, adaptive binary segmentation, morphological and topological joint denoising, and skeletonization reconstruction and flow function construction. The multi-directional enhancement process employs a directional filter bank to convolve the original crack image to highlight crack features in different orientations. The directional filter bank contains multiple filter kernels in different directions, typically 8 to 16 directions, each designed using a Gabor function or a directional Gaussian derivative function. The expression for the directional response function is: ( , )= ( , ) ( , ) in, ( , () represents the grayscale value of the original image. ( , ) is the direction An enhanced filter kernel is used; the maximum envelope is taken for the response in all directions to obtain the enhanced image.

[0031] For each pixel (x, y), calculate its response value in all directions, and then take the maximum envelope to obtain the enhanced image. ( , )=max( This processing method effectively enhances the characteristics of cracks in any orientation, avoiding the directional sensitivity caused by single-directional filters. The design parameters of the directional filter bank include the number of filter directions, filter scale, and Gaussian envelope parameters. The selection of these parameters needs to be optimized and adjusted according to the resolution of the crack image and the crack width distribution characteristics.

[0032] An adaptive binary segmentation is achieved by constructing a local contrast-grayscale joint threshold function based on the enhanced image. The expression for this function is: ( , )= + * + *

[0033] in, This represents the average gray level of a local window. For local standard deviation, For gradient terms, and This is the adjustment factor. The size of the local window is typically chosen to be 3 to 5 times the crack width to ensure the representativeness of the statistical characteristics. Adjustment factor and The value of is determined based on the overall image contrast and noise level. Through preliminary experiments, the Otsu method is used to calibrate typical crack images, yielding the desired result. The value range is 2.0 to 5.0. The value range is 0.5 to 2.0. Pixel values ​​less than or equal to the threshold are identified as crack regions, while pixel values ​​higher than the threshold are identified as background regions. Adaptive binary segmentation automatically determines the segmentation threshold based on the statistical characteristics of the local regions of each pixel. Specifically, this threshold function integrates the local gray-level mean (reflecting base brightness), local standard deviation (reflecting texture complexity), and gradient term (reflecting edge strength). In regions with low contrast between cracks and the background, the values ​​of the local standard deviation and gradient term are small, making the threshold mainly determined by the local gray-level mean and at a low level, thus effectively reducing the risk of missed crack detection and improving the detection rate. In contrast, in regions with high contrast, the contributions of the local standard deviation and gradient term increase, and the threshold increases accordingly, helping to suppress background noise and reduce false positives.

[0034] The morphological and topological joint denoising step, based on the binary segmentation results, employs a morphological and topological joint denoising method to screen candidate crack regions and eliminate various types of noise. The constructed crack candidate discrimination index includes thinness length: = 2 /

[0035] in, A is the target skeleton length, and A is the target projected area; local connectivity length. Lc Defined as the maximum connectivity length of the fracture target in an eight-connected domain; directional consistency Q is defined as the cosine similarity between the principal direction of the fracture target and the local average direction, retaining only those that simultaneously satisfy... > 0、 > 0、 > In areas with zero conditions, short-branch noise, isolated noise points, and non-fractured patches are eliminated. S 0、 L 0 、Q 0 is the preset threshold. The expression is:

[0036] in, For the first connected path One skeleton pixel, This represents the total number of pixels along the path.

[0037] The expression for directional consistency Q is:

[0038] in, The main direction angle is obtained by performing principal component analysis (PCA) on the coordinates of all pixels within the target area. This is the local average direction angle.

[0039] The morphological-topological joint denoising discrimination index establishes a criterion for distinguishing real fractures from various types of noise by comprehensively considering three topological features: the elongation, connectivity length, and directional consistency of the fracture target.

[0040] Experiments show that when the thin length threshold S 0 is set to 5, local connectivity length threshold. L 0 is set to 10 pixels, and the orientation consistency threshold is set to 0. Q The noise reduction effect is optimal when the value is set to 0.7.

[0041] In the skeletonization and flow function construction stages, morphological refinement of the preserved region is performed, and the fracture skeleton is extracted and parameters such as skeleton density, average aperture, maximum connected cluster ratio, and directional consistency are calculated. (Crack skeleton density) Defined as the total length of the skeleton per unit area, average opening Defined as the ratio of fracture area to skeleton length, the local connectivity index. Defined as the ratio of the area of ​​the largest connected cluster to the total area of ​​the fracture.

[0042] Based on these parameters, a flow guidance function for the fracture image is constructed, firstly by adjusting the fracture skeleton density. and average opening Perform dimensionless processing: Define dimensionless skeleton density * = / re Dimensionless average aperture * = / re In this context, the superscript * indicates dimensionless, and the subscript ref indicates a reference value. re The reference density (which can be the reciprocal of the model tendency length) is used. re Reference aperture (can be taken as the maximum crack aperture or 1 mm). Local connectivity index. It is a dimensionless quantity. The expression for the flow conduction function in the fracture image, constructed based on the aforementioned dimensionless parameters, is as follows: ( , )=1+ 1 * ( , )+ 2 * ( , )+ 3 ( , ) in, 1. 2. 3 represents a dimensionless correction coefficient used to balance the contribution weights of different parameters to conductivity. The value of the correction coefficient is determined through multivariate nonlinear regression analysis using indoor core seepage experimental data and fracture image parameters at corresponding locations. Preferably, in the absence of experimental conditions, empirical values ​​can be obtained based on sensitivity analysis, with a range of [value missing]. , , The subscript f represents the crack, F f It is a dimensionless number used to quantitatively characterize the amplification effect of fracture networks on local permeability, and its value is greater than or equal to 1.

[0043] The fracture image conductivity function establishes a quantitative relationship between fracture microscopic characteristics and macroscopic seepage capacity by linearly weighting three key parameters: fracture skeleton density, average aperture, and local connectivity index. The correction coefficient needs to be calibrated based on seepage experimental data from actual coal and rock samples.

[0044] In the inversion module, the continuous permeability field of the bent subsidence zone is inverted based on discrete element numerical simulation. This module includes three steps: discrete element numerical model establishment, permeability function construction, and parameter inversion.

[0045] The discrete element numerical model (DEM) is consistent with the physical similarity model, and simulations are performed at a 1:1 scale based on the actual site dimensions and physical properties. The model boundary conditions are set as follows: the bottom boundary is a stress boundary condition, applying vertical stress equivalent to the original strata; the two side boundaries are displacement boundary conditions, restricting horizontal displacement; and the top boundary is a free boundary. A parallel bond model is used to describe the interlayer bond strength, and a smooth joint model is used to describe the pre-defined weak surface directions. The coal seam mining simulation adopts the same step-by-step excavation scheme as the physical model to ensure the comparability of the numerical simulation results with the physical simulation results.

[0046] Numerical simulation was used to obtain continuous stress and displacement field data of the bent subsidence zone at different stages, including stress components, displacement components, and displacement rate components.

[0047] The permeability function includes a stress relief term, a displacement gradient term, and a displacement rate term, and its expression is:

[0048] in, The original stratum permeability, For average stress, For the original stress, For displacement field, For displacement gradient, For displacement rate, and These are displacement and velocity characteristic quantities, respectively. , , These are parameters to be determined. The permeability function couples multiple factors, including stress state, displacement gradient, and displacement rate, in an exponential form. The stress state reflects the enhancing effect of stress release caused by mining on permeability; permeability increases when the average stress is less than the original stress. The displacement gradient reflects the deformation gradient effect; a larger displacement gradient indicates more severe rock deformation, more developed microfractures, and a more significant increase in permeability. The displacement rate reflects the characteristics of time-dependent evolution; under long-term shutdown conditions, the displacement rate approaches zero, but the cumulative displacement is large. At this time, permeability evolution is mainly determined by cumulative deformation.

[0049] In the modeling module, the fracture-dominated seepage zone of the collapse zone and the fracture zone is coupled with the pore seepage zone of the tortuous subsidence zone. This module includes four steps: determining the permeability of the fracture-dominated seepage zone, determining the permeability of the pore seepage zone, zonal transition treatment, and constructing a unified permeability field.

[0050] For caving zones and fracture zones, based on the original strata permeability k Based on 0, the permeability of the fracture-dominant seepage zone is obtained by correcting the local permeability using the fracture image conductivity function. Its expression is: ( , )= 0 ( , ).

[0051] Due to the flow function of the crack image ( , The value of ) is usually greater than 1, therefore the permeability of the fracture-dominated seepage zone is higher than that of the original formation. For the tortuous subsidence zone, the continuous permeability field obtained by inversion is used. ( , , ) is used as the permeability of the pore seepage zone.

[0052] A smooth transition function is constructed at the zone interface to achieve a continuous transition of permeability between different zones. The expression for the smooth transition function is:

[0053] in, The distance to the zone boundary, 0 represents the transition control distance. This is the smoothing coefficient.

[0054] The smooth transition function adopts an S-shaped curve form, with its value gradually transitioning from 0 to 1 near the zonal boundary. This achieves a continuous transition in permeability between the fracture-dominated seepage zone and the pore seepage zone, effectively avoiding abrupt changes in permeability at the zonal interface and controlling the transition distance. d 0 is typically taken as 10% to 20% of the strip thickness, and the smoothness coefficient is... δ Control the steepness of the transition zone.

[0055] By continuously splicing the permeability fields of the caving zone, fracture zone, and tortuous subsidence zone according to the smooth transition function, a unified permeability field is obtained, the expression of which is: ( , , )=Φ ( , )+(1 Φ) ( , , ).

[0056] The unified permeability field approximates the permeability of the fracture-dominated seepage zone in the fracture-dominated seepage zone, and the permeability of the pore seepage zone in the pore seepage zone. In the zonal transition zone, it is weighted and fused according to the smooth transition function. The unified permeability field can also be defined in segments within different zones. That is, the permeability of the fracture-dominated seepage zone is used in the collapse zone and fracture zone, the permeability of the pore seepage zone is used in the tortuous subsidence zone, and the smooth transition function is used for correction near the zonal boundary.

[0057] A unified permeability field is incorporated into the governing equations for coalbed methane migration to establish a zonal migration model that reflects the combined effects of fracture conduction and pore diffusion. The governing equations comprehensively consider the convection-diffusion mechanism in the fracture network and the adsorption-desorption mechanism in the porous medium, achieving differentiated descriptions of migration mechanisms in different zones by introducing a zonal permeability field. In the caving and fracture zones, coalbed methane mainly moves along the fracture network via convection and diffusion, exhibiting high permeability and Darcy velocity. In the tortuous subsidence zone, coalbed methane mainly seeps and diffuses through pores, while adsorption-desorption of methane by the coal matrix also occurs. The zonal migration model is solved using a multiphysics coupling method, simultaneously considering the interactions of the gas flow field, solid deformation field, and concentration field.

[0058] In the output module, the migration paths and enrichment zones of coalbed methane (CBM) are output based on the zonal migration model. The dominant migration paths of CBM along highly connected fracture networks within the caving and fracture zones are analyzed. These dominant migration paths are determined by tracing high-permeability connected regions to identify the main flow channels of CBM. The slow diffusion and enrichment range of CBM in the permeability-enhancing zone of the tortuous subsidence zone are analyzed. Favorable enrichment areas are determined by calculating the CBM concentration distribution and accumulation. The asymmetric enrichment characteristics of CBM preferentially expanding towards higher elevations under inclined coal seam conditions are identified. These asymmetric enrichment characteristics consider the influence of gravity and zonal boundary morphology on CBM migration under inclined coal seam conditions, providing a basis for the layout of extraction areas.

[0059] To better understand this invention, we will take an abandoned inclined coal seam mine as an example to illustrate the specific application process of this invention.

[0060] The coal seam in this mine has a dip angle of 25 degrees, a thickness of 3.5 meters, and a mining depth of approximately 400 meters. The overlying strata are mainly interbedded sandstone and mudstone. First, the geometric parameters of the coal seam and the mechanical parameters of the overlying strata were determined through geological data collection. A physical similarity model was constructed and mining simulations were conducted. Images of overlying strata fractures were acquired during key mining stages.

[0061] Secondly, the fracture image is processed using the multi-directional enhancement, adaptive binary segmentation and morphological topology joint denoising method proposed in this invention. The fracture skeleton density, average aperture and local connectivity index of the collapse zone and fracture zone are extracted, and the fracture image flow function is constructed to obtain the spatial distribution of permeability in the fracture-dominated seepage zone.

[0062] Next, a discrete element numerical model consistent with the physical model is established, and numerical simulations are run to obtain stress and displacement field data of the bent subsidence zone. A permeability function is constructed and parameter inversion is performed to obtain the continuous distribution of permeability in the pore seepage zone. Then, a smooth transition function is used to process the zonal boundaries to obtain a unified permeability field, which is then imported into the coalbed methane migration control equation to establish a zonal migration model.

[0063] Finally, numerical simulations were performed using a zonal migration model to output the dominant migration paths and the distribution of enrichment zones. The results show that coalbed methane forms dominant migration paths along the high-level fracture network within the fracture zone, and enrichment zones form near the fracture zone boundary in the tortuous subsidence zone, exhibiting a clear asymmetric enrichment characteristic. Based on these analysis results, it is recommended to locate extraction wells on the high-level side of inclined coal seams to effectively develop coalbed methane resources in abandoned mines.

[0064] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. It will be apparent to those skilled in the art that the invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered illustrative and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the scope of the invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0065] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A method for characterizing the migration and zonation of overburden coalbed methane, wherein the method is executed by a system, characterized in that, The system includes at least a data connection acquisition module, a simulation module, a processing module, an inversion module, a modeling module, and an output module; The acquisition of basic parameters of coal seam and overlying rock includes at least geometric parameters and physical and mechanical parameters, and at the same time determines the coal seam mining sequence; The simulation module constructs a physical similarity model based on the principle of similarity to simulate the overburden evolution process after coal seam mining and cessation of mining, and collects overburden fracture images at different mining and cessation stages; The processing module identifies and denoises the acquired rock fracture images, extracts fracture skeleton density, average aperture, and local connectivity index, and constructs a fracture image flow function. The inversion module is based on the discrete element numerical model to invert the continuous permeability field of the bent subsidence zone, and constructs a permeability function including stress relief term, displacement gradient term and displacement rate term to obtain a permeability field with continuous spatial distribution characteristics. The modeling module couples the fracture-dominated seepage zone of the collapse zone and fracture zone with the pore seepage zone of the tortuous subsidence zone through a smooth transition function to establish a zonal transport model. The output module outputs the coalbed methane migration path and enrichment zone based on the zonal migration model.

2. The method according to claim 1, characterized in that, The processing module includes multi-directional enhancement processing, which uses a directional filter bank to perform convolution processing on the original crack image. The expression for the directional response function is: ( , )= ( , ) ( , ) in, ( , () represents the grayscale value of the original image. ( , ) is the direction The enhancement filter kernel is used; the maximum envelope is taken for the response in all directions to obtain the enhanced image: ( , )=max( ).

3. The method according to claim 2, characterized in that, Adaptive binary segmentation is achieved by constructing a local contrast-grayscale joint threshold function based on the enhanced image. The expression for the threshold is: ( , )= + + in, The average gray level of a local window. For local standard deviation, For gradient terms, and This is an adjustment coefficient; when the pixel value is below the threshold, it is determined to be a crack region.

4. The method according to claim 3, characterized in that, Based on the binary segmentation results, a morphological and topological joint denoising method is used to screen candidate crack regions, and a crack candidate discrimination index is constructed, including thinness length: = 2 / in, For the target skeleton length, The target projected area; simultaneously, the local connectivity length is calculated. and direction consistency Only retain those that meet the requirements. > 0、 > 0、 > The area of ​​0 is used to remove short-branch noise, isolated noise points, and non-cracked patches, among which S 0、 L 0、 Q 0 is the preset threshold.

5. The method according to claim 4, characterized in that, The preserved region is reconstructed using skeletonization. Fracture skeleton density, average aperture, maximum connected cluster ratio, and orientation consistency parameters are extracted. A fracture image flow function is then constructed, the expression of which is: ( , )=1+ 1 ( , )+ 2 ( , )+ 3 ( , ) in, The density of the fracture skeleton, For average opening, The local connectivity index.

1.

2. 3 is the correction factor; subscript f For cracks, F f The physical meaning of is a dimensionless number, used to quantitatively characterize the amplification effect of fracture networks on local permeability, and its value is greater than or equal to 1.

6. The method according to claim 1, characterized in that, The inversion module includes the establishment of a discrete element numerical model, which is consistent with the physical similarity model. The continuous stress field and displacement field data of the bending subsidence zone at different stages are obtained through numerical simulation.

7. The method according to claim 6, characterized in that, The expression for the permeability function is: in, The original stratum permeability, For average stress, For the original stress, For displacement field, For displacement gradient, For displacement rate, and These are displacement and velocity characteristic quantities, respectively. , , These are parameters to be determined.

8. The method according to claim 1, characterized in that, For the collapse zone and fracture zone, based on the original formation permeability k0, the local permeability is corrected using the fracture image conductivity function to obtain the permeability of the fracture-dominated seepage zone, expressed as: ( , )= 0 ( , )； For the tortuous subsidence zone, the continuous permeability field obtained by inversion is used. ( , , ) is used as the permeability of the pore seepage zone.

9. The method according to claim 1, characterized in that, A smooth transition function is constructed at the interface between the two zones. The expression of the smooth transition function is as follows: in, The distance to the zonal boundary, in the fracture zone region. d Negative values ​​are found in the curved subsidence zone region. d Positive value; at the transition boundary d=d 0, d 0 represents the transition control distance. d 0 is determined based on the unit grid size and band thickness, and is usually taken as >3 grid unit sizes; This is a smoothing coefficient that controls the steepness of the transition curve; the larger the value, the more pronounced the transition. Based on a smooth transition function, the permeability fields of the caving zone, fracture zone, and tortuous subsidence zone are continuously spliced ​​together to obtain a unified permeability field, the expression of which is: ( , , )=Φ ( , )+(1 F) ( , , ); By incorporating a unified permeability field into the coalbed methane migration control equation, a zonal migration model that can reflect the combined effects of fracture conduction and pore diffusion is established.

Images