Quantitative characterization method for depicting mining-induced fractures of goaf roof

Abstract

The invention discloses a quantitative characterization method for mining-induced fractures of a goaf roof, which belongs to the technical field of mine hydrogeology and water control engineering, and comprises the following steps of: obtaining the development height of a fracture guiding zone through network parallel electrical method CT (Computed Tomography) field monitoring, carrying out similar material physical simulation and carrying out graying and binarization processing on an image; the method comprises the following steps: extracting a void ratio discrete data sequence, establishing a coupling equation comprising an exponential attenuation base item, a multi-peak Gaussian item and a constant background value, introducing a segmented dynamic attenuation coefficient to describe stress dissipation differences of different rock stratums, continuously correcting initial strength at a lithology mutation interface, and performing parameter inversion by adopting a damping least square method. And outputting a characteristic equation for ponding evolution prediction and water guide channel identification by taking the root-mean-square error and the correlation coefficient as evaluation indexes. The method is high in fitting precision and clear in physical significance, and reliable parameter support can be provided for water control engineering design.

Description

Technical Field

[0001] This invention belongs to the field of mine hydrogeology and water control engineering technology, and specifically relates to a quantitative characterization method for delineating mining-induced fractures in the roof of goaf areas. Background Technology

[0002] During coal seam mining, the overlying strata undergo structural damage due to stress adjustment and gravity, forming a complex fracture system from bottom to top, including caving zones, fracture zones, and tortuosity zones. Among these, the roof fracture zone is the main water-conducting channel between the goaf and the overlying aquifer, and its development pattern directly determines the goaf's water storage capacity, hydrodynamic response, and risk of water inrush. Accurately describing the spatial distribution characteristics of fracture zones is of great significance for goaf water hazard prevention and control.

[0003] Current research on the structural characteristics of fracture zones mostly employs empirical formulas or unimodal function fitting methods, commonly including exponential decay functions, power functions, or Gaussian function models. While these models are applicable to describing the overall trend of fracture development, they are significantly insufficient in characterizing the distribution patterns of fractures under complex geological conditions.

[0004] The main problems with traditional models are as follows. First, the roof of a coal seam is usually composed of strata of different lithologies, and the stress distribution at the interface between soft and hard rocks can change abruptly, leading to multi-peak characteristics in fracture development. A single-peak function cannot accurately describe this local anomaly. Second, most existing models assume that the attenuation coefficient is constant, failing to reflect the dynamic attenuation law of fracture development in different strata with increasing height, and ignoring the differential influence of different lithologies on stress transmission and fracture response. Third, the parameters of many models lack actual physical meaning, being obtained only through mathematical fitting, and lack coupling verification with field monitoring data and physical simulation experiments, thus limiting the reliability and generalizability of the models. In addition, traditional methods often separate the overall attenuation trend and local peak characteristics when describing fracture development, making it difficult to form a unified mathematical expression.

[0005] These shortcomings result in significant errors in existing models when providing parameter inputs for predicting the evolution of water accumulation in goaf areas and for hydrodynamic analysis, affecting the scientific rigor and relevance of water control engineering design. Therefore, there is an urgent need to establish a quantitative descriptive equation that can simultaneously reflect the overall attenuation trend of fracture development and the response to local lithological abrupt changes. Furthermore, the accuracy and practicality of the model should be improved through multi-source data coupling verification, thereby achieving a refined characterization of fracture development and water accumulation processes in goaf areas. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a quantitative characterization method for the delineation of mining-induced fractures in the roof of goaf areas. This method determines the development height of the fracture zone through on-site monitoring using network-parallel electrical resistivity tomography (OTCT), extracts discrete data sequences of fracture porosity using similar material physics simulation, establishes an exponentially decaying corrected multi-peak Gaussian mixture equation, introduces a piecewise dynamic decay coefficient to describe the differences in stress dissipation rates among different rock strata, and continuously corrects the initial strength at lithological abrupt change interfaces. Parameter inversion optimization is performed using the damped least squares method, ultimately outputting a characteristic equation with clear physical meaning for predicting goaf water accumulation evolution and identifying water-conducting channels, providing reliable parameter support for water control engineering design.

