Automatic modeling method for multi-layered shale reservoirs based on layer identification

Abstract

The application discloses a multi-layered shale reservoir automatic modeling method based on bedding identification, specifically comprising the following steps: obtaining core samples from a shale reservoir, recording interval information, and obtaining core images; carrying out denoising processing on the core images and carrying out contrast enhancement; using a Canny edge detection algorithm to identify bedding edges of the core images; using a Hough line transformation algorithm to detect straight lines and identify bedding structure lines; judging whether the lines are straight by calculating the change of adjacent pixel points on the y-axis, calculating the depth of the lines according to the pixel coordinates of the lines; creating planes representing different beddings according to the depth values of each bedding, generating color of brightness and darkness, and further obtaining a shale multi-layered reservoir model. The automatic modeling method significantly improves the accuracy and efficiency of shale reservoir modeling, accurately identifies shale bedding characteristics by combining the bedding identification technology, and provides more detailed and reliable basic data for reservoir modeling.

Description

Technical Field

[0001] This invention belongs to the field of unconventional shale oil and gas exploration and development technology, specifically involving an automated modeling method for multi-layered shale reservoirs based on bedding identification. Background Technology

[0002] In the exploration and development of unconventional oil and gas resources, accurate reservoir characterization and modeling are crucial steps, directly impacting reservoir distribution, development potential assessment, and the effectiveness of engineering design. However, shale oil and gas resources differ from conventional oil and gas resources; their reservoirs typically exhibit complex heterogeneity and multi-layered structures, which limits the application of traditional reservoir modeling methods in shale reservoirs.

[0003] The heterogeneity of shale reservoirs is typically manifested in multi-scale, multi-type bedding structures. These bedding structures often significantly impact reservoir fluid occurrence, compressibility during fracturing, and fracture propagation morphology. Existing reservoir modeling methods usually rely on geophysical and conventional logging data. First, well information is used to establish lithofacies and bedding development patterns. Then, different characteristic parameters are assigned to these patterns to create a three-dimensional model. However, when faced with the complex bedding structures of shale reservoirs, these methods often suffer from insufficient accuracy in bedding identification and fine modeling, failing to accurately capture these microscopic bedding features, easily leading to model distortion and decreased accuracy. Furthermore, traditional methods rely heavily on manual intervention in practical applications; expert experience and judgment play a decisive role in model accuracy, resulting in significant subjectivity and uncertainty. How to better handle bedding features of different scales and types during bedding identification, seamlessly integrate bedding identification results with reservoir modeling, improve modeling efficiency while ensuring modeling accuracy, and ensure model reliability while reducing manual intervention, are the key and challenging aspects of multi-bedding shale modeling.

[0004] Chinese Patent (Application No.: 202310320605.8, Publication No.: CN116299738A, Publication Date: 2023-03-29) discloses a method and apparatus for identifying shale bedding based on well logging curves. This method interpolates well logging data points and extracts frequency data using discrete Fourier transform to determine the degree of shale bedding development. By setting the relationship between frequency ranges and the degree of bedding development, shale bedding can be accurately identified. While this method can identify bedding relatively accurately, it suffers from insufficient accuracy and inability to handle multi-scale and multi-type bedding structures when dealing with complex bedding structures. This is particularly problematic in practical applications of shale reservoir modeling, where it can easily lead to modeling errors and uncertainties.

[0005] Chinese Patent (Application No.: 202210221224.X, Publication No.: CN114359569A, Publication Date: 2022-03-09) discloses a method, apparatus, device, and storage medium for identifying rock bedding. This method, based on multiple sets of CT slice images, obtains the dual-energy index of the rock through scanning at different energy levels, thereby determining the rock's bedding structure. While this method demonstrates good accuracy in identifying multi-layered bedding rocks, it relies on the accuracy and scanning quality of the CT image data, and still faces certain challenges in identifying the details of complex shale reservoir structures, especially small-scale bedding.

[0006] In summary, the existing technologies have the following main drawbacks: 1) Data processing methods based on well logging curves are not accurate enough for identifying complex bedding structures and cannot effectively capture multi-scale and multi-type bedding features; 2) Existing bedding identification methods have high requirements for data quality and limited ability to identify small-scale micro-bedding, making them difficult to be widely used in actual reservoir modeling; 3) The application of existing methods in automated modeling processes is not yet mature and still requires significant human intervention, resulting in greater uncertainty and subjectivity in the model and making it difficult to achieve efficient and accurate reservoir modeling. Summary of the Invention

