A method for extracting fracture attitude parameters based on borehole images

Abstract

A method for extracting fracture attitude parameters based on borehole images, which relates to the field of geological resources and geological engineering technology, mainly includes the following steps: collecting and importing the original borehole images of the drilling area, enhancing and denoising the original borehole images, binarizing the images through dynamic region division, using the Canny operator that adaptively selects thresholds based on the gradient histogram of the borehole images to detect the edges of the fracture structural planes, fitting the fracture structural planes with sine curves, and extracting the fracture attitude parameters on the fitted fracture sine curves; based on engineering borehole images, the present invention can automatically identify the rock mass structural planes and extract the attitude parameters, replacing the previous manual intervention and operation, and improving the efficiency of analyzing geological conditions and groundwater migration in geological engineering.

Description

Technical Field

[0001] The present invention relates to the technical field of geological resources and geological engineering, and particularly relates to a method for extracting fracture occurrence parameters based on borehole images. Background Art

[0002] The structural plane is one of the basic information required in geological engineering exploration. Various fracture structural planes in rock masses are also important places and channels for groundwater storage and migration. Their occurrence information such as dip angle and dip direction will affect the flow rate and direction of groundwater. Correctly analyzing the occurrence information of regional structural planes is the basis for engineering geological activities.

[0003] With the development of digital panoramic borehole camera technology, extracting structural plane parameters through borehole wall images has become an important means. In borehole images, the interpretation of structural plane information usually relies on manual operation, which is cumbersome, time-consuming, and different people may have different understandings and interpretation results for the same borehole wall image. And in existing automatic recognition technologies, such as wavelet analysis, clustering algorithms, three-dimensional Hough transform and other methods, although they have certain effects in automatically recognizing structural planes, the edge detection effect of fracture structural planes in images is relatively fuzzy, affecting the accuracy of final parameter extraction. In order to extract the parameters of rock mass structural planes practically, conveniently and accurately, it is particularly important to develop a set of efficient methods to obtain key information of structural planes and realize the automatic extraction of structural plane parameters.

[0004] Since the trace of the fracture structural plane after cutting the borehole has the characteristics of a sine curve in the planar image of the borehole wall image, the prior art usually relies on extracting structural plane edge information and parameter space mapping theory to identify the fracture structural plane, so as to obtain parameter information such as the strike, dip direction and dip angle of the fracture structural plane. In this process, the Canny operator is often used. The Canny operator is a commonly used means for detecting image edges, mainly including several steps such as image smoothing and denoising, calculating gradient magnitude and direction, non-maximum suppression, double-threshold detection, and edge connection. Among them, for double-threshold detection, the traditional method is to set two thresholds T H , T L , points with gradient values larger than T H are strong edge points, smaller than T H but larger than T LThe large dots are weak edge points. When performing threshold processing, only strong edge points and weak edge points connected to strong edges will be retained to effectively identify edges and exclude non-edge points. However, there is no unified standard for selecting these two thresholds, and they are often set manually, lacking adaptability and robustness. When applied to borehole images, due to the complexity of the geological environment, borehole images often have many noise interferences, such as fractured rock strata, dim light, and interlayer cements, etc., resulting in the difficulty of manually setting threshold parameters to adapt to different borehole images, and the need to repeatedly try and adjust to correct the detection results. Summary of the Invention

[0005] In view of this, the present invention proposes a method for extracting fracture attitude parameters based on borehole images. When applying the Canny operator to detect the edges of fracture structural planes in borehole images, the present invention proposes an adaptive threshold selection based on the gradient histogram of borehole images, which can automatically find appropriate thresholds to distinguish edges and non-edges.

[0006] To solve the above-mentioned at least one technical problem, the technical solution provided by the present invention is a method for extracting fracture attitude parameters based on borehole images, including the following steps:

[0007] Step S1: Collect and import the original borehole images of the drilling area;

[0008] Step S2: Enhance and denoise the original borehole images;

[0009] Step S3: Binarize the image through dynamic region division;

[0010] Step S4: Use the Canny operator with adaptive threshold selection based on the gradient histogram of borehole images to detect the edges of fracture structural planes;