[0007] To achieve the above objectives, the present invention provides a quantitative characterization method for delineating mining-induced fractures in the roof of a goaf, comprising the following steps: Step S1: Obtain key hydrogeological parameters of the roof of the goaf, use network parallel electrical resistivity CT technology to dynamically monitor the resistivity of the coal mine working face, achieve high-resolution imaging of the underground structure by arranging multiple sets of electrode arrays, interpret the resistivity data to extract the development law of the fracture network, the range of the goaf, and the development of the roof fracture conduction height, and obtain the development height of the fracture zone.

[0008] Step S2: Conduct physical simulation of similar materials. Set up a similar model according to the development height of the fracture zone. Perform image grayscale and binarization processing on the simulation results. Extract porosity through gridded statistics to form a discrete data sequence of porosity with height. Identify the distribution height of fracture development peak and the critical interface height of lithological abrupt change.

[0009] Step S3: Construct a coupled equation. Based on the discrete data sequence of porosity with height, establish a fracture development characteristic equation including an exponentially decaying base term, multiple Gaussian superposition terms, and a constant background value. The exponentially decaying base term describes the overall decay trend of fracture development with height, the multiple Gaussian superposition terms describe the peak fracture density at different lithological interfaces, a piecewise dynamic decay coefficient is introduced to describe the change in stress dissipation rate in different rock layers, and the initial strength of the exponentially decaying base term is corrected at the critical interface height of lithological abrupt change to make the function continuous at the piecewise points.

[0010] Step S4, Parameter Inversion and Accuracy Evaluation: The coupled equations are fitted and analyzed based on a nonlinear least squares optimization algorithm. The fitted curves and parameter estimates are output. The residual distribution, root mean square error, and correlation coefficient are used as evaluation indicators to evaluate the model fitting effect. If the evaluation fails to meet the standard, the coupled equations are corrected and the parameter inversion is re-executed. When the evaluation meets the standard, the coupled equations are output for prediction of goaf water accumulation evolution and identification of water diversion channels.

[0011] Furthermore, step S1 also includes obtaining the groundwater recharge and drainage conditions, formation thickness, aquitard relationship, borehole data and monitoring well water level in the mining area, and verifying the resistivity data interpretation results using hydrogeological borehole field monitoring data.

[0012] Further, in step S2, the physical simulation of similar materials satisfies the conditions of geometric similarity, mechanical similarity, and dynamic similarity, wherein: the formula for calculating the geometric similarity coefficient is: In the formula: The geometric similarity coefficient; The geometric dimensions are those of the actual geological prototype; These are the geometric dimensions of similar models.

[0013] The formula for calculating the mechanical similarity coefficient is: In the formula: The mechanical similarity coefficient; The stress is based on the actual geological prototype; Stress for similar models; is the material density similarity coefficient.

[0014] The formula for calculating the time similarity coefficient is: .

[0015] Furthermore, in step S2, the image grayscale conversion uses a weighted average method, and the conversion formula is: In the formula: Grayscale value; , , These are the red, green, and blue channel values ​​of the pixel, respectively.

[0016] The formula for binarization is: In the formula: These are the pixel values ​​of a binary image, where 1 represents white and 0 represents black. These are the pixel values ​​of a grayscale image; The threshold value is set.

[0017] Furthermore, in step S2, the image processing program is used to extract partition information from the binarized image, calculate the porosity within each grid, identify the peak range caused by stress concentration at the interface between soft and hard rocks based on the vertical differentiation characteristics of the porosity, reveal the vertical differentiation characteristics and peak causes of the porosity, and thus determine the critical interface height of the lithological abrupt change.

[0018] Furthermore, in step S3, the basic expression of the coupling equation is: In the formula: For height Porosity at that location; The initial strength for fracture development; The attenuation coefficient reflects the stress dissipation rate; For the first The area of ​​a superimposed Gaussian peak; For the first The distribution height of each fracture peak; For the first The width of each peak crack distribution; Baseline porosity; denoted as the number of Gaussian superposition terms.

[0019] Further, in step S3, the expression for the segmented dynamic attenuation coefficient is: In the formula: The segmented attenuation coefficient; The height of the critical interface for the lithological abrupt change is given.

[0020] Further, in step S3, the exponentially decaying base term is corrected, and the end value of the previous segment is used as the initial intensity coefficient of the next segment at the segmentation point. The corrected final expression is: In the formula: The first paragraph, the second paragraph, and so on are respectively the first paragraph, the second paragraph, and so on. The function value of the exponentially decaying base term at the segmentation point.

