Surrounding rock fracture identification and three-dimensional model reconstruction method

Abstract

The invention relates to the technical field of coal mine tunnel surrounding rock detection, and particularly provides a surrounding rock crack identification and three-dimensional model reconstruction method. The method comprises the steps of collecting a data set for preprocessing, performing fracture recognition based on an improved U-Net network, and outputting a fracture binary image; based on the fracture binary image, sinusoidal fracture quantitative analysis modeling is carried out; performing roughness coefficient JRC calculation and sampling method sensitivity analysis according to modeling; according to JRC calculation and sensitivity analysis, multi-borehole fracture network three-dimensional reconstruction and hole patching optimization are carried out, and the method improves the accuracy and reliability of surrounding rock fracture recognition and improves the engineering applicability.

Description

Technical Field

[0001] This invention relates to the field of coal mine roadway surrounding rock detection technology, and in particular to a method for identifying surrounding rock fissures and reconstructing a three-dimensional model. Background Technology

[0002] In the complex environment of underground coal mining, the stability of the surrounding rock in the roadway is always the core link between production efficiency and operational safety. As a naturally occurring structural defect in the underground rock mass, the degree of development and geometry of the surrounding rock fissures directly determine the integrity and bearing capacity of the rock mass. When the fissures are closed or slightly open, they may only affect the local rock mass strength; however, as the mining depth increases and the ground stress is redistributed, the fissures are prone to expand and penetrate, thereby inducing a series of disasters: at best, they can cause local collapse of the roadway roof, resulting in damage to the support structure and production interruption; at worst, they can trigger large-scale water inrush (fissures act as water channels, connecting the aquifer and the roadway) or gas outburst (the instantaneous release of high-concentration gas accumulated in the fissures), posing a fatal threat to the lives of underground personnel and core production equipment. Summary of the Invention

[0003] In view of this, the present invention provides a method for identifying surrounding rock fissures and reconstructing a three-dimensional model, in order to improve the accuracy and reliability of identifying surrounding rock fissures and to improve engineering applicability.

[0004] In a first aspect, the present invention provides a method for identifying surrounding rock fractures and reconstructing a three-dimensional model, the method comprising: Step 1: Collect the dataset, preprocess it, and perform crack identification based on the improved U-Net network to output a binary crack map; Step 2: Based on the binary fracture map, perform quantitative analysis and modeling of sinusoidal fractures; Step 3: Based on the modeling, perform roughness coefficient JRC calculation and sampling method sensitivity analysis; Step 4: Based on JRC calculations and sensitivity analysis, perform three-dimensional reconstruction and hole filling optimization of the multi-hole fracture network.

[0005] Optionally, step 1 includes: By comparing image preprocessing with multi-method fusion and optimizing and training the network structure, pixel-level separation of cracks and interference in borehole images can be achieved. Step 11: Data acquisition and preprocessing. Obtain the dataset of borehole images and perform grayscale conversion, Gaussian denoising, and histogram equalization on the images. By random cropping, the images are randomly cropped into 244×244 pixel samples, followed by rotation and mosaic enhancement to expand the dataset samples. Step 12: Multi-method fusion comparison, designing 5 methods: a) Employing multi-method edge feature extraction plus voting mechanism; b) Employing multi-dimensional feature extraction and weighted fusion, combined with adaptive threshold for crack detection; c) Employing machine learning algorithms, based on multi-dimensional feature extraction and random forest algorithm to classify crack pixels; d) Employing multi-scale feature extraction and KMeans clustering algorithm to cluster crack pixels; e) Employing deep learning algorithms for segmentation, using improved U-Net, U2Net, and YOLOv11 for training respectively; the improved U-Net network includes an encoder, decoder, multi-scale fusion layer, and output layer; Step 13: Model training and evaluation. The results are displayed in a visual form. Select the crack extraction results. First, divide the dataset into training set, validation set and test set in a 7:2:1 ratio. Evaluate the results by Intersection over Union (IoU) and Mean Average Precision (mAP) and output the binary crack map.

[0006] Optionally, it includes: Step 111: Grayscale conversion. The drill image is in RGB format. A weighted average method Y=0.299R+0.587G+0.114B is used to convert it into a grayscale image. Gaussian denoising is performed using a 5×5 Gaussian kernel for convolution. Histogram equalization is achieved using the contrast-limited adaptive histogram equalization algorithm CLAHE. The image is divided into an 8×8 grid, and histogram equalization is performed on each grid separately, while limiting contrast enhancement. Random cropping is used to randomly crop the image into 244×244 pixel regions to simulate different imaging perspectives. Rotation is used to randomly rotate the samples from 0-360°. Mosaic enhancement is achieved by stitching four different samples together in a 2×2 pattern, with linear transitions at the edges. Step 112: After edge extraction using three edge extraction methods—multi-directional gradient Sobel edge detection, Canny edge detection, LoG edge detection combined with Gaussian blur and Laplacian operator edge detection—a voting method is then used. If at least two detection results indicate a crack, then the crack is determined to be a crack. Step 113: After extracting three features—Sobel edge detection, LBP texture detection, and Hessian—and combining them with an adaptive threshold, the weights of the three features are 50% for multi-scale, 30% for Sobel edge detection, and 20% for Hessian, forming multi-scale and multi-directional features. Finally, the Otsu automatic thresholding is used to automatically calculate the threshold with the minimum intra-class variance between the foreground and background to obtain the final detection result. Step 114: Supervised learning based on random forest sequentially involves feature extraction, pseudo-label generation, forest training, and confidence filtering. Feature extraction is diverse, extracting multi-dimensional features such as original intensity, gradient magnitude, LBP texture, Canny edges, and Hessian matrix features. Randomly selecting some features to train each tree aims to explore the discriminative power of different feature combinations. Pseudo-label generation uses mathematical and statistical methods to generate pseudo-labels, transforming unsupervised or weakly supervised data into labeled data, providing a training foundation for supervised learning of the random forest. Forest training parameter design includes n_estimators=50, indicating 50 trees to control the ensemble size; max_depth=8, indicating a maximum tree depth of 8, to limit the complexity of a single tree and balance bias and variance. Confidence evaluation based on gradient and standard deviation filters the random forest predictions; the gradient reflects the stability of the predictions, and the standard deviation reflects the discrepancy between predictions between trees. Finally, the classifier learns the mapping relationship between features and fissure or non-fissure labels to classify new data. Step 115: Multi-scale feature extraction. Comprehensive features are obtained through multiple window sizes and structural tensors. Basic statistical features include calculating the mean, variance, standard deviation, and contrast of each window to capture local gray-scale distribution patterns. Structural tensor features are obtained by constructing a Gaussian-smoothed structural tensor using the Sobel operator to calculate gradients, and then calculating the determinant, trace, and coherence to distinguish between edge and uniform regions. K-means coarse clustering (unsupervised) algorithm is used, setting the number of clusters K=3 to quickly classify multi-dimensional features into background, interference, and crack categories. The core is to minimize the sum of distances from samples within a cluster to the cluster center to achieve initial separation of target and non-target regions. DBSCAN refines the density-based clustering algorithm and optimizes the K-means results. Statistical feature selection calculates the intensity statistical features of the mean and variance of each cluster. By setting thresholds or statistical tests, clusters that meet the crack characteristics are identified. Finally, crack clusters are automatically identified based on the intensity statistical features, grouping pixels belonging to the same crack into one cluster to obtain the identification result.

[0007] Step 116: The U-Net network extracts multi-scale features by progressively downsampling the encoder on the left and progressively upsampling the decoder on the right to restore spatial resolution. It then fuses the high-resolution detail features of the encoder with the semantic features of the decoder using skip connections. U2Net is an improved model based on UNet, adding multiple branches and attention mechanisms to enhance segmentation performance. The U2-Net architecture is based on the RSU module, which is used in every En_x. The RSU module combines residual blocks and the U-Net symmetric encoder-decoder structure, preserving feature information through skip connections. Further improvements are made by introducing the ResNet residual module for multi-scale feature fusion. For the encoder part, there are a total of 4 downsampling modules. Each module contains 2 ResNet residual blocks and 1 2×2 max pooling. The ResNet residual block takes the input feature map as a 3×3 convolution, a batch normalization (BN) layer, a ReLU activation, a 3×3 convolution, and a BN layer. The shortcut connection directly superimposes the input feature map onto the output of the second BN layer. Its formula is F(x)=Conv2(ReLU(BN(Conv1(x)))) + x. The max pooling has a stride of 2 and doubles the number of channels to gradually expand the receptive field and capture global features of the rift. For the decoder part, there are a total of 4 upsampling modules. Each module contains a 2×2 transposed convolution and a feature concatenation layer. The transposed convolution restores the size of the feature map by padding with zeros and convolution operation. Its formula is Output = ConvTranspose2. The feature concatenation module concatenates the current feature map of the decoder with the feature map of the same scale of the encoder along the channel dimension. The number of channels is unified to 256 by 1×1 convolution, which preserves the high-resolution details of the encoder. For the multi-scale fusion layer, the sampling features of each layer of the fusion decoder are unified to 128 channels through 1×1 convolution, and the elements are added and then activated by ReLU: low-magnification downsampling features retain the details of small crack edges, while high-magnification downsampling features provide global context. For the output layer, a 1×1 convolution compresses 256-channel features into 1 channel, and a Sigmoid activation function outputs pixel-level probabilities. When the probability is >0.5, it is considered a crack. Then, model training and evaluation are performed. The cross-entropy loss function is selected. To address the class imbalance problem of low crack pixel proportion, a weight factor is introduced to ensure that the model pays attention to small cracks. The optimizer is set to Adam optimizer, which combines momentum gradient descent and adaptive learning rate. The initial learning rate is 0.001, which decays to 1 / 10 every 50 rounds, and the weight decays to 1e-5. The evaluation metric is Intersection over Union (IoU), which measures the degree of overlap between the predicted and real regions. The mean accuracy (mAP) is calculated within the IoU threshold of 0.5-0.95 to reflect the robustness of the model in recognizing cracks of different sizes. When the mean accuracy (mAP) is ≥92%, it is considered qualified. Finally, after multiple rounds of training to meet the requirements, a binary crack map is output.

[0008] Optionally, step 2 includes: By employing multi-sampling fusion and nonlinear optimization, parameterized fitting of sinusoidal cracks is achieved. First, a modified U-Net network is used to obtain a binary crack map, which is then preprocessed. Next, the coordinates are standardized, and multiple sampling methods are used to obtain the centerline. Finally, various parameters of the sinusoidal crack are derived. Based on the crack binary map, a technical process of coordinate standardization, geometric filtering, and multi-sampling fitting is constructed: a unified dimension is achieved through pixel-to-millimeter mapping, and non-sinusoidal interference is eliminated using the minimum bounding circle screening method. A weighted fusion sampling strategy is designed, incorporating equidistant sampling, adaptive sampling, and curvature keypoint sampling. This strategy, combined with the Levenberg-Marquardt algorithm and iterative reweighted least squares (IRLS), achieves the fitting of the sinusoidal model. Step 21: Crack preprocessing. The dataset is segmented using an improved U-Net network, and the segmentation results are preprocessed. The image coordinates are converted to physical dimensions through pixel-to-millimeter mapping. The crack skeleton is sampled at equal intervals and the mean value of the same horizontal coordinate is deduplicated. Step 22: Multi-mode sampling, equidistant sampling, taking points with a step size of 5 pixels along the center line x-direction to ensure uniform distribution; adaptive sampling, calculating curvature k, when k>0.1mm -1 At that time, the step size is halved to 2.5 pixels; curvature key points are sampled, the second derivative is calculated, and points with absolute values ​​> 70% quantile are retained; finally, weighted fusion is performed, and then duplicates are removed with the same horizontal coordinate to obtain the center line of the crack; Step 23: Parameter fitting. The sampled points are deduplicated and the mean is taken. Initial parameters for the sinusoidal model y = Rsin(2πx / P + β) + C are globally searched using the differential evolution algorithm. These are then refined using the weighted total least squares method, outputting the amplitude R, period P, phase β, and centerline position C. The parameters are then determined using the coefficient of determination R0. 2 Verification with root mean square error (RMSE).

[0009] Optionally, it includes: Step 211: Coordinate standardization. First, under the constraints of uniform aperture and unfolded aperture, calculate the pixel-to-millimeter conversion factor. Second, define the coordinate system, specifying the horizontal direction as circumferential length and the vertical direction as axial depth. Simultaneously, record the image size (W, H), resolution, period P, and retain the initial depth y0 of each image. Finally, apply a uniform dimension to sampling, fitting, and clustering to ensure that the parameters of different images have the same physical units. Step 212: Geometric Filtering and Coordinate Mapping. First, convert the image pixel coordinates to real physical millimeter coordinates. Set a length threshold Lth, calculate the real millimeter length L of each candidate line segment, retain line segments with L≥Lth, and remove noise or minor impurity line segments with L<Lth. Geometric Shape Constraints. Through two sub-constraints, filter line segments that conform to the sinusoidal crack shape. Slenderness Ratio Threshold Filtering. Set a slenderness ratio threshold Rth, and retain line segments with a length-to-width ratio Rslenderness≥Rth. Continuity Constraints. Calculate the line segment direction angle θ. For adjacent line segments with an endpoint distance <2mm, calculate the angle difference. θ, setting the angle difference threshold θ, to eliminate noise that leads to random directions; Step 213: Parameter estimation, estimation of the centerline position C. The upper and lower halves of a sine curve are symmetrical within a complete period, and their mean is approximately equal to the centerline position C. Let the dataset be... Then the mean The calculation formula is: The initial estimate of the center location is: ; The amplitude R is estimated, and the range of the sine curve is approximately equal to 2R. (Data set range) The calculation formula is: The initial estimate of the amplitude is: ; The period P is estimated, describing the length of one complete oscillation of the sine curve, while the frequency f is the number of oscillations per unit length. The two satisfy the reciprocal relationship: P = 1 / f. The Fast Fourier Transform (FFT) converts the time-domain data into the frequency-domain data, thereby extracting the periodic components of the data. The steps are as follows: First, to focus on the volatile aspects of the data, the data is processed to remove the mean: ; Secondly, the data after removing the mean... Perform an FFT to obtain a complex array in the frequency domain: ; in, Indicates frequency as The amplitude and phase of the components; frequency The relationship with array index k is as follows: ; in, Let x be the sampling frequency, i.e., the number of sampling points per unit length. If x is an equidistant sampling, and the interval between adjacent x is... ,but ; Next, we find the dominant frequency; in the frequency domain, the component with the largest amplitude corresponds to the dominant periodicity of the data. We then calculate the periodicity of each frequency. Corresponding amplitude The frequency corresponding to the maximum amplitude is the dominant frequency. ; Finally, based on the reciprocal relationship between period and frequency, the estimated value of the period is: ; Step 214: Equal-interval sampling, generating sampling points uniformly within a specified range; Adaptive sampling, dynamically adjusting the sampling density according to the curvature of the crack; Curvature key point sampling, selecting the key points with the most drastic morphological changes for sampling by calculating the second derivative; Step 215: Deduplication and Mean Taking. Remove duplicate points at the same x-coordinate and take the mean of the remaining y-values ​​to obtain the unique y-value corresponding to the current x-coordinate. Global Search and Local Refinement: First Stage: Differential Evolution (DEE), i.e., global search to find an approximate solution. As a population-based global optimization algorithm, it iteratively searches through three steps: mutation, crossover, and selection. Mutation: For each individual xi in the population, a mutation vector vi is generated using the difference vectors of other random individuals. Crossover: The mutation vector vi is crossed with the original individual xi to generate an experimental vector ui. Selection: The residuals of ui and xi are compared, and individuals with smaller residuals are retained for the next generation, ultimately obtaining globally approximate optimal candidate parameters. Second Stage: Weighted Overall Least Squares (WLS), i.e., local refinement to improve accuracy. Based on the global search for an approximate solution, further optimization is performed on the inlier set. The core is weighting, assigning weights to different data points. The goal is to minimize the weighted sum of squared residuals. The optimal parameters are solved by taking the partial derivative of the objective function, ultimately obtaining a high-precision parameter estimation result. After crack acquisition and preprocessing, initial parameter estimation, sampling and weighted fusion of three sampling methods, and deduplication of sampling points to obtain the centerline, a global search and local refinement are performed. The results are then evaluated and visualized, and the goodness-of-fit index, the coefficient of determination R, is calculated. 2 The root mean square error (RMSE) is used to assess the quality of the fit. For cracks with poor fit quality, the clustering and centerline extraction processes are re-examined and adjusted. The fitted sine curves are visualized and compared with the original cracks to visually verify the effect. The parameter estimates of each crack are summarized to form a quantitative analysis results table.

