A method for identifying shale fracture development patterns based on electrical imaging logging images

Abstract

The invention discloses a shale fracture development pattern recognition method based on electrical imaging logging images, comprising the following steps: step 1: acquiring electrical imaging and performing bilateral filtering processing; step 2: binarizing the electrical imaging, then extracting the occurrence information of structural fractures through improved Hough transform, and counting the number of structural fractures and the opening degree per meter; step 3: graying and sharpening the image filtered in step 1, and counting the number of bedding strips n per meter; step 4: calculating the sum of the opening degrees of structural fractures per meter according to the opening degree obtained in step 2, and dividing the sum of the opening degrees by the number of bedding strips n per section to obtain σ; the maximum value of σ multiplied by n is recorded as k; the shale fracture development pattern recognition method based on electrical imaging logging images disclosed by the invention can finely pick up the fracture information on the imaging and accurately identify the fracture development pattern of the section, and provides a new idea for imaging fracture evaluation.

Description

Technical Field

[0001] The invention relates to the technical field of oil and gas exploration, and in particular to a shale fracture development pattern recognition method based on electrical imaging logging images. Background Art

[0002] The development of shale fractures is related to the preservation and migration of shale gas, and is also one of the important factors affecting fracturing transformation and profitable development. However, shale is extremely heterogeneous, and the types and combination patterns of developed fractures are complex and diverse. Due to the size, it is difficult to completely match the characteristics of a single or certain fractures with the logging curve response. Therefore, summarizing and identifying the main development patterns and characteristics of fractures in complex fracture development sections can better provide guidance for fracturing development. At present, relatively clear and obvious structural fractures and bedding information can be extracted through electrical imaging logging images, but there is still a lot of unclear fracture information that cannot be displayed on the image, which brings great challenges to the evaluation and analysis of shale gas reservoir fractures.

[0003] At present, the identification and prediction of shale fracture sections are mainly through special logging methods such as array acoustic waves. A large number of geologists have also established fracture prediction models based on conventional logging data with the help of signal analysis and deep learning. Conventional logging and array acoustic wave data are difficult to extract information such as fracture type, quantity, and occurrence. Therefore, the evaluation of fractures by the above data and methods is limited to qualitative evaluation of the presence or absence of fractures and the degree of fracture development. Because electrical imaging logging images have the advantage of high resolution, there have been a number of methods for extracting fracture occurrence information through electrical imaging logging images. However, many scholars have not achieved obvious results in extracting fracture information on electrical imaging logging images with complex or unclear information, especially dense bedding information. Moreover, the research on fractures in imaging often focuses on image processing, ignoring the type and combination characteristics of fractures. However, the different fracture types, quantities, and development characteristics of reservoir sections have a great impact on later fracturing transformation.

[0004] Existing methods such as patent application number 2019110288605, a method for identifying and evaluating fine-grained sedimentary continental shale fractures, include the following steps: S1: extract logging attributes, including wavelet high frequency attribute WHF, wavelet singular attribute WSV and prediction filter error attribute PFE; S2: construct shale fracture evaluation sub-parameters; S3: adopt the coefficient of variation weighting method to construct shale fracture evaluation comprehensive sub-parameters; S4: construct shale fracture evaluation comprehensive parameters; this method realizes the identification of shale fractures in a single well and the evaluation of the degree of fracture development in the fracture section, solving the problem of identifying fracture sections using conventional logging data. However, it does not involve information such as the type and occurrence of fractures, and only qualitatively evaluates whether there are fractures in a certain section and the degree of fracture development, which has limited guiding significance for later fracturing transformation.

[0005] For example, the patent application number is 2019103754535, which is an intelligent calculation method for quantitative parameters of fractures in oil-based mud wells based on ultrasonic imaging. On the basis of analyzing the characteristics of ultrasonic imaging fracture response, the fracture response is highlighted through filtering processing to improve the signal-to-noise ratio. Then, the fracture trajectory is automatically picked up through the image intelligent recognition algorithm and the fracture length, density and inclination are calculated. Then, the fracture width is automatically explored and statistically calculated along the normal direction of the fracture trajectory, and the fracture surface ratio is calculated based on the fracture length, width and density parameters. The core roll scanning images are automatically processed by manual or computer to obtain the core fracture parameters. The core fracture parameters calculated by ultrasonic imaging are calibrated using core fracture parameters to improve the accuracy and reliability of logging evaluation of complex unconventional reservoirs, thereby providing a basis for formulating reasonable and efficient development plans. Although this method realizes the extraction of fracture density and occurrence information, it does not involve the combined characteristics of the development of different fracture types. Summary of the invention

