A method for reconstructing three-dimensional connectivity of fractures in surrounding rock in a well based on a CS-HAT algorithm

Abstract

The application discloses a kind of based on CS-HAT algorithm's underground surrounding rock fissure three-dimensional communication reconstruction method, belong to coal rock engineering monitoring and spatial structure reconstruction analysis technical field.This method takes binary fissure image as input, carries out automatic clustering and parameterized modeling to the fissure section of sinusoidal distribution in borehole imaging development diagram, obtains fissure geometric parameter;Fissure profile carries out physical scale conversion, sampling and roughness evaluation, constructs geometric roughness index and obtains JRC value;Two-dimensional fissure point is projected back to real cylindrical hole wall space, constructs three-dimensional point cloud and carries out plane fitting, in combination with spatial geometric relationship and roughness factor establishes fissure communication probability discrimination model, realizes the three-dimensional communication reconstruction of multiple borehole fissure network;Finally, the position and drilling depth of supplementary borehole are optimized.The application can improve the reliability of underground surrounding rock fissure communication discrimination and three-dimensional reconstruction accuracy, provide technical basis for underground roadway surrounding rock stability analysis and disaster prevention.

Description

Technical Field

[0001] This invention relates to the field of spatial structure reconstruction and analysis technology for coal and rock engineering monitoring, and in particular to a method for three-dimensional connectivity reconstruction of underground surrounding rock fractures based on the CS-HAT algorithm. Background Technology

[0002] As coal mining depths increase, the surrounding rock of underground roadways is subjected to high ground stress, high disturbance, and complex hydrogeological conditions for extended periods. This results in numerous well-hidden, fine-scale, and complexly distributed fracture structures within the rock. The spatial distribution, connectivity, and roughness characteristics of these fractures not only affect the stability of the roadway's surrounding rock but also, under the influence of support, construction, and changes in hydrological conditions, can easily trigger safety accidents such as roof collapse, water inrush, and gas outbursts. Therefore, three-dimensional connectivity analysis of underground rock fractures is of great significance for preventing and controlling related disasters.

[0003] Currently, borehole imaging technology can obtain two-dimensional distribution information of surrounding rock fractures in borehole unfolding diagrams and further acquire fracture contours and their geometric features. However, existing methods mostly remain at the level of two-dimensional fracture identification and local parameter extraction. There is a lack of unified and effective processing methods for automatic clustering of sinusoidal fractures, extraction of key geometric parameters, evaluation of fracture roughness, and spatial connectivity between fractures under multi-bore conditions, making it difficult to establish a reliable and stable three-dimensional fracture network model.

[0004] In addition, factors such as fracture roughness, geometric continuity, hole spacing, and spatial coverage can all affect the fracture connectivity analysis results in actual engineering. Relying solely on manual experience and simple geometric splicing methods can easily lead to large errors in connectivity judgment and insufficient accuracy in three-dimensional reconstruction. Under the condition of a limited number of holes, the existing borehole information does not cover high-risk areas completely, making it difficult to accurately reflect the true distribution of the fracture network.

[0005] To address existing shortcomings, this invention proposes a three-dimensional connectivity reconstruction method for downhole surrounding rock fractures based on the CS-HAT algorithm. By automatically clustering, evaluating roughness, performing three-dimensional backcasting, and analyzing connectivity probability of two-dimensional fracture images, the method achieves three-dimensional connectivity reconstruction of the downhole surrounding rock fracture network. Furthermore, it optimizes supplementary boreholes to improve the reliability of fracture connectivity identification and the accuracy of three-dimensional reconstruction. Summary of the Invention

[0006] To address the shortcomings of the existing technologies, this invention provides a method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm. The overall process of this model includes the following four steps.

[0007] Step 1: Using the binarized fracture image as input, perform automatic clustering and quantitative modeling of fracture structures that are distributed in a "sine-like" pattern in the borehole imaging unfolded image. Through processes such as connected region marking, composite distance merging, skeleton extraction, trunk smoothing, and sine curve fitting, key geometric parameters such as the amplitude, period, phase, and centerline position of the fracture are automatically obtained.

[0008] Furthermore, the specific steps for automatic clustering and quantitative modeling of "sinusoidal" cracks in step 1 above are as follows.

[0009] Step 1.1: Let the obtained binarized crack image be I. b (x,y), where (x,y) are the image pixel coordinates, x∈[0,W]. i -1],y∈[0,H i -1] When pixel value I b A fracture is defined when (x,y)=1. b When (x,y)=0, it is defined as a non-crack. The input binary crack image is initially segmented by the 8-connected region labeling algorithm, as shown in equation (1).

[0010] (1)

[0011] In the formula, C i This represents the nth connected region; This represents an 8-connected neighborhood.

[0012] Furthermore, considering that actual cracks may be divided into multiple disconnected regions due to noise or interruptions, two connected regions C are defined. i and C j The composite distance function between them is shown in equation (2).

[0013] (2)

[0014] In the formula, ||pq|| is the Euclidean distance between pixels p and q; , For region C i and C j The mean of the ordinate center; λ c This is a weighting coefficient used to balance the effects of spatial distance and longitudinal position differences.

[0015] Furthermore, when D c (C i C j )≤T m At that time, it is assumed that two regions belong to the same fracture, and after merging, the final set of fracture regions {F1, F2, ..., F} is obtained. m}

[0016] Step 1.2: For each clustered fracture region F k Calculate its distance transformation graph T k (x, y), its formula is shown in formula (3).

[0017] (3)

[0018] In the formula, 𝜕F k For region F k The set of boundary points.

[0019] Furthermore, after performing a distance transformation, the Zhang-Suen thinning algorithm is used to extract the skeleton line S of the fracture. k Its formula is shown in formula (4).

[0020] (4)

[0021] In the formula, ZS represents the Zhang-Suen refinement operator.

[0022] Furthermore, since the refined skeleton line may contain branches, the main trunk is extracted using the longest path algorithm, as shown in equation (5).

[0023] (5)

[0024] In the formula, LP represents the longest path extraction operator; P k This represents the main skeleton of the fissure.

[0025] Furthermore, the main skeleton line is subjected to moving average filtering to reduce local noise, as shown in equation (6).