[0010] Optionally, step 3 includes: The macroscopic orientation trend and microscopic roughness residuals of cracks are separated using smooth spline functions. The effects of three sampling methods—equal-interval sampling, equal-arc-length sampling, and multi-scale adaptive sampling—are compared. The influence of roughness and JRC value is investigated; various geometric features of the crack are extracted, decoupled, and multi-sampling is used to calculate the JRC value based on the sampling points, thus obtaining the roughness; Step 31: Contour extraction and trend separation. First, cracks are extracted and preprocessed using an improved U-Net network. Then, trend separation and residual extraction are performed. After reparameterizing the crack skeleton by arc length, a smooth spline or low-order polynomial is used to fit the trend line. The high-frequency rough residual is obtained by the difference between the original curve and the trend line. Step 32: Calculate roughness indices using multiple sampling methods, compare the differences between equal-interval sampling, equal-arc-length sampling, and multi-scale adaptive sampling, and employ denser sampling in high-curvature regions. Calculate the length-normalized first-order slope energy index based on the residual data. ; Step 33: JRC Inversion and Uncertainty Assessment. Introduce endpoint missing measurement masks to complete JRC mapping and sensitivity analysis. Apply empirical formulas to... The values ​​are converted into JRC exponents, and the Bootstrap method is used to quantify the impact of sampling interval, window size and detrending strength on the results, outputting the JRC mean and 95% confidence interval.

[0011] Optionally, it includes: Step 311: Feature decoupling separates the low-frequency shape and high-frequency roughness features of the crack using the smooth spline method. The core is to minimize the joint objective function of data fitting error and curve smoothness to obtain a smooth curve that reflects the overall trend. Then, the smooth curve is subtracted from the original data to obtain a residual sequence containing only high-frequency roughness information. Step 312: Select the appropriate sampling method for different scenarios and calculate the roughness index. Equal-distance sampling, sampling at equal intervals along the x-axis, focuses on lateral comparability, and calculates the distance between adjacent points based on the elevation difference and lateral step size. It is suitable for lateral feature comparison, and its expression is: ; Equal arc length sampling: Convert parameters to cumulative arc length. Under equal arc length sampling, the roughness index... The expression is: ; Multi-scale adaptive sampling: Combining multi-scale analysis with dynamic curvature adjustment of sampling density, dense sampling is performed in complex regions and sparse sampling is performed in smooth regions to balance feature capture accuracy and computational efficiency; the set of sampling points obtained by multi-scale adaptive sampling is represented as: ; Curvature-weighted optimization: Assigning high weights to regions with high curvature to correct... Indicators that enhance sensitivity to irregular features of fracture peaks; Step 313, JRC inversion utilization The mapping relationship with JRC's engineering experience will be used to map different sampling methods and scales. The value is converted into a dimensionless JRC value to achieve a quantitative conversion from geometric features to roughness coefficients; in engineering, empirical formulas are used to convert... Mapped to JRC, using a logarithmic calibration expression: ; Where JRC(w) represents the joint roughness coefficient on window w, which is a dimensionless empirical level; The calibration data is related to specimen dimensions, test line length, and unit consistency. This represents the roughness index value on window w; Step 314: Bootstrap uncertainty assessment calculates the standard error and confidence interval of the JRC by repeatedly resampling and propagating the discrete point selection, taking into account the influence of noise factors, quantifying the reliability of the results, and finally outputting the JRC value of the crack with confidence. The formula is as follows: ; in, This represents the data from the b-th resampled dataset. Indicates detrending, The computation chain of JRC, where B represents the number of resampling steps. This represents the mean obtained from multiple resampling operations, and SE represents the standard error. This represents the empirical quantile of the 1-α confidence interval.

[0012] Optionally, step 4 includes: Establish a technical framework for coordinate mapping, connectivity assessment, uncertainty optimization, 3D reconstruction and supplementary point generation: Based on the borehole opening or bottom coordinates, the two-dimensional fracture parameters are inverted into fracture patches in three-dimensional space through Lie algebra SE(3) pose transformation. The normal angle, spatial distance, fracture size and JRC difference are integrated to construct a multi-dimensional consistency index, which is mapped to the connectivity probability in the 0-1 interval by the sigmoid function. The breadth-first search algorithm is used to extract fracture clusters, and the representative plane of the fractures in the cluster is fitted by the least squares algorithm. By calculating the fracture density variance and connectivity probability entropy of each point in space, the high uncertainty region is located. Step 41: 3D coordinate mapping and fracture extraction. The borehole opening is taken as the origin (x=0, y=0, z=0), where x represents the parallel direction of the roadway, y represents the perpendicular direction of the roadway wall, and z represents the depth. A homogeneous transformation matrix T_BH is constructed through the borehole pose to map the local fracture parameters to the global coordinate system. Based on the borehole opening or bottom coordinates, the single-hole fracture parameters are transformed to the global coordinate system. The 3D plane corresponding to the fracture is inverted through the SE(3) pose transformation and clipped into finite patches according to the observation length. For complex fractures, Fourier series or spline surface extension models are used. The fractures are extracted by improving the U-Net network and fracture preprocessing is performed to filter out interference information in the fractures. Step 42: Connectivity analysis. Construct a multi-dimensional consistency index. At the geometric level, features such as normal angle, plane constant difference, and Hasdorf distance are used. At the texture level, JRC difference is introduced, and weights are designed to weight the two features. The comprehensive score is mapped to the connectivity probability in the 0-1 interval through logistic regression. A weighted undirected graph is constructed and connected gap clusters are extracted. Within the cluster, the total least squares joint fitting is used to fit the representative plane. Step 43: Uncertainty and Hole Filling Optimization. By calculating the crack density variance and connectivity probability entropy of each point in space, high uncertainty regions are located, and candidate hole filling regions are selected in combination with engineering constraints. Hole filling position optimization involves traversing spatial coordinates and using the information gain maximization algorithm to evaluate the uncertainty reduction benefits of each candidate point. The top three optimal positions are selected, and finally, the crack network, connectivity probability, and hole filling scheme are presented through three-dimensional visualization, forming a closed loop of identification-modeling-decision.

[0013] Optionally, it includes: Step 411: Data Integration and Coordinate Unification. Establish the rigid body transformation relationship between the borehole local coordinate system and the global coordinate system to unify different borehole data into the same spatial system and ensure consistency in cross-hole analysis; at the same time, integrate the core data of sinusoidal model parameters, JRC, and borehole radius to provide a unified quantitative basis for 3D reconstruction. Step 412: 2D to 3D Inversion and Patch Construction. Based on the sinusoidal model parameters of the borehole unfolding diagram, the normal vector and plane equation of the fracture plane are inversely derived and mapped to the global coordinate system through pose transformation. Combined with the effective observation length of the borehole, a finite fracture patch is cut out and the metadata of the center and scale is recorded to complete the transformation from 2D measurement to 3D morphology. Step 413: Cross-hole connectivity determination. Using geometric and texture multi-dimensional indicators, the connectivity probability of different borehole fractures is quantified. Geometric consistency: Calculate the normal angle, planar position difference, and spatial distance, and combine attitude similarity and spatial coplanarity to determine the fit between the fracture's spatial orientation and position. Texture corroboration: Introduce JRC difference to supplement roughness consistency. Perform probability mapping, obtain a comprehensive score through weighted fusion of multi-dimensional indicators, and then convert it into connectivity probability using the Sigmoid function to form a symmetric probability matrix. Obtain a comprehensive similarity score by weighted fusion of geometric and roughness indicators. Comprehensive similarity... The formula is: ; in, As weight, This is a scaling constant used to normalize the various indicators. For roughness inversion of the corresponding intersecting window of the facet; Logistic regression mapping of connectivity probability converts the comprehensive similarity score into connectivity probability, facilitating threshold determination and graph model calculation; the sigmoid function is used to achieve a monotonic mapping from score to probability, obtaining the probability that two fractures belong to the same continuous fracture volume. ; in, α represents the probability that patches i and j belong to the same continuous fracture body, and b are calibration coefficients determined through training data. Step 414: Probabilistic Graph Construction and Cluster Extraction. Using fracture patches as nodes and connectivity probabilities as edge weights, a weighted undirected graph is constructed. A probability threshold is set to filter valid edges. Interrelated fracture clusters are obtained through connected component extraction or community detection algorithms. The overall plane or surface of the cluster is then fitted using the global least squares method to characterize its geometric features. The connectivity relationships of the fracture patches are transformed into a graph structure. Using connected component extraction or community detection algorithms in graph theory, the set of closely related fracture patches, i.e., fracture clusters, is identified. ; Where G is the crack diagram. For a set of nodes, Let it be the set of edges; Step 415: Uncertainty Propagation and Aggregation. The uncertainty of a single parameter is propagated to the position and orientation parameters of the 3D patch through the Jacobian matrix. At the same time, the covariance of parameters within the cluster is aggregated to quantify the uncertainty of the fractured mass cluster and the connectivity probability. Step 416: Hole Filling Optimization. Under the engineering reachability domain and construction constraints, achieve efficient hole filling selection; define the constraint range, define the reachability domain in combination with construction conditions, and generate a set of candidate holes that are feasible for the project based on the existing borehole distribution; lock high-value regions, integrate cluster plane uncertainty and connectivity probability, and identify high-uncertainty-high-value hole filling hotspots; select the optimization objective, using the expected information gain of candidate holes as the objective function, and screen the optimal hole filling combination.

[0014] The technical solution provided by this invention includes a method that involves collecting and preprocessing a dataset, identifying fractures based on an improved U-Net network, and outputting a binary fracture map; performing quantitative analysis and modeling of sinusoidal fractures based on the binary fracture map; calculating the roughness coefficient JRC and performing sensitivity analysis of the sampling method based on the model; and performing three-dimensional reconstruction and borehole optimization of the multi-hole fracture network based on the JRC calculation and sensitivity analysis. This method improves the accuracy and reliability of fracture identification in surrounding rock and enhances its engineering applicability. Attached Figure Description

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

[0016] Figure 1 A flowchart of the method for identifying surrounding rock fractures and reconstructing a three-dimensional model provided in an embodiment of the present invention; Figure 2 This is a diagram of the improved U-Net network structure provided in an embodiment of the present invention; Figure 3 This is a detection effect diagram of the traditional modeling method provided in the embodiments of the present invention; Figure 4 The diagram shows the detection effect of the adaptive threshold method provided in the embodiment of the present invention. Figure 5 The detection effect diagram of the machine learning method provided in the embodiment of the present invention; Figure 6 This is a diagram showing the detection effect of a traditional deep learning method provided in an embodiment of the present invention. Figure 7 The images provided in this embodiment of the invention show the detection results of YOLOv11, U2Net, and the improved version Unet, where (a) is the detection result of YOLOv11; (b) is the detection result of U2Net; and (c) is the detection result of the improved version Unet. Figure 8 The following are multiple sets of visualization effect diagrams for crack detection provided in the embodiments of the present invention, wherein (a) is an image feature enhancement effect diagram; (b) is a visualization binary image of crack detection; and (c) is an overlay diagram of the detection effect and the original image. Figure 9a The following is a visualization of the fitted result provided in the embodiments of the present invention, wherein (a) is the original image; (b) is the fracture clustering result; (c) is the fracture skeleton extraction result; and (d) is the superimposed image of the original image and the fitted curve. Figure 9b Another visualization of the fitted result provided in this embodiment of the invention, wherein (a) is the original image; (b) is the fracture clustering result; (c) is the fracture skeleton extraction result; and (d) is the superimposed image of the original image and the fitted curve. Figure 9c Another visualization of the fitted result provided in this embodiment of the invention, wherein (a) is the original image; (b) is the fracture clustering result; (c) is the fracture skeleton extraction result; and (d) is the superimposed image of the original image and the fitted curve. Figure 10The images provided in this embodiment of the invention are a crack detection and contour extraction effect diagram, wherein (a) is a crack detection effect diagram; and (b) is a contour extraction effect diagram. Figure 11a This is a visualization of the equally spaced sampling results provided in an embodiment of the present invention; Figure 11b This is a visualization of the equal arc length sampling results provided in an embodiment of the present invention; Figure 11c A visualization of multi-scale sampling results provided in an embodiment of the present invention; Figure 12a A segmented JRC analysis diagram provided for an embodiment of the present invention; Figure 12b This is a JRC comparison chart showing different sampling effects provided in embodiments of the present invention; Figure 13 A comparison chart of JRC values ​​for different sampling methods provided in embodiments of the present invention; Figure 14 This is a rendering of the 3D connectivity reconstruction provided in an embodiment of the present invention; Figure 15 Supplementary borehole location diagram provided for embodiments of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] It should be understood that the described embodiments are merely some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0019] The terminology used in the embodiments of this invention is for the purpose of describing particular embodiments only and is not intended to limit the invention. The singular forms “a,” “the,” and “the” used in the embodiments of this invention are also intended to include the plural forms unless the context clearly indicates otherwise.

[0020] It should be understood that the term "and / or" used in this article 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.

[0021] Depending on the context, the word "if" as used here can be interpreted as "when," "when," "in response to determination," or "in response to detection." Similarly, depending on the context, the phrase "if determination" or "if detection (of the stated condition or event)" can be interpreted as "when determination," "in response to determination," "when detection (of the stated condition or event)," or "in response to detection (of the stated condition or event)."

[0022] Figure 1 The flowchart of the surrounding rock fracture identification and three-dimensional model reconstruction method provided in the embodiments of the present invention is as follows: Figure 1 As shown, the method includes: Step 1: Collect the dataset, preprocess it, and perform crack identification based on the improved U-Net network to output a binary crack map.

[0023] To address the issue that interference factors such as natural rock texture, borehole wall mud contamination, and drill bit marks in borehole imaging unfolded images are similar to the visual features of fractures, a three-level technical comparison system of traditional methods, machine learning, and deep learning is designed. Traditional methods initially extract fractures through a voting mechanism of multi-directional gradient detection, Canny edge detection, and LoG edge detection (at least two methods must agree). An adaptive thresholding method integrates Sobel edge features, LBP texture features, and Hessian linear structure features, combined with the Otsu automatic thresholding algorithm to optimize segmentation. The machine learning method is based on a random forest classifier and uses pseudo-label generation technology (constructing labeled data using the results of traditional methods) to solve the problem of scarce training data. The traditional deep learning method uses multi-scale feature extraction (Gaussian filtering, LoG response, gradient magnitude) and K-Means+DBSCAN clustering to achieve unsupervised recognition. Finally, an improved U-Net model is proposed, replacing some convolutional layers with residual modules to reduce gradient vanishing, integrating features from multiple decoding layers to improve fine-grained recognition capabilities, and combining image segmentation (244×244 pixels) and mosaic data augmentation strategies to improve the model's generalization ability. Comparative experiments show that the improved U-Net achieves an average accuracy (mAP) of 92.6%, which is 8 percentage points higher than U2Net. It can effectively suppress various interferences and completely preserve the details of the crack contour.