[0006] In view of the problems existing in the prior art, the present invention provides a shale fracture development pattern recognition method based on electrical imaging logging images.

[0007] The technical solution adopted by the present invention is:

[0008] A method for identifying shale fracture development patterns based on electrical imaging logging images comprises the following steps:

[0009] Step 1: Obtain electrical imaging and perform bilateral filtering;

[0010] Step 2: Binarize the electrical image, then extract the occurrence information of the structural fractures through the improved Hough transform, and count the number and opening of the structural fractures per meter;

[0011] Step 3: grayscale and sharpen the image filtered in step 1, and count the number of bedding strips per meter n;

[0012] Step 4: Calculate the sum of the openings of each meter of structural fractures based on the openings obtained in step 2. Divide the sum of the openings by the number of bedding strips in each section n to get σ; the maximum value of σ multiplied by n is recorded as k;

[0013] A rectangular coordinate system is constructed with n as the abscissa and σ as the ordinate; a curve σ=k / n is drawn in the rectangular coordinate system; the area formed by the abscissa, the ordinate and the curve of the rectangular coordinate system is divided into six areas corresponding to six fracture development modes; the sample points are projected into the rectangular coordinate system, and the fracture development mode is identified according to the projected position.

[0014] Furthermore, the bilateral filtering process in step 1 is as follows:

[0015]

[0016] Where: I is the input image, BF is the image after bilateral filtering, p is the center coordinate of the filter kernel window, q is the non-center coordinate of the filter kernel window, S is the spatial domain, W p is the normalization factor, is the spatial domain weight coefficient, is the range weight coefficient, I q is the pixel value corresponding to the q coordinate.

[0017] Furthermore, the process of extracting the occurrence information of the structural fractures in step 2 is as follows: extracting the occurrence information of the structural fractures by extracting the fixed periodic sinusoidal baseline through improved Hough transform;

[0018] The baseline is determined as follows:

[0019] S1: Divide the image into two equal parts, I1 and I2;

[0020] S2: Set a one-dimensional accumulator array ACC[y0] and initialize each element to 0; y0 is the baseline position of the sine curve;

[0021] S3: For any boundary point P(x p ,y p ), execute steps S4 to S7;