[0011] Step S5: Fit a sine curve to the fracture structural plane;

[0012] Step S6: Extract fracture attitude parameters on the fitted fracture sine curve.

[0013] The technical effects achieved by the present invention are as follows:

[0014] 1. Based on engineering borehole images, the present invention can automatically identify rock mass structural planes and extract attitude parameters, replacing the previous manual intervention and operation, and improving the efficiency of analyzing geological conditions and groundwater migration in geological engineering.

[0015] 2. Compared with traditional binarization methods, the present invention uses dynamic region division to achieve the binarization of borehole images, and will more accurately retain the information of fracture structural planes when applied to borehole images.

[0016] 3. When the present invention performs crack structural plane edge detection, it improves the conventional Canny operator and adopts an adaptive threshold selection based on the gradient histogram of the borehole image to replace the artificial trial-and-error threshold parameters, which can more conveniently and accurately detect the crack edges in the borehole image.

[0017] 4. The present invention utilizes the characteristic that the traces after the structural plane cuts the borehole show a sine curve when the image on the borehole wall is unfolded into a planar image, and directly uses the characteristics of the sine function to fit the crack structural plane. Compared with the conventional curve fitting, it improves the detection efficiency and accuracy.

[0018] 5. The method of the present invention runs fast, has good controllability during the parameter extraction process, can provide real-time feedback on the recognition effect, and the accuracy and precision of the recognition can be adjusted by parameters for subsequent data use. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following will briefly introduce the drawings required for the embodiments. It should be understood that the following drawings only show some embodiments of the present invention, so they should not be regarded as limiting the scope. For those of ordinary skill in the art, without creative efforts, other related drawings can also be obtained based on these drawings.

[0020] Figure 1 is the specific flowchart of the present invention;

[0021] Figure 2 is a schematic diagram of three cases where the inclined structural plane cuts the borehole when the borehole image of the present invention is unfolded into a planar image;

[0022] Figure 3 where a is the planar development diagram of the original borehole image in the present invention;

[0023] Figure 3 where b is the edge detection result diagram of the crack structural plane using the improved Canny operator in the present invention;

[0024] Figure 3 where c is the search detection result diagram of the sine curve baseline y0 in the present invention;

[0025] Figure 3 where d is the sine curve fitting effect diagram of the structural plane in the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0026] The following will further elaborate on the present invention in detail in combination with the embodiments and the drawings.

[0027] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Apparently, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the claimed invention, but merely represents selected embodiments of the present invention.

[0028] Example:

[0029] See Figure 1 , a method for extracting fracture occurrence parameters based on borehole images, comprising the following steps:

[0030] Step S1: Collect the original borehole images in the drilling area. The original borehole images include borehole images and planar unfolded diagrams after unfolding the borehole images, as shown in Figure 2 and Figure 3 a shown therein;

[0031] Step S2: Enhance and denoise the original borehole images. Among them, the enhancement and denoising method in Step S2 is adaptive median filtering. Different from conventional median filtering, this adaptive median filtering dynamically adjusts the window size according to the variance of the gray values of the pixels within the window. When the variance of the gray values within the window is large, it indicates that there may be noise points within the window. At this time, the filter increases the window size to include more pixel points to calculate the median, so as to better eliminate the noise. On the contrary, when the variance of the gray values within the window is small, it indicates that the main part within the window is the true features of the image. At this time, the filter reduces the window size to retain more image details;

[0032] Step S3: Binarize the image through dynamic region division;

[0033] The steps of binarizing the image through the dynamic region division described in Step S3 are as follows:

[0034] 1) Calculate the gray histogram and local variance of the image. Regions with large local variance often contain more details or target boundaries. Judge whether there are multi-peaks, the positions and widths of the peaks, etc. information, and use the valleys between the peaks of the gray histogram as the threshold;

[0035] 2) Start small-window scanning from the upper left corner of the image. Scan with a small window (such as 3x3 or 5x5 pixels). When the gray-scale difference of the pixels within the window is less than the threshold and the pixels form a connected region, divide this region into an initial region. Scan the entire image in sequence to complete the division. When the gray-scale difference is greater than the threshold, it can be considered that the region division straddles the background and boundaries of different targets, and the scanning window needs to be adjusted until the gray-scale difference is less than the threshold.

