A method for depicting and quantitatively representing mining cracks of a goaf roof

Abstract

The application discloses a kind of goaf roof mining fracture delineation quantitative characterization method, belong to mine hydrogeology and water prevention and control engineering technical field, by network parallel electric method CT field monitoring obtains the development height of guide crack zone, carries out similar material physical simulation and carries out gray scale and binary processing to image, extracts void fraction discrete data sequence, establishes the coupling equation including exponential decay base item, multi-peak gaussian item and constant background value, introduces subsection dynamic attenuation coefficient to describe the difference of stress dissipation in different rock strata, continuously revises initial strength at lithology mutation interface, adopts damping least square method to carry out parameter inversion, with root mean square error and correlation coefficient as evaluation index, output characteristic equation for water accumulation evolution prediction and water channel identification.The application has high fitting precision and clear physical meaning, and can provide reliable parameter support for water prevention and control engineering design.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of mine hydrogeology and water prevention and control engineering, and particularly relates to a method for depicting and quantitatively representing mining-induced fractures of a goaf roof. BACKGROUND

[0002] In the process of coal mining, the overburden strata are structurally damaged by stress adjustment and gravity, and a complex fracture system such as a caving zone, a fracture zone and a bending zone is formed from bottom to top. Among them, the roof fracture zone is the main water-conducting channel between the goaf and the overlying aquifer, and its development law directly determines the water storage capacity, hydrodynamic response and water inrush risk of the goaf. Accurate description of the spatial distribution characteristics of the fracture zone is of great significance for water disaster prevention and control of the goaf.

[0003] At present, the research on the structural characteristics of the fracture zone mostly adopts empirical formula or single-peak function fitting method, and common ones include exponential decay function, power function or Gaussian function model. These models have certain applicability in describing the overall trend of fracture development, but have obvious shortcomings in describing the fracture distribution law under complex geological conditions.

[0004] The main problems of the traditional model are reflected in several aspects. First, the coal roof is usually composed of different lithological strata, and the stress distribution will change abruptly at the interface of soft and hard rocks, resulting in the multi-peak characteristics of fracture development, while the single-peak function cannot accurately describe this local anomaly. Second, most existing models assume that the decay coefficient is constant, and fail to reflect the dynamic decay law of fracture development with increasing height in different rock layers, ignoring the different influences of different lithologies on stress transmission and fracture response. Third, the parameters of many models lack practical physical meaning, and are obtained only by mathematical fitting, lacking coupling verification with field monitoring data and physical simulation tests, which limits the reliability and generalizability of the model. In addition, the traditional method often separates the overall decay trend and local peak characteristics when describing the fracture development, and it is difficult to form a unified mathematical expression.

[0005] These shortcomings make the existing model have a large error when providing parameter input for goaf water accumulation evolution prediction and hydrodynamic analysis, affecting the scientificity and pertinence of water prevention and control engineering design. Therefore, it is urgent to establish a quantitative description equation that can reflect the overall decay trend and local lithology mutation response of fracture development at the same time, and improve the accuracy and practicality of the model through multi-source data coupling verification, so as to realize the fine characterization of goaf fracture development and water accumulation evolution process. SUMMARY

[0006] The purpose of the present application is to overcome the deficiencies in the prior art, provide a goaf roof mining fracture sketch quantitative characterization method, determine the development height of the crack guide zone through network parallel electric method CT field monitoring, combine similar material physical simulation to extract crack porosity discrete data sequence, establish an exponential decay correction multi-peak Gaussian mixed equation, introduce a segmented dynamic attenuation coefficient to describe the difference in stress dissipation rate in different rock layers, and continuously correct the initial strength at the lithology mutation interface, use the damped least squares method for parameter inversion optimization, and finally output the characteristic equation with clear physical meaning for goaf water accumulation evolution prediction and water channel identification, providing reliable parameter support for water prevention and control engineering design.

[0007] To achieve the above purpose, the present application provides a goaf roof mining fracture sketch quantitative characterization method, comprising 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, extracting the development law of the fracture network, the range of the goaf and the development of the roof fracture guide height, and obtaining the development height of the crack guide zone.

[0008] Step S2, carrying out similar material physical simulation, setting a similar model according to the development height of the crack guide zone, carrying out image gray scale and binary processing on the simulation results, extracting the porosity through gridding statistics, forming the discrete data sequence of the porosity with the height, and identifying the distribution height of the fracture development peak value and the critical interface height of the lithology mutation.