[0007] The purpose of this invention is to provide an automated modeling method for multi-layered shale reservoirs based on bedding identification, which significantly improves the accuracy and efficiency of shale reservoir modeling.

[0008] The technical solution adopted in this invention is an automated modeling method for multi-layered shale reservoirs based on bedding identification, which is implemented according to the following steps: Step 1: Obtain core samples from shale reservoirs, record the stratigraphic information of the core samples, including the top and bottom depths, and obtain core images; Step 2: Denoise the core image and enhance its contrast. Step 3: Use the Canny edge detection algorithm to identify the bedding edges of the core image and identify the edge pixels in the core image; Step 4: Use the Hough line transform algorithm to detect straight lines and identify the layered structure lines in the image; Step 5: Determine whether the line is straight by calculating the changes of adjacent pixels on the y-axis, calculate the depth of the line based on the pixel coordinates, and save the layer depth value to an Excel file; Step 6: Obtain the depth value of each bedding layer from the Excel file, create planes representing different bedding layers, and generate colors of varying brightness to obtain a multi-bedding shale reservoir model.

[0009] The invention is further characterized in that, Step 3 specifically involves: Step 3.1: Convert the core image to a grayscale image to eliminate color information; The grayscale conversion formula is shown in equation (1): (1); in, , , These represent the pixel values ​​of the red, green, and blue channels in the core image, respectively. Step 3.2: Use the Sobel operator to calculate the grayscale image in the horizontal direction. and vertical direction The changes in are used to calculate the gradient of the image; Step 3.3: Apply the non-maximum suppression method to retain only local gradient maxima as potential edges; Step 3.4: Determine the final edge through dual thresholding and hysteresis edge tracking.

[0010] In step 3.2, the Sobel operator performs convolution operations on the image using the following convolution kernel; Horizontal Sobel core, S x As shown in equation (2): (2); Vertical Sobel core, S y As shown in equation (3): (3); For each pixel, the gradient in the horizontal and vertical directions is calculated using a convolution operation, as shown in equation (4): (4); in, i, j For each pixel in a grayscale image; gradient magnitude The gradient direction ▽ represents the edge intensity of a pixel in the image, and the gradient direction ▽ represents the orientation of the edge. The calculation method is shown in Equation (5) and Equation (6): (5); (6).

[0011] Step 3.3 involves the following specific steps: Based on the gradient direction ▽, the gradient magnitude of each pixel in the image is compared with the neighboring pixels along the gradient direction. If the gradient magnitude of the current pixel is greater than the pixel values ​​in both of its neighboring directions, the value of the pixel is retained; otherwise, its value is set to 0.

[0012] Step 3.4 involves the following specific steps: Set two thresholds: high threshold and low threshold If the gradient magnitude of the pixel If the pixel's gradient magnitude is not high, it is considered a strong edge and marked as an edge; if the pixel's gradient magnitude is low, it is considered a strong edge. If the pixel's gradient magnitude is not an edge, it is considered not to belong to the edge and is marked as non-edge; if the pixel's gradient magnitude is... If the pixel is weak, it will be considered a weak edge and will only be preserved as an edge if it is connected to a strong edge; otherwise, it will be suppressed. The process of hysteresis edge tracking is as follows: For each weak edge pixel, check its neighboring pixels. If there is a strong edge pixel in the neighborhood of the pixel, mark it as an edge; if there is no strong edge pixel in the neighborhood, remove it; for weak edge pixels marked as edges, continue to check its neighboring pixels; if there are more weak edge pixels in the neighborhood and they are connected to the marked strong edges, mark these weak edge pixels as edges as well; continue propagation until all weak edge pixels have been checked. Ultimately, all strong edge pixels are preserved; only weak edge pixels connected to strong edges are preserved, while other weak edges are suppressed; the final edge image contains all strong edges as well as the weak edges connected to them.

[0013] Step 4 specifically involves: Step 4.1: Take the edge detection result from step 3 as input to obtain a binary image, where the edge pixel value is 1 and the non-edge pixel value is 0. Step 4.2, Select the angle range and distance range and initialize an accumulator. To store the cumulative value of each parameter combination; the angle range is The distance range is ; Step 4.3: For each edge point in the image, the Hough transform is used to convert it into a curve in the parameter space. Finally, the peak value in the parameter space is found by accumulating these curves, thereby identifying the straight line. Step 4.4, will Combined votes are added to the accumulator; each element in the accumulator This indicates the number of votes for this parameter combination, i.e., how many edge points support this line; Step 4.5: Search for local peaks in the accumulator; the parameter combinations corresponding to the local peaks. The method involves finding the straight line in the image, identifying the maximum value in the accumulator, and then using the resulting polar coordinate parameters. Transform back to the image coordinate system to obtain the corresponding straight line in the image.