[0024] In this embodiment of the invention, step 1 includes: By comparing image preprocessing and multi-method fusion, and optimizing and training the network structure, pixel-level separation of cracks and interference in borehole images is achieved, providing a high-quality data foundation for subsequent quantitative analysis. This invention designs five methods to solve this problem, and finally, the method with the best effect is obtained through comparison. Step 11: Data acquisition and preprocessing. Obtain the dataset of borehole images and perform grayscale conversion, Gaussian denoising, and histogram equalization on the images. By random cropping, the images are randomly cropped into 244×244 pixel samples, followed by rotation and mosaic enhancement to expand the dataset samples. Step 12: Multi-method fusion and comparison, designing 5 methods: a) Employing multi-method edge feature extraction plus a voting mechanism; b) Employing multi-dimensional feature extraction and weighted fusion, combined with adaptive thresholding for crack detection; c) Employing machine learning algorithms, based on multi-dimensional feature extraction and random forest algorithms to classify crack pixels; d) Employing multi-scale feature extraction and KMeans clustering algorithms to cluster crack pixels, ultimately achieving crack detection; e) Employing deep learning algorithms for segmentation, using improved U-Net, U2Net, and YOLOv11 for training respectively; The improved U-Net network consists of an encoder (4 downsampling modules, each containing 2 ResNet residual blocks (3×3 convolution + BN + ReLU, with shortcut connections skipping 2 convolutional layers) + 1 2×2 max pooling (stride 2), with 64, 128, 256, and 512 channels), a decoder (4 upsampling modules, 2×2 transposed convolution (stride 2) + concatenation of features at the same scale as the encoder (1×1 convolution for unified channels)), a multi-scale fusion layer (enhancing small gap response), and an output layer (1×1 convolution + Sigmoid activation). Step 13: Model training and evaluation. The results are visualized, and the fracture extraction results are selected. Deep learning methods have more objective evaluation criteria. First, the dataset is divided into training, validation, and test sets in a 7:2:1 ratio. The dataset is evaluated by Intersection over Union (IoU) (IoU = predicted ∩ label / predicted ∪ label) and mean accuracy (mAP). The resulting binary fracture map is output (mAP ≥ 92%, IoU ≥ 0.85).

[0025] In this embodiment of the invention, it includes: Step 111: Grayscale Conversion. The borehole image is in RGB format, and color information offers no benefit for fracture identification (the color difference between the fracture and the background is not significant). A weighted average method, Y=0.299R+0.587G+0.114B, is used to convert it to a grayscale image. This preserves the brightness information sensitive to the human eye while reducing the amount of data by 2 / 3, improving computational efficiency. Gaussian Denoising. During borehole imaging, sensor noise and downhole light fluctuations cause high-frequency noise in the image (manifested as isolated bright spots / dark spots). A 5×5 Gaussian kernel (σ=1.5) is used for convolution. The Gaussian function smooths high-frequency noise while preserving... To address the issue of concentrated grayscale distribution in images due to uneven downhole illumination (e.g., localized over-darkness / over-brightness), the CLAHE algorithm, an adaptive histogram equalization algorithm with limited contrast, was employed. The image was divided into an 8×8 grid, and histogram equalization was performed on each grid individually, while limiting contrast enhancement (clipLimit=2.0) to avoid noise amplification. This improved the grayscale difference between the fracture and the background from 10-20 grayscale levels to 30-50 grayscale levels. To address the problem of insufficient sample size and uneven distribution in the original dataset, the dataset was expanded to four times the original size using the following methods: random cropping (images were randomly cropped into 244×244 pixel regions to simulate different imaging perspectives); random rotation (samples were randomly rotated 0-360° to address the diversity of fracture orientations); and mosaic enhancement (four different samples were stitched together in a 2×2 pattern with linear transitions at the edges to enhance the model's robustness to incomplete fracture areas). Step 112: After edge extraction using three edge extraction methods—multi-directional gradient Sobel edge detection, Canny edge detection, and LoG edge detection combined with Gaussian blur and Laplacian operator edge detection—a voting method is used. A crack is identified by combining the results of at least two detection methods. These three methods complement each other in terms of directional sensitivity, noise resistance and continuity, and adaptability to special morphologies: Sobel edge detection quickly captures multi-directional edge intensity, Canny edge detection accurately extracts continuous edges, and LoG edge detection supplements the detection of blob-like structures. By combining criteria (e.g., a crack is identified only if at least two methods agree), false detections (e.g., noise is misjudged) and false negatives (e.g., special-morphology cracks are not identified) caused by a single method can be significantly reduced, ultimately yielding a more reliable preliminary crack region.

[0026] Step 113: After extracting three features—Sobel edge detection, LBP texture detection, and Hessian—and combining them with an adaptive threshold, the weights of the three features are 50% for multi-scale, 30% for Sobel edge detection, and 20% for Hessian, forming multi-scale and multi-directional features. Finally, the Otsu automatic thresholding is used to automatically calculate the threshold with the minimum intra-class variance between the foreground and background to obtain the final detection result. Step 114: Supervised learning based on random forest sequentially involves feature extraction, pseudo-label generation, forest training, and confidence filtering. Each step is deeply adapted to the characteristics of random forest. Feature extraction is diverse, extracting multi-dimensional features such as original intensity, gradient magnitude, LBP texture, Canny edges, and Hessian matrix features to provide sufficient material for random forest features. Randomly selected features are used to train each tree to explore the discriminative power of different feature combinations. Pseudo-label generation uses mathematical and statistical methods to generate pseudo-labels, transforming unsupervised or weakly supervised data into labeled data, providing a training foundation for supervised learning of random forest (random forests inherently rely on labeled data). Forest training parameter design, n_estimator... `s=50` indicates that the number of trees is 50, used to control the scale of the ensemble. More trees usually improve performance but increase computation. `max_depth=8` indicates that the maximum depth of the trees is 8, used to limit the complexity of a single tree and balance bias and variance (too large a depth can lead to overfitting, while too small a depth can lead to underfitting). Confidence assessment based on gradient and standard deviation filters the results of the random forest predictions. The gradient reflects the stability of the predictions, and the standard deviation reflects the divergence between predictions between trees. Filtering reduces the impact of incorrect predictions. Finally, the classifier learns the mapping relationship between features and fissure or non-fissure labels to classify new data. The confidence of the classification results is evaluated based on gradient and standard deviation, filtering out predictions with low confidence and retaining those with high confidence to obtain the final recognition result.

[0027] Step 115: The core of traditional deep learning methods is multi-scale feature extraction + clustering optimization. It achieves gap recognition through four steps: full feature coverage, coarse classification, fine-tuning, and target selection. Multi-scale feature extraction uses multiple window sizes (5×5, 7×7, 9×9, 11×11) and structural tensors to obtain comprehensive features; basic statistical features are calculated by measuring the mean, variance, standard deviation, and contrast of each window to capture local gray-level distribution patterns; structural tensor features are constructed by using the Sobel operator to calculate gradients, constructing a Gaussian-smoothed structural tensor, calculating determinant, trace, and coherence to distinguish between edges (high coherence) and uniform regions (low coherence); the K-means coarse clustering unsupervised clustering algorithm, with K=3 clusters, quickly divides the multi-dimensional features into... The system categorizes data into three types: background, interference, and cracks. The core principle is to minimize the sum of distances from samples within a cluster to the cluster center to initially separate the target from non-target data. DBSCAN refines the density-based clustering algorithm (parameters eps=0.5, min_samples=5) and optimizes the K-means results. Statistical feature filtering calculates the intensity statistical features of the mean and variance of each cluster. By setting thresholds or statistical tests, clusters that meet the crack characteristics are identified. Finally, based on the intensity statistical features, crack clusters are automatically identified, and pixels belonging to the same crack are grouped into one cluster to obtain the identification result.

[0028] Step 116: The U-Net network is a fully convolutional neural network with an encoder-decoder structure, initially designed for biomedical image segmentation. The U-Net network extracts multi-scale features through progressive downsampling by the encoder on the left and progressive upsampling by the decoder on the right to restore spatial resolution. It also uses skip connections to fuse the high-resolution detail features of the encoder with the semantic features of the decoder. U2Net is an improved model based on UNet, adding multiple branches and attention mechanisms to enhance segmentation performance. The U2-Net architecture is based on the RSU (ReSidual U-block) module. U2-Net uses an RSU module in each En_x. The RSU module combines residual blocks and the U-Net symmetric encoder-decoder structure, preserving feature information through skip connections and avoiding the gradient vanishing problem in deep networks. The RSU module is mainly used to improve edge detection capabilities in image segmentation tasks, learning multi-scale features without significantly increasing computational cost. However, the traditional U-Net network is prone to gradient vanishing during deep training and has a weak response to small gaps (width < 2mm). Therefore, this invention improves upon U-Net by introducing a ResNet residual module and multi-scale feature fusion, such as... Figure 2 As shown: For the encoder part, there are a total of 4 downsampling modules. Each module contains 2 ResNet residual blocks and 1 2×2 max pooling (stride 2). The ResNet residual block takes the input feature map as a 3×3 convolution (the number of convolution kernels is the same as the number of input channels), a batch normalization (BN) layer (to accelerate training convergence and prevent overfitting), ReLU activation, 3×3 convolution, and BN layer. The shortcut connection directly superimposes the input feature map to the output of the second BN layer. Its formula is F(x)=Conv2(ReLU(BN(Conv1(x)))) + x, which solves the gradient vanishing problem through residual propagation. The max pooling has a stride of 2 and doubles the number of channels (64, 128, 256, 512), gradually expanding the receptive field (from 3×3 to 24×24) to capture global features of the gap. For the decoder section, there are four upsampling modules. Each module contains a 2×2 transposed convolution (stride 2, channel count halved) and a feature concatenation layer. The transposed convolution restores the feature map size (e.g., from 512 channels to 256 channels, feature map size × 2) through zero padding and convolution operations. Its output formula is Output = ConvTranspose2(Input, kernel=2×2, stride=2, padding=0), avoiding feature blurring caused by traditional interpolation. The feature concatenation module concatenates the current feature map of the decoder with the feature map of the same scale of the encoder along the channel dimension (e.g., 256-channel decoder features + 256-channel encoder features equals 512 channels). The number of channels is unified to 256 through 1×1 convolution, preserving the high-resolution details of the encoder (e.g., crack edges). For the multi-scale fusion layer, the sampled features of each layer of the fusion decoder are unified to 128 channels through 1×1 convolution, and the elements are summed and then activated by ReLU. This layer can enhance the feature response of small cracks: low-magnification downsampled features preserve the edge details of small cracks, while high-magnification downsampled features provide global context. This is equivalent to combining the advantages of convolution and transform. After fusion, the feature values ​​of small cracks (width <2mm) are improved by 30%-50%. For the output layer, a 1×1 convolution compresses 256-channel features into 1 channel. The Sigmoid activation function outputs pixel-level probabilities (0-1), and a probability > 0.5 is considered a crack. Then, model training and evaluation are performed. The cross-entropy loss function is chosen. To address the class imbalance problem of low crack pixel proportion (typically <10%), a weight factor is introduced to ensure the model's attention to small cracks. The optimizer is set to the Adam optimizer, combining momentum gradient descent and adaptive learning rate, with an initial learning rate of 0.001. Balancing convergence speed and accuracy, the weights are reduced to 1 / 10 every 50 rounds (for fine-tuning in later stages), with a weight reduction of 1e-5 (to suppress overfitting). The evaluation metric is the Intersection over Union (IoU), which measures the degree of overlap between the predicted and real regions. The mean accuracy (mAP) is calculated within an IoU threshold of 0.5-0.95 (step size 0.05) to reflect the model's robustness in identifying cracks of different sizes. A mean accuracy (mAP) ≥ 92% is considered acceptable. Finally, after multiple rounds of training to meet the metric requirements, a binary crack map is output.

[0029] Step 2: Based on the binary fracture map, perform quantitative analysis and modeling of sinusoidal fractures.

[0030] In this embodiment of the invention, step 2 includes: By employing multi-sampling fusion and nonlinear optimization, parameterized fitting of sinusoidal fractures is achieved. First, an improved U-Net network is used to obtain a binary fracture map, which is then preprocessed. Next, the coordinates are standardized, and multiple sampling methods are used to obtain the centerline. Finally, various parameters of the sinusoidal fracture are derived. Based on the fracture binary map, a technical process of coordinate standardization, geometric filtering, and multi-sampling fitting is constructed: Dimensions are unified through pixel-to-millimeter mapping (with borehole circumference as the circumferential reference), and the minimum bounding circle screening method is used (the ratio of fracture area to minimum bounding circle area > 0). 35) Remove non-sinusoidal interference (such as irregular drilling marks); design a weighted fusion sampling strategy that includes equal-interval sampling (weight 0.3 to ensure data uniformity), adaptive sampling (weight 0.5 for denser sampling in high-curvature regions), and curvature key point sampling (weight 0.2 to focus on points with drastic morphological changes). Combine this with the Levenberg-Marquardt algorithm and iterative reweighted least squares IRLS to fit the sinusoidal model (amplitude R, period P, phase β, centerline position C); analysis of the data in the dataset shows that the coefficient of determination (R²) is high. 2 The mean is greater than 0.95.

[0031] Step 21: Crack preprocessing. The improved U-Net network is used to segment the dataset. The segmentation results are preprocessed to filter out most irrelevant factors. The image coordinates are converted to physical dimensions through pixel-millimeter mapping. The crack skeleton is sampled at equal intervals and the mean value of the same horizontal coordinate is deduplicated to reduce redundant point interference. Step 22: Multi-mode sampling, equidistant sampling, taking points with a step size of 5 pixels along the center line x-direction to ensure uniform distribution; adaptive sampling, calculating curvature k, when k>0.1mm -1 At that time, the step size is halved to 2.5 pixels; curvature key points are sampled, the second derivative is calculated, and points with absolute values ​​> 70% quantile are retained; finally, weighted fusion is performed, and then duplicates are removed with the same horizontal coordinate to obtain the center line of the crack; Step 23: Parameter fitting. The sampled points are deduplicated and the mean is taken. A differential evolution algorithm (population 50, 100 iterations, mutation factor F=0.5, crossover probability CR=0.7) is used to globally search for the initial parameters of the sinusoidal model y=Rsin(2πx / P+β)+C (R∈[1,50]mm, P∈[50,200]mm, β∈[-π,π], C∈[0,1200]mm). This is then refined using a weighted overall least squares method, outputting the amplitude R, period P, phase β, and centerline position C. The parameters are then determined using the coefficient of determination R0. 2 Verification with root mean square error (RMSE).