[0036] 3) Use the watershed algorithm to perform segmentation processing on the divided regions, and perform binary processing on the segmented regions. For regions with relatively uniform gray-scale distribution, determine the threshold according to the valley value of the histogram. For regions with non-uniform gray-scale distribution within the region, determine the threshold according to the local mean and standard deviation. Assign pixels with gray-scale values greater than or equal to the region threshold to 255, and assign pixels with gray-scale values less than the region threshold to 0.

[0037] Among them, the threshold of the region with uniform gray-scale distribution after segmentation processing is the valley value between the peaks of the gray-scale histogram, and the threshold of the region with non-uniform gray-scale distribution after segmentation processing is determined according to the local mean and standard deviation.

[0038] Among them, the threshold of the region with non-uniform gray-scale distribution after segmentation processing is shown in Equation (1):

[0039] T = μ + kσ (1)

[0040] In Equation (1), T is the threshold of the region with non-uniform gray-scale distribution after segmentation processing, μ is the local mean of the region gray-scale, σ is the standard deviation of the region gray-scale, and k is an empirical coefficient with a value between 0.5 and 2. The specific value is obtained through testing of the drilling map and is usually taken as 1.

[0041] Step S4: Use the Canny operator that adaptively selects the threshold based on the gradient histogram of the drilling image to detect the edges of the fracture structural plane;

[0042] The basic principle is to calculate the amplitude and direction of the image gradient, suppress non-maximum values, leave the points with the largest region gradient, and suppress other points to zero, thereby refining the edges and making the fracture edges clearer.

[0043] Combined with Figure 3 in b, the steps to use the Canny operator that adaptively selects the threshold based on the gradient histogram of the drilling image to detect the edges of the fracture structural plane are as follows:

[0044] 1) Set the image gradient amplitude range to n intervals of [0, M max , where the width of any interval is shown in Equation (2):

[0045] ΔM = M max / n (2)

[0046] In formula (2), ΔM is the width of any interval, M max is the maximum amplitude value, and n is the number of intervals;

[0047] 2) Count the number of pixels H(k) in each interval, where k = 1, 2, …, n. Here, H(k) represents the number of pixels whose gradient amplitude falls within the k-th interval, and construct a gradient amplitude histogram;

[0048] 3) Determine the valleys in the histogram. The low-amplitude gradients are mainly generated by noise, and the valleys represent regions with relatively few pixels. Therefore, calculate the first-order difference of the gradient amplitude histogram according to the content shown in formula (3):

[0049] D(k) = H(k + 1) - H(k) (3)

[0050] In formula (3), D(k) is the first-order difference in the k-th region of the gradient amplitude histogram, H(k + 1) is the number of pixels in the adjacent (k + 1)-th region, and H(k) is the number of pixels in the k-th interval;

[0051] After that, calculate the second-order difference according to the content shown in formula (4) based on the first-order difference:

[0052] D2(k) = D(k + 1) - D(k) (4)

[0053] In formula (4), D2(k) is the second-order difference in the k-th region of the gradient amplitude histogram, and D(k + 1) is the first-order difference in the adjacent (k + 1)-th region.

[0054] When D(k) < 0 and D2(k) > 0, it is considered a valley. At the same time, compare the magnitudes of D2(k) corresponding to multiple obtained valleys, and select the position with the largest D2(k) as the most likely valley;

[0055] 4) Take the gradient amplitude corresponding to the valley as the low threshold, and determine the high threshold according to the distance from the valley to the high-amplitude part and the trend of gradient amplitude change. The part on the right side of the valley is mainly the gradient amplitude related to edges. Determine the high threshold according to the distance from the valley to the high-amplitude part and the trend of gradient amplitude change;