[0009] Step S3, constructing a coupled equation, based on the discrete data sequence of the porosity 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 value at different lithology interfaces, a segmented dynamic attenuation coefficient is introduced to describe the change in stress dissipation rate in different rock layers, and the initial strength of the exponential decay base term is corrected at the critical interface height of the lithology mutation, so that the function is continuous at the segmented point.

[0010] Step S4, parameter inversion and precision evaluation, based on the nonlinear least squares optimization algorithm, fitting analysis is carried out on the coupled 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 coupled equation is corrected and the parameter inversion is re-executed, and when the evaluation is up to standard, the coupled equation is output for goaf water accumulation evolution prediction and water channel identification.

[0011] Further, in the step S1, the supplement, runoff and discharge 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 obtained, and the hydrogeological drilling field monitoring data are used to verify the interpretation result of the resistivity data.

[0012] Further, in the step S2, the physical simulation of the similar material satisfies the geometric similarity, the mechanical similarity and the dynamic similarity conditions, wherein, the calculation formula of the geometric similarity coefficient is: ; in the formula, is the geometric similarity coefficient; is the geometric size of the real geological prototype; is the geometric size of the similar model.

[0013] The calculation formula of the mechanical similarity coefficient is: ; in the formula, is the mechanical similarity coefficient; is the stress of the real geological prototype; is the stress of the similar model; is the material density similarity coefficient.

[0014] The calculation formula of the time similarity coefficient is: .

[0015] Further, in the step S2, the image gray scale is converted by using the weighted average method, and the conversion formula is: ; in the formula, is the gray scale value; , , are respectively the red, green and blue channel values of the pixel.

[0016] The formula of the binary processing is: ; in the formula, is the pixel value of the binary image, 1 represents white, and 0 represents black; is the pixel value of the gray scale image; is the set threshold value.

[0017] Further, in the step S2, the image processing program is used to extract the partition information of the binary image, to calculate the void ratio in each grid, to identify the peak interval caused by the stress concentration of the soft and hard rock interface according to the vertical differentiation characteristics of the void ratio, to reveal the vertical differentiation characteristics of the void ratio and the peak value causes, and to determine the lithology mutation critical interface height.

[0018] Further, in the step S3, the basic expression of the coupling equation is: ; in the formula, is the void ratio at the height ; is the initial intensity of the fracture development; is the attenuation coefficient, reflecting the stress dissipation rate; is the area of the first Gaussian superposition peak; is the distribution height of the first fracture peak value; is the width of the first peak fracture distribution; is the baseline porosity; is the number of Gaussian superposition terms.

[0019] Further, in the step S3, the expression of the segmented dynamic attenuation coefficient is: ; in the formula: is the segmented attenuation coefficient; is the critical interface height of the lithology mutation.

[0020] Further, in the step S3, the exponential decay base term is modified, and the end value of the previous segment is taken as the initial intensity coefficient of the next segment at the segmentation point. The final modified expression is: ; in the formula: are the function values of the exponential decay base term of the first segment, the second segment to the segment at the segmentation point, respectively.

[0021] Further, in the step S4, the nonlinear least squares optimization algorithm adopts a damped least squares method to optimize the parameters, so as to minimize the error between the predicted value and the porosity discrete data sequence, perform fitting analysis, and obtain optimal characteristic parameters. The damped least squares method is preferably a Levenberg-Marquardt algorithm.

[0022] Further, in the 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, it is determined that the model fitting effect evaluation meets the standard. If the evaluation does not meet the standard, the background value of the coupling equation is corrected or the continuity is corrected, and the parameter inversion is re-executed until the evaluation meets the standard, and then the goaf fracture development water accumulation evolution characteristic equation for predicting the water channel is output.

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

[0024] (1) The present application establishes an exponential decay corrected multi-peak Gaussian mixed model by using network parallel electrical method CT technology for resistivity dynamic monitoring and combining with similar material physical simulation, the model describes the overall decay trend of fracture development with height through the exponential decay base item, accurately captures the fracture density peak value at the soft and hard rock interface through multiple Gaussian superposition items, reflects the change of stress dissipation rate in different rock layers by introducing a segmented dynamic decay coefficient, and corrects the initial strength at the critical interface of lithology mutation to maintain the continuity of the function, thereby realizing the unified mathematical description of the fracture heterogeneity and multi-peak characteristics. This method overcomes the defects of traditional single-peak function that cannot describe the stress mutation phenomenon of lithology interface, so that the model meets the hydrogeological constraint of the decay trend approaching the background value with height on the macro level, and accurately describes the local peak response on the micro level. Compared with traditional models such as polynomial function, simple exponential function, power function or single-peak Gaussian function, the present application significantly reduces the systematic error and fitting residual, and improves the interpretation rate and prediction accuracy of the heterogeneous fracture field.