[0032] In this embodiment of the invention, it includes: Step 211: Coordinate Standardization. The original unit of the borehole image is pixels, while engineering analysis (such as 3D inversion and cross-hole comparison) requires actual physical units (millimeters). Without a unified unit of measurement, 100 pixels in different images may correspond to different actual lengths, rendering the parameters meaningless. First, under the constraints of a unified borehole diameter and unfolded aperture, calculate the pixel-to-millimeters conversion factor. For example, if 244 pixels in the circumferential direction of the image corresponds to an actual borehole circumference of 94.2 mm, then the conversion factor is... mm / pixel; secondly, coordinate conventions are defined, specifying the horizontal direction as the circumferential length (corresponding to the circumferential direction after the cylindrical hole wall is unfolded) and the vertical direction as the axial depth. At the same time, the image size (W, H), resolution (dpi), period P, and the initial depth y0 of each image are recorded; finally, a unified dimension is applied to sampling, fitting, and clustering to ensure that the parameters (such as sinusoidal amplitude and period) of different images have the same physical units. Step 212, Geometric Filtering (Preprocessing): Geometric filtering is a key screening step for sinusoidal crack detection. Its core is millimeter coordinate mapping + dual rule constraints, which removes noise from candidate line segments and retains high-quality line segments that meet the characteristics of large slenderness ratio and continuous direction, paving the way for subsequent sinusoidal fitting. The coordinate mapping process first converts the image pixel coordinates into real physical millimeter coordinates, ensuring that the calculation of geometric parameters such as length and angle conforms to the actual physical scale, providing accurate data support for subsequent screening. Length threshold filtering removes excessively short line segments that do not conform to physical characteristics, avoiding noise interference. A length threshold Lth is set (usually 5mm~10mm, which can be lowered to 3mm in small crack scenarios). The actual millimeter length L of each candidate line segment is calculated, retaining line segments with L≥Lth and removing noisy or minor impurity line segments with L<Lth. Geometric shape constraints use two sub-constraints to screen line segments that conform to the sinusoidal crack shape. Slenderness ratio threshold filtering sets a slenderness ratio threshold Rth (usually ≥5), retaining line segments with a length-to-width ratio Rslenderness≥Rth and removing short, thick impurity, and scratched line segments. directional continuity constraints calculate the line segment directional angle θ (the angle with the positive x-axis). For adjacent line segments with an endpoint distance <2mm, the angle difference is calculated. θ (for angles greater than 90°, supplementary angle correction is used), setting an angle difference threshold. θ (10°~25°) is used to eliminate noise with chaotic directions. Step 213: Parameter estimation. Nonlinear optimization is sensitive to initial values. If the initial values ​​deviate too far from the true solution, it is easy to get trapped in local optima. Obtaining initial estimates close to the true values ​​through statistics (mean, range) and frequency domain analysis (FFT) can significantly improve fitting efficiency and accuracy. Estimation of the centerline position C: Within a complete period, the areas of the upper and lower halves of the sine curve are symmetrical, and their mean is approximately equal to the centerline position C; let the dataset be... Then the mean The calculation formula is: The initial estimate of the center location is: ; The amplitude R is estimated, and the range of the sine curve (the difference between the maximum and minimum values) is approximately equal to 2R (because the amplitude is the distance from the center line to the crest / trough, and the total distance from the crest to the trough is 2R). (Data set range) The calculation formula is: The initial estimate of the amplitude is: ; The period P is estimated, describing the length of one complete oscillation of the sine curve, while the frequency f is the number of oscillations per unit length. The two satisfy the reciprocal relationship: P = 1 / f. The Fast Fourier Transform (FFT) converts the time-domain data (y-values ​​varying with x) into the frequency-domain data, thereby extracting the periodic components of the data. The steps are as follows: First, to focus on the volatile parts of the data (eliminating the influence of the center line), the data is processed to remove the mean: ; Secondly, the data after removing the mean... Perform an FFT to obtain a complex array in the frequency domain: ; in, Indicates frequency as The amplitude and phase of the components; frequency The relationship with array index k is as follows: ; in, Let x be the sampling frequency, i.e., the number of sampling points per unit length. If x is an equidistant sampling, and the interval between adjacent x is... ,but ; Next, we find the dominant frequency; in the frequency domain, the component with the largest amplitude corresponds to the dominant periodicity of the data. We then calculate the periodicity of each frequency. Corresponding amplitude The frequency corresponding to the maximum amplitude is the dominant frequency. ; Finally, based on the reciprocal relationship between period and frequency, the estimated value of the period is: ; Step 214: Equal-interval sampling generates sampling points uniformly within a specified range. This is suitable for crack scenarios with uniform morphology and low noise, ensuring the uniformity of data distribution and often serving as a basic reference for initial analysis or comparative experiments. Adaptive sampling dynamically adjusts the sampling density based on the crack's curvature, sampling densely in areas with high curvature (complex morphology, obvious transitions) and sparsely in areas with low curvature (straight). This is suitable for scenarios with complex crack morphology, accurately restoring key features while reducing redundant computation. Curvature key point sampling uses the second derivative to select key points with the most dramatic morphological changes for sampling. This is suitable for scenarios with limited computational resources but requiring accurate capture of key crack morphology, capturing core features with fewer samples. Step 215: Deduplication and Mean Taking. When fitting a sinusoidal model to the data, the same x-coordinate may correspond to multiple y-values. This is usually caused by imaging noise (such as jitter at the edge of a crack due to imaging problems) or discrete noise. The deduplication and mean taking operation is designed to address this situation. It removes duplicate points under the same x-coordinate and takes the average of the remaining multiple y-values ​​to obtain the unique y-value corresponding to the current x-coordinate. This can eliminate the influence of local thick boundaries (jatter at the edge of a crack due to imaging noise) and discrete noise, making the data smoother and more suitable for fitting a sinusoidal model. Global Search and Local Refinement. The sinusoidal model is a nonlinear model containing sinusoidal terms. Directly using ordinary least squares can easily lead to local optima (missing the global optimal solution). However, differential evolution has a strong global search capability, which can traverse the feasible region of parameters to find an approximate optimal solution. Then, it is handed over to weighted least squares for local refinement, which can ensure the global optimality of the solution and improve the accuracy of parameter estimation. The first stage is differential evolution, which is a global search to find an approximate solution. As a population-based global optimization algorithm, it iteratively searches through three steps: mutation, crossover, and selection. Mutation involves generating a mutation vector vi for each individual xi in the population using the difference vectors of other random individuals. Crossover involves crossing the mutation vector vi with the original individual xi to generate an experimental vector ui. Selection involves comparing the residuals of ui and xi, retaining individuals with smaller residuals for the next generation, and finally obtaining globally approximate optimal candidate parameters. The second stage is weighted total least squares, which is a local refinement to improve accuracy. Based on the global search for an approximate solution, it further optimizes the inlier set (data points with good fitting effect and low noise). The core is weighting, assigning weights to different data points (higher weights for points with small residuals and lower weights for noisy points) to reduce the interference of noise on parameter estimation. The goal is to minimize the weighted sum of squared residuals and solve for the optimal parameters by taking the partial derivative of the objective function, ultimately obtaining high-precision parameter estimation results. After crack acquisition and preprocessing, initial parameter estimation, sampling and weighted fusion of three sampling methods, and deduplication of sampling points to obtain the centerline, a global search and local refinement are performed. The results are then evaluated and visualized, and the goodness-of-fit index, the coefficient of determination R, is calculated. 2The root mean square error (RMSE) is used to assess the quality of the fit. For fractures with poor fit quality, the clustering and centerline extraction processes are re-examined and adjusted. The fitted sine curve is visualized and compared with the original fracture to visually verify the effect. The parameter estimates of each fracture are summarized to form a quantitative analysis result table. Based on these parameters, the inclination angle between the fracture plane and the borehole axis, the direction of the fracture in three-dimensional space, etc., can also be calculated.

[0033] Step 3: Based on the modeling, perform roughness coefficient JRC calculation and sampling method sensitivity analysis.

[0034] In this embodiment of the invention, step 3 includes: The macroscopic orientation trend and microscopic roughness residuals of cracks are separated using smooth spline functions. The effects of three sampling methods—equal-interval sampling, equal-arc-length sampling, and multi-scale adaptive sampling—are compared. The study investigated the influence of the first-order slope energy index (length-normalized) and the JRC value. Various geometric features of the fracture were extracted and decoupled. Multiple sampling was employed, and the JRC value was calculated based on the sampling points. Since the JRC formula clearly involves calculations related to fracture scatter points, different sampling methods yield different results. A better sampling method naturally leads to more accurate JRC results and reduces uncertainty, ultimately yielding the roughness. Results showed that equal-interval sampling tends to overestimate the roughness of steep sections (JRC deviation reaches 20%). Equal-arc-length sampling, through arc-length reparameterization, improves recognition accuracy by 15% in high-curvature regions. Multi-scale adaptive sampling balances computational efficiency and recognition accuracy (JRC confidence interval is reduced by 12%). Bootstrap resampling (1000 times) quantified the impact of sampling interval and window size on the uncertainty of the results, revealing a weak positive correlation between fracture area and JRC (correlation coefficient 0.32). A 5-segment segmented calculation strategy can control the mean error of JRC within 5%.

[0035] Step 31: Contour extraction and trend separation. First, cracks are extracted and preprocessed using an improved U-Net network. Then, trend separation and residual extraction are performed. After reparameterizing the crack skeleton by arc length, the trend line is fitted using a smooth spline or low-order polynomial. The high-frequency roughness residual is obtained by the difference between the original curve and the trend line, eliminating the interference of the change in direction on the roughness. Step 32: Calculate roughness indices using multiple sampling methods, compare the differences between equal-interval sampling, equal-arc-length sampling, and multi-scale adaptive sampling, and employ denser sampling in high-curvature regions (e.g., curvature > 0.1 mm). -1 (Sampling interval halved), and first-order slope energy index with length normalized based on residual data. ; Step 33: JRC Inversion and Uncertainty Assessment. An endpoint missing measurement mask is introduced to ensure consistent effective calculation lengths for different cracks. JRC mapping and sensitivity analysis are completed, and empirical formulas are applied to... The values ​​are converted into JRC indices, and the Bootstrap method is used to quantify the impact of sampling interval, window size and detrending intensity on the results. The JRC mean and 95% confidence interval are output, providing texture feature basis for fracture connectivity assessment.

[0036] In this embodiment of the invention, it includes: Step 311, Feature Decoupling (Detrending): The low-frequency shape (macroscopic profile) and high-frequency roughness (local undulation) features of the crack are separated by the smooth spline method. The core is to minimize the joint objective function of data fitting error and curve smoothness to obtain a smooth curve that reflects the overall trend. Then, the smooth curve is subtracted from the original data to obtain a residual sequence containing only high-frequency roughness information, which provides a clean data carrier for subsequent JRC calculations. Step 312: Select the appropriate sampling method for different scenarios and calculate the roughness index. Equal-distance sampling, sampling at equal intervals along the x-axis, focuses on lateral comparability, and calculates the distance between adjacent points based on the elevation difference and lateral step size. It is suitable for lateral feature comparison, and its expression is: ; Equal arc length sampling: Convert parameters to cumulative arc length to ensure consistent sampling density across steep and gentle sections, avoiding underestimation of undulations in steep sections and improving efficiency. Geometric fairness of computation, roughness index under equal arc length sampling The expression is: ; Multi-scale adaptive sampling: Combining multi-scale analysis with dynamic curvature adjustment of sampling density, dense sampling is performed in complex regions and sparse sampling is performed in smooth regions to balance feature capture accuracy and computational efficiency; the set of sampling points obtained by multi-scale adaptive sampling is represented as: ; Curvature-weighted optimization: Assigning high weights to regions with high curvature (severe bending) to correct... Indicators that enhance sensitivity to irregular features of fracture peaks; Step 313, JRC inversion utilization The mapping relationship with JRC's engineering experience (logarithmic calibration formula) will be used to map different sampling methods and scales. The value is converted into a dimensionless JRC value to achieve a quantitative conversion from geometric features to roughness coefficients; in engineering, empirical formulas are used to convert... Mapped to JRC, using a logarithmic calibration expression: ; Where JRC(w) represents the joint roughness coefficient on window w, which is a dimensionless empirical level; The calibration data is related to specimen dimensions, test line length, and unit consistency. This represents the roughness index value on window w; Step 314: Bootstrap uncertainty assessment calculates the standard error and confidence interval of the JRC by repeatedly resampling and propagating the discrete point selection, taking into account the influence of noise factors, quantifying the reliability of the results, and finally outputting the JRC value of the crack with confidence. The formula is as follows: ; in, This represents the data from the b-th resampled dataset. Indicates detrending, The computation chain of JRC, where B represents the number of resampling steps. This represents the mean obtained from multiple resampling operations, and SE represents the standard error. This represents the empirical quantile of the 1-α confidence interval.

[0037] Step 4: Based on JRC calculations and sensitivity analysis, perform three-dimensional reconstruction and hole filling optimization of the multi-hole fracture network.

[0038] In this embodiment of the invention, step 4 includes: Establish a technical framework for coordinate mapping, connectivity assessment, uncertainty optimization, 3D reconstruction and supplementary point generation: Based on the borehole opening or bottom coordinates, the two-dimensional fracture parameters (sine model, JRC value) are inverted into fracture patches in three-dimensional space through Lie algebra SE(3) pose transformation (rotation matrix + translation vector). The normal angle (weight 25%), spatial distance (weight 25%), fracture size (weight 20%) and JRC difference (weight 30%) are integrated to construct a multi-dimensional consistency index, which is mapped to the connectivity probability in the 0-1 interval by the sigmoid function. The breadth-first search algorithm is used to extract fracture clusters, and the representative plane of fractures in the cluster is fitted by the least squares algorithm. By calculating the fracture density variance and connectivity probability entropy of each point in space, high uncertainty areas (such as the depth below 5000mm and the borehole gap area) are located, providing accurate geological basis for roadway support design and water inrush prevention. Step 41: 3D coordinate mapping and fracture extraction. The borehole opening is taken as the origin (x=0, y=0, z=0), where x represents the parallel direction of the roadway, y represents the perpendicular direction of the roadway wall, and z represents the depth. A homogeneous transformation matrix T_BH is constructed using the borehole pose (borehole (x0, y0, z0), dip angle θ, azimuth angle φ) to map the local fracture parameters to the global coordinate system. Based on the borehole opening or bottom coordinates, the single-hole fracture parameters (sine curve, JRC) are transformed to the global coordinate system. The 3D plane corresponding to the fracture is inverted through the SE(3) pose transformation and clipped into finite patches according to the observation length. For complex fractures, Fourier series or spline surface extension models are used to improve the accuracy of 3D representation. The fracture is extracted by improving the U-Net network and the fracture preprocessing is performed to filter out interference information in the fracture. Step 42: Connectivity analysis. Construct a multi-dimensional consistency index. At the geometric level, features such as normal angle, plane constant difference, and Hasdorf distance are used. At the texture level, JRC difference is introduced, and weights are designed to weight the two features. The comprehensive score is mapped to the connectivity probability in the 0-1 interval through logistic regression. A weighted undirected graph is constructed and connected gap clusters are extracted. Within the cluster, the total least squares joint fitting is used to fit the representative plane. Step 43: Uncertainty and Hole Filling Optimization. By calculating the fracture density variance and connectivity probability entropy of each point in space, high uncertainty areas (entropy value > 0.6) are located. Candidate hole filling areas are selected by combining engineering constraints (such as roadway reachability and drilling depth < 100m). Hole filling location optimization involves traversing spatial coordinates and using the information gain maximization algorithm to evaluate the uncertainty reduction benefits of each candidate point. The top three optimal locations are selected. Finally, the fracture network, connectivity probability, and hole filling scheme are presented through three-dimensional visualization, forming a closed loop of identification-modeling-decision.

[0039] In this embodiment of the invention, it includes: The core of 3D reconstruction and hole repair optimization of multi-hole fracture networks is coordinate unification, 3D inversion, connectivity determination, and hole repair optimization. Through full-process quantification and probabilistic modeling, accurate reconstruction of fracture networks and efficient hole repair suggestions are achieved.

