Fracture occurrence parameter extraction method based on borehole image

Abstract

A fissure occurrence parameter extraction method based on a borehole image relates to the technical field of geological resources and geological engineering, and mainly comprises the following steps: collecting and importing an original borehole image of a drilling area, enhancing and denoising the original borehole image, binarizing the image through dynamic area division, and obtaining a fissure occurrence parameter extraction result; detecting the edge of the fracture structural plane by using a Canny operator for carrying out adaptive threshold selection based on a borehole image gradient histogram, carrying out sine curve fitting on the fracture structural plane, and extracting fracture occurrence parameters on a fitted fracture sine curve; based on the engineering drilling image, the rock mass structural surface can be automatically identified, occurrence parameters can be extracted, previous manual intervention and operation are replaced, and the efficiency of analyzing geological conditions and underground water migration conditions in geological engineering is improved.

Description

Technical Field

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

[0002] Structural surface is one of the basic information that needs to be obtained in geological engineering surveys. Various fractured structural surfaces in rock masses are also important places and channels for groundwater storage and migration. Their dip and inclination information will affect the flow rate and direction of groundwater. Correct analysis of regional structural surface information is the basis for engineering geological activities.

[0003] With the development of digital panoramic borehole camera technology, extracting structural surface parameters through borehole wall images has become an important means. In borehole images, the interpretation of structural surface information is usually done manually, which is cumbersome and time-consuming, and different people may have different understandings and interpretations of the same borehole wall image. In addition, in existing automatic recognition technologies, such as wavelet analysis, clustering algorithm, three-dimensional Hough transform and other methods, although they have certain effects on automatic recognition of structural surfaces, the edge detection effect of fracture structural surfaces in images is relatively fuzzy, affecting the accuracy of final parameter extraction. In order to extract the parameters of rock mass structural surfaces in a practical, convenient and accurate manner, it is particularly important to develop an efficient method to obtain key information of structural surfaces and realize the automatic extraction of structural surface parameters.

[0004] Since the traces of the fracture structure surface after cutting and drilling are unfolded in the plane image of the hole wall, the existing technology usually relies on extracting the edge information of the structure surface and the parameter space mapping theory to identify the fracture structure surface, so as to obtain the parameter information such as the direction, inclination and inclination of the fracture structure surface. In this process, the Canny operator is often needed. The Canny operator is a commonly used method for detecting image edges, which mainly includes image smoothing and denoising, calculating gradient amplitude and direction, non-maximum suppression, dual threshold detection and edge connection. Among them, the dual threshold detection, the traditional approach is to set two high and low thresholds T H , T L , the gradient value is T H The larger point is a strong edge point, which is larger than T H Smaller than T LLarge points 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 the selection of these two thresholds. They are often set manually and are not adaptive and robust. When applied to borehole images, due to the complexity of the geological environment, borehole images often have a lot of noise interference, such as broken rock layers, dim light, and interlayer cementation. As a result, it is difficult for manually set threshold parameters to adapt to different borehole images, and continuous trial and error is required to repeatedly correct the detection results. Summary of the invention

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

[0006] In order to solve at least one of the above technical problems, the technical solution provided by the present invention is a method for extracting fracture occurrence parameters based on borehole images, comprising the following steps:

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

[0008] Step S2: enhancing and denoising the original drilling image;

[0009] Step S3: binarizing the image by dynamic region segmentation;

[0010] 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;

[0011] Step S5: performing sinusoidal curve fitting on the fracture structure surface;

[0012] Step S6: extracting fracture occurrence parameters from the fitted fracture sinusoidal curve.

[0013] The technical effects of the present invention are:

[0014] 1. Based on engineering drilling images, the present invention can automatically identify rock structure surfaces and extract occurrence parameters, replacing previous manual intervention and operation, and improving the efficiency of analyzing geological conditions and groundwater migration in geological engineering.

[0015] 2. Compared with the traditional binarization method, the present invention uses dynamic area division to realize the binarization of the borehole image, which will more accurately retain the fracture structure surface information when applied to the borehole image.

[0016] 3. When detecting the edge of a crack structure surface, the present invention 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 edge in the borehole image.