[0056] Among them, the method for determining the high threshold is to detect the magnitude of the second-order difference along the direction of increasing gradient amplitude and judge its front and back fluctuation amplitudes. As the amplitude increases, the value of the second-order difference gradually decreases. When its value fluctuates less than 0.5 with the previous and next values, it means that the growth rate of the gradient amplitude significantly slows down. Take the gradient amplitude corresponding to this position as the high threshold.

[0057] 5) Mark the pixels with gradient magnitude greater than the high threshold as strong edge pixel points, mark the pixels with gradient magnitude less than the low threshold as non-edge pixel points, and mark the pixels between the two as weak edge pixel points. According to the strong edge points, connect the weak edge points adjacent to the strong edge points to construct the edge contour of the complete fracture structural plane.

[0058] Step S5: Fit a sine curve to the fracture structural plane;

[0059] Combined with Figure 3 c in Figure 3 and d in

[0060] Map the image coordinate space to the parameter space. In this process, curves in the original image (such as straight lines, circles, etc.) are mapped to points or curves in the parameter space (the space composed of slope and intercept), that is, the dual property between points and lines in two different coordinate systems. The fracture trajectory line satisfies Equation (4):

[0061]

[0062] In Equation (4), A is the amplitude of the sine function; y0 is the baseline position of the sine function; is the initial phase of the sine function; ω = 2π / T; T is the period, specifically the width of the image.

[0063] Among them, the position of the baseline y0 is detected by means of a voting mechanism. The basic process is to create a one-dimensional accumulator array and assign all elements in the array to 0; traverse the boundary points in the image and all points on the vertical scan line where they are located, and try to pair these points; after finding the paired points, calculate the ordinate of the midpoint of these two paired points and add 1 to the value of the corresponding accumulator array element; continue the whole process until all point pairs in the image are visited and processed; finally, find the maximum value in the accumulator array, and the position corresponding to this maximum value is the baseline position of the sine curve.

[0064] After determining the baseline position of the sine curve, continue to extract the amplitude and initial phase through parameter space transformation. According to Equation (4) deformed to obtain Equation (5):

[0065]

[0066] Among them, The range of is from 0 to T; the range of A is from 0 to A Max ;

[0067] Create a two-dimensional accumulator array and assign an initial value of 0; for the detected y0 above, traverse each boundary point P(x, y) in the image; for For each value (from 0 to T), A can be calculated through formula (5). If A is between 0 and A Max If it is within this range, then the count of the corresponding element in the accumulator array is incremented by 1. Repeat the above steps until all boundary pixel points have been processed, and find the peak position in the accumulator array. Repeat the above steps until all y0 points have been processed;

[0068] Step S6: Extract the fracture occurrence parameters from the fitted fracture sine curve;

[0069] The fracture occurrence parameters include the dip angle and dip direction of the fracture structural plane. The vertical distance H (i.e., twice the amplitude A of the sine curve) between the peak and trough of the fracture sine curve obtained through step S5 and the drilling diameter D are used to calculate the dip angle β through the arctangent function, as shown in formula (6) specifically:

[0070]

[0071] The dip direction can be directly read from the azimuth corresponding to the lowest point of the trace (i.e., the trough of the sine curve).

[0072] In the description of the present invention, it should be noted that the orientation or positional relationship indicated by the terms "upper", "lower", "front", "rear", "left", "right", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings. It is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and should not be construed as a limitation to the present invention.

[0073] The above is only a preferred specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed in the embodiments of the present invention should be covered by the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims

1. A method for extracting fracture occurrence parameters based on borehole images, characterized in that: The following steps are involved: Step S1: Collect and import the original borehole images of the drilling area; Step S2: enhancing and denoising the original drilling image; Step S3: binarizing the image by dynamic region segmentation; Step S4: using the Canny operator with adaptive threshold selection based on the gradient histogram of the drilling image to detect the edge of the fracture structure surface; Step S5: performing sinusoidal curve fitting on the fracture structure surface; Step S6: extracting fracture occurrence parameters from the fitted fracture sinusoidal curve; The step of detecting the edge of the fracture structure surface by using the Canny operator for adaptive threshold selection based on the gradient histogram of the drilling image in step S4 is as follows: 1) Set the image gradient amplitude range to [0, M max ], where the width of any interval is as shown in formula (2): ΔM=M max / n (2) In formula (2), ΔM is the width of any interval, M max is the maximum amplitude, n is the number of intervals; 2) Count the number of pixels in each interval and construct a gradient amplitude histogram; 3) Determine the trough in the histogram, where the trough is determined by calculating the first-order difference of the gradient amplitude histogram according to formula (3): D(k)=H(k+1)-H(k) (3) In formula (3), D(k) is the first-order difference in the kth region of the gradient magnitude histogram, H(k+1) is the number of pixels in the adjacent k+1th region, and H(k) is the number of pixels in the kth region; Then, the second-order difference is calculated based on the first-order difference according to formula (4): D2(k)=D(k+1)-D(k) (4) In formula (4), D2(k) is the second-order difference in the kth region of the gradient amplitude histogram, and D(k+1) is the first-order difference in the adjacent k+1th region; The position where the second-order difference D2(k) is the largest is taken as the trough; 4) The gradient amplitude corresponding to the trough is used as the low threshold, and the high threshold is determined according to the distance from the trough to the high amplitude part and the gradient amplitude change trend; 5) Pixels with gradient amplitude greater than the high threshold are marked as strong edge pixels, pixels with gradient amplitude less than the low threshold are marked as non-edge pixels, and pixels between the two are marked as weak edge pixels. The weak edge points adjacent to the strong edge points are connected according to the strong edge points to construct the edge contour of the complete fracture structure surface; The method for determining the high threshold is to detect the size of the second-order difference along the direction of increasing gradient amplitude and determine its fluctuation amplitude before and after. When the fluctuation between its value and the two values ​​before and after is less than 0.5, the gradient amplitude corresponding to this position is used as the high threshold.

2. The method for extracting fracture occurrence parameters based on borehole images according to claim 1, characterized in that: The enhancement denoising method described in step S2 is adaptive median filtering.

3. The method for extracting fracture occurrence parameters based on borehole images according to claim 1, characterized in that: The steps of binarizing the image by dynamic area division in step S3 are: 1) Calculate the grayscale histogram and local variance of the image, and use the valley value between the peaks of the grayscale histogram as the threshold; 2) Start scanning a small window from the upper left corner of the image. When the grayscale difference of the pixels in the window is less than the threshold and the pixels form a connected area, the area is divided into an initial area, and the entire image is scanned in sequence to complete the division; 3) Use the watershed algorithm to segment the divided area, and perform binarization on the segmented area. Pixels with grayscale values ​​greater than or equal to the regional threshold are assigned a value of 255, and pixels with grayscale values ​​less than the regional threshold are assigned a value of 0; The threshold of the region with uniform grayscale distribution after segmentation is the valley value between the peaks of the grayscale histogram, and the threshold of the region with uneven grayscale distribution after segmentation is determined according to the local mean and standard deviation of the regional grayscale.

4. The method for extracting fracture occurrence parameters based on borehole images according to claim 3, characterized in that: The window size of the small window scanning is 3×3 or 5×5 pixels.

5. The method for extracting fracture occurrence parameters based on borehole images according to claim 3, characterized in that: The threshold of the area with uneven grayscale distribution after the segmentation process is shown in formula (1): T = μ + kσ （1） In formula (1), T is the threshold of the area with uneven grayscale distribution after segmentation processing, μ is the local mean of regional grayscale, σ is the regional grayscale standard deviation, k It is an empirical coefficient with a value between 0.5-2.

6. The method for extracting fracture occurrence parameters based on borehole images according to claim 1, characterized in that: The steps of performing sinusoidal curve fitting on the fracture structure surface in step S5 are: Mapping the image coordinate space to the parameter space, the crack trajectory satisfies equation (5): （5） In formula (5), A is the amplitude of the sine function; y0 is the baseline position of the sine function; is the initial phase of the sine function; ; T is the period, specifically the width of the image.

7. The method for extracting fracture occurrence parameters based on borehole images according to claim 1, characterized in that: The fracture occurrence parameters extracted in step S6 include the inclination and dip of the fracture structural surface.

Images