[0014] In step 5, the straightness of the line is determined by calculating the changes in adjacent pixels along the y-axis. Specifically: Suppose the current line contains multiple pixels. Then for two adjacent pixels and Calculate the difference between its y-axis coordinates, as shown in equation (8); (8); Set a threshold This is used to determine if the height difference is too large; if If the two points are on the same horizontal line, then they are considered to belong to the same horizontal line; otherwise, they do not belong to the same horizontal line. In the entire detected straight line, if the height difference between adjacent points is less than the threshold, the straight line is considered to be straight; otherwise, the straight line is not straight. The depth of the line is calculated based on its pixel coordinate y1. The depth calculation formula is shown in equation (9). Depth = Top depth + ((y1 / image height) * (bottom depth - top depth)) (9); During the depth data calculation process, a depth difference threshold is set. If the depth difference between two adjacent bedding planes is less than the threshold, they are considered to belong to the same bedding plane, and their average value is calculated as a single bedding plane.

[0015] The beneficial effects of this invention are: This invention employs image processing and data analysis techniques to accurately identify bedding features in core images and construct detailed three-dimensional reservoir models based on these features. The automated process reduces manual intervention, providing an automated modeling method for multi-layered shale reservoirs that combines bedding identification. This method can meet the needs of reservoir models with complex multi-layered bedding structures, significantly improving the accuracy and efficiency of shale reservoir modeling and ensuring the accuracy and practicality of the models. In summary, the method of this invention not only improves the accuracy and efficiency of shale reservoir modeling but also adapts to the modeling needs of complex geological structures, providing reliable technical support for the exploration and development of shale oil and gas resources. Attached Figure Description

[0016] Figure 1 This is a flowchart of the automated modeling method for multi-layered shale reservoirs based on bedding identification according to the present invention. Figure 2 Core images obtained by scanning the core using a high-resolution imaging device; Figure 3The image shows the result of Gaussian denoising and contrast enhancement using an adaptive histogram equalization algorithm on the core image. Figure 4 This is the result of edge detection applied to an image that has undergone denoising and contrast enhancement using the Canny edge detection algorithm. Figure 5 This is the result of using the Hough line transform algorithm to detect straight lines in a core image. Figure 6 An Excel spreadsheet obtained after calculating the bedding depth and adjusting for similar bedding depths; Figure 7 The image shows the results of building a realistic multi-layered shale reservoir model based on depth values ​​and visualizing it using PyVista. Detailed Implementation

[0017] The present invention will now be described in detail with reference to specific embodiments and accompanying drawings.

[0018] Example 1 This invention presents an automated modeling method for multi-layered shale reservoirs based on bedding identification. Specifically, the method involves: obtaining core samples from the shale reservoir, recording the bedding information of the core samples (including top and bottom depths), and obtaining a core image; denoising the core image and enhancing its contrast; using the Canny edge detection algorithm to identify bedding edges in the core image; using the Hough line transform algorithm to detect straight lines and identify bedding structure lines in the image; determining the straightness of lines by calculating the changes in adjacent pixels along the y-axis, calculating the depth of the lines based on their pixel coordinates, and saving the bedding depth values ​​to an Excel file; obtaining the depth value of each bedding layer from the Excel file, creating planes representing different bedding layers, and generating colors to represent brightness and darkness, thereby obtaining a multi-layered shale reservoir model.

[0019] Example 2 This invention provides an automated modeling method for multi-layered shale reservoirs based on bedding identification, which is implemented according to the following steps: Step 1: Obtain core samples from shale reservoirs, record the stratigraphic information of the core samples, including the top and bottom depths, and scan them using high-resolution imaging equipment to obtain high-quality core images; Step 2: Denoise the core image using Gaussian convolution and enhance its contrast using the Adaptive Histogram Equalization (CLAHE) algorithm to further improve image quality; When denoising core images, Gaussian denoising, also known as Gaussian blurring, is a commonly used image processing technique to reduce noise while preserving edge information as much as possible. In practice, Gaussian denoising involves convolving the image with a Gaussian kernel matrix. The Gaussian kernel matrix is ​​a weight matrix generated based on a Gaussian distribution function, with the highest weight at the center and gradually decreasing outwards. For each image pixel, its new value is the weighted average of that pixel and its neighboring pixels, where the weights are determined by the Gaussian kernel matrix.