[0040] Step 411: Data Integration and Coordinate Unification. Establish the rigid body transformation relationship between the borehole local coordinate system and the global coordinate system to unify different borehole data into the same spatial system and ensure consistency in cross-hole analysis; at the same time, integrate the core data of sinusoidal model parameters, JRC, and borehole radius to provide a unified quantitative basis for 3D reconstruction. Step 412: 2D to 3D Inversion and Patch Construction. Based on the sinusoidal model parameters of the borehole unfolding diagram, the normal vector and plane equation of the fracture plane are inversely derived and mapped to the global coordinate system through pose transformation. Combined with the effective observation length of the borehole, a finite fracture patch is cut out (to avoid computational redundancy), and the metadata of the center and scale is recorded to complete the transformation from 2D measurement to 3D morphology. Step 413: Cross-hole connectivity determination. Using geometric and texture multi-dimensional indicators, the connectivity probability of different borehole fractures is quantified. Geometric consistency: Calculate the normal angle, planar position difference, and spatial distance, combining similarity in attitude (dip angle, dip direction) and spatial coplanarity to determine the fit between the fracture's spatial orientation and location. Texture corroboration: Introduce JRC difference to supplement roughness consistency. Perform probability mapping, obtain a comprehensive score through weighted fusion of multi-dimensional indicators, and then convert it into connectivity probability using the Sigmoid function to form a symmetric probability matrix. Obtain a comprehensive similarity score by weighted fusion of geometric and roughness indicators. Comprehensive similarity... The formula is: ; in, As weight, This is a scaling constant used to normalize the various indicators. For roughness inversion of the corresponding intersecting window of the facet; Logistic regression mapping of connectivity probability converts the comprehensive similarity score into connectivity probability, facilitating threshold determination and graph model calculation; the sigmoid function is used to achieve a monotonic mapping from score to probability, obtaining the probability that two fractures belong to the same continuous fracture volume. ; in, α represents the probability that patches i and j belong to the same continuous fracture body, and b are calibration coefficients determined through training data. Step 414: Probabilistic Graph Construction and Cluster Extraction. Using fracture patches as nodes and connectivity probabilities as edge weights, a weighted undirected graph is constructed. A probability threshold is set to filter valid edges. Interrelated fracture clusters are obtained through connected component extraction or community detection algorithms. The overall plane or surface of the cluster is then fitted using the global least squares method to characterize its geometric features. The connectivity relationships of the fracture patches are transformed into a graph structure. Using connected component extraction or community detection algorithms in graph theory, the set of closely related fracture patches, i.e., fracture clusters, is identified. ; Where G is the crack diagram. It is a set of nodes (faces). The set of edges (with a connectivity probability greater than a threshold) Clusters can be extracted using connected components, spectral clustering, or modularity maximization. Step 415: Uncertainty Propagation and Aggregation. The uncertainty of a single parameter (such as the parameters of a sinusoidal model) is propagated to the position and orientation parameters of the 3D patch through the Jacobian matrix. At the same time, the covariance of parameters within the cluster is aggregated to quantify the uncertainty of the fracture cluster and the connectivity probability, thereby improving the reliability of the model. Step 416: Hole Filling Optimization. Under the constraints of engineering reachability and construction, efficient hole filling selection is achieved. The constraint range is defined, and the reachability is defined in combination with construction conditions. Based on the existing borehole distribution, a set of candidate holes that are feasible for the project is generated. High-value areas are identified by integrating cluster plane uncertainty and connectivity probability to identify high-uncertainty-high-value hole filling hotspots. The selection objective is to use the expected information gain of candidate holes as the objective function to select the optimal hole filling combination, maximizing the improvement of understanding of fracture connectivity under the premise of limited resources.

[0041] This invention employs a comprehensive technical framework encompassing data preprocessing, intelligent identification, parameter fitting, roughness quantification, 3D reconstruction, and optimization decision-making. This framework enables pixel-level separation of fractures in borehole imaging images, quantitative characterization of sinusoidal and complex fractures, spatial correlation analysis of multi-bore fracture networks, and monitoring optimization. Specifically, it includes: improving fracture identification accuracy based on an improved U-Net network (integrating ResNet residual modules and multi-scale features); jointly optimizing sinusoidal model parameters using differential evolution and weighted least squares methods; reliably calculating the JRC value of complex fractures by combining trend separation and Bootstrap uncertainty assessment; and constructing a connectivity probability model through weighted fusion of geometric and texture features, combined with particle swarm optimization to identify high-uncertainty regions and select optimal borehole locations. This invention effectively solves the problems of low identification accuracy, large parameter fitting deviations, and insufficient reliability of 3D reconstruction in traditional methods. It provides a scientific basis for coal mine roadway surrounding rock stability assessment, disaster early warning, and support scheme optimization, meeting the demands of intelligent coal mining for high-precision, comprehensive, and real-time surrounding rock monitoring. The method of this invention is applicable to the quantitative analysis, three-dimensional network construction and monitoring optimization of rock fissures in underground coal mines. It can be widely used in scenarios such as rock stability assessment, early warning of water inrush / gas outburst disasters, and support scheme design in intelligent coal mining, providing key technical support for safe coal mine production.

[0042] I. Data Preparation and Preprocessing Specifications for this Invention: To verify the feasibility and effectiveness of the present invention, boreholes in a coal mine roadway (buried at a depth of 500-700m, with sandstone as the lithology) were used as the experimental object. Image data of the boreholes were collected, and finally, the imaging data of 6 boreholes were used as the verification object to evaluate and visualize each task.

[0043] Data format: Borehole imaging images are in JPG format, with borehole parameters recorded synchronously (diameter 94.25mm, imaging resolution 0.1mm / pixel, borehole coordinates, inclination angle, and azimuth angle); Preprocessing operations: Grayscale conversion, using a weighted average method to convert RGB images into grayscale images, with a grayscale range of 0-255; Gaussian denoising, 5×5 kernel, σ=1.5, to eliminate sensor noise; Adaptive histogram equalization (CLAHE), clipLimit=2.0, 8×8 grid, to enhance the contrast between the crack and the background; Data augmentation: Cropped to 244×244 pixel samples, the dataset was expanded through random cropping (15% range), rotation (0-360°), and mosaic enhancement.

[0044] II. Detailed operational examples of each step in this invention: First, it is clear that the core of the problem is to classify fractures, and the model should take into account the influence of interference factors such as rock texture, drilling marks, and mud contamination. Therefore, this invention designs five methods to solve the problem, and finally, through comparison, the method with the best effect is determined. The five methods include traditional modeling methods, adaptive thresholding methods, machine learning methods, traditional deep learning methods, and deep learning methods. The first step for all methods is to unify coordinates and preprocess image data. The task data is read, and the horizontal expansion is defined as the circumference of a cylinder, with the vertical axis representing depth. Then, the data is preprocessed. The original data contains noise, uneven lighting, and other interference, so data augmentation is performed first. The original image is converted to grayscale, because in computer vision tasks, color channels do not significantly affect the model's accuracy; instead, they provide many meaningless features. Gaussian denoising removes most of the noise from the image, and histogram equalization improves contrast, which reduces the impact of lighting on the data to some extent and makes cracks more visible. Finally, MedianBlur is used to suppress impurities such as mud.

[0045] Traditional modeling methods for crack extraction combine gradient detection (edge ​​detection) and morphological operations, employing multi-directional gradient detection, Canny edge detection, and LoG edge detection to detect image edge information. Since cracks inevitably exhibit significant edge differences, different edge detection methods can focus on different features. Multi-directional gradient detection has strong directional characteristics, addressing the challenge of capturing complex crack orientations that are difficult to capture with single-directional gradients. Canny edge detection offers good noise resistance and, through non-maximum suppression, retains only local maxima along the gradient direction, refining wide edges to single-pixel widths and accurately locating the crack contour, thus resolving edge blurring and excessive width issues that may occur with other gradient detection methods. A dual-threshold mechanism is employed: a high threshold identifies strong edges, a low threshold marks weak edges, and only weak edges connected to strong edges are retained. This connects broken crack edges into a continuous contour, making it particularly suitable for crack edge breaks caused by uneven illumination and local blurring. LoG edge detection is sensitive to blob-like structures. The LoG operator essentially performs a Laplacian operation on a Gaussian-smoothed image. The Laplacian operator is sensitive to abrupt changes in grayscale (edges) and extreme grayscale regions (such as the center of a blob), making it particularly suitable for detecting cracks in circular, elliptical, and other blob-like structures (such as small-area hole-like cracks), thus compensating for the weak response of Sobel and Canny to such curved structures. Furthermore, compared to the Sobel operator, LoG has a more uniform response to edges and is less affected by edge direction, enabling it to more stably identify crack edges in complex textured backgrounds. This invention further employs a voting method, meaning that if two of the three methods identify a crack, then it is indeed a crack. The three methods complement each other in three dimensions: directional sensitivity, noise resistance and continuity, and adaptability to special morphologies: Sobel is used to quickly capture multi-directional edge intensity, Canny is used to accurately extract continuous edges, and LoG is used to supplement the detection of blob-like structures. By combining criteria (such as requiring at least two methods to agree before classifying a crack), false detections (e.g., noise being misjudged) and false negatives (e.g., special crack morphologies not being identified) caused by a single method can be significantly reduced, ultimately yielding a more reliable preliminary crack region. Figure 3 As shown, the traditional modeling method does not perform well. Even if it can detect some cracks, it fails to detect many other non-crack parts. This may be because the image is too complex, and the edge features in many places are similar to the cracks.

[0046] The adaptive thresholding method first enhances the difference between cracks and the background through edge detection and texture analysis. Sobel edge detection, LBP texture detection, and Hessian are used to process the image. Sobel edge detection extracts edge features, LBP texture analysis extracts texture features based on the difference between texture and background, and Hessian distinguishes between linear (cracks), planar (background), and corner points to highlight the crack structure. The weights of these three methods are set as follows: multi-scale response 50%, Sobel edge detection 30%, and current response 20%, thus forming multi-scale, multi-directional features. These features are then normalized to prevent any single feature from excessively affecting other feature representations. Finally, Otsu's automatic thresholding method is used to automatically calculate the threshold that minimizes the intra-class variance between the foreground and background. Morphological operations are then performed to connect the scattered contours as much as possible before outputting the detection results, such as... Figure 4 As shown, the extraction effect is still not good. Although it detects more details compared to traditional modeling methods, it also identifies more non-crack areas.

[0047] The core principle of machine learning methods is supervised learning based on random forest classifiers. Random forests consist of multiple decision trees. By learning from a large number of features, a classification model is built to classify pixels or regions as cracked or non-cracked. First, pseudo-labels are generated using the results of mathematical and statistical modeling methods to provide training data for subsequent supervised learning. Because there is insufficient manually labeled data, pseudo-labels that approximate true labels are generated using the results of two existing methods to address the lack of training data. Then, the random forest algorithm is trained using the `sklearn.ensemble.RandomForestClassifier` library in Python, with settings `n_estimators=60` (the number of decision trees; more trees generally result in stronger generalization ability, but also higher computational cost) and `max_depth=10` (the maximum depth of the decision trees, controlling tree complexity and preventing overfitting). The classifier learns the mapping relationship between features and crack / non-crack labels, thus gaining the ability to classify new data. The confidence of the classification results is evaluated based on gradients and standard deviations. Predictions with low confidence are filtered out, while those with high confidence are retained to further improve the accuracy of crack identification. Figure 5 As shown, machine learning methods, by employing a supervised learning strategy, make the detected details more accurate and the overall impression better, but they also detect many non-cracked parts.

[0048] The core principle of traditional deep learning methods is based on multi-scale feature extraction and cluster analysis. Multi-scale feature extraction can capture image features at different scales, providing a more comprehensive description of cracks; cluster analysis groups pixels with similar features into one class, distinguishing between cracks, background, and interference. First, Gaussian filtering (smoothing the image with varying degrees of smoothness at different scales, capturing structures of different sizes) is calculated using different scales (1.0, 1.5, 2.0, 2.5). Multi-scale features can cover various crack features from small to large, making the feature representation more comprehensive. Next, the K-means clustering algorithm is used, implemented using the KMeans class from Python's scikit-learn library, setting n_clusters=3 to divide pixels into three classes: background, interference, and cracks. K-means iteratively optimizes, clustering pixels with similar features into one class, using unsupervised learning to initially delineate possible crack regions. Then, DBSCAN refinement is performed. Using the scikit-learn library's DBSCAN, eps=0.5 (neighborhood radius, controlling neighborhood size) and min_samples=5 (minimum number of samples required for core points) are set to refine the K-means clustering results. DBSCAN can discover clusters of arbitrary shapes and is insensitive to noise, allowing for more accurate segmentation of crack regions and removal of interference and noise points. Finally, crack clusters are automatically identified based on intensity statistical features, grouping pixels belonging to the same crack into a single cluster to obtain the identification result. For example... Figure 6 As shown, although the results lack many details of impurities, they still result in a loss of details of the cracks themselves. Therefore, multi-scale feature extraction and cluster analysis are not the most ideal methods.

[0049] Through experiments and testing using the methods described above, it is evident that traditional methods are ineffective for this task. For the automatic identification of crack pixels in borehole imaging unfolded images, a modern deep learning method is considered: a convolutional neural network (CNN) that performs semantic classification of pixels, namely the U-Net segmentation network. Many optimized versions of this network are used in medical segmentation tasks, such as cardiovascular and cerebrovascular segmentation, and coronary artery segmentation. These tasks often involve data that is small in volume and difficult to obtain (due to the limitations of medical data and the scarcity of open-source data), poor image quality (most are angiography and X-ray images, whose quality is affected by many factors and is far inferior to natural images), and minimal foreground representation (due to vascular data, blood vessels only occupy a small portion of the image, with the majority being background). Considering these factors, the U-Net series is well-suited for handling such tasks. U-Net consists of an encoder and a decoder. Its essence is to map to a deep space to extract high-dimensional features of the image. The features of a few dimensions may correspond to most of the key features of the crack, so that the model can learn these deep features and make the model better at classification. At the same time, the model's loss function can be backpropagated based on the loss of the classification result and the true result to adjust the model to learn in a direction that is closer to the true value, so that the model achieves the best performance.

[0050] Therefore, this invention employs the U-Net series of algorithms for final recognition, optimizing it by improving the original U-Net. A residual network, ResNet, is used to replace part of the encoder and decoder models, storing the data after each decoding layer and concatenating the feature layers before outputting. This achieves a hybrid sampling effect. The reason is that the feature layers are relatively shallow at the beginning of the network, so ordinary convolutions can be used to extract surface adjustments and local features first. Then, when the network has deeper layers, the residual module in ResNet replaces the ordinary convolutions, preventing feature loss due to excessive depth. Furthermore, the feature fusion technique at each layer allows the model to extract global and multi-scale features, achieving comprehensive feature extraction and better model performance. Finally, to demonstrate the model's advantages, this invention compares it with YOLOv11 and U2Net. YOLOv11 is currently an excellent object detection and semantic segmentation algorithm, while U2Net is a leading algorithm for saliency segmentation, specifically designed for binary classification tasks, making it well-suited for the task presented in this paper.

[0051] Because the amount of data is relatively small, it needs to be expanded. Firstly, since the data collected in this invention is in a 244*1350 pixel format, with a small width and a large length, the invention can first segment each data point into 244*244 pixel images. This expands the data and facilitates training. Then, data augmentation methods such as random cropping, rotation, contrast adjustment, and mosaic data augmentation are applied to further expand the data.

[0052] In embodiments of the present invention, such as Figure 7 As shown in (a), (b), and (c), the optimized UNet performs best. While YOLOv11 has good confidence and labels the cracked areas, it also detects some non-cracked areas. U2Net, on the other hand, detects the cracked areas very well, although some impurities and disconnections are present. The optimized UNet, however, completely identifies the crack outline with very little impurity. Furthermore, the improved UNet has a map accuracy of 92.6%, which is about 8 percentage points higher than that of U2Net. Therefore, the improved UNet in this invention performs best, and subsequent tasks in this invention will be based on this model for detection.