[0025] (2) The present application obtains key hydrogeological parameters such as the recharge, runoff and discharge conditions of underground water in the mining area, stratum thickness, water-resisting layer relationship, drilling data and monitoring well water level, and verifies the resistivity data interpretation results by using the hydrogeological drilling site monitoring data, to establish 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 a background value with actual geological significance for subsequent parameter inversion, avoids the randomness of model parameters, ensures the consistency of the constructed equation with the actual geological conditions, and makes the model have stronger physical basis and engineering applicability.

[0026] (3) The present application establishes a similar material physical simulation that meets the conditions of geometric similarity, mechanical similarity and dynamic similarity, and uses weighted average method for image gray scaling and threshold-based binaryzation processing, extracts the void ratio formation with height through grid statistics, calculates the void ratio in each grid by extracting the partition information of the binaryzation image using an image processing program, and reveals the vertical differentiation characteristics and peak value causes of the void ratio. This quantitative image analysis method converts the physical simulation test results into discrete data that can be used for mathematical modeling, provides a clear physical basis for the setting of multi-peak items, makes the determination of fracture development peak value distribution height and lithology mutation critical interface height have traceability and repeatability, greatly improves the scientificity and standardization of model construction, and is convenient for popularization and application under similar geological conditions.

[0027] (4) The nonlinear parameter optimization of the coupling equation is carried out by adopting the damping least square method combining the gradient descent method and the Gauss-Newton method, and a strict precision evaluation system is established by taking the residual distribution, the root mean square error and the correlation coefficient as evaluation indexes, and the model evaluation is determined to be up to standard when the root mean square error is significantly lower than the data fluctuation range and the correlation coefficient is greater than 0.95. The method ensures the global optimality and numerical stability of the fitting parameters, avoids the phenomena of overfitting or underfitting of the model through the quantitative evaluation standard, guarantees the prediction reliability of the output equation, and enables the constructed fracture development and waterlogging evolution characteristic equation to provide high-precision and high-confidence parameter support for the prediction of goaf water evolution, the identification of water conducting channels and the design of water prevention and control engineering, and significantly improves the scientific decision-making level and engineering safety guarantee capability of the mine water disaster prevention. BRIEF DESCRIPTION OF DRAWINGS

[0028] The accompanying drawings, which constitute a part of this specification, illustrate embodiments of the present application and serve to explain the principles of the present application.

[0029] The present application can be understood more readily by reference to the following detailed description, when considered in connection with the attached drawings, wherein:

[0030] Figure 1 A technical flow chart of a goaf roof mining fracture sketch quantitative characterization method provided for the embodiments of the present application.

[0031] Figure 2 A coal seam roof fracture zone development characteristic monitoring result schematic diagram.

[0032] Figure 3 A similar model profile.

[0033] Figure 4 A similar simulation result binary processing image schematic diagram.

[0034] Figure 5 A binary image partition identification result schematic diagram.

[0035] Figure 6 A roof fracture zone porosity with height change measured value and fitting curve schematic diagram. DETAILED DESCRIPTION

[0036] The technical solutions in the embodiments of the present application will be described clearly and completely below with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments. The following description of at least one example embodiment is merely illustrative in nature and not intended to further limit the present application or its application or uses. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative efforts fall within the scope of the present application. Unless otherwise specified, the relative arrangement, expressions and values of components and steps set forth in these embodiments do not limit the scope of the present application. The technology, methods and devices known to those of ordinary skill in the relevant art can not be discussed in detail, but should be considered as part of the specification in appropriate cases. In all examples shown and discussed herein, any specific value should be interpreted as merely exemplary, rather than limiting. Therefore, other examples of the example embodiments can have different values.

[0037] The term "and / or" herein merely describes an association relationship of associated objects, which means that there can be three relationships, for example, A and / or B can mean that A exists alone, A and B exist together, and B exists alone. In addition, the character " / " herein generally means that the associated objects before and after the " / " have an "or" relationship.