[0020] Contrast Limited Adaptive Histogram Equalization (CLAHE) is a technique for enhancing image contrast, particularly suitable for processing images with insufficient local contrast. Compared to traditional histogram equalization, CLAHE can enhance contrast in different local regions of the image, avoiding the problems of over-enhancement or loss of detail that may occur with global histogram equalization.

[0021] Step 3: Use the Canny edge detection algorithm to perform bedding edge recognition on the core image, accurately identifying edge pixels in the core image; specifically: Step 3.1: Convert the core image into a grayscale image to eliminate color information and simplify calculations; The RGB image is converted to a grayscale image using the grayscale conversion formula, as shown in equation (1): (1); in, , , These represent the pixel values ​​of the red, green, and blue channels in the core image, respectively. Step 3.2, use the Sobel operator to calculate the grayscale image in the horizontal direction ( ) and vertical direction ( The changes in ) are used to calculate the gradient of the image; The Sobel operator performs convolution operations on the image using the following convolution kernels; Horizontal Sobel core, S x As shown in equation (2): (2); Vertical Sobel core, S y As shown in equation (3): (3); For each pixel, the gradient in the horizontal and vertical directions is calculated using a convolution operation, as shown in equation (4): (4); in, i, j For each pixel in a grayscale image; gradient magnitude The gradient direction ▽ represents the edge intensity of a pixel in the image, and the gradient direction ▽ represents the orientation of the edge. The calculation method is shown in Equation (5) and Equation (6): (5); (6); Step 3.3: Apply the non-maximum suppression method, retaining only local gradient maxima as potential edges. The specific steps are as follows: Based on the gradient direction ▽, the gradient magnitude of each pixel in the image is compared with the neighboring pixels along the gradient direction. If the gradient magnitude of the current pixel is greater than the pixel values ​​in both of its neighboring directions, the value of the pixel is retained; otherwise, its value is set to 0. Step 3.4: Determine the final edge through dual thresholding and hysteresis edge tracking. The specific steps are as follows: Set two thresholds: high threshold and low threshold If the gradient magnitude of the pixel If the pixel's gradient magnitude is not high, it is considered a strong edge and marked as an edge; if the pixel's gradient magnitude is low, it is considered a strong edge. If the pixel's gradient magnitude is not an edge, it is considered not to belong to the edge and is marked as non-edge; if the pixel's gradient magnitude is... If the pixel is weak, it will be considered a weak edge and will only be preserved as an edge if it is connected to a strong edge; otherwise, it will be suppressed. Hysteresis edge tracking is a process used to connect weak and strong edges; its purpose is to determine which weak edges should be retained as edges and which should be suppressed. For each weak edge pixel, its neighborhood (including the top, bottom, left, right, and eight diagonal neighboring pixels) is checked. If a strong edge pixel exists in the neighborhood of this pixel, it is marked as an edge; if no strong edge pixel exists in the neighborhood, it is removed. For a weak edge pixel marked as an edge, its neighboring pixels are checked again. If there are more weak edge pixels in the neighborhood that are connected to the already marked strong edges, these weak edge pixels are also marked as edges. This process continues until all weak edge pixels have been checked. Ultimately, all strong edge pixels are preserved; only weak edge pixels connected to strong edges are preserved, while other weak edges are suppressed; the final edge image contains all strong edges as well as the weak edges connected to them. Step 4: Based on the edge detection results, use the Hough line transform algorithm to detect straight lines and identify the layered structure lines in the image; Step 4.1: Take the edge detection result from step 3 as input to obtain a binary image, where the edge pixel value is 1 and the non-edge pixel value is 0. Step 4.2, select a suitable angle range and distance range and initialize an accumulator. (Using a two-dimensional array) to store the cumulative value of each parameter combination. The angle range is... The distance range is .

[0022] Step 4.3: For each edge point in the image, the Hough transform is used to convert it into a curve in the parameter space. Finally, the peak value in the parameter space is found by accumulating these curves, thereby identifying the straight line.