[0021] Furthermore, in step S4, the nonlinear least squares optimization algorithm uses damped least squares to optimize parameters, and performs fitting analysis with the goal of minimizing the error between the predicted value and the discrete data sequence of the porosity to obtain the optimal feature parameters. The damped least squares method is preferably the Levenberg-Marquardt algorithm.

[0022] Further, in step S4, when the root mean square error is significantly lower than the data fluctuation range and the correlation coefficient is greater than 0.95, the model fitting effect is deemed to have met the evaluation criteria; if the evaluation does not meet the criteria, the coupling equation is corrected for background values ​​or continuity, and the parameter inversion is re-executed until the evaluation meets the criteria, at which point the characteristic equation for predicting the development of fissures and water accumulation in the goaf is output.

[0023] Compared with the prior art, the present invention has at least the following beneficial effects.

[0024] (1) This invention establishes an exponentially decaying corrected multi-peak Gaussian mixture model by using network parallel electrical resistivity CT technology for dynamic resistivity monitoring and combining it with similar material physics simulation. This model describes the overall decay trend of fracture development with height through exponentially decaying base terms, accurately captures the peak fracture density at the interface between soft and hard rocks through multiple Gaussian superposition terms, reflects the change in stress dissipation rate in different rock layers by introducing piecewise dynamic decay coefficients, and corrects the initial strength at the critical interface of lithological abrupt change to maintain the continuity of the function, thereby achieving a unified mathematical description of the heterogeneity and multi-peak characteristics of fractures. This method overcomes the defect that traditional single-peak functions cannot describe the stress abrupt change phenomenon at the lithological interface, enabling the model to meet the hydrogeological constraints of decaying with height approaching the background value on a macroscopic level, and accurately characterize the local peak response on a microscopic level. Compared with traditional models such as polynomial functions, simple exponential functions, power functions, or single-peak Gaussian functions, this invention significantly reduces systematic errors and fitting residuals, and improves the explanatory power and prediction accuracy of heterogeneous fracture fields.

[0025] (2) This invention obtains key hydrogeological parameters such as groundwater recharge and drainage conditions, stratum thickness, aquitard relationships, borehole data, and monitoring well water levels in the mining area, and verifies the resistivity data interpretation results using hydrogeological borehole field monitoring data, establishing a dual constraint mechanism of field monitoring and physical simulation. This multi-source data cross-validation method effectively improves the reliability of fracture development law identification, provides background values ​​with practical geological significance for subsequent parameter inversion, avoids the arbitrariness of model parameters, ensures the consistency between the constructed equations and actual geological conditions, and makes the model have a stronger physical basis and engineering applicability.

[0026] (3) This invention establishes a physical simulation of similar materials that meets the conditions of geometric similarity, mechanical similarity, and dynamic similarity. It uses a weighted average method for image grayscale conversion and threshold-based binarization. By extracting porosity through gridded statistics, a discrete data sequence with height is formed. The image processing program is used to extract partition information from the binarized image and calculate the porosity within each grid, revealing the vertical differentiation characteristics and peak formation of porosity. This quantitative image analysis method transforms the physical simulation test results into discrete data that can be used for mathematical modeling, providing a clear physical basis for setting multi-peak terms. It makes the distribution height of fracture development peaks and the determination of critical interface heights of lithological abrupt changes traceable and repeatable, greatly improving the scientific nature and standardization of model construction, and facilitating its application under similar geological conditions.

[0027] (4) This invention optimizes the nonlinear parameters of the coupled equations using a damped least squares method combining gradient descent and Gauss-Newton methods. A rigorous accuracy evaluation system is established using residual distribution, root mean square error, and correlation coefficient as evaluation indicators. The model is deemed to have met the evaluation criteria when the root mean square error is significantly lower than the data fluctuation range and the correlation coefficient is greater than 0.95. This method ensures the global optimality and numerical stability of the fitted parameters. The quantitative evaluation criteria avoid overfitting or underfitting of the model, guaranteeing the predictive reliability of the output equations. This enables the constructed characteristic equations of fracture development and water accumulation evolution to provide high-precision, high-confidence parameter support for predicting water accumulation evolution in goaf areas, identifying water channels, and designing water control engineering projects. This significantly improves the scientific decision-making level and engineering safety assurance capabilities for mine water hazard prevention. Attached Figure Description