[0026] (6)

[0027] In the formula, r s The radius of the sliding window; The main skeleton lines after smoothing.

[0028] Step 1.3: Use the mathematical expression of the "sinusoidal" fracture to fit the sine curve, as shown in equation (7).

[0029] (7)

[0030] In the formula, A S T represents the amplitude; S For periodicity; φ s For phase; Y S This is the position of the center line.

[0031] Furthermore, let the set of centerline points be {(xi ,y i )|i=1,2,…,N p}, define the fitting error function, which is shown in equation (8).

[0032] (8)

[0033] In the formula, E s N is the fitting error function; p This represents the number of points along the center line.

[0034] Furthermore, nonlinear optimization is performed using the Levenberg-Marquardt algorithm for parameter optimization, as shown in equation (9).

[0035] (9)

[0036] In the formula, θ S J is a parameter vector; s For Jacobian matrices; e t μ is the residual vector; s is the damping factor; I is the identity matrix.

[0037] Step 1.4: To verify the reliability and stability of the established clustering and sinusoidal fitting models, a quantitative effectiveness evaluation is conducted. Evaluation indicators are established based on detection rate, goodness of fit, root mean square error, stability, and robustness. Furthermore, a comprehensive performance scoring function is constructed to comprehensively measure the effectiveness of the model.

[0038] Furthermore, the detection rate is defined as the proportion of the number of components correctly identified as cracks to the total number of candidate components, as shown in equation (10).

[0039] (10)

[0040] In the formula, N c To satisfy the aspect ratio and coverage threshold, the number of crack sample components, N l This represents the number of connected components whose area exceeds the threshold.

[0041] Furthermore, the fitting effect is measured. Given a set of crack pixels {(x i ,y i )} N i=1 and the corresponding fitted values Define the goodness of fit R 2 The root mean square error (RMSE) is given by equation (11).

[0042] (11)

[0043] In the formula, These are the fitted values; This is the mean of the ordinates of all fracture points.

[0044] Furthermore, stability and robustness indices are introduced to measure the consistency of the model's performance and its ability to resist interference under different samples, as shown in equation (12).

[0045] (12)

[0046] In the formula, d This represents the standard deviation of the detection rate across different samples; and r S represents the mean and standard deviation of the root mean square error, respectively; b For stability indicators; R b This is a stability indicator.

[0047] Furthermore, a unified scoring function is introduced for comprehensive judgment, as shown in equation (13).

[0048] (13)

[0049] In the formula, The mean detection rate for all samples. For the average goodness of fit, Pr(Q) f >τ q ) represents the proportion of high-quality fits; τ q Goodness-of-fit threshold; Q s This is a comprehensive performance score.

[0050] Step 2: Based on the clustering and parameterization results, the physical scale of the crack profile is converted and sampled, and the crack roughness is calculated and evaluated. The pixels are uniformly converted into millimeters to ensure that the data of different cracks and different attachments can be compared under the same physical scale. The formula is shown in Equation (14).

[0051] (14)

[0052] In the formula, W i H is the image width. i L represents the image height. c L is the perimeter of the unfolded diagram. d The depth corresponding to the unfolded diagram; s x ,s y These represent the actual lengths of each pixel in the horizontal and vertical directions, respectively.

[0053] Furthermore, based on pixel conversion to obtain physical coordinates, the distance between adjacent points is calculated and accumulated segment by segment to obtain the total arc length, as shown in equation (15).

[0054] (15)

[0055] In the formula, x i ,y i Let i be the physical coordinates of the i-th crack point; The arc length between two adjacent points; Let be the cumulative arc length from the starting point to the j-th point.

[0056] Furthermore, a curvature adaptive method is introduced, which takes more points at turning points and fewer points on straight segments. The curvature magnitude is approximated by the second-order difference as the sampling weight, and combined with the arc length formula, a curvature-weighted adaptive sampling method is formed, as shown in equation (16).

[0057] (16)

[0058] In the formula, P i Let be the physical coordinates of the i-th point; к i This is an approximate value of the curvature at that point.

[0059] Furthermore, we need to establish a unified geometric roughness index, Z, to measure the overall volatility of the curve. r The essence of the index is the "mean square slope". It takes the average of the slopes of discrete points and takes the square root to obtain a value that represents the degree of undulation of the curve. It is directly related to the sampling method in the previous step and serves as the core input for the final JRC calculation. Its formula is shown in equation (17).

[0060] (17)

[0061] In the formula, N s ∆y is the number of sampling points; i The longitudinal difference between adjacent points; ∆x i Z represents the lateral difference between adjacent points. r This is the mean square slope index.

[0062] Furthermore, due to Z r It is a geometric quantity, but the JRC value is actually used in engineering practice, so Z needs to be converted using empirical formulas. r Convert to JRC, as shown in equation (18).

[0063] (18)

[0064] In the formula, b J c Jis the empirical fitting coefficient; JRC is the crack roughness coefficient.

[0065] Furthermore, the actual area of ​​the crack is calculated and standardized. The true area of ​​the crack is obtained by converting the pixel area with the pixel-to-millimeter ratio. Then, it is divided by the unfolded length to obtain the area per unit length. Combined with the JRC results, the effect of the crack area on roughness evaluation is analyzed, as shown in equation (19).

[0066] (19)

[0067] In the formula, N a A represents the total number of pixels in the crack; r L represents the actual area of ​​the crack. x A represents the horizontal length of the unfolded diagram. r / L x The area of ​​the crack per unit length.

[0068] Step 3: Project the two-dimensional fracture pixels in the borehole unfolding diagram back onto the actual cylindrical borehole wall space to construct a three-dimensional point cloud of fractures and perform planar fitting. Based on this, integrate indicators such as normal angle, geometric proximity, longitudinal continuity, borehole distance and roughness factor to establish a fracture connectivity probability discrimination model, thereby realizing the connectivity analysis and three-dimensional reconstruction of the multi-bore fracture network.

[0069] Furthermore, the specific steps of the fracture connectivity probability analysis and three-dimensional back projection in step 3 above are as follows.

[0070] Step 3.1: To restore the planar coordinates to the real space, the actual radius of the borehole must first be obtained, as shown in equation (20).

[0071] (20)