[0053] In embodiments of the present invention, such as Figure 8 As shown in (a), (b), and (c), each set of images includes an image feature enhancement effect, a visualization binary image of crack detection, and an overlay image of the detection effect and the original image. It can be seen that whether it is the first image, which has clear and few cracks but many interfering impurities (such as texture, mud, etc.), the second image, which contains many cracks, or the third image, which is relatively blurry and the cracks are not clear, the recognition effect is very good, and they basically fit the cracks in the original image.

[0054] This paper proposes a multi-model comparison and automatic fitting algorithm for sinusoidal fracture fitting in coal mine roadway surrounding rock fracture analysis. The core of this algorithm is to extract geometric feature parameters from the binary fracture image identified in step 1 and construct a mathematical model to characterize the spatial distribution of the fracture. Since a sinusoidal fracture is the intersection line of a planar fracture and a cylindrical borehole, the position and orientation of the fracture plane in three-dimensional space can be inferred by fitting a sinusoidal curve. The parameters of the sinusoidal curve (amplitude R, period P, phase, and centerline position C) are directly related to the geometric characteristics of the fracture. In practice, the fracture pixels need to be clustered first, then the centerline extracted, and finally a sinusoidal model fitted. Challenges such as image noise and fracture deformation must be addressed, and multiple fractures must be distinguished and the fitting error evaluated.

[0055] First, the coordinates must be confirmed, as this involves physical properties. Since the data in this invention is an image based on pixels, the pixel data and actual coordinates must be processed first. Then, the model from step 1 is used to identify the data and obtain a binary image of the crack. Because the second task requires analysis of sinusoidal cracks, the identified cracks need to be processed. The identification results from the model in step 1 may incorrectly identify some non-sinusoidal cracks and some small interferences; measures must be taken to remove these potential interferences.

[0056] Data preprocessing and initial clustering: First, the binary image of the crack from step 1 is loaded, and preprocessing such as removing small-area noise and filling small holes is performed. This reduces interference in subsequent clustering and improves accuracy. Basic features such as area, perimeter, and centroid of each region are calculated, and regions with excessively small areas are filtered out. For regions with abnormal shapes, the minimum enclosing circle method is used for filtering. This is because the sinusoidal cracks required in this task are characterized by their narrow and small area, while many interference points, such as misidentified drilling marks and non-sinusoidal cracks, are mostly irregular elliptical. Therefore, by calculating the center point and area of ​​each connected region and the area of ​​its minimum enclosing circle, and using the ratio of the original area to the minimum enclosing circle, a threshold greater than 0.35 is considered a non-sinusoidal crack, thus optimizing the overall detection effect to better meet the requirements of this invention.

[0057] Then, the spatial distance and morphological similarity between each connected region are calculated to construct a similarity matrix. The regions are then merged using a hierarchical clustering algorithm, taking into account spatial distance, axial position differences, and morphological similarity to ensure that regions with the same crack are correctly aggregated.

[0058] After clustering, the centerlines of the fractures need to be extracted. For each clustered fracture region, a distance transformation algorithm is applied to obtain a distance map, and then a morphological thinning algorithm (Zhang-Suen or Guo-Hall algorithm) is used to extract the skeleton lines. Since the skeleton lines may have branches, a pruning algorithm is used to remove short branches, and the main trunk is obtained through a graph theory-based longest path algorithm. Spline interpolation or moving average filtering is used to smooth the main trunk lines and reduce jagged noise. For fractured fractures, a line segment connection algorithm based on directional consistency and distance constraints is used to bridge the fractured parts to ensure the continuity and integrity of the centerline.

[0059] Simultaneously, since the centerline needs to be extracted, sampling along candidate line segments is required. This invention employs a combination of three sampling methods for final sampling. First, equal-interval sampling, as the most basic sampling method, ensures that sampling points are evenly distributed circumferentially without bias. This is suitable for scenarios with relatively uniform crack morphology and low noise. Next, adaptive sampling addresses the need for denser sampling at complex crack morphology sections (such as bends and bifurcations) to accurately reconstruct their shape, as these sections often contain more crucial feature information. Conversely, straight sections have less feature information, and sparse sampling reduces computational load. This sampling method optimizes sampling efficiency and improves fitting accuracy for complex cracks while ensuring the capture of core crack features. Finally, curvature key point sampling is used. Curvature key points are where crack morphology changes most drastically, concentrating the core feature information of the crack. Focusing sampling on these key points allows for the reconstruction of the crack's key morphology with the fewest sampling points. This is suitable for scenarios with limited computational resources but requiring accurate capture of crack features, reducing computational costs while preserving as many important crack morphological features as possible. A (x, y) point set is generated (x is the circumferential coordinate, y is the depth coordinate). Simultaneously, for multiple y values ​​corresponding to the same x-coordinate, duplicate points are removed and the average is taken to eliminate the influence of local thick boundaries (such as jitter caused by imaging noise at the edge of a crack) and discrete noise. Finally, the point sets obtained from the three sampling methods are weighted and fused. Point set A is obtained by sampling at equal intervals and assigned a weight (w_A = 0.3) (basic distribution, ensuring uniformity); point set B is obtained by adaptive sampling and assigned a weight (w_B = 0.5) (key segments, high accuracy); and point set C is obtained by sampling curvature key points and assigned a weight (w_C = 0.2) (core features, strong representativeness).

[0060] Parameter initialization and nonlinear fitting: To improve the convergence speed and stability of the optimization algorithm, reasonable initial values ​​need to be provided for the sinusoidal model parameters. The initial value of amplitude R is taken as half the longitudinal span of the centerline; the period P is set according to the known circumference of the borehole (approximately 94.25 mm); the phase is estimated based on the starting position of the centerline; the centerline position C is taken as the average ordinate of the centerline. Simultaneously, the feasible region of the sinusoidal model parameters is adaptively set according to the circumferential range (lateral span) and depth range (longitudinal span) of the data. Amplitude R: upper limit is half the depth range, lower limit is a small value (e.g., 1 mm, to avoid fitting a straight segment as a sine). Period P: upper limit is the total circumferential length, lower limit is 1 / 2 of the borehole circumference (consistent with the prior knowledge of matching the sinusoidal period with the developed diameter). Phase: range is (0, 2...). (Physical period of phase).

[0061] The Levenberg-Marquardt algorithm is used for nonlinear least squares fitting, iteratively adjusting parameters to minimize the error between the model prediction and the actual centerline. A weighting function (such as a double-weighting method or Huber weighting function) is introduced to assign lower weights to outliers, improving fit robustness. Iterative optimization and outlier handling employ an iterative reweighted least squares (IRLS) strategy. In each iteration, the weights of points are adjusted based on the residual size (the larger the residual, the lower the weight), and the fitting is re-fitted after updating the weights until the parameters stabilize or the maximum number of iterations is reached.

[0062] Finally, the results are evaluated and visualized, and the goodness-of-fit index (Rfit) is calculated. 2 The fitting quality is evaluated using RMSE (Real-Time Sequence of Analysis). For fractures with poor fitting quality, the clustering and centerline extraction processes are re-examined and adjusted. The fitted sine curve is visualized and compared with the original fracture to visually verify the effect. The parameter estimates of each fracture are summarized to form a quantitative analysis results table. Based on these parameters, the dip angle of the fracture plane to the borehole axis and the direction of the fracture in three-dimensional space can also be calculated, providing information for subsequent three-dimensional reconstruction.

[0063] In embodiments of the present invention, such as Figure 9a , Figure 9b and Figure 9c As shown in (a), (b), (c), and (d), each image includes four visualizations: the original image (a), the crack clustering result (b), the crack skeleton extraction result (c), and the overlay of the original image and the fitted curve (d), allowing for a comprehensive observation of the fitting effect. The figures demonstrate that the model can fit a sine curve well, whether for cases with many cracks like the first image or for cases with lower resolution like the third image, with an average R0. 2 With a success rate of over 95%, the results are excellent.

[0064] This paper proposes an adaptive multi-dimensional crack roughness calculation algorithm for complex cracks. The core objective is to calculate the JRC value to represent the crack complexity, while also exploring the impact of different sampling methods on the JRC value. The JRC formula is related to the geometric characteristics of the crack and essentially reflects the degree of undulation in the contour line; the more drastic the change, the larger the value, and the higher the corresponding JRC value. Therefore, the solution to this problem still requires extracting various geometric features of the crack and then further calculating the JRC value based on the sampling points for different sampling methods. Because the JRC formula clearly involves calculations related to crack scatter points, different sampling methods will yield different results. A better sampling method will naturally bring more accurate JRC results, ultimately yielding the roughness.

[0065] Therefore, the first step of this task is still to use the algorithm in step 1 to extract the cracks from each image in the dataset. Then, since this task still needs to analyze sinusoidal cracks, it is necessary to use the same preprocessing method as in step 2 to filter out various interference factors and obtain a clean binary crack image.

[0066] Because this task still requires a lot of physical information about the cracks and the mapping of the centerline, after extracting the cracks, it is necessary to map the pixels to their coordinate values. We define the horizontal direction as the circumferential length and the vertical direction as the axial depth, and convert the pixels to millimeters using image resolution to ensure that the coordinates of different images are physically comparable.

[0067] Three sampling methods were employed and compared: equidistant sampling, equal arc length sampling, and adaptive multi-scale sampling. In equidistant sampling (taking points at equal intervals along the x-direction), if the fracture has steep sections (e.g., large local slopes), equidistant x-sampling leads to sparse point distribution in these sections (because steep sections change rapidly in the x-direction), thus underestimating the calculated slope. The actual slope of a steep section is large, but due to fewer points, the calculated average slope is smaller. Arc length, on the other hand, is the actual length along the fracture centerline. Equal arc length sampling takes points at equal intervals along the arc length of the centerline. This ensures a uniform number of points within the same arc length interval, regardless of whether the fracture is steep or gentle. Because steep sections change rapidly in the x-direction, the number of points is consistent with that of gentle sections with the same arc length, allowing for a denser capture of the steep section's morphology. Multi-scale adaptive sampling combines the advantages of equidistant sampling and another curvature-weighted sampling method. It first analyzes the contour line at multiple scales and then dynamically adjusts the sampling density based on the complexity of local features. Because the roughness of a crack is related to the observation scale, a small window will capture fine details of the surface (high local roughness), while a large window will smooth out fine details. Therefore, multi-scale windows are needed for sampling at different scales.

[0068] The purpose of trend separation is to isolate macro-level trends and focus on micro-level residual trend separation and residual construction. The macro-level orientation of the crack (low-frequency component) is extracted using a smoothed trend term. The residual sequence (high-frequency component) is obtained by subtracting the trend term from the original profile. The residuals carry the local roughness fluctuations. This step is meaningful if the goal is to analyze micro-level roughness because the residuals can eliminate the interference of macro-level trends on roughness assessment (for example, a crack that is generally sloping but locally smooth will have small residuals and a small JRC after trend separation, which aligns with the intuition that local smoothness indicates low roughness).

[0069] After sampling using three sampling methods and separating the trends, the JRC (Joint Rating Scale) can be calculated. First, the JRC is calculated based on the sampling points and the formula, which is the square root of the average of the squares of the slopes between adjacent points. Then, the JRC is derived from this. After calculating the JRC, uncertainty assessment and robustness verification are performed to quantify the reliability of the results. Resampling with replacement is performed on the crack point series using Bootstrap resampling to generate multiple sets of virtual point series. Link recalculation: For each set of virtual point series, the trend separation, roughness, and JRC processes are repeated to obtain multiple sets of JRC results. The mean, dispersion, and confidence interval of the JRC are calculated to quantify the impact of the sampling method on the JRC. Finally, considering the diverse shapes of the detected cracks, a segmented approach to calculating the JRC is the most reliable and precise method. Therefore, the JRC is calculated in five segments, and the mean is calculated (because considering that approximately two pixels represent 1mm on average, choosing 5 is a more reliable approach to balance accuracy and precision). Finally, the roughness of the crack is obtained. Finally, the sampling effect and roughness of each crack are visualized.

[0070] In embodiments of the present invention, such as Figure 10 As shown in (a) and (b) in the figure, the overall detection effect and contour extraction effect of the crack in a randomly selected data image are very ideal; the visualization of the crack sampling results is shown in the figure. Figure 11a , 11b As shown in Figure 11c, clearly different methods produce different results in areas of rapid or smooth change. Equal arc length sampling produces more compact details in areas with large arc variations, while multi-scale sampling combines the advantages of both methods, taking into account both local and global sampling effects; for example... Figure 12a and Figure 12b As shown, the JRC for each segment and the JRC calculated by different sampling methods are visualized. It is clear that for this case of drastic change, the other three sampling results, except for the segmented sampling results, all yielded higher JRCs.

[0071] In embodiments of the present invention, such as Figure 13 As shown, the core function of the JRC value comparison chart of different sampling methods is to quantify the systematic bias of the sampling methods. Equal-interval sampling tends to overestimate JRC values, while multi-scale or piecewise averaging tends to be conservative. Equal-interval sampling, by taking points at fixed pixel intervals, overcaptures high-frequency fluctuations in the fracture locality, such as sharp protrusions and small depressions, leading to an overestimation of the calculated mean squared slope (Z2 parameter), and consequently, an overestimation of the JRC value. Multi-scale sampling, through multi-scale curvature fusion, balances local details with the overall trend, while piecewise averaging weakens the influence of extreme values ​​through multiple segmented averages, thus yielding more robust results. Finally, three images were selected as key case studies, not randomly, but because these three images cover typical fracture morphology types: Figure 13Medium 3-1-1 (medium roughness, simple shape): As a benchmark case, the results of different sampling methods are consistent, which verifies the rule that the sampling method has little impact when the shape is simple. This is because the local undulations of simple shapes (such as gentle curves) are uniform, and the difference in the calculated slope distribution is small regardless of which method is used to select points. Figure 13 3-2-1 (extremely high roughness, complex morphology): As an extreme case, the sampling methods differ significantly, which confirms that the sampling method is crucial when the morphology is complex. In complex morphologies (such as frequent turns and violent undulations), equal-interval sampling will over-amplify local details, while multi-scale sampling can filter out noise-level fluctuations and retain the true roughness features. Therefore, the results of the two methods have a large deviation. Figure 13 3-3-1 (Multiple fractures mixed, with large differences in roughness): As a mixed case, it includes both low roughness (fracture 2) and high roughness (fractures 1, 3, 4), which verifies both simple morphology and consistent method (fracture 2) and complex morphology and large differences in method (fractures 1, 3, 4), further strengthening the necessity of classification assessment.

[0072] Therefore, the following conclusions can be drawn: the difference in sampling methods is highly correlated with the fracture morphology. The calculation of JRC depends on the slope of the contour points, and the sampling method determines which points are selected. For fractures with simple morphology (such as gentle curves): any method of sampling results in small slope changes between adjacent points, so the difference in Z2 (mean of squared slopes) is small, and the JRC results are relatively similar. For fractures with complex morphology (such as turbulent undulations): equal-interval sampling may just take points on the steep section, resulting in an overestimated Z2. Multi-scale sampling balances the steep and gentle sections through multi-scale smoothing, making the Z2 more reasonable, and therefore the JRC results differ greatly.

[0073] In summary, differentiated assessment strategies should be selected based on the type of crack. Although the multi-scale sampling method is generally superior to the first two methods and is more stable, the merits of a method depend on the applicable scenario. There is no absolutely optimal sampling method. Therefore, in practice, the three methods will be analyzed separately and the best method will be selected in real time.

[0074] This project involves 3D reconstruction and hole filling optimization, constructing a connectivity model based on a weighted calculation of geometric and roughness features. The task requires first performing 3D reconstruction to demonstrate the connectivity between six pillars, then identifying high-uncertainty regions in space based on this connectivity, and finally providing suggestions for optimal hole placement.