[0023] The Hough transform represents a straight line in polar coordinates, as shown in equation (7): (7); in, This is the distance offset from the line (the perpendicular distance from the origin to the line). The angle between the line and the horizontal axis. and These are the coordinates of a point in the image. In polar coordinates, each edge point in the image corresponds to a curve.

[0024] Step 4.4, take these Combine the votes into the accumulator; specifically, for each Edge points, calculate a series of The value corresponds to a line, and the count value of these line parameter combinations in the accumulator is increased. Each element in the accumulator... This represents the "votes" for this parameter combination, i.e., how many edge points support the line; Step 4.5: Search for local peaks in the accumulator; the parameter combinations corresponding to the local peaks. In an image, a straight line can be found by identifying the maximum value (or local maximum) in the accumulator. The resulting polar coordinate parameters... Transform back to the image coordinate system to obtain the corresponding straight line in the image; Step 5: For the line detection results obtained in Step 4, determine whether the line is straight by calculating the changes of adjacent pixels on the y-axis (i.e., the vertical direction). Suppose the current line contains multiple pixels. Then for two adjacent pixels and Calculate the difference between its y-axis coordinates, as shown in equation (8); (8); Set a threshold This is used to determine if the height difference is too large; if If the two points are on the same horizontal line, then they are considered to belong to the same horizontal line; otherwise, they do not belong to the same horizontal line. In the entire detected straight line, if the height difference between adjacent points is less than the threshold, the straight line is considered to be straight; otherwise, the straight line is not straight. threshold This can be determined experimentally; for example, if the longitudinal resolution of the core image is high, a smaller threshold (e.g., 2-5 pixels) can be selected; if the resolution is low, a slightly larger threshold (e.g., 5-10 pixels) can be selected.

[0025] Calculate the depth of the lines based on their pixel coordinates and save the layer depth values ​​to an Excel file; Step 6: Obtain the depth value of each bedding layer from the Excel file, use the PyVista library to create multiple planes representing different bedding layers, and generate random brightness colors for each depth. Visualize the data in the 3D model to obtain the shale multi-bedding reservoir model.

[0026] The PyVista library is a Python library for 3D plotting and processing mesh data, built on top of VTK (Visualization Toolkit). PyVista simplifies the complex interface of VTK, providing a more intuitive and easy-to-use API, and is widely used in the field of geological modeling.

[0027] Example 3 Unlike Example 2, in step 5, the depth of the line is calculated based on the pixel coordinates (y1), and the depth calculation formula is shown in Equation (9). Depth = Top depth + ((y1 / image height) * (bottom depth - top depth)) (9); During the depth data calculation process, a depth difference threshold is set. This threshold is used to determine whether the depth difference between two adjacent bedding planes is small. If the depth difference between two adjacent bedding planes is less than this threshold, they are considered to belong to the same bedding plane, and their average value is calculated as a single bedding plane.

[0028] Example 4 This invention provides an automated modeling method for multi-layered shale reservoirs based on bedding identification, significantly improving the accuracy and efficiency of shale reservoir modeling. By combining bedding identification technology, it accurately identifies shale bedding characteristics, providing more detailed and reliable basic data for reservoir modeling. Utilizing computer algorithms and automation technology, bedding identification is seamlessly integrated with the reservoir modeling process, reducing manual intervention and improving work efficiency. By achieving these objectives, this invention provides a novel modeling technique for the exploration and development of shale oil and gas resources, contributing to improved scientific rigor and reliability of resource evaluation.