[0028] The accompanying drawings, which form part of this specification, illustrate embodiments of the invention and, together with the specification, serve to explain the principles of the invention.

[0029] The invention will become clearer with reference to the accompanying drawings and the following detailed description.

[0030] Figure 1 A flowchart illustrating a quantitative characterization method for delineating mining-induced fractures in the roof of a goaf, provided in an embodiment of the present invention.

[0031] Figure 2 This is a schematic diagram of the monitoring results of the development characteristics of the fracture zone in the roof of a coal seam.

[0032] Figure 3 This is a cross-sectional view of a similar model.

[0033] Figure 4 This is a schematic diagram of the binarized image processing of similar simulation results.

[0034] Figure 5 This is a schematic diagram of the binarized image partitioning recognition results.

[0035] Figure 6 This is a schematic diagram showing the measured values ​​and fitted curves of the porosity of the fracture zone in the roof as a function of height. Detailed Implementation

[0036] 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. The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the present invention or its application or use. 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. Unless otherwise specifically stated, the relative arrangement, expressions, and values ​​of components and steps set forth in these embodiments do not limit the scope of the present invention. Techniques, methods, and devices known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and devices should be considered part of the specification. In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0037] In this article, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0038] like Figure 1 As shown, this invention provides a quantitative characterization method for delineating mining-induced fractures in the roof of a goaf, comprising the following steps: Step S1, obtaining key hydrogeological parameters of the goaf roof. Networked parallel electrical resistivity CT technology is used to dynamically monitor the resistivity of the coal mine working face. Multiple electrode arrays are deployed to achieve high-resolution imaging of the underground structure. The resistivity data is interpreted to extract the development pattern of the fracture network, the extent of the goaf, and the development of the roof fracture conduction height, thus obtaining the development height of the fracture zone. Simultaneously, the groundwater recharge and drainage conditions, stratum thickness, aquitard relationships, borehole data, and monitoring well water levels in the mining area are obtained. The resistivity data interpretation results are verified using on-site monitoring data from hydrogeological boreholes.

[0039] Step S2: Conduct physical simulation of similar materials. Based on the fracture development height obtained in Step S1, a similar model is set up, specifically including the following sub-steps: Step S21: Construct a physical model that satisfies the conditions of geometric similarity, mechanical similarity, and dynamic similarity. Geometric similarity means that the shape of the model is consistent with the prototype, and the dimensions of each part are scaled according to a uniform ratio.

[0040] Geometric similarity coefficient The calculation method is shown in equation (1): In the formula: The geometric similarity coefficient; The geometric dimensions are those of the actual geological prototype; These are the geometric dimensions of the similar model. Mechanical similarity requires that stress, strain, and material properties be adjusted proportionally.

[0041] Mechanical similarity coefficient Related to the geometric similarity coefficient and the material density similarity coefficient, the calculation formula is Equation (2): In the formula: The mechanical similarity coefficient; The stress is based on the actual geological prototype; Stress for similar models; This is the material density similarity coefficient. Dynamic similarity requires ensuring the proportional relationships of relevant physical quantities such as time, fluid parameters, and permeability, so that the dynamic processes such as displacement and deformation in the similar material model are consistent with the principles of the geological prototype.

[0042] Time similarity coefficient Related to the geometric similarity coefficient, its calculation formula is as shown in equation (3): .

[0043] Step S22: Perform image grayscale and binarization processing on the simulation results. The commonly used grayscale method is the weighted average method. Given that the human eye is most sensitive to green and least sensitive to blue, in order to preserve the brightness contrast information of rock texture and fissures to the greatest extent in grayscale conversion, this invention selects weighting coefficients that conform to the visual characteristics of the human eye. The conversion formula is as shown in equation (4): In the formula: Grayscale value; , , These represent the red, green, and blue channel values ​​of the pixel, respectively. Using these specific coefficients effectively suppresses background noise and improves the accuracy of crack edge recognition in subsequent binarization processing. It is worth noting that while these coefficients are not unique in the field of image grayscale conversion, they are the optimal coefficients for the application scenario of this invention. Selecting these coefficients enhances the repeatability and technical certainty of the method of this invention.