[0072] In the formula, L C R is the circumference of the hole. h This represents the actual radius of the borehole.

[0073] Step 3.2: Establish a pixel-to-millimeter mapping, and make the width and height of the image correspond to the actual perimeter and depth, so that the pixel crack points can be accurately converted into millimeter coordinates. After combining with the hole radius calculated in the previous step, spatial restoration is achieved, as shown in equation (21).

[0074] (twenty one)

[0075] In the formula, L s The actual depth corresponding to the vertical direction of each unfolded diagram; s x ,s y These represent the millimeter values ​​corresponding to the horizontal and vertical pixels, respectively.

[0076] Furthermore, the crack point is projected onto the cylindrical hole wall, as shown in equation (22).

[0077] (twenty two)

[0078] In the formula, x m ,z m The x coordinates are the millimeter coordinates after the unfolded diagram is converted. h ,y h R represents the position of the hole center in the planar layout. h (X,Y,Z) represents the radius of the hole; (X,Y,Z) represents the coordinates of the point in three-dimensional space.

[0079] Step 3.3: Use the least squares method to fit the crack point cloud to obtain the plane normal vector and constant term, as shown in equation (23).

[0080] (twenty three)

[0081] In the formula, n p The plane normal vector; x=(X,Y,Z) are the coordinates of the point; d p This is a plane constant term.

[0082] Furthermore, the RMSE and fitting formula are used to verify whether the fracture segment can be regarded as the same plane, as shown in equation (24).

[0083] (twenty four)

[0084] In the formula, P i N represents the coordinates of the fracture point; p n is the number of points; p d p These are the normal vector and constant term of the fitted plane, respectively.

[0085] Step 3.4: Introduce a roughness index to prevent the problem of reduced penetration probability due to excessive tortuosity, as shown in equation (25).

[0086] (25)

[0087] In the formula, z i z is the longitudinal coordinate of the fracture point; l,i Z represents the fitted linear trend; r This is the mean square of fluctuation index.

[0088] Step 3.5: Combine the weighted sum with the penalty factor to unify the aforementioned indicators to a probability value of 0–1, and obtain a comprehensive score for whether the two pores are connected, as shown in equation (26).

[0089] (26)

[0090] In the formula, S n For normality consistency scoring; S p Scoring for consistency in planar distance; S v For longitudinal continuity scoring; S h J represents the weight of the distance between holes; p ν is the roughness penalty factor; n ,ν p ,ν ν P represents the weighting coefficient. c The connectivity probability is scored.

[0091] Step 4: To further improve the completeness of fracture network identification and the reliability of connectivity discrimination, uncertainty analysis is performed on the study area and supplementary borehole optimization is carried out. By constructing comprehensive indicators such as coverage field, planar proximity and inter-hole distance, the current insufficient borehole area is determined, and the optimal supplementary borehole location and drilling depth are recommended using greedy point selection and dispersion constraint strategies.

[0092] Furthermore, the specific steps for uncertainty region analysis and supplementary borehole optimization in step 4 above are as follows.

[0093] Step 4.1: Use the Gaussian kernel function to spatially superimpose the center points of the existing fracture segments to obtain a coverage field. The larger the value, the denser the surrounding existing fracture information. The formula is shown in Equation (27).

[0094] (27)

[0095] In the formula, (X,Y) are the coordinates of a point on the planar grid; (x,Y) are the coordinates of a point on the grid. c,i ,y c,i ) represents the coordinates of the center point of the fracture segment; σ g Gaussian kernel scale; G min G max These are the minimum and maximum values ​​under normalization, respectively; ε g To prevent decimals with a denominator of zero.

[0096] Furthermore, the minimum distance from the point to the plane is calculated to measure the proximity of the hole-filling point to these planes, for the calculation of the uncertainty field, as shown in equation (28).

[0097] (28)

[0098] In the formula, (x,y,z) are the coordinates of the candidate points; (n xi ,n yi ,n zi ) represents the unit normal vector component of the i-th plane; d iThis is a plane constant term.

[0099] Furthermore, the shortest distance from the candidate point to all the arranged orifices is calculated, as shown in equation (29).

[0100] (29)

[0101] In the formula, (x,y) are the planar coordinates of the candidate point; (x j ,y j ) represents the plane coordinates of the j-th borehole.

[0102] Step 4.2: Combine the relevant indicators of the rationality of the hole filling point into a unified uncertainty index, as shown in equation (30).

[0103] (30)

[0104] In the formula, δ h (x,y) is the shortest distance from the point to the existing hole; δ p (x,y,z) is the distance from the point to the nearest crack plane; G c (x,y) represents the coverage field; q1,q2,q3 are weight coefficients.

[0105] Furthermore, a distance penalty is introduced between points to ensure a more uniform distribution of recommended boreholes, as shown in equation (31).

[0106] (31)

[0107] In the formula, Q represents the planar coordinates of the candidate point; P is the selected borehole point; ||Q−P|| is the Euclidean distance between the two points; R d To influence the scale of the scope; This represents the mean value of the workable area U.

[0108] Furthermore, an objective function is constructed to filter points. In each round of selection, the candidate points are scored, and the one with the highest score is selected as the new hole-filling point. This is closely combined with the penalty formula in the previous step, as shown in equation (32).

[0109] (32)

[0110] In the formula, U c (x,y) represents the uncertainty value of the candidate point; П(x,y) is the distance penalty function; F o (x,y) represents the overall objective value.

[0111] Step 4.3: Find the optimal depth within the longitudinal z range, as shown in equation (33).

[0112] (33)

[0113] In the formula, z ∗ (x,y) represents the optimal depth corresponding to the candidate point; Z c h represents the depth of the search window center. z The width of the window is half its width; L b Recommended drilling length.

[0114] The beneficial effects of adopting the above technical solution are as follows: This invention provides a method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm. This method, based on the obtained two-dimensional fracture results, can achieve automatic clustering of fracture segments, parameter extraction, quantitative evaluation of roughness, and three-dimensional connectivity reconstruction of multi-bore fracture networks, and complete supplementary borehole optimization, thereby improving the reliability of fracture connectivity discrimination and the accuracy of three-dimensional reconstruction. Attached Figure Description