[0022] S4: For the vertical scan line x in I2 = x p Any edge point P′(x p ′,y p '), match P and P' into a point pair, and execute steps S5 to S7;

[0023] S5: Calculate the ordinate y of the end point of point pair P′ c ;

[0024] S6: with y c The value of the corresponding accumulator array element is increased by 1, ACC[y0]=ACC[y0]+1;

[0025] S7: loop until all point pairs in the image are visited;

[0026] S8: The peak position of the accumulator array is the baseline position of the sine line.

[0027] Furthermore, the sinusoidal curve is as follows:

[0028] y0=yA sin(ωx-β)

[0029] Where: x, y are the coordinates of the pixel points in the corresponding image space, A is the amplitude of the sine line, ω is the angular velocity, β is the initial phase, and y0 is the baseline position of the sine curve.

[0030] Furthermore, the crack opening calculation method is as follows:

[0031]

[0032] Where: w i is the number of pixels occupied by the crack in the vertical direction, m is the number of equally spaced positions, To obtain the average number of pixels covered by the crack in the longitudinal direction at m locations at equal intervals, w d The actual height of a single pixel.

[0033] Furthermore, the grayscale and sharpening process is as follows:

[0034] Grayscale is processed using the rgb2gray function in MATLAB; sharpening is processed using the sobel operator and imfilter function in MATLAB.

[0035] Furthermore, the six fracture development modes in step 4 are: mixed development of fractures and bedding, medium-density bedding, dense bedding, strong structural fracture and weak bedding, large-scale structural fracture, and dense section.

[0036] The technical solution adopted by the present invention is:

[0037] (1) The fracture pattern recognition method of the present invention can accurately pick up the fracture information on the image and accurately identify the fracture development pattern of the section, which provides a new idea for imaging fracture evaluation and a basis for optimizing the subsequent fracturing scheme;

[0038] (2) The present invention can achieve noise reduction and retain edge information through bilateral filtering, making the extraction of crack information more accurate;

[0039] (3) The present invention highlights the layer information on the image through image grayscale and sharpening processing, making human-computer interactive statistical layer information more convenient and accurate;

[0040] (4) The present invention uses a 1-meter segment as a sample point for research, focusing on and highlighting the pattern and characteristics of fracture development in this segment, and reducing the impact of errors in fracture information extraction. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 It is a schematic diagram of the process of the present invention.

[0042] Figure 2 This is a schematic diagram of the principle of extracting fixed-period sinusoidal lines using the improved Hough transform in Embodiment 1 of the present invention.

[0043] Figure 3The principle of extracting fixed periodic sinusoidal lines by improved Hough transform in Embodiment 1 of the present invention is applied to a simple image ( Figure 3 a structural fracture) and complex images ( Figure 3 b The effect diagram on the surface (with more bedding information).

[0044] Figure 4 These are the effect diagrams of the electrical imaging logging image (a), bilateral filtering (b), grayscale image sharpening (c) and binary image (d) in Example 1 of the present invention.

[0045] Figure 5 Schematic diagram of the fracture development pattern type in Example 1 of the present invention.

[0046] Figure 6 This is a conceptual diagram of the regional division of different fracture development modes in the σ and n coordinate systems in Example 1 of the present invention.

[0047] Figure 7 The sample points are projected into the σ and n coordinate systems in Example 1 of the present invention to determine the dividing boundary diagram of different fracture development modes.

[0048] Figure 8 This is a diagram showing the effect of applying the deep shale fracture development pattern recognition method in Example 1 of the present invention. DETAILED DESCRIPTION

[0049] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0050] like Figure 1 As shown, the electrical imaging image is obtained and bilateral filtering is performed on the electrical imaging image. The bilateral filtering process is as follows:

[0051]

[0052] Where: I is the input image, BF is the image after bilateral filtering, p is the center coordinate of the filter kernel window, q is the non-center coordinate of the filter kernel window, S is the spatial domain, W p is the normalization factor, is the spatial domain weight coefficient, is the range weight coefficient, I q is the pixel value corresponding to the q coordinate.

[0053] When the window has edge cases, The weight becomes larger, so that the edge information is retained. When the pixel values ​​in the window are smooth, the weight becomes smaller. The area is de-noised. The electro-imaging image and the image after bilateral filtering are shown in Figure 4 As shown in b.

[0054] Step 2: Binarize the electrical image, then extract the occurrence information of the structural fractures through the improved Hough transform, and count the number and opening of the structural fractures per meter;

[0055] The improved Hough transform changes the straight line equation into a sine curve equation through the point-line duality of the image space parameter space. The sine curve is:

[0056] y0=yA sin(ωx-β)

[0057] Where: x, y are the pixels in the corresponding image space, A is the amplitude of the sine line, ω is the angular velocity, β is the initial phase, y0 is the baseline position of the sine curve, and the baseline is determined as follows:

[0058] The baseline is determined as follows:

[0059] S1: Divide the image into two equal parts, I1 and I2;

[0060] S2: Set a one-dimensional accumulator array ACC[y0] and initialize each element to 0; y0 is the baseline position of the sine curve;

[0061] S3: For any boundary point P(x p ,y p ), execute steps S4 to S7;

[0062] S4: For the vertical scan line x in I2 = x p Any edge point P′(x p ′,y p '), match P and P' into a point pair, and execute steps S5 to S7;

[0063] S5: Calculate the ordinate y of the end point of point pair P′ c ;

[0064] S6: with y c The value of the corresponding accumulator array element is increased by 1, ACC[y0]=ACC[y0]+1;

[0065] S7: loop until all point pairs in the image are visited;

[0066] S8: The peak position of the accumulator array is the baseline position of the sine line.

[0067] The above process can quickly and effectively identify the baseline of the sine line of the structural fracture on the electrical imaging logging image ( Figure 3 a), but for images with more layers, this method cannot effectively extract the layer information (such as Figure 3 b).

[0068] The calculation formula of the opening of structural fracture is as follows:

[0069]

[0070] Where: w i is the number of pixels occupied by the crack in the vertical direction, m is the number of equally spaced positions, To obtain the average number of pixels covered by the crack in the longitudinal direction at m locations at equal intervals, w d is the actual height of a single pixel. In this embodiment, m is 8.

[0071] Step 3: grayscale and sharpen the image filtered in step 1, and count the number of bedding strips per meter n through human-computer interaction;

[0072] For unclear and complex images with many layers, grayscale and sharpening can highlight the layer information (such as Figure 4 c). It is more accurate to count the number of bedding strips per meter through human-computer interaction.

[0073] Grayscale is processed using the rgb2gray function in MATLAB;

[0074] The function is called as follows:

[0075] A_gray=rgb2gray(BF)

[0076] Where: BF is the image after bilateral filtering, and A_gray is the grayscale image.

[0077] The sharpening is processed using the sobel operator and imfilter function in MATLAB. The sharpening process is as follows:

[0078] f_sobel=fspecial(′sobel′)

[0079] B = imfilter(A_gray,f_sobel)

[0080] A_sobel=A_gray-B

[0081] Where: f_sobel is the generated filter function, A_gray is the grayscale image, B is the filtered information, and A_sobel is the grayscale image that highlights the details after sharpening.

[0082] The sum of the openings of the structural fractures per meter is calculated based on the openings obtained in step 2, and the sum of the openings is divided by the number of bedding strips in each section n to obtain σ; the maximum value of σ multiplied by n is recorded as k;

[0083] A rectangular coordinate system is constructed with n as the abscissa and σ as the ordinate; a curve σ=k / n is drawn in the rectangular coordinate system; the area formed by the abscissa, the ordinate and the curve of the rectangular coordinate system is divided into six areas corresponding to six fracture development modes; the sample points are projected into the rectangular coordinate system, and the fracture development mode is identified according to the projected position.

[0084] Through a large number of surveys and core and imaging observations, combined with the characteristics of structural fractures and bedding development and combination, deep shale fractures are divided into six categories according to their development patterns (such as Figure 5 As shown in the figure, they are: mixed development type of structural fractures and bedding (type I), medium-density bedding type (type II), dense bedding type (type III), strong structural fracture and weak bedding type (type IV), large-scale structural fracture type (type V), and dense segment (type VI). The specific process of the method for dividing the above six fracture development patterns is as follows:

[0085] S11: Divide the electrical imaging logging image into images of 1m each segment, and calculate the sum of the structural fracture apertures of each segment (1m) using the method of steps 1 to 3 (if there are j structural fractures in the 1m segment, the sum of the structural fracture apertures of the segment is ) and the number of bedding strips n in each segment.

[0086] S12: Calculate the value of each structural fracture opening and divide it by the number of bedding strips (σ):

[0087]

[0088] S13: Calculate the maximum value of each segment σ multiplied by n, denoted as k;

[0089] S14: According to the characteristics of σ and n for each of the six fracture development modes, a rectangular coordinate system with σ as the ordinate and n as the abscissa is established, and the region (n>0, σ>0, σ=k / n) is divided into six sub-regions (such as Figure 6 As shown). Each area corresponds to a fracture development mode (the six fracture development modes are: mixed development mode of structural fractures and bedding (type I), medium-density bedding (type II), density bedding (type III), strong structural fracture and weak bedding (type IV), large-scale structural fracture (type V), and dense section (type VI)). Among them, the n and σ values ​​of type I are relatively large, the n value of type II is relatively large, the σ value is relatively small, the n value of type III is the largest, and the σ value is very small, the n value of type IV is relatively small, and the σ value is very large, the n value of type V is the smallest, and the σ value is the largest, and the n value and σ value of type VI are both the smallest when the fractures are underdeveloped.

[0090] S15: Take each segment as a sample point and project it into a rectangular coordinate system with σ as the ordinate and n as the abscissa (e.g. Figure 7As shown in Table 1). The projection of each sample point should be within the limits of the σ=k / n curve (n>0,σ>0). Based on the characteristics that each of the six fracture development modes has its own σ and n, the limits of σ and n for each fracture development mode are determined by combining the electrical imaging logging images (as shown in Table 1). The imaging fracture development mode recognition method is established.

[0091] Table 1. Boundaries of different fracture development patterns

[0092]

[0093] The shale fracture development pattern recognition method based on electrical imaging logging images disclosed in the present invention can accurately pick up the fracture information on the image and accurately identify the fracture development pattern of the section, providing a new idea for imaging fracture evaluation; the method of the present invention can achieve noise reduction and retain edge information through bilateral filtering, making the extraction of fracture information more accurate. By processing the image grayscale and sharpening, the bedding information on the image is highlighted, making the human-computer interactive statistics of bedding information more convenient and accurate. Taking the 1m section as a sample point for research, the pattern and characteristics of the fracture development in this section are focused on and highlighted, which reduces the influence of errors in the extraction of complex fracture and bedding information, and provides a basis for the optimization of the later fracturing scheme.

Claims

1. A method for identifying shale fracture development patterns based on electrical imaging logging images, characterized in that: The following steps are involved: Step 1: Obtain electrical imaging and perform bilateral filtering; Step 2: Binarize the electrical image, and then extract the occurrence information of the structural fractures through the improved Hough transform, and count the number and opening of the structural fractures per meter; The process of extracting the occurrence information of the structural fractures is as follows: extract the occurrence information of the structural fractures by extracting the fixed period sine line baseline through the improved Hough transform; The baseline is determined as follows: S1: Divide the image into two equal parts, I1 and I2; S2: Set a one-dimensional accumulator array ACC[y0] and initialize each element to 0; y0 is the baseline position of the sine curve; S3: For any boundary point P(x p ,y p ), execute steps S4 to S7; S4: For the vertical scan line x in I2 = x p Any edge point P′(x p ′,y p '), match P and P' into a point pair, and execute steps S5 to S7; S5: Calculate the ordinate y of the end point of point pair P′ c ; S6: with y c The value of the corresponding accumulator array element is increased by 1, ACC[y0]=ACC[y0]+1; S7: loop until all point pairs in the image are visited; S8: The peak position of the accumulator array is the baseline position of the sine line; Step 3: grayscale and sharpen the image filtered in step 1, and count the number of bedding strips per meter n; Step 4: Calculate the sum of the openings of structural fractures per meter based on the openings obtained in step 2. Divide the sum of the openings by the number of bedding strips in each section n to get σ; the maximum value of σ multiplied by n is recorded as k; A rectangular coordinate system is constructed with n as the abscissa and σ as the ordinate; a curve σ=k / n is drawn in the rectangular coordinate system; the area formed by the abscissa, the ordinate and the curve of the rectangular coordinate system is divided into six areas corresponding to six fracture development modes; the sample points are projected into the rectangular coordinate system, and the fracture development mode is identified according to the projected position.

2. The method for identifying shale fracture development patterns based on electrical imaging logging images according to claim 1, characterized in that: The bilateral filtering process in step 1 is as follows: Where: I is the input image, BF is the image after bilateral filtering, p is the center coordinate of the filter kernel window, q is the non-center coordinate of the filter kernel window, S is the spatial domain, W p is the normalization factor, is the spatial domain weight coefficient, is the range weight coefficient, I q is the pixel value corresponding to the q coordinate.

3. The method for identifying shale fracture development patterns based on electrical imaging logging images according to claim 1, characterized in that: The sine curve is as follows: y0=y-Asin(ωx-β) Where: x, y are the coordinates of the pixel points in the corresponding image space, A is the amplitude of the sine line, ω is the angular velocity, β is the initial phase, and y0 is the baseline position of the sine curve.

4. The method for identifying shale fracture development patterns based on electrical imaging logging images according to claim 1, characterized in that: The crack opening calculation method is as follows: Where: w i is the number of pixels occupied by the crack in the vertical direction, m is the number of equally spaced positions, To obtain the average number of pixels covered by the crack in the longitudinal direction at m locations at equal intervals, w d The actual height of a single pixel.

5. The method for identifying shale fracture development patterns based on electrical imaging logging images according to claim 1, characterized in that: The grayscale and sharpening process is as follows: Grayscale is processed using the rgb2gray function in MATLAB; sharpening is processed using the sobel operator and imfilter function in MATLAB.

6. The method for identifying shale fracture development patterns based on electrical imaging logging images according to claim 1, characterized in that: The six fracture development modes in step 4 are: fracture and bedding mixed development type, medium-density bedding type, dense bedding type, strong structural fracture and weak bedding type, large-scale structural fracture type, and dense section.

Images