[0075] The first step still requires crack detection on different sections of different pillars at each borehole location. Therefore, the model from step 1 is used for preliminary detection, followed by filtering out interference factors using the image processing methods from the previous two questions. Then, a spatial coordinate system is established according to the description in the question. After obtaining complete crack data, it is necessary to solve for the various physical features in the two-dimensional plane and then transform them into three-dimensional features and the location of cracks in three-dimensional space. Therefore, three-dimensional parameter inversion and coordinate unification (from two-dimensional to global three-dimensional) are required. The core task is to integrate the second stage of the multi-stage iterative algorithm with the core steps of cross-borehole analysis.

[0076] Based on the two-dimensional fracture parameters and sinusoidal model results, the three-dimensional orientation of each fracture is inferred. The borehole circumference is verified by the fracture period, the inclination angle between the fracture plane and the borehole axis is calculated by the amplitude, and the position of the fracture plane in the local coordinate system is determined by the phase and centerline position. Finally, the normal vector and spatial position of the fracture plane are obtained. Then, combined with the effective observation length of the borehole, the infinitely extended fracture plane is trimmed into finite three-dimensional patches. The center position, scale (range along depth and circumference), normal vector direction, and JRC value of the patches are recorded. Then, it is checked whether there are spatial conflicts between fracture patches of different boreholes (such as patches that overlap but are not connected). Parameters that exceed the physical reasonable range are corrected (such as the normal vector direction of fracture planes with an inclination angle greater than 90 degrees needs to be recalculated).

[0077] To visualize the 3D patches, connectivity calculation and probabilistic network construction (from patches to connected clusters) are performed first. Connectivity indices are calculated by enumerating all fracture patch pairs across boreholes and calculating consistency indices from both geometric and textural dimensions. The geometric dimension includes the similarity of the normal angle and computational direction (measuring differences in plane orientation; smaller angles indicate higher consistency), patch spacing and computational space distance (measuring spatial distance; smaller distances suggest greater connectivity), and fracture size (larger sizes increase the likelihood of connectivity with other fractures). The textural dimension uses the JRC difference (measuring roughness consistency; smaller differences suggest the same fracture). All indices are then normalized to the 0-1 range.

[0078] Connectivity probability calculation and network construction: Consistency score is calculated with geometric weight of 0.7 and texture weight of 0.3. Among the set weights, directional similarity accounts for 25%, spatial distance accounts for 25%, and gap size accounts for 20%. The score is converted into connectivity probability using the sigmoid function. A weighted undirected graph is constructed with gap patches as nodes and connectivity probability as edge weights. An initial threshold is set (0.6-0.8 for ground connectivity, and above 0.8 for strong connectivity). Edges with a probability greater than or equal to the threshold are retained as valid connectivity relationships.

[0079] A breadth-first search is used to group interconnected patches into initial clusters. The representative plane of each cluster is fitted by total least squares, the residuals from the patches in the cluster to the representative plane are calculated, and outlier patches with excessive residuals are removed. The clusters are then rated for reliability based on the number of patches in the cluster, the fitted residuals, and the average connectivity probability.

[0080] Finally, end-to-end uncertainty calculations are performed, propagating uncertainty layer by layer from local parameters to the global model. Based on the fitting error of the two-dimensional parameters, the error range of the three-dimensional patch parameters (normal vector, center position) is calculated using derivative relationships. The uncertainty of the plane representing the cluster is aggregated according to the patch reliability weights to obtain the confidence interval of the cluster position and the standard deviation of the connectivity probability. Uniform sampling is performed on a grid within the analysis domain, and the uncertainty value of each sampling point is calculated to generate an uncertainty distribution map (the darker the color, the higher the uncertainty). Finally, based on the uncertainty distribution, high uncertainty regions are identified, and combined with engineering constraints, the optimal location for supplementary boreholes is determined.

[0081] In embodiments of the present invention, such as Figure 14 As shown, the effect of 3D connectivity reconstruction of the fractures in the six columns of the selected well is demonstrated. Based on the connectivity model solving algorithm constructed above, the fracture information in the six boreholes is accurately extracted, the three-dimensional spatial parameters of the fractures are calculated, the connectivity probability of fractures between adjacent boreholes is evaluated, high uncertainty areas are identified, and finally the optimal location of the supplementary borehole is determined, forming a complete analysis loop.

[0082] The fracture parameter statistics and spatial distribution characteristics clearly show that 49 major fractures were identified in the six boreholes. The spatial distribution of these fractures exhibits significant heterogeneity, with the following specific characteristics: In terms of depth distribution, there are fewer fractures in deeper and very shallower areas. This distribution pattern suggests that the central part of the rock mass in the study area may contain more structural planes or weak zones. These areas have relatively poor rock mass integrity and are key areas requiring attention during engineering.

[0083] In embodiments of the present invention, such as Figure 15As shown, the closer the color is to cool tones such as red and purple, the higher the spatial uncertainty of the area; the closer the color is to warm tones such as yellow and orange, the lower the spatial uncertainty of the area. The overall uncertainty in the shallow rock mass is at a moderate level. The color distribution shows a certain degree of color gradation, indicating that the shallow rock mass is affected by surface weathering and unloading, resulting in a certain complexity in fracture development. However, due to its proximity to the surface, there may have been relatively more indirect observation methods in the past, and existing borehole data also provides some coverage of this area, preventing the uncertainty from reaching a very high level, but it is not very low either, remaining in an intermediate transitional state. The spatial differences in uncertainty in the central region are more significant. Areas closer to existing boreholes have colors that lean towards low uncertainty (such as yellow), indicating that these areas have relatively good information about fracture distribution and connectivity due to existing borehole data, resulting in lower uncertainty. Conversely, areas farther from existing boreholes have colors that lean towards high uncertainty (such as red and purple), indicating a significant increase in uncertainty due to a lack of direct observation data, a large cognitive bias in the fracture network, and other factors. The central rock mass is influenced by factors such as tectonic stress, resulting in fractures with specific orientations and connectivity patterns. Far from the borehole, existing data struggles to accurately characterize these features, leading to increased uncertainty. Overall, uncertainty is high in deeper areas. Color-coded regions show higher uncertainty. This is due, in part, to the greater difficulty of exploring deep rock masses, resulting in reduced data density from existing borehole observations and insufficient direct observation. Furthermore, the deeper geological environment is more complex, potentially exhibiting different fracture formation mechanisms compared to the shallower and central parts. Existing models based on shallower and central data are difficult to accurately extrapolate to deeper areas, further increasing uncertainty and reflecting a significant rise in spatial uncertainty at depth. The optimal supplementary points obtained by the model are (1000, 1500), (2000, 1500), and (1500, 500). In fact, according to the uncertainty heat map and the fracture connectivity diagram of each layer, the first two points are located in the middle area of ​​the six pillars. Moreover, there is a high degree of uncertainty in the vicinity of these two areas at each depth. Therefore, it is the most suitable and cost-effective choice to build new pillars at these two points. The selection of the last point may be due to the high uncertainty in the vicinity of this point in the deeper area.

[0084] The technical solution provided by this invention includes a method that involves collecting and preprocessing a dataset, identifying fractures based on an improved U-Net network, and outputting a binary fracture map; performing quantitative analysis and modeling of sinusoidal fractures based on the binary fracture map; calculating the roughness coefficient JRC and performing sensitivity analysis of the sampling method based on the model; and performing three-dimensional reconstruction and borehole optimization of the multi-hole fracture network based on the JRC calculation and sensitivity analysis. This method improves the accuracy and reliability of fracture identification in surrounding rock and enhances its engineering applicability.

[0085] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for surrounding rock fracture identification and three-dimensional model reconstruction, characterized in that, The method comprises: Step 1, collecting data sets for pretreatment, and identifying cracks based on an improved U-Net network to output a crack binary graph; Step 2, based on the crack binary graph, performing sinusoidal crack quantitative analysis modeling; Step 3, according to the modeling, performing roughness coefficient JRC calculation and sampling method sensitivity analysis; Step 4, according to the JRC calculation and sensitivity analysis, performing multi-drill hole crack network three-dimensional reconstruction and hole optimization.

2. The method of claim 1, wherein, The step 1 comprises: Through image preprocessing, multi-method fusion comparison, network structure optimization and training, pixel-level separation of cracks and interference in the drill hole image is realized; Step 11, data acquisition and pretreatment, obtaining the data set of the drill hole image, performing gray scale, Gaussian denoising and histogram equalization on the image; through random cropping, the image is randomly cropped into a 244x244 pixel sample, and then rotation and mosaic enhancement are performed to expand the data set sample; Step 12, multi-method fusion comparison, five methods are designed: a, multi-method edge feature extraction and voting mechanism; b, multi-dimensional feature extraction and weighted fusion combined with adaptive threshold for crack detection; c, using a machine learning algorithm, based on multi-dimensional feature extraction and random forest algorithm to classify crack pixels; d, using multi-scale feature extraction and KMeans clustering algorithm to cluster crack pixels; e, using a deep learning algorithm for segmentation, using improved U-Net, U2Net and YOLOv11 respectively for training; wherein the improved U-Net network comprises an encoder, a decoder, a multi-scale fusion layer and an output layer; Step 13, model training and evaluation, the results are displayed in the form of result visualization, and the crack extraction result is selected; first, the data set is divided into training set, validation set and test set according to the ratio of 7:2:1, the intersection over union IoU and the average precision mean mAP are evaluated, and the crack binary graph is output.

3. The method of claim 2, wherein, It comprises: Step 111, gray scale, the drill hole image is in RGB format, and the weighted average method Y=0.299R+0.587G+0.114B is used to convert it into a gray scale image; Gaussian denoising, 5x5 Gaussian kernel convolution operation is adopted; histogram equalization, using the contrast limited adaptive histogram equalization CLAHE algorithm, the image is divided into 8x8 grids, and histogram equalization is performed on each grid separately, and the contrast is limited to improve; through random cropping, the image is randomly cropped into a 244x244 pixel area to simulate different imaging angles; through rotation, the sample is randomly rotated by 0-360 degrees; mosaic enhancement, four different samples are spliced according to 2x2, and the edges are linearly transitioned; Step 112, through multi-directional gradient Sobel edge detection, Canny edge detection, LoG edge detection combined with Gaussian blur and Laplacian operator edge detection, then using voting method, through combination judgment, that is, at least two detection results are cracks, then it is judged as a crack; Step 113, after Sobel edge detection, LBP texture detection and Hessian feature extraction, combined with adaptive threshold, the weights of the three are 50% for multi-scale, 30% for Soble edge detection and 20% for Hessian, forming multi-scale and multi-directional features, and finally through Otsu automatic threshold, the threshold value of the minimum intra-class variance of foreground and background is automatically calculated to obtain the final detection result; Step 114, supervised learning based on random forest successively experiences feature extraction, pseudo-label generation, forest training and confidence filtering; the diversity of feature extraction extracts multi-dimensional features of original intensity, gradient amplitude, LBP texture, Canny edge and Hessian matrix feature, and randomly selects part of the features to train each tree to mine the discriminative nature of different feature combinations; pseudo-label generation generates pseudo-labels through mathematical and statistical methods, converts unsupervised or weakly supervised data into labeled data, and provides a training basis for supervised learning of random forest; the parameter design of forest training, n_estimators=50 represents the number of trees as 50, which is used to control the size of the ensemble; max_depth=8 represents the maximum depth of the tree as 8, which is used to limit the complexity of a single tree and balance bias and variance; confidence evaluation based on gradient and standard deviation, screening the results predicted by random forest, gradient is used to reflect the stability of prediction, and standard deviation is used to reflect the divergence between trees; finally, the classifier learns the mapping relationship between features and crack or non-crack labels to classify new data; Step 115, multi-scale feature extraction obtains comprehensive features through multi-window size and structure tensor; Basic statistical features, calculating the mean, variance, standard deviation and contrast of each window, are used to capture the local gray distribution rule; structure tensor features, through gradient calculation by Sobel operator, construct structure tensor after Gaussian smoothing, calculate determinant, trace and coherence, and distinguish edge and uniform area; K-means coarse clustering unsupervised clustering algorithm, setting the clustering number K=3, quickly divides multi-dimensional features into three categories of background, interference and crack, the core is to minimize the sum of distances from samples in the cluster to the cluster center to complete the preliminary separation of target and non-target; DBSCAN refinement based on density-based clustering algorithm, optimizing the K-means result; Statistical feature selection calculates the intensity statistical features of the mean and variance of each cluster, and locks the cluster meeting the crack feature through threshold setting or statistical test; Finally, based on the intensity statistical feature, the crack cluster is automatically identified, the pixels belonging to the same crack are aggregated into a cluster, and the recognition result is obtained. Step 116, the left side of the U-Net network extracts multi-scale features through the encoder step by step down-sampling, and the right side restores the spatial resolution through the decoder step by step up-sampling, and the high-resolution detail features of the encoder are fused with the semantic features of the decoder through the jump connection; U2Net is an improved model based on UNet, which improves the segmentation effect by increasing multiple branches and attention mechanism; the architecture of U2-Net is based on RSU module, and RSU module is used in each En_x of U2-Net, which combines residual block and U-Net symmetric encoder-decoder structure to retain feature information through jump connection; Through the introduction of ResNet residual module and multi-scale feature fusion for improvement: For the encoder part, there are 4 down-sampling modules, each module contains 2 ResNet residual blocks and 1 2x2 max pooling, wherein the ResNet residual block, the input feature map is subjected to 3x3 convolution, batch normalization BN layer, ReLU activation, 3x3 convolution, BN layer; shortcut connection directly stacks the input feature map to the output end of the second BN layer, and its formula is F(x)=Conv2(ReLU(BN(Conv1(x))))+x; Max pooling, step 2, channel number doubled, gradually expanding the receptive field, used to capture global features of the fracture; For the decoder part, there are 4 up-sampling modules, each module contains 1 2x2 transpose convolution and 1 feature splicing layer; transpose convolution recovers the feature map size through zero padding and convolution operation, and its formula is Output=ConvTranspose2; Feature Splicing module splices the current feature map of the decoder with the same size feature map of the encoder along the channel dimension, and unifies the channel number to 256 through 1x1 convolution, and retains the high-resolution details of the encoder; For the multi-scale fusion layer part, the sampling features of each layer of the decoder are fused, the channel number is unified to 128 through 1x1 convolution, and the elements are added after ReLU activation: low down-sampling features retain small fracture edge details, and high down-sampling features provide global context; For the output layer part, 1x1 convolution compresses 256-channel features to 1 channel, and Sigmoid activation function outputs pixel-level probability, and when the probability>0.5, it is judged as a fracture; then the model is trained and evaluated, the loss function selects cross-entropy loss function, and the weight factor is introduced to ensure the attention of the model to small fractures; the optimizer setting selects Adam optimizer, which combines momentum gradient descent and adaptive learning rate, the initial learning rate is 0.001, it is attenuated to 1 / 10 every 50 rounds, and the weight attenuation is 1e-5; the evaluation index is intersection over union IoU, which measures the overlap degree of the predicted region and the real region, and the average precision mean mAP is calculated under the IoU threshold of 0.5-0.95, which is used to reflect the recognition robustness of the model to different size fractures, and when the average precision mean mAP is greater than or equal to 92%, it is qualified; finally, the fracture binary graph is output after multiple rounds of training to meet the index requirements.