[0029] Example 5 According to the method of the present invention, based on core samples from example shale reservoirs, multi-layered shale reservoir modeling combined with bedding identification is carried out, specifically as follows: (1) Core samples were obtained from the example shale reservoir, and the cores were scanned using a high-resolution imaging device to acquire their image data. The stratigraphic information of the core sample was also recorded: top depth 450.3 m, bottom depth 451.5 m, as shown below. Figure 2 As shown; (2) Denoising the core image: A 5x5 Gaussian kernel was used for convolution to denoise the core image. Adaptive Histogram Equalization (CLAHE) was used for contrast enhancement, converting the core image from the BGR color space to the LAB color space and separating the L, A, and B channels. Then, the L channel was enhanced for contrast, with a contrast limit parameter of 2.0 and a grid size of 8x8. Finally, the enhanced L channel was merged with the original A and B channels, and the image was converted back to the BGR color space to improve the contrast and clarity of the core image. Figure 3 As shown; (3) The Canny edge detection algorithm was used to identify the layered edges of the image after denoising and contrast enhancement. The low threshold of the Canny edge detection parameters was set to 50, and the high threshold was set to 150 to determine the edge pixels in the core image, such as... Figure 4 As shown; (4) Based on the edge detection results, the Hough line transform algorithm is used to detect straight lines in the core image. The Hough line transform distance resolution ρ is set to 1, the angular resolution θ to π / 180, the threshold to 20, the minimum line length to 100, and the maximum line gap to 20. Through the Hough line transform, the bedding structure lines in the core image can be identified, such as… Figure 5 As shown; (5) For each detected line, if its vertical height difference is less than 10 pixels (i.e., approximately horizontal), the line is drawn on the image, and its depth is calculated based on the line's y-coordinate. The depth calculation formula is: Depth = 450.3 + ((y1 / image height) * (451.5 - 450.3)), where 450.3 and 451.5 are predetermined ranges of layer depth values. Next, small differences in the depth data are detected and adjusted. If the difference between two adjacent depth values ​​is less than a threshold of 0.01m, their average is calculated as the adjusted depth. Finally, the adjusted depth values ​​are converted into a data frame and saved as an Excel file, such as... Figure 6 As shown; (6) A realistic shale multi-layered bedding reservoir model was built based on depth values ​​and visualized using PyVista. First, a depth array of bedding layers was defined, and a color with random brightness was generated for each depth to facilitate visualization in the 3D model. Next, multiple planes representing different bedding layers were created using the PyVista library and added to the drawing object for 3D display. Then, by setting the camera position and angle of the drawing object, different viewing angles were achieved, and the model intuitively displayed the structure of the shale multi-layered bedding reservoir, such as... Figure 7 As shown, this model can be used to carry out subsequent reservoir stimulation and analysis work, as well as to explore the influence of bedding on reservoir characteristics.

[0030] Example 6 This invention provides an automated modeling method for multi-layered shale reservoirs based on bedding identification, specifically as follows: Core sample acquisition and image acquisition: Multiple core samples were acquired in the target study area of ​​the shale reservoir; subsequently, the core samples were scanned using a high-resolution scanner to obtain high-quality core image data; the scanned images clearly show the texture and structural information inside the core. Image denoising: Gaussian denoising technology was applied to the acquired core images. A 5×5 Gaussian kernel matrix was used, with the central element having a weight of 0.3. The weights of other elements were determined according to the Gaussian distribution function and gradually decreased outwards. The image was then convolved using this Gaussian kernel matrix to calculate the new value for each pixel. This operation effectively reduced noise generated by scanning and other factors in the image while preserving the edge information of the core image relatively well.

[0031] Contrast Enhancement: The denoised image is enhanced using an Adaptive Histogram Equalization (CLAHE) algorithm. The image is divided into 4×4 local regions for processing. Within these local regions, previously blurred textures and color differences become clearer after processing, making details in the image easier to distinguish. This effectively avoids the loss of local details or over-enhancement that may occur with global histogram equalization, further improving image quality and facilitating accurate identification of subsequent layer edges.

[0032] Canny edge detection calculates image gradients by first converting the core image to grayscale using a grayscale conversion formula. Then, the Sobel operator is used to calculate the changes in the grayscale image in the horizontal (x) and vertical (y) directions to determine edge intensity and direction, thereby calculating the gradient magnitude and direction.

[0033] Non-maximum suppression: Based on the calculated gradient direction, the gradient magnitude of each pixel is compared with its neighboring pixels along the gradient direction. For example, within a certain region, for pixel P, if its gradient magnitude is greater than the pixel values ​​in its two neighboring directions along the gradient direction, then the value of point P is retained; otherwise, its value is set to 0. This step filters out potential edge pixels with large gradient magnitudes, reducing false edge detection.

[0034] Dual thresholding and hysteresis edge tracking: A high threshold of 80 and a low threshold of 30 are set. For a pixel in the image, if its gradient magnitude is greater than or equal to 80, it is considered a strong edge and marked as an edge; if it is less than or equal to 30, it is considered not an edge and marked as a non-edge; if it is between 30 and 80, it is considered a weak edge. Only when a weak edge is connected to a strong edge will it be retained as an edge; otherwise, it will be suppressed. For example, in a local region of the image, there are some weak edge pixels. By checking their neighboring pixels, if there are strong edge pixels in the neighborhood, these weak edge pixels are marked as edges; otherwise, they are removed. After hysteresis edge tracking processing, the true edges in the image are finally determined, and the edge pixels in the core image are accurately identified. These edge pixels provide important basis for the subsequent detection of bedding structure lines.