[0115] Figure 1 This is a flowchart of the automatic clustering of cracks and extraction of geometric parameters in this invention.

[0116] Figure 2 This is a diagram showing the results of automatic clustering and numbering of image gaps in various embodiments of the present invention.

[0117] Figure 3 This is a flowchart of the crack profile sampling and roughness evaluation steps of the present invention.

[0118] Figure 4 This is a comparison chart of JRC results under different numbers of sampling points in an embodiment of the present invention.

[0119] Figure 5 This is a flowchart of the three-dimensional reconstruction steps of multi-hole fractures in this invention.

[0120] Figure 6 This is a planar diagram showing the crack penetration result in an embodiment of the present invention.

[0121] Figure 7 This is a diagram showing the optimization results of hole filling in an embodiment of the present invention. Detailed Implementation

[0122] The specific embodiments of the present invention will be described in further detail below with reference to the accompanying drawings and examples. The following examples are used to illustrate the present invention, but are not intended to limit the scope of the present invention.

[0123] This embodiment takes the three-dimensional connectivity reconstruction of fractures in downhole surrounding rock drilling imaging as an example, according to... Figure 1 , Figure 3 and Figure 5 The process shown sequentially completes automatic fracture clustering, parameter extraction, fracture roughness evaluation, three-dimensional connectivity reconstruction of fractures in multiple boreholes, and optimization of supplementary boreholes. The specific details are as follows.

[0124] Step 1: Using the binarized fracture image as input, perform automatic clustering and quantitative modeling of fracture structures that are distributed in a "sine-like" pattern in the borehole imaging unfolded image. Through processes such as connected region marking, composite distance merging, skeleton extraction, trunk smoothing, and sine curve fitting, key geometric parameters such as the amplitude, period, phase, and centerline position of the fracture are automatically obtained.

[0125] Furthermore, the specific steps for automatic clustering and quantitative modeling of "sinusoidal" cracks in step 1 above are as follows.

[0126] Step 1.1: Let the obtained binarized crack image be I. b (x,y), where (x,y) are the image pixel coordinates, x∈[0,W]. i -1],y∈[0,H i -1] When pixel value I b A fracture is defined when (x,y)=1. b When (x,y)=0, it is defined as a non-crack. The input binary crack image is initially segmented by the 8-connected region labeling algorithm, as shown in equation (1).

[0127] (1)

[0128] In the formula, C i This represents the nth connected region; This represents an 8-connected neighborhood.

[0129] Furthermore, considering that actual cracks may be divided into multiple disconnected regions due to noise or interruptions, two connected regions C are defined. i and C j The composite distance function between them is shown in equation (2).

[0130] (2)

[0131] In the formula, ||pq|| is the Euclidean distance between pixels p and q; , For region C i and C j The mean of the ordinate center; λ c This is a weighting coefficient used to balance the effects of spatial distance and longitudinal position differences.

[0132] Furthermore, when D c (C i C j )≤T m At that time, it is assumed that two regions belong to the same fracture, and after merging, the final set of fracture regions {F1, F2, ..., F} is obtained.m}

[0133] Step 1.2: For each clustered fracture region F k Calculate its distance transformation graph T k (x, y), its formula is shown in formula (3).

[0134] (3)

[0135] In the formula, 𝜕F k For region F k The set of boundary points.

[0136] Furthermore, after performing a distance transformation, the Zhang-Suen thinning algorithm is used to extract the skeleton line S of the fracture. k Its formula is shown in formula (4).

[0137] (4)

[0138] In the formula, ZS represents the Zhang-Suen refinement operator.

[0139] Furthermore, since the refined skeleton line may contain branches, the main trunk is extracted using the longest path algorithm, as shown in equation (5).

[0140] (5)

[0141] In the formula, LP represents the longest path extraction operator; P k This represents the main skeleton of the fissure.

[0142] Furthermore, the main skeleton line is subjected to moving average filtering to reduce local noise, as shown in equation (6).

[0143] (6)

[0144] In the formula, r s The radius of the sliding window; The main skeleton lines after smoothing.

[0145] Step 1.3: Use the mathematical expression of the "sinusoidal" fracture to fit the sine curve, as shown in equation (7).

[0146] (7)

[0147] In the formula, A S T represents the amplitude; S For periodicity; φ s For phase; Y S This is the position of the center line.

[0148] Furthermore, let the set of centerline points be {(x i ,y i )|i=1,2,…,N p}, define the fitting error function, which is shown in equation (8).

[0149] (8)

[0150] In the formula, E s N is the fitting error function; p This represents the number of points along the center line.

[0151] Furthermore, nonlinear optimization is performed using the Levenberg-Marquardt algorithm for parameter optimization, as shown in equation (9).

[0152] (9)

[0153] In the formula, θ S J is a parameter vector; s For Jacobian matrices; e t μ is the residual vector; s is the damping factor; I is the identity matrix.

[0154] Step 1.4: To verify the reliability and stability of the established clustering and sinusoidal fitting models, a quantitative effectiveness evaluation is conducted. Evaluation indicators are established based on detection rate, goodness of fit, root mean square error, stability, and robustness. Furthermore, a comprehensive performance scoring function is constructed to comprehensively measure the effectiveness of the model.

[0155] Furthermore, the detection rate is defined as the proportion of the number of components correctly identified as cracks to the total number of candidate components, as shown in equation (10).

[0156] (10)

[0157] In the formula, N c To satisfy the aspect ratio and coverage threshold, the number of crack sample components, N l This represents the number of connected components whose area exceeds the threshold.

[0158] Furthermore, the fitting effect is measured. Given a set of crack pixels {(x i ,y i )} N i=1 and the corresponding fitted values Define the goodness of fit R 2 The root mean square error (RMSE) is given by equation (11).

[0159] (11)

[0160] In the formula, These are the fitted values; This is the mean of the ordinates of all fracture points.

[0161] Furthermore, stability and robustness indices are introduced to measure the consistency of the model's performance and its ability to resist interference under different samples, as shown in equation (12).

[0162] (12)

[0163] In the formula, d This represents the standard deviation of the detection rate across different samples; and r S represents the mean and standard deviation of the root mean square error, respectively; b For stability indicators; R b This is a stability indicator.