[0044] The binarization formula is shown in equation (5): In the formula: These are the pixel values ​​of a binary image, where 1 represents white and 0 represents black. These are the pixel values ​​of a grayscale image; The threshold value is set.

[0045] Step S23: Form a discrete data sequence of porosity with height. Porosity is extracted through gridded statistics. Image processing software is used to extract partition information from the binarized image, and the porosity within each grid is calculated. Based on the vertical differentiation characteristics of porosity, the peak interval caused by stress concentration at the interface between soft and hard rocks is identified, revealing the vertical differentiation characteristics and peak causes of porosity. This determines the critical interface height of lithological abrupt change and identifies the distribution height of fracture development peaks.

[0046] Step S3: Constructing the coupling equation. Based on the discrete data sequence of porosity with height obtained in Step S2, a characteristic equation for fracture development is established, which includes the following sub-steps: Step S31: Establishing the basic descriptive equation. An exponentially decaying base term is used to describe the overall decay trend of fracture development with height, and multiple Gaussian superposition terms are used to describe the peak fracture density at different lithological interfaces, combined with a constant background value. The basic expression of the coupling equation is Equation (6): In the formula: For height Porosity at that location; It is an exponentially decaying term; For Gaussian terms; The initial strength for fracture development; The attenuation coefficient reflects the stress dissipation rate; For the first The area of ​​a superimposed Gaussian peak; For the first The distribution height of each fracture peak; For the first The width of each peak crack distribution; Baseline porosity; denoted as the number of Gaussian superposition terms.

[0047] Step S32: Introduce the piecewise dynamic attenuation coefficient. To accurately describe the changes in stress dissipation rate in different rock strata, the piecewise dynamic attenuation coefficient is defined as shown in equation (7): In the formula: The segmented attenuation coefficient; This represents the critical interface height for lithological abrupt changes.

[0048] Step S33: Perform continuity correction. Correct the initial strength of the exponentially decaying base term at the critical interface height of lithological abrupt change to ensure the function is continuous at the segmentation point. When correcting the exponentially decaying base term, use the end value of the previous segment as the initial strength coefficient of the next segment at the segmentation point. The final corrected expression is Equation (8): In the formula: The first paragraph, the second paragraph, and so on are respectively the first paragraph, the second paragraph, and so on. The function value of the exponentially decaying base term at the segmentation point.

[0049] Step S4, Parameter Inversion and Accuracy Evaluation. This includes the following sub-steps: Step S41, Parameter Optimization Inversion. The coupled equations are fitted using a nonlinear least squares optimization algorithm. The algorithm employs damped least squares for parameter optimization, with the Levenberg-Marquardt algorithm being preferred. The fitting analysis aims to minimize the error between the predicted value and the discrete data sequence of porosity, obtaining the optimal characteristic parameters and outputting the fitted curve and parameter estimates.

[0050] Step S42: Model Accuracy Evaluation. The model fitting effect is evaluated using residual distribution, root mean square error, and correlation coefficient as evaluation indicators. When the root mean square error is significantly lower than the data fluctuation range and the correlation coefficient is greater than 0.95, the model fitting effect is deemed satisfactory. If the evaluation is not satisfactory, the coupling equation is corrected for background values ​​or continuity, and parameter inversion is re-executed until the evaluation is satisfactory. Then, the characteristic equation for the development and water accumulation evolution of goaf fractures in the goaf is output for predicting the water-conducting channel.

[0051] Specifically, the root mean square error (RMSE) must be less than 55% of the minimum value of the measured discrete sequence of porosity to be considered significantly below the data fluctuation range. The physical meaning of this standard is that, due to the extremely strong heterogeneity of the roof crack distribution in the goaf, the background noise of the measured data is relatively large. If the RMSE can be controlled within approximately half of the minimum ripple in the measured data, it indicates that the equation not only captures the macroscopic attenuation trend but also accurately identifies the characteristics of minute local crack development. When the RMSE meets the standard and the correlation coefficient is... If the model evaluation is deemed satisfactory, the process is to return to step S3 to correct the background value C of the coupling equation or adjust the continuity at the segmentation points until the above convergence conditions are met and the characteristic equation is output.