[0017] 4. The present invention utilizes the traces left after the structural surface is cut and drilled, and the hole wall image is unfolded into a plane image to show the characteristics of a sine curve. The characteristics of the sine function are directly used to fit the crack structural surface, which improves the detection efficiency and accuracy compared with conventional curve fitting.

[0018] 5. The method of the present invention has a fast running speed, good controllability during parameter extraction, and can provide real-time feedback on recognition effects. The precision and accuracy of recognition can be controlled by adjusting parameters to facilitate the use of subsequent data. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without paying creative work.

[0020] Figure 1 It is a specific flow chart of the present invention;

[0021] Figure 2 Schematic diagrams of three situations of cutting drilling holes on an inclined structural surface when the drilling image of the present invention is unfolded into a plane image;

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

[0023] Figure 3 b is a diagram showing the result of detecting the edge of the crack structure surface using the improved Canny operator in the present invention;

[0024] Figure 3 c is a graph showing the search and detection results of the sinusoidal baseline y0 in the present invention;

[0025] Figure 3 Figure d is a sinusoidal curve fitting effect diagram of the structural surface in the present invention. DETAILED DESCRIPTION

[0026] The present invention will be further described in detail below in conjunction with the embodiments and drawings.

[0027] In order to make the purpose, 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 in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection 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 invention claimed for protection, but merely represents selected embodiments of the present invention.

[0028] Example:

[0029] See also Figure 1 A method for extracting fracture occurrence parameters based on borehole images comprises the following steps:

[0030] Step S1: Collect the original borehole image of the drilling area. The original borehole image includes the borehole image and the plane expansion image after the borehole image is expanded. Figure 2 and Figure 3 As shown in a;

[0031] Step S2: Enhance and denoise the original borehole image, wherein the enhancement and denoising method described in step S2 is an adaptive median filter. Different from the conventional median filter, the adaptive median filter dynamically adjusts the size of the window according to the variance of the grayscale values ​​of the pixels in the window. When the variance of the grayscale value in the window is large, it means that the window may contain noise points. At this time, the filter will increase the size of the window to include more pixels to calculate the median, thereby better eliminating the noise. On the contrary, when the variance of the grayscale value in the window is small, it means that the window is mainly the real features of the image. At this time, the filter will reduce the size of the window to retain more image details.

[0032] Step S3: binarizing the image by dynamic region segmentation;

[0033] The steps of binarizing the image by dynamic area division in step S3 are:

[0034] 1) Calculate the grayscale histogram and local variance of the image. Areas with large local variance often contain more details or target boundaries. Determine whether there are multiple peaks, the location and width of the peaks, and use the valley value between the peaks of the grayscale histogram as the threshold;

[0035] 2) Start scanning with a small window from the upper left corner of the image. Scan with a smaller window (such as 3x3 or 5x5 pixels). 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. Scan the entire image in sequence to complete the division. When the grayscale difference is greater than the threshold, it can be considered that the area division spans the background and boundaries of different targets, and the scanning window needs to be adjusted until the grayscale difference is less than the threshold.

[0036] 3) Use the watershed algorithm to segment the divided area, and perform binarization on the segmented area. For areas with relatively uniform grayscale distribution, the threshold is determined according to the valley value of the histogram. For areas with uneven grayscale distribution, the threshold is determined according to the local mean and standard deviation. 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.

[0037] The threshold of the area with uniform grayscale distribution after segmentation is the valley value between the peaks of the grayscale histogram, and the threshold of the area with uneven grayscale distribution after segmentation is determined according to the local mean and standard deviation.

[0038] Among them, the threshold of the area with uneven grayscale distribution after segmentation is shown in formula (1):

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

[0040] In formula (1), T is the threshold of the area with uneven grayscale distribution after segmentation, μ is the local mean of the grayscale of the area, σ is the standard deviation of the grayscale of the area, and k is an empirical coefficient with a value between 0.5 and 2. The specific value is obtained through drilling diagram testing and is usually 1.

[0041] 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;

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

[0043] Combination Figure 3 In b, the steps of detecting the edge of the fracture structure surface using the Canny operator with adaptive threshold selection based on the gradient histogram of the drilling image are as follows:

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

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

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

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