[0164] Furthermore, a unified scoring function is introduced for comprehensive judgment, as shown in equation (13).

[0165] (13)

[0166] In the formula, The mean detection rate for all samples. For the average goodness of fit, Pr(Q) f >τ q ) represents the proportion of high-quality fits; τ q Goodness-of-fit threshold; Q s A comprehensive performance score is given.

[0167] Figure 2 This is an image showing the results of automatic clustering and numbering of cracks in each image in this embodiment, such as... Figure 2 As shown, after connected region marking, composite distance merging, and skeleton backbone extraction, different fracture segments can be automatically distinguished and assigned independent numbers, providing a basis for subsequent sinusoidal fitting and geometric parameter extraction.

[0168] Step 2: Based on the clustering and parameterization results, perform physical scale conversion and sampling on the crack profile, and calculate and evaluate the crack roughness. The steps for crack profile sampling and roughness evaluation are as follows: Figure 3 As shown, pixels are uniformly converted to millimeters to ensure that data from different cracks and attachments can be compared on the same physical scale, as shown in equation (14).

[0169] (14)

[0170] In the formula, Wi H is the image width. i L represents the image height. c L is the perimeter of the unfolded diagram. d The depth corresponding to the unfolded diagram; s x ,s y These represent the actual lengths of each pixel in the horizontal and vertical directions, respectively.

[0171] Furthermore, based on pixel conversion to obtain physical coordinates, the distance between adjacent points is calculated and accumulated segment by segment to obtain the total arc length, as shown in equation (15).

[0172] (15)

[0173] In the formula, x i ,y i Let i be the physical coordinates of the i-th crack point; The arc length between two adjacent points; Let be the cumulative arc length from the starting point to the j-th point.

[0174] Furthermore, a curvature adaptive method is introduced, which takes more points at turning points and fewer points on straight segments. The curvature magnitude is approximated by the second-order difference as the sampling weight, and combined with the arc length formula, a curvature-weighted adaptive sampling method is formed, as shown in equation (16).

[0175] (16)

[0176] In the formula, P i Let be the physical coordinates of the i-th point; к i This is an approximate value of the curvature at that point.

[0177] Furthermore, we need to establish a unified geometric roughness index, Z, to measure the overall volatility of the curve. r The essence of the index is the "mean square slope". It takes the average of the slopes of discrete points and takes the square root to obtain a value that represents the degree of undulation of the curve. It is directly related to the sampling method in the previous step and serves as the core input for the final JRC calculation. Its formula is shown in equation (17).

[0178] (17)

[0179] In the formula, N s ∆y is the number of sampling points; i The longitudinal difference between adjacent points; ∆x i Z represents the lateral difference between adjacent points. r This is the mean square slope index.

[0180] Furthermore, due to Z rIt is a geometric quantity, but the JRC value is actually used in engineering practice, so Z needs to be converted using empirical formulas. r Convert to JRC, as shown in equation (18).

[0181] (18)

[0182] In the formula, b J c J is the empirical fitting coefficient; JRC is the crack roughness coefficient.

[0183] Furthermore, the actual area of ​​the crack is calculated and standardized. The true area of ​​the crack is obtained by converting the pixel area with the pixel-to-millimeter ratio. Then, it is divided by the unfolded length to obtain the area per unit length. Combined with the JRC results, the effect of the crack area on roughness evaluation is analyzed, as shown in equation (19).

[0184] (19)

[0185] In the formula, N a A represents the total number of pixels in the crack; r L represents the actual area of ​​the crack. x A represents the horizontal length of the unfolded diagram. r / L x The area of ​​the crack per unit length.

[0186] Figure 4 This is a comparison chart of JRC results under different numbers of sampling points in this embodiment, as shown below. Figure 4 As shown, under different sampling density conditions, the crack roughness coefficient JRC varies with the sampling method and the number of sampling points, exhibiting different distribution patterns. This can be used to compare the impact of different sampling strategies on the roughness evaluation results.

[0187] Step 3: Project the 2D fracture pixels from the borehole unfolded diagram back onto the actual cylindrical borehole wall space to construct a 3D fracture point cloud and perform planar fitting. Based on this, and considering indices such as normal angle, geometric proximity, longitudinal continuity, borehole distance, and roughness factor, establish a fracture connectivity probability discrimination model. This enables connectivity analysis and 3D reconstruction of the multi-bore fracture network. The process flow for 3D connectivity reconstruction of multi-bore fractures is as follows: Figure 5 As shown.

[0188] Furthermore, the specific steps of the fracture connectivity probability analysis and three-dimensional back projection in step 3 above are as follows.

[0189] Step 3.1: To restore the planar coordinates to the real space, the actual radius of the borehole must first be obtained, as shown in equation (20).

[0190] (20)

[0191] In the formula, L C R is the circumference of the hole. h This represents the actual radius of the borehole.

[0192] Step 3.2: Establish a pixel-to-millimeter mapping, and make the width and height of the image correspond to the actual perimeter and depth, so that the pixel crack points can be accurately converted into millimeter coordinates. After combining with the hole radius calculated in the previous step, spatial restoration is achieved, as shown in equation (21).

[0193] (twenty one)

[0194] In the formula, L s The actual depth corresponding to the vertical direction of each unfolded diagram; s x ,s y These represent the millimeter values ​​corresponding to the horizontal and vertical pixels, respectively.

[0195] Furthermore, the crack point is projected onto the cylindrical hole wall, as shown in equation (22).

[0196] (twenty two)

[0197] In the formula, x m ,z m The x coordinates are the millimeter coordinates after the unfolded diagram is converted. h ,y h R represents the position of the hole center in the planar layout. h (X,Y,Z) represents the radius of the hole; (X,Y,Z) represents the coordinates of the point in three-dimensional space.

[0198] Step 3.3: Use the least squares method to fit the crack point cloud to obtain the plane normal vector and constant term, as shown in equation (23).

[0199] (twenty three)

[0200] In the formula, n p The plane normal vector; x=(X,Y,Z) are the coordinates of the point; d p This is a plane constant term.

[0201] Furthermore, the RMSE and fitting formula are used to verify whether the fracture segment can be regarded as the same plane, as shown in equation (24).