[0052] As an example, the technical solution of the present invention will be described in detail below with reference to a specific example of the 7201 working face of a certain mine.

[0053] Step S1: Obtain key hydrogeological parameters of the goaf roof. In this embodiment, network parallel electrical resistivity CT technology is used to dynamically monitor the resistivity of the 7201 working face. By interpreting the resistivity data, the development pattern of the fracture network, the extent of the goaf, and the development of fracture conductivity in the roof are initially extracted. Figure 2 As shown, Figure 2 The resistivity monitoring results show the characteristics of fracture zone development in the coal seam roof. Among them, Figure 2 The left sub-figure (a) is a cloud map of the original background resistivity distribution before coal seam mining. At this stage, the resistivity distribution of the roof rock strata is relatively uniform and the value is relatively high, indicating that the rock strata are relatively intact, the original fractures are very few and the water accumulation conditions are not obvious. Figure 2The right-hand subplot (b) is a cloud map showing the evolution response of roof resistivity after coal seam mining. Monitoring reveals a significant low-resistivity anomaly zone in the roof strata, exhibiting a "zonal" characteristic. This indicates that as the coal seam is mined, the roof strata are fractured due to mining activities, and a fracture network gradually develops, leading to a local decrease in resistivity. Figure 2 The resistivity gradient variation characteristics in (b) were delineated, and the development height of the roof collapse zone of the 7201 working face was finally determined to be 13m, and the development height of the fracture zone was 40m. Based on the actual mining thickness of 3.6m, the fracture-to-mining ratio was 11.1:1.

[0054] Step S2: Conduct physical simulations of similar materials. For example... Figure 3 As shown, based on the fracture zone height of approximately 40m, the geometric similarity coefficient is set to 100:1. The strike length of the 7201 working face in a certain mine is approximately 1070m, and the dip length is approximately 160m. The dimensions of the test platform for this simulation experiment are length × height × width = 300cm × 200cm × 30cm. The model is designed with a total length of 300cm, a total width of 30cm, a partial height of 110cm, and a total height of 200cm. A portion of the Quaternary loose strata is included at the top of the model. Due to size limitations, the actual thickness of the Quaternary strata (approximately 250m on average) cannot be completely covered. Therefore, under experimental conditions, an additional compensation load of 18.9 kPa is applied to the top to simulate actual conditions. Simultaneously, 50cm protective coal pillars are reserved at each end of the model, with a total mining length of 200cm. Based on the fact that the lithology of the working face roof is mostly fine sandstone, medium sandstone, mudstone, and sandy mudstone, river sand was selected as the aggregate for this simulation, and gypsum, kaolin, and calcium carbonate were used as the cementing materials. After mixing and compaction, the mixture served as the material for the model strata, with a unit weight of approximately 2.45 × 10⁴ N / m³. 3 The similarity coefficient of the model strength, including compressive and tensile strength, is 490:3. The specific strength parameters are shown in Table 1.

[0055] Table 1

[0056]

[0057] The similarity model was generalized into 12 main strata, with river sand used to simulate the actual geological environment, and loads applied to the topmost rock strata to reproduce the real working conditions as closely as possible. Table 2 shows the stratigraphic generalization and material composition of the similarity model.

[0058] Table 2

[0059]

[0060] Based on the parameters, overall structure, and actual mining conditions of the similar model, after determining the laying scheme for each layer, the raw materials were mixed according to the design ratio, and an appropriate amount of water was added and stirred evenly to prepare a similar material. This material was then laid into a similar physical model to simulate the roof fracture development under the mining conditions of the No. 7 coal seam. According to actual conditions, the mine's average daily mining volume is approximately 4 m³ / d. The simulated mining started 50 cm from the left side and ended 250 cm from the right side, advancing an average of 4 cm per hour. Continuous mining was adopted to synchronize with the actual mining process as much as possible. Simulation results show that the caving zone of the No. 7 coal seam has a development height of approximately 12 m, and the fracture zone is approximately 44 m. Assuming a mining thickness of 4 m (4 cm in the model), the fracture-to-mining ratio is 11:1. Finally, the mining thickness: caving zone: fracture zone ratio under the mining conditions of the No. 7 coal seam in a certain coal mine was determined to be 1:3:11. Figure 4 As shown, the simulation results are processed by image grayscale and binarization, with the black areas in the figure representing cracks.