[0038] As shown in Figure 1 , the present application provides a method for depicting and quantitatively characterizing mining-induced fractures of a goaf roof, comprising the following steps: step S1, obtaining key hydrogeological parameters of the goaf roof. Network parallel electric method CT technology is used to dynamically monitor the resistivity of the coal mine working face, high-resolution imaging of the underground structure is achieved by arranging multiple electrode arrays, the development law of the fracture network, the range of the goaf and the development of the roof fracture guide high are extracted by interpreting the resistivity data, and the development height of the guide fracture zone is obtained. At the same time, the recharge-discharge 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 obtained, and the resistivity data interpretation results are verified by using the hydrogeological drilling site monitoring data.

[0039] Step S2, carrying out physical simulation of similar materials. The development height of the guide fracture zone obtained in step S1 is used to set up a similar model, which specifically comprises the following steps: step S21, constructing a physical model that meets 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 unified proportion.

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

[0041] Mechanical similarity coefficient Related to the geometric similarity coefficient and the material density similarity coefficient, the calculation formula is formula (2): ; In the formula: is the mechanical similarity coefficient; is the stress of the real geological prototype; is the stress of the similar model; is the material density similarity coefficient. Dynamics similarity requires that the proportional relationship of time, fluid parameters and permeability and other related physical quantities be ensured, so that the displacement, deformation and other dynamic processes in the similar material model are unified with the principle of the geological prototype.

[0042] Time similarity coefficient Related to the geometric similarity coefficient, the calculation formula is formula (3): .

[0043] Step S22, image gray scale and binary processing are performed on the simulation result. The commonly used gray scale method is the weighted average method. In view of the fact that the human eye is most sensitive to green and least sensitive to blue, in order to maximize the brightness contrast information of the rock layer texture and fissure in the gray scale conversion, the application selects a weighted coefficient that conforms to the visual characteristics of the human eye, and the conversion formula is formula (4): ; In the formula: is the gray value; , , respectively, the red, green and blue channel values of the pixel. Using this specific coefficient can effectively suppress image background noise and improve the accuracy of fissure edge recognition in subsequent binary processing. It should be noted that the coefficient of this formula is not the only coefficient in the field of image gray scale, but it is the optimal coefficient adapted to the application scenario of the application. Selecting this coefficient strengthens the repeatability and technical certainty of the method of the application.

[0044] The binary formula is formula (5): ; In the formula: is the pixel value of the binary image, 1 represents white, and 0 represents black; is the pixel value of the gray scale image; is the set threshold value.

[0045] Step S23, forming a discrete data sequence of porosity with height. The porosity is extracted by gridding statistics, the partition information of the binary image is extracted by using an image processing program, the porosity in each grid is calculated, the peak interval caused by stress concentration of the soft and hard rock interface is identified according to the vertical differentiation characteristics of the porosity, the vertical differentiation characteristics and the peak value cause are revealed, so as to determine the critical interface height of lithology mutation, and identify the distribution height of the fracture development peak value.

[0046] Step S3, constructing a coupling equation. Based on the discrete data sequence of porosity with height obtained in step S2, a fracture development characteristic equation is established, which specifically includes the following steps: step S31, establishing a basic description equation. An exponential decay base term is used to describe the overall decay trend of the fracture development with height, a plurality of Gaussian superposition terms are used to describe the fracture density peak values at different lithology interfaces, and a constant background value is combined. The basic expression of the coupling equation is formula (6): ; in the formula: is the porosity at the height ; is the exponential decay term; is the Gaussian term; 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 porosity; is the number of Gaussian superposition terms.

[0047] Step S32, introducing a segmented dynamic decay coefficient. In order to accurately describe the change of stress dissipation rate in different rock layers, a segmented dynamic decay coefficient is defined as formula (7): ; in the formula: is the segmented decay coefficient; is the critical interface height of lithology mutation.

[0048] Step S33, performing continuity correction. The initial intensity of the exponential decay base term is corrected at the critical interface height of lithology mutation, so that the function is continuous at the segmentation point. When the exponential decay base term is corrected, 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 correction is formula (8): ; in the formula: are the function values of the exponential decay base term of the first segment, the second segment to the th segment at the segmentation point, respectively.