4. The method of claim 3, wherein, The step 2 comprises: The parameterized fitting of sinusoidal cracks is realized by multi-sampling fusion and nonlinear optimization. Firstly, the improved U-Net network is used to obtain the crack binary graph and perform preprocessing, then the coordinates are standardized, and then the center line is obtained by sampling in multiple ways, and finally the various parameters of the sinusoidal crack are obtained; based on the crack binary graph, the technical process of coordinate standardization, geometric filtering and multi-sampling fitting is constructed: through pixel-millimeter mapping to unify the dimension, the minimum enclosing circle screening method is used to remove non-sinusoidal interference; a weighted fusion sampling strategy is designed for equidistant sampling, adaptive sampling and curvature key point sampling, and the fitting of the sinusoidal model is realized by combining the Levenberg-Marquardt algorithm and the iterative reweighted least squares IRLS; Step 21, crack preprocessing, the improved U-Net network is used to segment the data set, and the segmented results are preliminarily preprocessed; the image coordinates are converted into physical dimension through pixel-millimeter mapping, the skeleton of the crack is equidistantly sampled, and the same horizontal coordinate is de-duplicated and averaged; Step 22, multi-mode sampling, equidistant sampling, take points along the center line x direction 5 pixel step length to ensure uniform distribution; adaptive sampling, calculate the curvature k, when k>0.1mm -1 , the step length is halved to 2.5 pixels; curvature key point sampling, calculate the second derivative, keep the points whose absolute value is greater than 70% quantile; finally, weighted fusion is carried out, and then the same horizontal coordinate is de-duplicated to obtain the center line of the fracture; Step 23, parameter fitting, de-duplication and averaging of sampling points, global search of initial parameters of the sine model y = Rsin(2πx / P + β) + C through differential evolution algorithm; combined with weighted total least squares, output amplitude R, period P, phase β, center line position C, and verified by determination coefficient R 2 and root mean square error RMSE.

5. The method of claim 4, wherein, It includes: Step 211, coordinate standardization, first, calculate the pixel-millimeter conversion factor under the constraint of unified aperture and developed aperture; second, coordinate convention, specify the horizontal direction as the circumferential length and the vertical direction as the axial depth, record the image size (W, H), resolution, period P, and keep the starting depth y0 of each image; finally, apply unified dimension to sampling, fitting and clustering to ensure that the parameters of different images have the same physical units; Step 212, geometric filtering, coordinate mapping is first used to convert image pixel coordinates into real physical millimeter coordinates; Set a length threshold Lth, calculate the real millimeter length L of each candidate line segment, and keep the line segments with L≥Lth and remove the noise or small impurity line segments with L<Lth; Geometric shape constraint selects line segments that meet the shape of sinusoidal cracks through two sub-constraints; set the slenderness ratio threshold Rth to retain line segments with a length-to-width ratio Rslenderness≥Rth; To the continuity constraint, the line segment direction angle θ is calculated, and for the adjacent line segments with the distance between the endpoints < 2mm, the angle difference is calculated θ, set the angle difference threshold θ, remove the noise with chaotic direction Step 213, parameter estimation, estimation of the center line position C, the mean value of the areas of the upper half and the lower half of the sinusoidal curve in one complete period is approximately equal to the center line position C; let the data set be Then the mean value The calculation formula is: The initial estimation of the center position is: ; The estimation of the amplitude R is that the range of the sinusoidal curve is approximately equal to 2R, and the range of the data set is The formula for calculating the range is: ; An initial estimate of the amplitude is: ; Period P estimation, period P is used to describe the length of a complete wave of a sinusoidal curve, and frequency f is the number of fluctuations per unit length, which satisfies the inverse relationship P=1 / f; fast Fourier transform FFT converts time domain data into frequency domain data to extract the periodic component of the data, and the steps are as follows: First, in order to focus on the fluctuation part of the data, the data is de-meaned: ； Secondly, the data after mean removal is subjected to FFT to obtain a complex array in frequency domain: ； wherein denotes the amplitude and phase of the component with frequency ; the frequency is related to the array index k by: ； wherein, is the sampling frequency, i.e. the number of sampling points per unit length, if x is equidistantly sampled and the interval between adjacent x is then ; Then find the dominant frequency, in the frequency domain, the component with the largest amplitude corresponds to the dominant periodicity of the data; calculate the frequency corresponding to the amplitude , find the frequency corresponding to the maximum amplitude, which is the dominant frequency ; Finally, according to the inverse relationship between period and frequency, the estimated value of the period is: ； Step 214, equidistant sampling, uniformly generate sampling points in a specified range; adaptive sampling, dynamically adjust the sampling density according to the curvature of the crack; curvature key point sampling, select the most dramatic key points by calculating the second derivative for sampling; Step 215, de-duplication and mean value, remove repeated points under the same abscissa x, take the mean value of the remaining multiple y values to obtain a unique y value corresponding to the current abscissa x; global search and local refinement, first stage: differential evolution, that is, global search for approximate solution, as a global optimization algorithm based on population, through three steps of iteration search of mutation, crossover and selection; mutation, for each individual xi in the population, use the difference vector of other random individuals to generate a mutation vector vi; crossover, cross the mutation vector vi with the original individual xi to generate a test vector ui; selection, compare the residuals of ui and xi, and retain the individual with smaller residual to enter the next generation, and finally obtain the globally approximate optimal candidate parameters; the second stage: weighted total least squares method, that is, local refinement for accuracy, on the basis of global search for approximate solution, further optimization is carried out for the inner point set, and the core is weighting, that is, different data points are assigned weights; the target is to minimize the weighted sum of squares of residuals, and the optimal parameters are solved by taking the partial derivative of the objective function, and finally the high-precision parameter estimation result is obtained; After the fracture acquisition and preprocessing, initial parameter estimation, three sampling methods sampling and weighted fusion, then the sampling point is removed to get the center line, finally the global search and local refinement, result evaluation and visualization calculation, goodness of fit index R 2 With the root mean square error RMSE, the fitting quality is evaluated, and for the fractures with poor fitting quality, the clustering and center line extraction process is rechecked and adjusted; the fitted sinusoidal curve is compared with the original fracture, and the effect is intuitively tested; the parameter estimation values of each fracture are summarized to form a quantitative analysis result table.

6. The method of claim 5, wherein, The step 3 comprises: The macroscopic trend and the microcosmic rough residual of the fracture are separated by using the smooth spline function. The influences of the equidistance sampling, equi-arc length sampling and multi-scale adaptive sampling on the JRC value and JRC value are compared. The various geometric features of the fracture are extracted, decoupled, and sampled. The JRC value is calculated according to the sampling points, and the roughness is finally obtained. The macroscopic trend and the microcosmic rough residual of the fracture are separated by using the smooth spline function. The influences of the equidistance sampling, equi-arc length sampling and multi-scale adaptive sampling on the JRC value and JRC value are compared. The various geometric features of the fracture are extracted, decoupled, and sampled. The JRC value is calculated according to the sampling points, and the roughness is finally obtained. Step 31, contour extraction and trend separation, first, crack extraction and pretreatment are carried out, the improved U-Net network is used for crack extraction, and then trend separation and residual extraction are carried out; after the crack skeleton is reparameterized by arc length, the trend line is fitted by using a smoothing spline or a low-order polynomial, and the high-frequency rough residual is obtained by the difference between the original curve and the trend line; Step 32, multi-mode sampling calculates roughness index, compares the differences between equal-interval sampling, equal-arc-length sampling and multi-scale adaptive sampling, adopts encryption sampling in high-curvature area, and calculates length-normalized first-order slope energy index based on residual data ; Step 33, JRC inversion and uncertainty evaluation Introduce endpoint missing mask, complete JRC mapping and sensitivity analysis, convert the value to JRC index according to empirical formula, use Bootstrap method to quantify the influence of sampling interval, window size and detrending strength on the results, output the mean and 95% confidence interval of JRC. Step 33, JRC inversion and uncertainty evaluation Introduce endpoint missing mask, complete JRC mapping and sensitivity analysis, convert the value to JRC index according to empirical formula, use Bootstrap method to quantify the influence of sampling interval, window size and detrending strength on the results, output the mean and 95% confidence interval of JRC.

7. The method of claim 6, wherein, Comprise: Step 311, feature decoupling separates the low-frequency shape and high-frequency rough features of the crack by a smoothing spline method; The core is to minimize the joint objective function of data fitting error and curve smoothness to obtain a smooth curve reflecting the overall trend, and then subtract the smooth curve from the original data to obtain a residual sequence containing only high-frequency rough information; Step 312, select the adaptive sampling method for different scenes, calculate the roughness index ; equidistant sampling, sampling by equal interval on x axis, focusing on horizontal comparability, calculated by the height difference of adjacent points and horizontal step , suitable for horizontal feature comparison, its expression is: ； Equal arc length sampling: Convert the parameters to cumulative arc length, the expression of roughness index under equal arc length sampling is: ​ ； Multi-scale adaptive sampling: combining multi-scale analysis and curvature to dynamically adjust the sampling density, dense sampling in complex areas and sparse sampling in smooth areas to balance feature capture accuracy and computational efficiency; the sampling point set obtained by multi-scale adaptive sampling is represented as: ； Curvature weighted optimization: Assign high weight to areas with high curvature, modify Indicators, enhance sensitivity to irregular features of the crack tip Step 313, JRC inversion utilization The mapping relationship with JRC's engineering experience will be used to map different sampling methods and scales. The values ​​are converted into dimensionless JRC values ​​to achieve a quantitative conversion from geometric features to roughness coefficients; In engineering, the formula is mapped to JRC by an empirical formula using a logarithmic calibration expression: ； where JRC(w) represents joint roughness coefficient in window w, which is a dimensionless empirical rating; represents calibration data, which is related to specimen size, gage length, and unit consistency; represents roughness index value in window w; Step 314, Bootstrap uncertainty evaluation, through multiple resampling propagation of discrete point selection and the influence of noise factors, the standard error and confidence interval of JRC are calculated, the result reliability is quantified, and finally the JRC value with reliability is output, and its formula is: ； where, denotes the bth resampled data, denotes detrending, , the chain of JRCs, B denotes the number of resamplings, denotes the mean of the multiple resamplings, SE denotes the standard error, denotes the empirical quantile of the 1 -a confidence interval.

8. The method of claim 7, wherein, The step 4 comprises: Establish a technical framework of coordinate mapping, connectivity evaluation, uncertainty optimization, three-dimensional (3D) reconstruction and supplementary point generation: based on the coordinates of the borehole opening or the bottom, the two-dimensional crack parameter is inversely transformed into the crack surface patch in the three-dimensional space through Lie algebra SE(3) pose transformation, the multidimensional consistency index is constructed by fusing the normal angle, spatial distance, crack size and JRC difference, and the connectivity probability in the 0-1 interval is mapped through the sigmoid function; the breadth-first search algorithm is used to extract the crack cluster, the least square algorithm is used to fit the representative plane of the crack in the cluster, and the high-uncertainty area is located by calculating the crack density variance and connectivity probability entropy of each point in the space; Step 41, three-dimensional coordinate mapping and fracture extraction, taking the borehole opening as the origin (x=0, y=0, z=0), x represents the parallel tunnel direction, y represents the vertical tunnel wall, and z represents the depth; a homogeneous transformation matrix T_BH is constructed through the borehole pose to map the local fracture parameters to the global, the single-hole fracture parameters are converted to the global coordinate system based on the opening or bottom coordinates, the corresponding three-dimensional plane of the fracture is inversely calculated through SE(3) pose transformation, and the fracture is cut into a limited face sheet according to the observed length; for complex fractures, Fourier series or spline surface expansion models are used; the fracture is extracted by improving the U-Net network, and the fracture is preprocessed to filter out interference information in the fracture; Step 42, connectivity analysis, construct multi-dimensional consistency index, geometric level adopts normal angle, plane constant difference and Hausdorff distance and other characteristics, texture level introduces JRC difference value, and designs weight to de-weight two characteristics; the comprehensive score is mapped to the connectivity probability in the interval of 0-1 through logistic regression, a weighted undirected graph is constructed, and connected fracture clusters are extracted, and the total least squares method is used to jointly fit the representative plane in the cluster; Step 43, uncertainty and hole optimization, the fracture density variance and connectivity probability entropy of each point in the space are calculated to locate the high uncertainty area, and the candidate hole area is selected combined with engineering constraints; hole position optimization, traverse the space coordinates, use the information gain maximization algorithm to evaluate the uncertainty reduction benefit of each candidate point, select the top three optimal positions, finally present the fracture network, connectivity probability and hole scheme through three-dimensional visualization, form the closed loop of identification-modeling-decision.

9. The method of claim 8, wherein, Including: Step 411, data docking and coordinate unification, establish the rigid transformation relationship between the local coordinate system and the global coordinate system of the borehole, unify different borehole data to the same spatial system, and ensure the consistency of cross-hole analysis; at the same time, integrate the core data of the sine model parameters, JRC and hole radius, and provide a unified quantitative basis for three-dimensional reconstruction; Step 412, 2D to 3D inversion and face sheet construction, based on the sine model parameters of the borehole development graph, the normal vector and plane equation of the fracture plane are inversely calculated, and are mapped to the global coordinate system through pose transformation; combined with the effective observation length of the borehole, the limited fracture face sheet is cut out, and the metadata of the center and scale are recorded, and the transformation from two-dimensional measurement to three-dimensional morphology is completed; Step 413: Cross-hole connectivity determination. Using geometric and texture multi-dimensional indicators, the connectivity probability of different borehole fractures is quantified. Geometric consistency: Calculate the normal angle, planar position difference, and spatial distance, and combine attitude similarity and spatial coplanarity to determine the fit between the fracture's spatial orientation and position. Texture corroboration: Introduce JRC difference to supplement roughness consistency. Perform probability mapping, obtain a comprehensive score through weighted fusion of multi-dimensional indicators, and then convert it into connectivity probability using the Sigmoid function to form a symmetric probability matrix. Obtain a comprehensive similarity score by weighted fusion of geometric and roughness indicators. Comprehensive similarity... The formula is: ； wherein, is a weight, is a scale constant for normalizing each index, is the roughness inversion of the corresponding facet intersection window; Logistic regression mapping of connectivity probability, converting the comprehensive similarity score into connectivity probability for threshold judgment and graph model calculation; Sigmoid function is used to realize the monotone mapping from scoring to probability, and the probability that two pieces of fractures belong to the same continuous fracture body is obtained: ； wherein, is the probability that patches i and j belong to the same continuous fracture body, and a and b are calibration coefficients determined by the training data; Step 414, probability graph construction and cluster extraction, taking the fracture face sheet as the node and the connectivity probability as the edge weight, a weighted undirected graph is constructed; set the probability threshold to filter the effective edges, and get the interrelated fracture body clusters through the connected component extraction or community discovery algorithm, and then use the total least squares method to fit the overall plane or surface of the cluster to represent its geometric characteristics; the connectivity relationship of the fracture face sheet is converted into the structure of the graph, and the connected component extraction or community discovery algorithm in graph theory is used to identify the closely related fracture face sheet set, i.e. the fracture body cluster: ； wherein G is a graph, is a set of nodes, is a set of edges; Step 415, uncertainty propagation and aggregation propagates the single-parameter uncertainty to the position and pose parameters of the three-dimensional facet through the Jacobian matrix, and aggregates the intra-cluster parameter covariance to quantify the uncertainty of the fracture cluster and the connectivity probability; Step 416, hole optimization realizes efficient hole optimization under the engineering reachable domain and construction constraints; the constraint range is drawn, the reachable domain is defined combined with the construction conditions, the candidate hole set that is feasible in engineering is generated based on the existing drilling distribution; the high-value area is locked, the cluster plane uncertainty and the connectivity probability are fused to identify the high-uncertainty-high-value hole hot area; the optimization target is selected, and the expected information gain of the candidate hole is taken as the objective function to select the optimal hole combination.

Images