[0048] 3) Determine the trough in the histogram. The low-amplitude gradient is mainly generated by noise. The trough represents the area with relatively few pixels. Therefore, the first-order difference of the gradient amplitude histogram is calculated according to formula (3):

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

[0050] 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 interval;

[0051] Then, the second-order difference is calculated based on the first-order difference according to formula (4):

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

[0053] 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.

[0054] When D(k)<0 and D2(k)>0, it is considered to be a trough. At the same time, the sizes of D2(k) corresponding to multiple troughs are compared, and the position with the largest D2(k) is selected as the most likely trough.

[0055] 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. The part to the right of the trough is mainly the edge-related gradient amplitude, and the high threshold is determined according to the distance from the trough to the high amplitude part and the gradient amplitude change trend;

[0056] Among them, 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. As the amplitude increases, the value of the second-order difference gradually decreases. When the fluctuation between its value and the previous and next values ​​is less than 0.5, it means that the growth rate of the gradient amplitude has slowed down significantly, and the gradient amplitude corresponding to this position is used as the high threshold.

[0057] 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 strong edge points are connected to the weak edge points adjacent to the strong edge points to construct the edge contour of the complete fracture structure surface.

[0058] Step S5: performing sinusoidal curve fitting on the fracture structure surface;

[0059] Combination Figure 3 Medium c and Figure 3 In step d, the steps of performing sinusoidal curve fitting on the fracture structure surface in step S5 are:

[0060] The image coordinate space is mapped to the parameter space. In this process, the 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 duality between points and lines in two different coordinate systems. The crack trajectory satisfies formula (4):

[0061]

[0062] In formula (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] The position of the baseline y0 is detected with the help 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 vertical coordinates of the midpoints of the two paired points and add 1 to the corresponding accumulator array element value; continue the whole process until all point pairs in the image have been 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 sinusoidal curve, the amplitude and initial phase are extracted through parameter space transformation. According to the transformation of formula (4), formula (5) is obtained:

[0065]

[0066] in, The range of is 0 to T; the range of A is 0 to A Max ;

[0067] Create a two-dimensional accumulator array and assign an initial value of 0; for the detected y0, traverse each boundary point P(x,y) in the image; for For each value of (0 to T), A can be calculated by formula (5). Max If the number of pixels between the boundary pixels and the peak value is greater than y0, the count of the corresponding element in the accumulator array is increased by 1. Repeat the above steps until all boundary pixels are processed, and find the peak position in the accumulator array. Repeat the above steps until all y0 points are processed.

[0068] Step S6: extracting fracture occurrence parameters from the fitted fracture sinusoidal curve;

[0069] The fracture occurrence parameters include the inclination and dip of the fracture structural surface. The inclination angle β is obtained by the inverse tangent function of the vertical distance H between the crest and the trough of the fracture sine curve obtained in step S5 (i.e., twice the amplitude A of the sine curve) and the borehole diameter D, as shown in formula (6):

[0070]

[0071] The inclination can be directly read from the direction 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 pointed out that the terms "upper", "lower", "front", "back", "left", "right", "top", "bottom", "inside" and "outside" etc. indicating directions or positional relationships are based on the directions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and should not be understood as a limitation on 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 a person skilled in the art within the technical scope disclosed in the embodiments of the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on 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.

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, perform binarization on the segmented area, assign a value of 255 to pixels whose grayscale value is greater than or equal to the regional threshold, and assign a value of 0 to pixels whose grayscale value is less than the regional threshold; 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 region 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 region with uneven grayscale distribution after segmentation, μ is the local mean of the regional grayscale, σ is the standard deviation of the regional grayscale, and k is an empirical coefficient with a value between 0.5 and 2.

6. The method for extracting fracture occurrence parameters based on borehole images according to claim 1, characterized in that: The steps 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 described in step S4 are: 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, wherein 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 interval; 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 strong edge points are connected to the weak edge points adjacent to the strong edge points to construct the edge contour of the complete fracture structure surface.

7. The method for extracting fracture occurrence parameters based on borehole images according to claim 6, characterized in that: 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.

8. 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): 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; ω = 2π / T; T is the period, specifically the width of the image.

9. 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