[0061] like Figure 5 As shown, porosity is extracted through gridded statistics, and partition information is extracted from the binarized image using an image processing program, revealing the vertical differentiation and peak causes of porosity. Figure 5 In this context, Avg represents the average porosity, and Cv represents the coefficient of variation. Specifically, Figure 5 The Avg value in the rightmost column represents the average porosity of all grids within the corresponding height layer. The change of this value with the layer height directly reflects the overall attenuation trend and local peak characteristics of fracture development in the vertical direction. Figure 5 The Avg value in the bottom row represents the average porosity of all grids at the corresponding horizontal position, reflecting the intensity of fracture distribution along the strike. The Cv value is calculated by dividing the standard deviation by the average value and is used to quantify the spatial heterogeneity of fracture distribution; a larger Cv value indicates a higher degree of non-uniformity in fracture development within that layer or region. Based on the above statistical indicators, the vertical differentiation and peak formation of porosity are revealed.

[0062] Step S3: Construct the coupling equation. In the data obtained in this embodiment, there are two peak values ​​for fracture development, with distribution heights of 10.95m and 25.55m, respectively. Assuming that the strata far from the coal seam that have not been disturbed by mining are aquitards, their porosity background value is very small and can be ignored, i.e., close to 0. The critical interface heights for lithological abrupt changes are 7m and 24m, respectively. Substituting the above known data into the descriptive equation, we obtain the specific equation for this embodiment as equation (9): .

[0063] Step S4: Parameter inversion and accuracy evaluation. For example... Figure 6 As shown, perform a fitting analysis and output the fitted curve and parameter estimates.

[0064] The estimated parameter values ​​for this embodiment are shown in Table 3.

[0065] Table 3

[0066]

[0067] The residual distribution of the fitted model results is shown in Table 4.

[0068] Table 4

[0069]

[0070] Residual analysis shows that the maximum absolute deviation between the model's predicted values ​​and the measured values ​​in this embodiment is 0.42%, and the minimum deviation is only 0.04%. The residuals are evenly distributed and there is no systematic shift. The specific calculation process is as follows: First, the measured values ​​at each height monitoring point are calculated. with fitted value absolute deviation To measure local error, after checking the data in Table 4, the height... The absolute deviation reaches its maximum value. At high altitude The absolute deviation is the minimum value Secondly, calculate the root mean square error. To measure overall accuracy, the sum of the squared residuals at each point in Table 4 is calculated, and the total is... Substitute the sample number Calculation .

[0071] Root mean square error The coefficient of determination (0.1912% < 0.35% × 55%) is significantly lower than the data fluctuation range of 0.35% to 5.24%, indicating that the model has high predictive accuracy. The coefficient of determination of the exponentially decaying Gaussian mixture model reaches [value missing]. This indicates that the fitting model can explain 95.9% of the data variation patterns, verifying its high characterization ability for fracture development patterns.

[0072] The foregoing description of the exemplary implementation of the quantitative characterization method for mining-induced fractures in the goaf roof proposed by the present invention has been described in detail with reference to preferred embodiments. However, those skilled in the art will understand that various modifications and alterations can be made to the above specific embodiments without departing from the concept of the present invention, and various combinations can be made to the various technical features and structures proposed by the present invention without exceeding the protection scope of the present invention, which is determined by the appended claims.