[0035] Hough Line Transform Line Detection: The edge detection results are used as input to generate a binary image, where edge pixels have a value of 1 and non-edge pixels have a value of 0. An angle range of [-30°, 30°] and a distance range of [0, 80°] are selected, and an accumulator is initialized to store the accumulated value for each parameter combination. For each edge point in the image, the Hough transform is used to convert it into a curve in the parameter space. A series of lines corresponding to ρ values ​​are calculated, and the parameter combinations of these lines are voted into the accumulator. By searching for local peaks in the accumulator, the lines in the image are identified. Transforming these lines back into the image coordinate system yields a layered structure line. In this way, multiple layered structure lines in the core image are successfully identified.

[0036] Straight line filtering: For the detected lines, calculate the change in the y-axis of adjacent pixels. Set a threshold of 2 pixels. For adjacent pixels A and B on a given line, if their y-axis coordinate difference is less than or equal to 2 pixels, then A and B are considered to belong to the same straight line; otherwise, they are not. After filtering, straight lines with small vertical height differences are obtained.

[0037] Depth Calculation: The depth of a line is calculated based on its pixel coordinates. Given a core sample with a top depth of 800 meters, a bottom depth of 900 meters, and an image height of 500 pixels, for a line with pixel coordinate y1 = 150 pixels, the depth calculation formula is: Depth = 800 + ((150 / 500) * (900 - 800)) = 830 meters. During depth calculation, if the depth values ​​of two adjacent bedding planes are small, for example, 830 meters and 832 meters respectively (the difference is less than a set threshold), their average value of 831 meters is calculated and considered as one bedding plane. All bedding plane depth values ​​are saved in an Excel file for subsequent data processing and model building.

[0038] 3D Model Visualization: Depth values ​​for each bedding layer were obtained from an Excel file, and multiple planes representing different bedding layers were created using the PyVista library. This method allows for the visualization of multi-bedding shale reservoirs within a 3D model, resulting in an intuitive model of multi-bedding shale reservoirs. This model clearly presents the multi-bedding structure of shale reservoirs and can be used for geological analysis, oil and gas resource assessment, and other research, contributing to a deeper understanding of the internal structure and characteristics of shale reservoirs.

[0039] This invention achieves accurate modeling and automated processing of actual shale reservoirs by organically combining bedding identification technology and reservoir modeling technology. The automated modeling method of this invention identifies shale bedding characteristics through a series of steps and constructs a multi-layered bedding reservoir model based on these characteristics. This method effectively improves the efficiency and accuracy of reservoir modeling while maintaining the high fidelity of actual core samples, reduces the uncertainty of human intervention, and can be widely applied to the exploration and development of shale oil and gas resources.

Claims

1. An automated modeling method for multi-layered shale reservoirs based on bedding identification, characterized in that, The specific steps are as follows: Step 1: Obtain core samples from shale reservoirs, record the stratigraphic information of the core samples, including the top and bottom depths, and obtain core images; Step 2: Denoise the core image and enhance its contrast. Step 3: Use the Canny edge detection algorithm to identify the bedding edges of the core image and identify the edge pixels in the core image; Step 4: Use the Hough line transform algorithm to detect straight lines and identify the layered structure lines in the image; Step 5: Determine whether the line is straight by calculating the changes of adjacent pixels on the y-axis, calculate the depth of the line based on the pixel coordinates, and save the layer depth value to an Excel file; Step 6: Obtain the depth value of each bedding layer from the Excel file, create planes representing different bedding layers, and generate colors of varying brightness to obtain a multi-bedding shale reservoir model.

2. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 1, characterized in that, Step 3 specifically involves: Step 3.1: Convert the core image to a grayscale image to eliminate color information; Step 3.2: Use the Sobel operator to calculate the grayscale image in the horizontal direction. and vertical direction The changes in are used to calculate the gradient of the image; Step 3.3: Apply the non-maximum suppression method to retain only local gradient maxima as potential edges; Step 3.4: Determine the final edge through dual thresholding and hysteresis edge tracking.

3. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 2, characterized in that, In step 3.1, the grayscale conversion formula is shown in equation (1): （1）； in, , , These represent the pixel values ​​of the red, green, and blue channels in the core image, respectively.

4. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 2, characterized in that, In step 3.2, the Sobel operator performs convolution operations on the image using the following convolution kernel; Horizontal Sobel core, S x As shown in equation (2): （2）； Vertical Sobel core, S y As shown in equation (3): （3）； For each pixel, the gradient in the horizontal and vertical directions is calculated using a convolution operation, as shown in equation (4): （4）； in, i, j For each pixel in a grayscale image; gradient magnitude The gradient direction ▽ represents the edge intensity of a pixel in the image, and the gradient direction ▽ represents the orientation of the edge. The calculation method is shown in Equation (5) and Equation (6): （5）； （6）。 5. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 4, characterized in that, In step 3.3, the specific steps are as follows: Based on the gradient direction ▽, the gradient magnitude of each pixel in the image is compared with the neighboring pixels along the gradient direction. If the gradient magnitude of the current pixel is greater than the pixel values ​​in both of its neighboring directions, the value of the pixel is retained; otherwise, its value is set to 0.

6. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 5, characterized in that, In step 3.4, the specific steps are as follows: Set two thresholds: high threshold and low threshold ; If the gradient magnitude of the pixel If the pixel is considered a strong edge, it is marked as an edge. If the gradient magnitude of the pixel If the pixel does not belong to the edge, it is considered not to be an edge and is marked as non-edge; If the gradient magnitude of the pixel If the pixel is weak, it will be considered a weak edge and will only be preserved as an edge if it is connected to a strong edge; otherwise, it will be suppressed. The process of hysteresis edge tracking is as follows: For each weak edge pixel, check its neighboring pixels. If there is a strong edge pixel in the neighborhood of the pixel, mark it as an edge; if there is no strong edge pixel in the neighborhood, remove it; for weak edge pixels marked as edges, continue to check its neighboring pixels; if there are more weak edge pixels in the neighborhood and they are connected to the marked strong edges, mark these weak edge pixels as edges as well; continue propagation until all weak edge pixels have been checked. Ultimately, all strong edge pixels are preserved; only weak edge pixels connected to strong edges are preserved, while other weak edges are suppressed; the final edge image contains all strong edges as well as the weak edges connected to them.

7. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 6, characterized in that, Step 4 specifically involves: Step 4.1: Take the edge detection result from step 3 as input to obtain a binary image, where the edge pixel value is 1 and the non-edge pixel value is 0. Step 4.2, Select the angle range and distance range and initialize an accumulator. To store the cumulative value of each parameter combination; the angle range is The distance range is ; Step 4.3: For each edge point in the image, the Hough transform is used to convert it into a curve in the parameter space. Finally, the peak value in the parameter space is found by accumulating these curves, thereby identifying the straight line. Step 4.4, will Combined votes are added to the accumulator; each element in the accumulator This indicates the number of votes for this parameter combination, i.e., how many edge points support this line; Step 4.5: Search for local peaks in the accumulator; the parameter combinations corresponding to the local peaks. The method involves finding the straight line in the image, identifying the maximum value in the accumulator, and then using the resulting polar coordinate parameters. Transform back to the image coordinate system to obtain the corresponding straight line in the image.

8. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 7, characterized in that, In step 5, the straightness of the line is determined by calculating the changes of adjacent pixels on the y-axis. Specifically: Suppose the current line contains multiple pixels. Then for two adjacent pixels and Calculate the difference between its y-axis coordinates, as shown in equation (8); （8）； Set a threshold This is used to determine if the height difference is too large; if If the two points are on the same horizontal line, then they are considered to belong to the same horizontal line; otherwise, they do not belong to the same horizontal line. In the entire detected straight line, if the height difference between adjacent points is less than the threshold, the straight line is considered to be straight; otherwise, the straight line is not straight. The depth is calculated based on the pixel coordinate y1 of the straight line and the layer depth value is saved as an Excel file.

9. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 8, characterized in that, The formula for depth calculation is shown in equation (9); Depth = Top depth + ((y1 / image height) * (bottom depth - top depth)) (9); During the depth data calculation process, a depth difference threshold is set. If the depth difference between two adjacent bedding planes is less than the threshold, they are considered to belong to the same bedding plane, and their average value is calculated as a single bedding plane.

10. The automated modeling method for multi-layered shale reservoirs based on bedding identification as described in claim 1, characterized in that, In step 6, the depth value of each bedding layer is obtained from the Excel file, multiple planes representing different bedding layers are created using the PyVista library, and a color with random brightness is generated for each depth. The results are then visualized in the 3D model to obtain a multi-bedding shale reservoir model.

Images