[0049] Step S4, parameter inversion and precision evaluation. Specifically comprising the following steps: Step S41, parameter optimization inversion. Based on the nonlinear least squares optimization algorithm, the fitting analysis is performed on the coupling equation, the nonlinear least squares optimization algorithm adopts the damped least squares method for parameter optimization, the Levenberg-Marquardt algorithm is preferably selected to minimize the error between the predicted value and the void ratio discrete data sequence as the target of the fitting analysis, and the optimal characteristic parameter is obtained, and the fitting curve and the parameter estimation value are output.

[0050] Step S42, model precision evaluation. The residual distribution, root mean square error and correlation coefficient are used as evaluation indexes to evaluate the model fitting effect. When the root mean square error is significantly 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 of the coupling equation is corrected or the continuity is corrected, and the parameter inversion is re-executed until the evaluation meets the standard, and the equation for predicting the water evolution characteristics of the goaf fracture of the water channel is output.

[0051] Specifically, the root mean square error RMSE value needs to be less than 55% of the minimum value of the measured void ratio discrete sequence, so as to be determined as significantly lower than the data fluctuation range. The physical meaning of this standard is that: since the roof fracture distribution of the goaf has a strong heterogeneity, the background noise of the measured data is large. If the root mean square error RMSE can be controlled within about half of the minimum value of the measured data, that is, the minimum value, it means that the equation not only captures the macroscopic attenuation trend, but also accurately identifies the small local fracture development characteristics. When the root mean square error RMSE meets the standard and the correlation coefficient , it is determined that the model evaluation meets the standard; if the evaluation does not meet the standard, the background value C of the coupling equation needs to be corrected or the continuity at the segmentation point needs to be adjusted, until the characteristic equation is output after the convergence condition is met.

[0052] As an embodiment, the technical solutions of the present application will be described in detail in combination with a specific example of 7201 working face of a certain mine.

[0053] Step S1, obtaining key hydrogeological parameters of the goaf roof. In this embodiment, the network parallel electric method CT technology is used to dynamically monitor the resistivity of 7201 working face, and the development law of the fracture network, the range of the goaf and the roof fracture high development are preliminarily extracted through the interpretation of the resistivity data. As shown in Figure 2 , Figure 2 shows the comparison results of resistivity monitoring of the development characteristics of the coal seam roof fracture zone. Among them, Figure 2 The left subgraph (a) is the original background resistivity distribution cloud chart before the coal seam is mined. At this stage, the resistivity distribution of the roof rock layer is relatively uniform and the value is high, indicating that the rock layer has good integrity, the original fracture development is very small and the water accumulation condition is not obvious. Figure 2The right subgraph (b) is the roof resistivity evolution response cloud chart after coal seam mining. The monitoring shows that the roof rock layer resistivity appears obvious low resistivity anomaly area and presents "zoning" characteristics, indicating that with the coal seam mining, the roof rock layer is affected by mining and broken, the fracture network gradually develops and leads to local resistivity drop. According to the variation characteristics of the resistivity gradient in (b), the development height of the roof caving zone of the 7201 working face is finally determined to be 13 m, and the development height of the guide crack zone is 40 m. According to the actual mining thickness of 3.6 m, the crack mining ratio is 11.1:1. Figure 2

[0054] Step S2, carry out similar material physical simulation. As shown in Figure 3 , according to the guide crack zone height of about 40 m, the geometric similarity coefficient is set to 100:1. The strike length of the 7201 working face of a certain mine is about 1070 m, and the inclined length is about 160 m. The size of the test bed for this simulation test is length x height x width = 300 cm x 200 cm x 30 cm. The full length of the model is 300 cm, the full width is 30 cm, the partial height is 110 cm, and the total height is 200 cm. A part of the Quaternary loose layer section is set at the top of the model. Due to size limitations, the actual thickness of the Quaternary system, which is about 250 m on average, cannot be completely covered. Therefore, under the allowable conditions of the test, an additional compensation load of 18.9 KPa is applied at the top to simulate the actual conditions. At the same time, 50 cm of protective coal pillars are reserved at both ends of the model, and the total length of the mining is 200 cm. According to the fact that the roof rock of the working face is mainly fine sandstone, medium sandstone, mudstone and sandy mudstone, river sand is used as the aggregate, gypsum, kaolin and calcium carbonate are used as the cementing material, and the mixed and compacted material is used as the material of the model stratum, with a bulk density of about 2.45 x 104 N / m 3 , and the model strength including compression and tension similarity coefficient is 490:3. The specific strength parameters are shown in Table 1.

[0055] Table 1

[0056]