[0202] (twenty four)

[0203] In the formula, P i N represents the coordinates of the fracture point; p n is the number of points; p dp These are the normal vector and constant term of the fitted plane, respectively.

[0204] Step 3.4: Introduce a roughness index to prevent the problem of reduced penetration probability due to excessive tortuosity, as shown in equation (25).

[0205] (25)

[0206] In the formula, z i z is the longitudinal coordinate of the fracture point; l,i Z represents the fitted linear trend; r This is the mean square of fluctuation index.

[0207] Step 3.5: Combine the weighted sum with the penalty factor to unify the aforementioned indicators to a probability value of 0–1, and obtain a comprehensive score for whether the two pores are connected, as shown in equation (26).

[0208] (26)

[0209] In the formula, S n For normality consistency scoring; S p Scoring for consistency in planar distance; S v For longitudinal continuity scoring; S h J represents the weight of the distance between holes; p ν is the roughness penalty factor; n ,ν p ,ν ν P represents the weighting coefficient. c The connectivity probability is scored.

[0210] Figure 6 This is a planar diagram showing the crack penetration result in this embodiment, as shown below. Figure 6 As shown, after the two-dimensional fracture points are projected back onto the borehole wall of the borehole cylinder, the connection relationship of the fractures in three-dimensional space and the fitting plane results can be obtained, realizing the spatial connectivity analysis of the multi-bore fracture network.

[0211] Step 4: To further improve the completeness of fracture network identification and the reliability of connectivity discrimination, uncertainty analysis is performed on the study area and supplementary borehole optimization is carried out. By constructing comprehensive indicators such as coverage field, planar proximity and inter-hole distance, the current insufficient borehole area is determined, and the optimal supplementary borehole location and drilling depth are recommended using greedy point selection and dispersion constraint strategies.

[0212] Furthermore, the specific steps for uncertainty region analysis and supplementary borehole optimization in step 4 above are as follows.

[0213] Step 4.1: Use the Gaussian kernel function to spatially superimpose the center points of the existing fracture segments to obtain a coverage field. The larger the value, the denser the surrounding existing fracture information. The formula is shown in Equation (27).

[0214] (27)

[0215] In the formula, (X,Y) are the coordinates of a point on the planar grid; (x,Y) are the coordinates of a point on the grid. c,i ,y c,i ) represents the coordinates of the center point of the fracture segment; σ g Gaussian kernel scale; G min G max These are the minimum and maximum values ​​under normalization, respectively; ε g To prevent decimals with a denominator of zero.

[0216] Furthermore, the minimum distance from the point to the plane is calculated to measure the proximity of the hole-filling point to these planes, for the calculation of the uncertainty field, as shown in equation (28).

[0217] (28)

[0218] In the formula, (x,y,z) are the coordinates of the candidate points; (n xi ,n yi ,n zi ) represents the unit normal vector component of the i-th plane; d i This is a plane constant term.

[0219] Furthermore, the shortest distance from the candidate point to all the arranged orifices is calculated, as shown in equation (29).

[0220] (29)

[0221] In the formula, (x,y) are the planar coordinates of the candidate point; (x j ,y j ) represents the plane coordinates of the j-th borehole.

[0222] Step 4.2: Combine the relevant indicators of the rationality of the hole filling point into a unified uncertainty index, as shown in equation (30).

[0223] (30)

[0224] In the formula, δ h (x,y) is the shortest distance from the point to the existing hole; δ p (x,y,z) is the distance from the point to the nearest crack plane; G c (x,y) represents the coverage field; q1,q2,q3 are weight coefficients.

[0225] Furthermore, a distance penalty is introduced between points to ensure a more uniform distribution of recommended boreholes, as shown in equation (31).

[0226] (31)

[0227] In the formula, Q represents the planar coordinates of the candidate point; P is the selected borehole point; ||Q−P|| is the Euclidean distance between the two points; R d To influence the scale of the scope; This represents the mean value of the workable area U.

[0228] Furthermore, an objective function is constructed to filter points. In each round of selection, the candidate points are scored, and the one with the highest score is selected as the new hole-filling point. This is closely combined with the penalty formula in the previous step, as shown in equation (32).

[0229] (32)

[0230] In the formula, U c (x,y) represents the uncertainty value of the candidate point; П(x,y) is the distance penalty function; F o (x,y) represents the overall objective value.

[0231] Step 4.3: Find the optimal depth within the longitudinal z range, as shown in equation (33).

[0232] (33)

[0233] In the formula, z ∗ (x,y) represents the optimal depth corresponding to the candidate point; Z c h represents the depth of the search window center. z The width of the window is half its width; L b Recommended drilling length.

[0234] Figure 7 To supplement this embodiment, a drilling optimization result diagram is provided, such as... Figure 7 As shown, after constructing an uncertainty field and integrating factors such as coverage, planar proximity, and hole-to-hole distance, recommended supplementary borehole locations and drilling depths can be obtained, which can improve the completeness of fracture network identification and the reliability of connectivity discrimination.

[0235] Although the present invention has been disclosed above with reference to preferred embodiments, these are not intended to limit the invention. Any person skilled in the art can make various changes or modifications without departing from the spirit and scope of the invention. Therefore, the scope of protection of the present invention should be defined by the scope of the claims of this application.

Claims