[0073] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for depicting and quantitatively characterizing mining-induced fractures in a goaf roof, characterized in that, The method comprises the following steps: Step S1, obtaining key hydrogeological parameters of the goaf roof, using network parallel electric method CT technology to dynamically monitor the resistivity of the coal mine working face, realizing high-resolution imaging of the underground structure by arranging multiple electrode arrays, interpreting the resistivity data to extract the development law of the fracture network, the range of the goaf, and the high development of the roof fracture, obtaining the development height of the fracture zone; Step S2, carrying out similar material physical simulation, setting up a similar model according to the development height of the fracture zone, carrying out image gray scale and binary processing on the simulation results, extracting the void ratio through gridding statistics to form a discrete data sequence of the void ratio with the height, and identifying the distribution height of the fracture development peak value and the critical interface height of the lithology mutation; Step S3, constructing a coupling equation, based on the discrete data sequence of the void ratio with the height, establishing a fracture development characteristic equation including an exponential decay base term, multiple Gaussian superposition terms and a constant background value, wherein the exponential decay base term describes the overall decay trend of the fracture development with the height, the multiple Gaussian superposition terms describe the fracture density peak values at different lithology interfaces, a segmented dynamic decay coefficient is introduced to describe the change of the stress dissipation rate in different rock layers, and the initial intensity of the exponential decay base term is modified at the critical interface height of the lithology mutation, so that the function is continuous at the segmentation point; Step S4, parameter inversion and precision evaluation, based on a nonlinear least squares optimization algorithm, fitting analysis is performed on the coupling equation, the fitting curve and the parameter estimate value are output, the residual distribution, the root mean square error and the correlation coefficient are used as evaluation indexes to evaluate the fitting effect of the model, if the evaluation is not up to standard, the coupling equation is modified and the parameter inversion is re-executed, and when the evaluation is up to standard, the coupling equation is output for goaf water accumulation evolution prediction and water channel identification.

2. The method of claim 1, wherein, In the step S1, the recharge, discharge and drainage conditions of the underground water in the mining area, the stratum thickness, the water-resisting layer relationship, the drilling data and the monitoring well water level are also obtained, and the resistivity data interpretation results are verified by using the hydrogeological drilling site monitoring data.

3. The method of claim 1, wherein, In the step S2, the similar material physical simulation meets the geometric similarity, mechanical similarity and dynamic similarity conditions, wherein: The calculation formula of the geometric similarity coefficient is: ； wherein: is a geometric similarity factor; is a geometric dimension of the real geological prototype; is a geometric dimension of the similar model; The calculation formula of the mechanical similarity coefficient is: ； wherein: is a mechanical similarity factor; is the stress of the real geological prototype; is the stress of the similar model; is a material density similarity factor; The calculation formula of the time similarity coefficient is: 。 4. The method of claim 1, wherein, In the step S2, the image gray scale adopts a weighted average method, and the conversion formula is: ； In the formula: is a gray value; , , are red, green, and blue channel values of the pixel, respectively. The formula of the binary processing is: ； wherein: is a binary image pixel value, 1 represents white color, and 0 represents black color; is a pixel value of a grayscale image; is a set threshold value.

5. The method of claim 1, wherein, In the step S2, the binary image is divided into zones by using an image processing program, the void ratio in each grid is calculated, the peak value interval caused by stress concentration at the soft and hard rock interface is identified according to the vertical differentiation characteristics of the void ratio, the vertical differentiation characteristics of the void ratio and the peak value causes are revealed, and thus the critical interface height of the lithology mutation is determined.

6. The method of claim 1, wherein, In the step S3, the basic expression of the coupling equation is: ； where: is the height of the voids at the location; is the initial intensity of fracture development; is the decay coefficient, reflecting the stress dissipation rate; is the area of the th Gaussian superposition peak; is the distribution height of the th fracture peak value; is the width of the th peak fracture distribution; is the baseline void ratio; is the number of Gaussian superposition terms.

7. The method of claim 6, wherein, In the step S3, the expression of the segmented dynamic decay coefficient is: ； In the formula: is the segmental attenuation coefficient; is the critical interface height of the lithology mutation.

8. The method of claim 7, wherein, In the step S3, the exponential decay base term is modified, the end value of the previous segment is taken as the initial intensity coefficient of the next segment at the segmentation point, and the final expression after modification is: ； wherein: are the first segment, the second segment to the the function value of the exponentially decaying base term at the segment point.

9. The method of claim 1, wherein, In the step S4, the nonlinear least square optimization algorithm adopts a damped least square method to optimize parameters, performs fitting analysis aiming at minimizing the error between the predicted value and the porosity discrete data sequence, and obtains optimal characteristic parameters. The damped least square method selects a Levenberg-Marquardt algorithm.

10. The method of claim 1, wherein, In the step S4, when the root mean square error is lower than the data fluctuation range and the correlation coefficient is greater than 0.95, it is determined that the model fitting effect evaluation meets the standard. If the evaluation does not meet the standard, the background value correction or continuity correction is performed on the coupling equation, and the parameter inversion is re-executed until the evaluation meets the standard, and then the equation for predicting the water evolution characteristics of the goaf fracture in the water channel is output.

Images