[0057] This similar model is generalized to 12 main layer groups, in which river sand is used to simulate the occurrence environment of the actual stratum, and a load is applied to the topmost rock layer to restore the actual working conditions as much as possible. Table 2 shows the stratum generalization and material proportioning of the similar model.

[0058] Table 2

[0059]

[0060] ​According to the parameters, overall structure and actual mining conditions of the similar model, the laying scheme of each layer is determined, the raw materials are mixed in the designed proportion, a proper amount of water is added and stirred uniformly to prepare the similar material and lay the similar physical model to simulate the roof fracture development law under the mining conditions of the 7 coal. According to the actual situation, the average daily mining of the mine is about 4m / d, the simulation mining starts from the left side 50cm and ends at the right side 250cm, the average hourly advance is 4cm, and the continuous mining mode is adopted to synchronize the actual mining process as far as possible. The simulation results show that the caving zone development height of the 7 coal is about 12m, the guide crack zone is about 44m, and the crack mining ratio is 11:1, taking the mining thickness of 4m in the model as an example. The final determination of the mining thickness of the 7 coal under the mining conditions of the coal mine is 1:3:11, as shown in the following table. Figure 4 As shown in the following table, the simulation results are subjected to image gray scale and binary processing, and the black area in the figure is the fracture.

[0061] As shown in the following table, the simulation results are subjected to image gray scale and binary processing, and the black area in the figure is the fracture. Figure 5 As shown in the following table, the simulation results are subjected to image gray scale and binary processing, and the black area in the figure is the fracture. Figure 5 In the table, Avg represents the average void ratio, and Cv represents the coefficient of variation. Specifically, Figure 5 The Avg value of the rightmost column in the table represents the average value of all grid void ratios in the corresponding height layer, and this value directly reflects the overall attenuation trend and local peak characteristics of the fracture development in the vertical direction; Figure 5 The Avg value of the bottom row in the table represents the average value of all grid void ratios in the corresponding horizontal position, which reflects the distribution intensity of the fracture in the strike direction. The Cv value is calculated by dividing the standard deviation by the average value, which is used to quantify the spatial heterogeneity of the fracture distribution. The larger the Cv value, the higher the non-uniformity of the fracture development in the layer or region. Based on the above statistical indicators, the vertical differentiation and peak causes of the void ratio are revealed.

[0062] Step S3, coupling equation is constructed. In the data obtained in the embodiment, there are two fracture development peaks, and the distribution heights are 10.95m and 25.55m respectively. It is assumed that the strata far away from the coal seam are not damaged by mining, and the void ratio background value is very small and can be ignored, that is, it is considered to be close to 0, that is, The critical interface heights of lithology mutation are 7m and 24m respectively. The specific equation of the embodiment can be obtained by bringing the above known data into the description equation as formula (9): .

[0063] Step S4, parameter inversion and precision evaluation. As shown in the following table, the fitting analysis is performed, and the fitting curve and parameter estimation value are output. Figure 6

[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 merely describes the preferred embodiments of the present application, and it should be pointed out that, for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present application, and these improvements and modifications should also be considered as the protection scope of the present application.

Claims

1. A quantitative characterization method for delineating mining-induced fractures in the roof of a goaf, characterized in that, Includes 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, and 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. 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. 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. 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.

2. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, Step S1 further includes obtaining the groundwater recharge and drainage conditions, stratum 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.

3. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, 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; For the geometric dimensions of similar models; 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; The material density similarity coefficient; The formula for calculating the time similarity coefficient is: 。 4. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, In step S2, the image grayscale conversion uses a weighted average method, and the conversion formula is as follows: ； In the formula: Grayscale value; , , These are the red, green, and blue channel values ​​of the pixel, respectively; 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.

5. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, In step S2, the image processing program is used to extract partition information from the binarized image, calculate the porosity in 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.

6. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, 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 ​​the superimposed Gaussian peaks; 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.

7. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 6, characterized in that, 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.

8. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 7, characterized in that, 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 final corrected 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.

9. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, In step S4, the nonlinear least squares optimization algorithm uses damped least squares to optimize parameters. The goal is to minimize the error between the predicted value and the discrete data sequence of porosity with height to perform fitting analysis and obtain the optimal feature parameters. The damped least squares method is the Levenberg-Marquardt algorithm.

10. The method for quantitative characterization of mining-induced fractures in the roof of a goaf according to claim 1, characterized in that, In step S4, when the root mean square error is 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 fails to 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.

Images