1. A method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm, characterized in that, The overall process of this model includes the following four steps; Step 1: Using the binarized fracture image as input, perform automatic clustering and quantitative modeling of fracture structures that are distributed in a "sinusoidal" pattern in the borehole imaging unfolded image. Through processes such as connected region marking, composite distance merging, skeleton extraction, trunk smoothing and sinusoidal curve fitting, key geometric parameters such as the amplitude, period, phase and centerline position of the fracture are automatically obtained. Step 2: Based on the clustering and parameterization results, the physical scale of the crack profile is converted and sampled, and the crack roughness is calculated for evaluation. The pixels are uniformly converted to millimeters to ensure that the data of different cracks and different attachments can be compared under the same physical scale. The formula is shown in formula (14). （14） In the formula, W i H is the image width. i L represents the image height. c L is the perimeter of the unfolded diagram. d The depth corresponding to the unfolded diagram; s x ,s y These represent the actual lengths of each pixel in the horizontal and vertical directions, respectively. Furthermore, based on pixel conversion to obtain physical coordinates, the distance between adjacent points is calculated and accumulated segment by segment to obtain the total arc length, as shown in equation (15). （15） In the formula, x i ,y i Let i be the physical coordinates of the i-th crack point; The arc length between two adjacent points; This represents the cumulative arc length from the starting point to the j-th point; Furthermore, a curvature adaptive method is introduced, which takes more points at the turning points and fewer points on the straight sections. The curvature magnitude is approximated by the second-order difference as the sampling weight, and combined with the arc length formula, a curvature-weighted adaptive sampling method is formed, as shown in equation (16). （16） In the formula, P i Let be the physical coordinates of the i-th point; к i This is an approximate value of the curvature at that point; Furthermore, we need to establish a unified geometric roughness index, Z, to measure the overall volatility of the curve. r The essence of the index is "mean square slope". It takes the average of the slopes of discrete points, takes the square root, and obtains a value representing the degree of undulation of the curve. It is directly related to the sampling method in the previous step and serves as the core input for the final JRC calculation. Its formula is shown in formula (17). （17） In the formula, N s ∆y is the number of sampling points; i The longitudinal difference between adjacent points; ∆x i Z represents the lateral difference between adjacent points. r The mean square slope index; Furthermore, due to Z r It is a geometric quantity, but the JRC value is actually used in engineering practice, so Z needs to be converted using empirical formulas. r Converted to JRC, its formula is shown in equation (18); （18） In the formula, b J c J The coefficients are empirical fitting coefficients; JRC is the crack roughness coefficient. Furthermore, the actual area of ​​the crack is calculated and standardized. The true area of ​​the crack is obtained by converting the pixel area with the pixel-to-millimeter ratio. Then, it is divided by the unfolded length to obtain the area per unit length. Combined with the JRC results, the effect of the crack area on roughness evaluation is analyzed. The formula is shown in formula (19). （19） In the formula, N a A represents the total number of pixels in the crack; r L represents the actual area of ​​the crack. x A represents the horizontal length of the unfolded diagram. r / L x The area of ​​the crack per unit length; Step 3: Project the two-dimensional fracture pixels in the borehole unfolding diagram back onto the actual cylindrical borehole wall space to construct a three-dimensional point cloud of fractures and perform plane fitting. Based on this, integrate indicators such as normal angle, geometric proximity, longitudinal continuity, borehole distance and roughness factor to establish a fracture connectivity probability discrimination model, thereby realizing the connectivity analysis and three-dimensional reconstruction of the multi-bore fracture network. Step 4: To further improve the completeness of fracture network identification and the reliability of connectivity discrimination, uncertainty analysis is performed on the study area and supplementary borehole optimization is carried out. By constructing comprehensive indicators such as coverage field, planar proximity and inter-hole distance, the current insufficient borehole area is determined, and the optimal supplementary borehole location and drilling depth are recommended using greedy point selection and dispersion constraint strategies.

2. The method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm according to claim 1, characterized in that, The specific steps for automatic clustering and quantitative modeling of "sinusoidal" fractures in step 1 above are as follows; Step 1.1: Let the obtained binarized crack image be I. b (x,y), where (x,y) are the image pixel coordinates, x∈[0,W]. i -1],y∈[0,H i -1] When pixel value I b A fracture is defined when (x,y)=1. b When (x,y)=0, it is defined as a non-crack; the input binary crack image is initially segmented by the 8-connected region labeling algorithm, as shown in equation (1); （1） In the formula, C i This represents the nth connected region; Indicates an 8-connected neighborhood; Furthermore, considering that actual cracks may be divided into multiple disconnected regions due to noise or interruptions, two connected regions C are defined. i and C j The composite distance function between them is shown in equation (2); （2） In the formula, ||pq|| is the Euclidean distance between pixels p and q; , For region C i and C j The mean of the center of the ordinate; λ c These are weighting coefficients used to balance the effects of spatial distance and longitudinal position differences; Furthermore, when D c (C i C j )≤T m At that time, it is assumed that two regions belong to the same fracture, and after merging, the final set of fracture regions {F1, F2, ..., F} is obtained. m }; Step 1.2: For each clustered fracture region F k Calculate its distance transformation graph T k (x, y), its formula is shown in formula (3); （3） In the formula, 𝜕F k For region F k The set of boundary points; Furthermore, after performing a distance transformation, the Zhang-Suen thinning algorithm is used to extract the skeleton line S of the fracture. k Its formula is shown in formula (4); （4） In the formula, ZS represents the Zhang-Suen refinement operator; Furthermore, since the refined skeleton line may contain branches, the main trunk is extracted by the longest path algorithm, as shown in equation (5). （5） In the formula, LP represents the longest path extraction operator; P k Indicates the main framework of the fracture; Furthermore, the main skeleton line is subjected to moving average filtering to reduce local noise, as shown in equation (6). （6） In the formula, r s The radius of the sliding window; The smoothed main skeleton lines; Step 1.3: Use the mathematical expression of "sinusoidal" fracture to fit a sine curve, as shown in equation (7); （7） In the formula, A S T represents the amplitude; S For periodicity; φ s For phase; Y S The centerline position; Furthermore, let the set of centerline points be {(x i ,y i )|i=1,2,…,N p }, define the fitting error function, as shown in equation (8); （8） In the formula, E s N is the fitting error function; p The number of points along the center line; Furthermore, nonlinear optimization is performed using the Levenberg-Marquardt algorithm for parameter optimization, as shown in equation (9). （9） In the formula, θ S J is a parameter vector; s For Jacobian matrices; e t μ is the residual vector; s I is the damping factor; I is the identity matrix; Step 1.4: To verify the reliability and stability of the established clustering and sinusoidal fitting models, a quantitative effectiveness evaluation is conducted. Evaluation indicators are established from aspects such as detection rate, goodness of fit, root mean square error, stability and robustness. Furthermore, a comprehensive performance scoring function is constructed to comprehensively measure the effectiveness of the model. Furthermore, the detection rate is defined as the proportion of the number of components correctly identified as cracks to the total number of candidate components, as shown in equation (10). （10） In the formula, N c To satisfy the aspect ratio and coverage threshold, the number of crack sample components, N l The number of connected components whose area exceeds the threshold; Furthermore, the fit is measured; Given a set of crack pixels {(x i ,y i )} N i=1 and the corresponding fitted values Define the goodness of fit R 2 The root mean square error (RMSE) is given by equation (11). （11） In the formula, These are the fitted values; This represents the mean of the ordinates of all fracture points; Furthermore, stability and robustness indices are introduced to measure the consistency of the model’s performance and its ability to resist interference under different samples, as shown in equation (12). （12） In the formula, d This represents the standard deviation of the detection rate across different samples; and r S represents the mean and standard deviation of the root mean square error, respectively; b For stability indicators; R b For stability indicators; Furthermore, a unified scoring function is introduced for comprehensive judgment, as shown in equation (13); （13） In the formula, The mean detection rate for all samples. For the average goodness of fit, Pr(Q) f >τ q ) represents the proportion of high-quality fits; τ q Goodness-of-fit threshold; Q s This is a comprehensive performance score.

3. The method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm according to claim 1, characterized in that, The specific steps of the fracture connectivity probability analysis and three-dimensional back projection in step 3 above are as follows; Step 3.1: To restore the planar coordinates to the real space, the actual radius of the borehole must first be obtained, as shown in equation (20); （20） In the formula, L C R is the circumference of the hole. h The actual radius of the borehole; Step 3.2: Establish a pixel-to-millimeter mapping, and make the width and height of the image correspond to the actual perimeter and depth, so that the pixel crack points can be accurately converted into millimeter coordinates. After combining with the hole radius calculated in the previous step, spatial restoration is achieved, as shown in equation (21). （21） In the formula, L s The actual depth corresponding to the vertical direction of each unfolded diagram; s x ,s y These represent the millimeter values ​​corresponding to the horizontal and vertical pixels, respectively. Furthermore, the crack point is projected onto the cylindrical hole wall, as shown in equation (22); （22） In the formula, x m ,z m The x coordinates are the millimeter coordinates after the unfolded diagram is converted. h ,y h R represents the position of the hole center in the planar layout. h Let (X,Y,Z) be the radius of the hole; let (X,Y,Z) be the coordinates of the point in three-dimensional space. Step 3.3: Fit the crack point cloud using the least squares method to obtain the plane normal vector and constant term, as shown in equation (23); （23） In the formula, n p The plane normal vector; x=(X,Y,Z) are the coordinates of the point; d p For plane constants; Furthermore, the RMSE and fitting formula are used to verify whether the fracture segment can be regarded as the same plane, as shown in equation (24). （24） In the formula, P i N represents the coordinates of the fracture point; p n is the number of points; p d p These are the normal vector and constant term of the fitted plane, respectively; Step 3.4: Introduce a roughness index to prevent the reduced possibility of penetration due to excessive tortuosity, as shown in equation (25); （25） In the formula, z i z is the longitudinal coordinate of the fracture point; l,i Z represents the fitted linear trend; r The mean square of fluctuations; Step 3.5: Combine the weighted sum with the penalty factor to unify the aforementioned indicators to a probability value of 0–1, and obtain a comprehensive score for whether the two pores are connected, as shown in equation (26). （26） In the formula, S n For normality consistency scoring; S p Scoring for consistency in planar distance; S v For longitudinal continuity scoring; S h J represents the weight of the distance between holes; p ν is the roughness penalty factor; n ,ν p ,ν ν P represents the weighting coefficient. c The connectivity probability is scored.

4. The method for reconstructing three-dimensional connectivity of fractures in downhole surrounding rock based on the CS-HAT algorithm according to claim 1, characterized in that, The specific steps for uncertainty region analysis and supplementary borehole optimization in step 4 above are as follows; Step 4.1: Use the Gaussian kernel function to spatially superimpose the center points of the existing fracture segments to obtain a coverage field. The larger the value, the denser the surrounding existing fracture information. The formula is shown in formula (27). （27） In the formula, (X,Y) are the coordinates of a point on the planar grid; (x,Y) are the coordinates of a point on the grid. c,i ,y c,i ) represents the coordinates of the center point of the fracture segment; σ g Gaussian kernel scale; G min G max These are the minimum and maximum values ​​under normalization, respectively; ε g To prevent decimals with a denominator of zero; Furthermore, the minimum distance from the point to the plane is calculated to measure the proximity of the hole-filling point to these planes, for the calculation of the uncertainty field, as shown in equation (28); （28） In the formula, (x,y,z) are the coordinates of the candidate points; (n xi ,n yi ,n zi ) represents the unit normal vector component of the i-th plane; d i For plane constants; Furthermore, the shortest distance from the candidate point to all the arranged orifices is calculated, as shown in equation (29); （29） In the formula, (x,y) are the planar coordinates of the candidate point; (x j ,y j () represents the plane coordinates of the j-th borehole; Step 4.2: Combine the relevant indicators of the rationality of the hole filling point into a unified uncertainty index, the formula of which is shown in formula (30); （30） In the formula, δ h (x,y) is the shortest distance from the point to the existing hole; δ p (x,y,z) is the distance from the point to the nearest crack plane; G c (x,y) represents the coverage field; q1,q2,q3 are weight coefficients; Furthermore, a distance penalty between points is introduced to ensure a more uniform distribution of recommended boreholes, as shown in equation (31). （31） In the formula, Q represents the planar coordinates of the candidate point; P is the selected borehole point; ||Q−P|| is the Euclidean distance between the two points; R d To influence the scale of the scope; The mean value of the workable area U; Furthermore, an objective function is constructed to filter points. In each round of selection, the candidate points are scored, and the one with the highest score is selected as the new hole-filling point. This is closely combined with the penalty formula in the previous step, as shown in equation (32). （32） In the formula, U c (x,y) represents the uncertainty value of the candidate point; П(x,y) is the distance penalty function; F o (x,y) represents the overall objective value; Step 4.3: Find the optimal depth within the longitudinal z range, as shown in equation (33); （33） In the formula, z ∗ (x,y) represents the optimal depth corresponding to the candidate point; Z c h represents the depth of the search window center. z The width of the window is half its width; L b Recommended drilling length.

Images