Method for reconstructing bedding structure rock finite element model based on rock core digital image

Abstract

The invention provides a bedding structure rock finite element model reconstruction method based on a rock core digital image, and belongs to the technical field of digital rock cores, and the method comprises the following steps: S1, obtaining a gray level image of a rock sample; s2, calculating a transverse mean value of the analysis area to obtain a longitudinal gray level distribution curve; s3, applying moving average filtering and smooth filtering to obtain a smooth gray curve; s4, first-order difference is carried out on the smooth gray level curve, the position with the absolute difference value larger than a threshold value is taken as an initial stripe boundary candidate, and starting and stopping pixels, the width and the average gray level of each stripe are determined; s5, mapping each clustering label interval into a physical interval based on FEM grid partition; and S6, constructing a sample finite element model, and generating a reconstructed finite element model with real partition information. According to the method, the macrostructure features of the rock can be efficiently recognized and extracted by directly utilizing the pictures and combining gray profile analysis, and efficient and low-cost reconstruction of the finite element model of the bedding structure rock is achieved.

Description

Technical Field

[0001] This invention belongs to the field of digital core technology, specifically relating to a method for reconstructing a finite element model of layered rock structure based on digital core images. Background Technology

[0002] For typical layered rocks such as siliceous banded dolomite, their mechanical properties are significantly affected by the heterogeneity of the internal structure and the distribution characteristics of the bands. Existing studies have mostly focused on homogeneous or micro-uniform models, while the overall and local mechanical response laws of layered rocks such as siliceous banded rocks are still unclear. This limits the inherent safety and efficient development of related underground engineering projects (such as oil and gas drilling, reservoir fracturing, and underground space utilization). Therefore, it is urgent to develop mechanical modeling and numerical analysis methods that can reflect the actual layered structure characteristics.

[0003] Traditional 3D digital core reconstruction techniques are mostly based on high-resolution X-ray CT or FIB-SEM to acquire 2D or 3D grayscale images of rock samples, and then use specific algorithms to achieve 3D reconstruction. While these methods can restore pores, fractures, and mineral distribution at the microscale, the details they reconstruct are mainly concentrated at the micrometer and below scale. They have limited ability to restore complex layering, banding, and transition zones of layered rocks at the macroscale (such as millimeter-level samples). At the same time, the experimental equipment is expensive, the data processing volume is large, and the model reconstruction cycle is long, resulting in high overall application costs and making it difficult to meet the needs of batch modeling and efficient simulation in engineering. Summary of the Invention

[0004] To address the problems existing in the prior art, a method for reconstructing a finite element model of layered rock based on core digital images is provided.

[0005] The technical solution adopted by this invention to solve its technical problem is: This technical solution proposes a method for reconstructing a finite element model of layered rock structure based on digital core images, including the following steps: S1: Acquire digital images of the rock sample and convert the digital images into grayscale images; S2: Using the center line of the grayscale image as the axis of symmetry, extract a vertical rectangular strip as the analysis area, calculate the horizontal mean of the analysis area, and obtain the vertical grayscale distribution curve; S3: Apply moving average filtering to eliminate high-frequency jitter in the vertical grayscale distribution curve, and then apply smoothing filtering to further denoise, resulting in a smooth grayscale curve. S4: Perform first-order difference on the smooth grayscale curve, take the positions where the absolute difference is greater than the threshold as the initial strip boundary candidates, apply the minimum interval and minimum width filtering rules to filter, and determine the start and end pixels, width and average grayscale of each strip. S5: The smooth gray-scale curve is partitioned based on the K-means clustering algorithm to obtain several clustering label intervals. Each clustering label interval corresponds to a physical region type. Based on the rock sample height and unit size, each clustering label interval is mapped to a physical region based on FEM grid partitioning. S6: Based on the rock sample size, construct a finite element model of the sample. Based on the mesh mapping algorithm, renumber each element in the finite element model of the sample to finally generate a reconstructed finite element model with real partition information.

[0006] Preferably, in step S2, the analysis area includes lithological segments that represent the typical bedding orientation and vertical variation of cross-banding. The average horizontal pixel value is calculated along each row of the extracted analysis area to obtain a vertical grayscale profile, including: Let the grayscale matrix of the original image be... I ( i , j After selecting the central vertical band range, the average grayscale value of each row is expressed by the formula: (1); In the formula, This represents the set of horizontal pixels of the selected vertical band. I Represents the grayscale matrix of the original image. i =1,…, H For vertical pixel rows, j =1,…, W Horizontal pixel column, y i Indicates the first i The average gray value of the row. N j This indicates the number of pixels horizontally.

[0007] Preferably, in step S3, a first-level moving average is performed on the longitudinal grayscale distribution curve, and the moving average uses a symmetrical window to obtain a first-level smooth curve. The formula is expressed as: (2); In the formula, x i Indicates the first i The original grayscale values ​​of each pixel row. x i+j Indicates the current line i Centered on, offset to the left and right k The grayscale value corresponding to the row m This represents the half-length parameter of the sliding window. m The value range is from 2 to 5, 2 m +1 indicates the length of the symmetrical window used in the moving average. Indicates the firsti The output value is the row moving average of the pixels.

[0008] Preferably, the first-level smooth curve The input is smoothed using a Savitzky–Golay filter with second-order smoothing, employing a local window of odd length and a polynomial order. p By performing local least-squares polynomial fitting at each center location, a smooth grayscale curve is obtained, expressed by the formula: (3); In the formula, c k The selected window length is 2 M +1 and polynomial order p The weighted combination weights of the center points are pre-calculated using least squares fitting, representing the weighted combination weights of the center points under this local fitting. y i SG This represents the grayscale value after two levels of smoothing. M Indicates half the window width. k Indicates the relative offset within the window. Indicates the first i + k The grayscale values ​​of each pixel row after first-level smoothing.

[0009] Preferably, in step S4, the smooth grayscale curve undergoes first-order difference processing, and the difference formula is expressed as: d i = y i - y i-1 , i =2,……, N (4); In the formula, y i Indicates the first i Smooth grayscale values ​​of each pixel row, N represents Section length, d i Indicates the first i The difference result of each pixel row y i-1 Indicates the first i- Smoothed grayscale values ​​for a single pixel row; The adaptive threshold is calculated by taking a proportion of the maximum absolute value of the difference, as expressed by the formula: (5); In the formula, T Indicates an adaptive threshold. α This represents the adjustment coefficient. αTypical values ​​range from 0.05 to 0.2, used to distinguish significant mutations from background fluctuations.

[0010] Preferred, all satisfying │ d i The position of │>T is used as an initial candidate point for the strip boundary. Minimum interval and minimum width filtering rules are applied to eliminate erroneous detections where the distance is less than the set minimum separation pixel number or the corresponding strip width is shorter than a set value. Based on the final boundary set, the start and end pixels of each strip are constructed. s k , e k ] Calculate its width w k = e k - s k +1, calculate the average gray level of this interval, expressed by the formula: (6); In the formula, s k Indicates the first k The starting pixel row number of each strip. e k Indicates the first k The row number of the terminating pixel of each stripe. w k Indicates the first k The width of each strip Indicates the first k The average gray value of each strip.

[0011] Preferably, in step S5, the smooth grayscale curve after smoothing and noise reduction is partitioned using the K-means clustering algorithm to obtain several clustering label intervals. Each clustering label interval corresponds to a physical region type. Based on the clustering results, the position of clustering label change is detected in sequence to determine the upper and lower boundaries of pixels in each region. Let the total pixel height of the rock sample in the image be... H p The actual physical height is H The conversion factor from pixel to physical height is λ = H / H p The starting pixel of each interval s k With terminating pixel e k Mapped to actual height range λ [( s k -1) λ ,(e k -1)], combined with the preset element height in the finite element model h elem Alignment is performed so that each image partition matches one or more cell layers of the FEM mesh, thereby achieving automatic mapping of grayscale image partitions to the finite element mesh and assignment of physical region attributes.

[0012] Preferably, in step S6, the node coordinates and element connection relationships of the finite element model of the sample are read, and the first... i The coordinates of the eight nodes of each unit in the longitudinal direction are: z i,1 , z i,2 , ..., z i,8 The formula for calculating the geometric center height of this unit is expressed as: (7); In the formula, Z i,j Indicates the first i Unit 1 j The coordinates of each node in the vertical direction Z i (c) Indicates the first i The geometric center height of each finite element. j Indicates the index of the node within the cell. i Indicates the index of the finite element.

[0013] Preferably, the starting height of each partition is loaded from an external partition information file. a p Termination height b p and the corresponding partition type number T p To form an ordered triple ( a p , b p , T p ),in, p =1, ..., P This represents all partitions, for each unit. i Iterate through all partitions in sequence, for the first partition... p Each partition is checked to determine if the center height of the unit meets the following requirements: (8); If satisfied, then number the components of that unit. part i The value is assigned to the type number of the partition, which is represented as part i = T p And terminate the partition traversal of that unit; If satisfied Z i (c) = T p Then let part i = T p After all finite element elements have been processed, a reconstructed finite element model with real partition information is generated.

[0014] Compared with the prior art, the present invention has the following advantages: 1. This application uses macroscopic rock surface images captured by a digital camera. Through image grayscale feature analysis and partitioning algorithms, it can accurately identify various stripes and matrix regions. At the same time, combined with automated data processing and physical parameter mapping, it can efficiently deploy FEM grid cells. This can significantly reduce reconstruction costs while fully preserving the macroscopic partitioning characteristics and engineering-related structural details of layered rocks. Thus, it provides a low-cost, high-efficiency, and structurally realistic model foundation for numerical research on the mechanical behavior of layered rocks.

[0015] 2. This invention effectively overcomes the dependence of existing technologies on expensive CT / SEM equipment and microscale reconstruction. It can directly utilize ordinary photographs combined with grayscale profile analysis and automatic partitioning algorithms to efficiently identify and extract macroscopic structural features of layered rocks, such as siliceous banded dolomite. Through noise reduction processing and clustering partitioning, it significantly improves the accuracy of macroscopic structural partitioning such as bands and matrix, and can directly output physical regions and types suitable for FEM modeling, realizing automatic mapping from pixels to finite element models. This greatly reduces modeling costs and technical barriers, filling the gaps in existing methods in the macroscopic structure FEM reconstruction stage. Attached Figure Description

[0016] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which: Figure 1 This is the overall flowchart of the present invention; Figure 2 This is a schematic diagram of the central vertical analysis area of ​​the siliceous banded dolomite sample in Example 2; Figure 3 This is the original grayscale profile partitioning diagram of the siliceous banded dolomite sample in Example 2; Figure 4 This is a grayscale profile partitioning diagram of the silica-bearing dolomite sample after noise reduction in Example 2; Figure 5 This is a schematic diagram of the finite element model of the siliceous banded dolomite sample in Example 2; Figure 6 This is a schematic diagram of the central vertical analysis area of ​​the siltstone-mudstone interbedded sample in Example 3; Figure 7 This is the original grayscale profile zoning diagram of the siltstone-mudstone interbedded sample in Example 3; Figure 8 This is a noise-reduced grayscale profile partitioning diagram of the siltstone-mudstone interbedded sample in Example 3; Figure 9 This is a schematic diagram of the finite element model of the siltstone-mudstone interbedded sample in Example 3. Detailed Implementation

[0017] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0018] Example 1 like Figure 1 As shown, this embodiment proposes a method for reconstructing a finite element model of layered rock structure based on digital core images, including the following steps: S1: Take pictures of rock samples with a digital camera to obtain digital images of the rock samples, and convert the digital images into grayscale images for subsequent processing; S2: Using the center line of the grayscale image as the axis of symmetry, a vertical rectangular strip is extracted as the analysis area. The horizontal mean of the analysis area is calculated to obtain the vertical grayscale distribution curve, which reflects the grayscale characteristics of different lithological intervals. S3: Apply moving average filtering to eliminate high-frequency jitter in the vertical grayscale distribution curve, and then apply smoothing filtering to further denoise, resulting in a smooth grayscale curve that retains the main trend while suppressing mid-to-high frequency interference. S4: Perform first-order difference on the smooth grayscale curve, take the positions where the absolute difference is greater than the threshold as the initial strip boundary candidates, apply the minimum interval and minimum width filtering rules to filter out false detections that are too close or too narrow, and determine the start and end pixels, width and average grayscale of each strip. S5: The smooth gray-scale curve is partitioned based on the K-means clustering algorithm (the number of clusters can be 2 to 8, K-means++ is used for initialization, and multiple restarts are used to avoid local optima), resulting in several cluster label intervals. Each cluster label interval corresponds to a physical region type. Based on the rock sample height and unit size, each cluster label interval is mapped to a physical interval based on FEM grid partitioning. S6: Based on the rock sample size, construct a finite element model of the sample, save the information of each element and node of the sample, and renumber each element in the sample finite element model based on the mesh mapping algorithm. First, calculate the center height of each element, and then assign the region type number corresponding to the height interval into which it falls to the element as the component number. The boundary height is processed according to the next interval assignment rule to avoid ambiguity. Finally, a reconstructed finite element model with real partition information is generated.

[0019] In S2, the analysis area includes lithological segments that represent the typical bedding orientation and vertical variation of cross-bands. To further ensure the representativeness of this central band, an adaptive check is introduced in the selection and verification process of candidate bands. The average horizontal pixel value is calculated along each row of the extracted analysis area to obtain a vertical grayscale profile, including: Let the grayscale matrix of the original image be... I ( i , j After selecting the central vertical band range, the average grayscale value of each row is expressed by the formula: (1); In the formula, This represents the set of horizontal pixels of the selected vertical band. I Represents the grayscale matrix of the original image. i =1,…, H For vertical pixel rows, j =1,…, W Horizontal pixel column, y i Indicates the first i The average gray value of the row. N j This indicates the number of pixels horizontally.

[0020] In S3, a first-level moving average is performed on the longitudinal grayscale distribution curve to suppress high-frequency fluctuations. The moving average uses a symmetrical window to obtain a first-level smooth curve. The formula is expressed as: (2); In the formula, x i Indicates the first i The original grayscale values ​​of each pixel row. x i+j Indicates the current line i Centered on, offset to the left and right k The grayscale value corresponding to the row m This represents the half-length parameter of the sliding window. m The value range is from 2 to 5, 2 m+1 indicates the length of the symmetrical window used in the moving average, with the corresponding window length 2m+1 ranging from 5 to 11, in order to strike a balance between noise reduction and trend preservation. Indicates the first i The output value is the row moving average of the pixels.

[0021] The first-level smooth curve The input is smoothed using a Savitzky-Golay filter in two stages. This two-stage noise reduction process gradually weakens high-frequency and mid-frequency disturbances while preserving the main grayscale trend to the maximum extent, providing a stable and high signal-to-noise ratio input signal for subsequent strip boundary localization and lithological differentiation. An odd-length local window is used, with a polynomial order. p By performing local least-squares polynomial fitting at each center location, a smooth grayscale curve is obtained, expressed by the formula: (3); In the formula, c k The selected window length is 2 M +1 and polynomial order p (Selecting 2 to 4) Pre-calculated using least squares fitting, representing the weighted combination weight of the center point under this local fitting. y i SG This represents the grayscale value after two levels of smoothing. M Indicates half the window width. M The value ranges from 7 to 25, corresponding to a window size between 15 and 51. k Indicates the relative offset within the window. Indicates the first i + k The grayscale values ​​of each pixel row after first-level smoothing.

[0022] In S4, the smooth grayscale curve is subjected to first-order difference processing, and the difference formula is expressed as: d i = y i - y i-1 , i =2,……, N (4); In the formula, y i Indicates the first i Smooth grayscale values ​​of each pixel row, N represents Section length, d i Indicates the first i The difference result of each pixel row y i-1 Indicates the firsti- Smoothed grayscale values ​​for a single pixel row; The adaptive threshold is calculated by taking a proportion of the maximum absolute value of the difference, as expressed by the formula: (5); In the formula, T Indicates an adaptive threshold. α This represents the adjustment coefficient. α Typical values ​​range from 0.05 to 0.2, used to distinguish significant mutations from background fluctuations.

[0023] All conditions are met | d i The position of │>T serves as an initial candidate point for the strip boundary. Minimum interval and minimum width filtering rules are applied to eliminate errors where the distance is less than the set minimum separation pixel count, or the corresponding strip width is shorter than a set value. This ensures that the retained boundaries are physically coherent and representative. Based on the final boundary set, the starting and ending pixels of each strip are constructed. s k , e k ] Calculate its width w k = e k - s k +1, calculate the average gray level of this interval, expressed by the formula: (6); In the formula, s k Indicates the first k The starting pixel row number of each strip. e k Indicates the first k The row number of the terminating pixel of each stripe. w k Indicates the first k The width of each strip Indicates the first k The average gray value of each strip.

[0024] In S5, when the sample has a complex strip structure or transition region, the K-means clustering algorithm is used to partition the smooth gray-scale curve after smoothing and noise reduction to obtain several clustering label intervals. Each clustering label interval corresponds to a physical region type. Based on the clustering results, the position of clustering label change is detected in sequence to determine the upper and lower boundaries of pixels in each region. Let the total pixel height of the rock sample in the image be... H p The actual physical height is HThe conversion factor from pixel to physical height is λ = H / H p The starting pixel of each interval s k With terminating pixel e k Mapped to actual height range λ [( s k -1) λ ,( e k -1)], combined with the preset element height in the finite element model h elem Alignment is performed so that each image partition matches one or more cell layers of the FEM mesh, thereby achieving automatic mapping of grayscale image partitions to the finite element mesh and assignment of physical region attributes.

[0025] In S6, the nodal coordinates and element connection relationships of the finite element model of the specimen are read, and the first... i The coordinates of the eight nodes of each unit in the longitudinal direction are: z i,1 , z i,2 , ..., z i,8 The formula for calculating the geometric center height of this unit is expressed as: (7); In the formula, Z i,j Indicates the first i Unit 1 j The coordinates of each node in the vertical direction Z i (c) Indicates the first i The geometric center height of each finite element. j Indicates the index of the node within the cell. i Indicates the index of the finite element.

[0026] Load the starting height of each partition from the external partition information file. a p Termination height b p and the corresponding partition type number T p To form an ordered triple ( a p , b p , T p ),in, p=1, ..., P This represents all partitions, for each unit. i Iterate through all partitions in sequence, for the first partition... p Each partition is checked to determine if the center height of the unit meets the following requirements: (8); If satisfied, then number the components of that unit. part i The value is assigned to the type number of the partition, which is represented as part i = T p And terminate the partition traversal of that unit; To avoid missing cases where the height of the topmost point equals the upper bound of the last interval, if the following conditions are met... Z i (c) = T p Then explicitly let part i = T p After all finite element elements have been processed, a reconstructed finite element model with real partition information is generated.

[0027] Example 2 like Figures 2-5 As shown, this embodiment proposes a method for reconstructing a finite element model of bedding structure rocks based on digital core images. Taking a siliceous banded dolomite sample from a certain stratum as an example, the method includes the following steps: S21: Obtain digital images of the rock sample surface by taking pictures with a digital camera, and after importing, convert the images to grayscale mode for easy subsequent processing.

[0028] S22: Automatically select the vertical stripe in the center of the image as the analysis area, such as... Figure 2 As shown, Figure 2 This is a schematic diagram of the central vertical analysis area of ​​the siliceous banded dolomite sample in Example 2. The longitudinal grayscale distribution curve is obtained by calculating its lateral mean, reflecting the grayscale characteristics of different lithological intervals. First, in a laboratory environment, a high-resolution digital camera was used to vertically photograph the side surface of the siliceous banded dolomite sample to ensure uniform light source and no shadows. The relative position of the camera and the sample was fixed to obtain repeatable images. After the images were taken, the original color images were imported into MATLAB and converted into single-channel grayscale images to remove interference from color information and retain only the grayscale distribution.

[0029] S23: Apply moving average filtering to eliminate high-frequency jitter in the longitudinal grayscale distribution curve, and then apply smoothing filtering to further denoise, resulting in a smooth grayscale curve that retains the main trend while suppressing mid-to-high frequency interference, providing a high-quality signal for partition boundary identification.

[0030] like Figure 3 As shown, Figure 3 This is the original grayscale profile partitioning map of the siliceous banded dolomite sample in Example 2. First, two smoothing methods were applied to the original grayscale profile (a sequence of horizontal average values ​​in pixel rows): one is short window moving average, and the other is Savitzky-Golay filtering with a larger window.

[0031] The moving average method can remove high-frequency "spurt" noise, while the Savitzky-Golay filter preserves a more complete overall trend. Combining the two methods can suppress local random fluctuations to the greatest extent possible while ensuring the main trend of the curve is not distorted, providing a high-quality signal for subsequent boundary detection and partitioning. The processing result is as follows: Figure 4 As shown, Figure 4 This is a grayscale profile partitioning diagram of the silica-containing banded dolomite sample in Example 2 after noise reduction.

[0032] S24: Based on the average gray level of each band, an automatic threshold is set to initially classify all bands into siliceous bands (low gray level), siliceous dolomite (medium gray level), and dolomite matrix (high gray level), thus achieving initial classification of macroscopic lithological types; First, the average gray level of each curve segment is calculated according to the upper and lower boundaries of each band interval. The overall average gray level of all bands is used as the dividing scale, resulting in two automatic threshold lines corresponding to two gray level levels of 105.7 and 162.5, respectively. Then, the area of ​​the curve below 105.7 is colored blue, representing siliceous bands; the area with gray level between 105.7 and 162.5 is colored green, representing siliceous dolomite; and the area with gray level above 162.5 is colored red, representing dolomite matrix. This completes the preliminary differentiation of macroscopic lithological types.

[0033] S25: The multi-interval automatic segmentation step uses K-means clustering to divide the grayscale curve into three partitions, and combines physical height and unit size parameters to map the image pixel segments into the physical intervals required for finite element modeling, thereby realizing automatic partitioning from image to model. By combining the total height of the sample (5cm) with the pixel-to-centimeter ratio of the number of rows in the image, each pixel row is mapped to a physical height, and the cluster partition number is assigned to the corresponding element in the finite element model, thus automatically completing the "part number" reconstruction from the two-dimensional image to the three-dimensional FEM element.

[0034] S26: Construct a finite element model of the specimen that conforms to the actual specimen size, save the information of each element and node of the specimen, and use LS-PrePost to create a new cylindrical geometry with a diameter of 50 mm and a height of 50 mm in its "Geometry" module. This geometry is completely consistent with the actual core specimen size.

[0035] Subsequently, in the "Mesh" module, a three-dimensional solid eight-node element mesh was generated according to the set element height of 0.04 cm and circumferential mesh density to ensure that the width of each element is close to the engineering requirements. After the mesh was generated, all node numbers and their three-dimensional coordinates, as well as the correspondence between each element and its eight nodes, were exported and saved to a text file, thus constructing complete finite element model node information and element topology information. Based on the principle of the mesh mapping algorithm and combined with the regional partitioning information obtained in S25, each element in the finite element model of the sample is renumbered to realize the automatic allocation of the element component number. First, each physical height interval obtained in step S25 is matched with the height mapping relationship corresponding to each row of elements when meshing. Using the exported node and element information, the center height of each element is compared with the pre-divided regional interval to determine its partition type.

[0036] Specifically, the program iterates through all units in the model sequentially, and for the first... i Each element reads its center height value and searches the interval array for the partition index that satisfies "interval start height ≤ element center height < interval end height". Then, the component number corresponding to that element in LS-PrePost is set to the lithology number corresponding to that partition. The result is as follows: Figure 5 As shown, Figure 5 This is a schematic diagram of the finite element model of the siliceous banded dolomite sample in Example 2.

[0037] Example 3 like Figures 6-9 As shown, this embodiment proposes a method for reconstructing a finite element model of layered rock structure based on core digital images. Taking an interbedded siltstone and mudstone sample from a certain stratum as an example, the method includes the following steps: S31: Digital images of the surface of siltstone-mudstone interbedded samples are obtained by taking pictures with a digital camera. After importing, the images are uniformly converted to grayscale mode for easy subsequent processing.

[0038] S32: Automatically select the vertical stripe slightly to the left of the center of the image as the analysis region, such as... Figure 6 As shown, Figure 6 This is a schematic diagram of the central vertical analysis area of ​​the siltstone-mudstone interbedded sample in Example 3. The longitudinal grayscale distribution curve is obtained by calculating its lateral mean, reflecting the grayscale characteristics of different lithological intervals. First, in a laboratory environment, a high-resolution digital camera was used to vertically photograph the side surface of the siltstone-mudstone interbedded sample to ensure uniform light source and no shadows. The relative position of the camera and the sample was fixed to obtain repeatable images. After the images were taken, the original color images were imported into MATLAB and converted into single-channel grayscale images to remove interference from color information and retain only the grayscale distribution.

[0039] S33: Apply moving average filtering to eliminate high-frequency jitter in the longitudinal grayscale distribution curve, and then apply smoothing filtering to further denoise, resulting in a smooth grayscale curve that retains the main trend while suppressing mid-to-high frequency interference, providing a high-quality signal for partition boundary identification.

[0040] like Figure 7 As shown, Figure 7 This is the original grayscale profile partition map of the siltstone-mudstone interbedded sample in Example 3. First, two smoothing methods were applied to the original grayscale profile (a sequence of horizontal average values ​​in pixel rows): one is short window moving average, and the other is Savitzky-Golay filtering with a larger window.

[0041] The moving average method can remove high-frequency "spurt" noise, while the Savitzky-Golay filter preserves a more complete overall trend. Combining the two methods can suppress local random fluctuations to the greatest extent possible while ensuring the main trend of the curve is not distorted, providing a high-quality signal for subsequent boundary detection and partitioning. The processing result is as follows: Figure 8 As shown, Figure 8 This is a grayscale profile partitioning diagram of the siltstone-mudstone interbedded sample in Example 3 after noise reduction.

[0042] S34: Based on the average gray level of each strip, an automatic threshold is set to initially classify all strips into siltstone (low gray level) and mudstone (high gray level), thus achieving initial classification of macroscopic lithological types; First, the average gray level of each curve segment is calculated according to the upper and lower boundaries of each band interval. The average gray level of all bands is used as the dividing scale to obtain an automatic threshold line, which corresponds to a gray level of 71.3. Then, the area of ​​the curve below 71.3 is colored blue to represent siltstone, and the area of ​​gray level above 71.3 is colored green to represent mudstone. This completes the preliminary differentiation of macroscopic lithological types.

[0043] S35: The multi-interval automatic segmentation process uses K-means clustering to divide the grayscale curve into two partitions, and combines physical height and unit size parameters to map image pixel segments into physical intervals required for finite element modeling, thereby realizing automatic partitioning from image to model. By combining the total height of the sample (4cm) with the pixel-to-centimeter ratio of the number of rows in the image, each pixel row is mapped to a physical height, and the cluster partition number is assigned to the corresponding element in the finite element model, thus automatically completing the "part number" reconstruction from the two-dimensional image to the three-dimensional FEM element.

[0044] S36: Construct a finite element model of the specimen that matches the actual specimen size, save the information of each element and node of the specimen, and use LS-PrePost to create a new cylindrical geometry with a diameter of 25 mm and a height of 40 mm in its "Geometry" module. This geometry is completely consistent with the actual core specimen size.

[0045] Subsequently, in the "Mesh" module, a three-dimensional solid eight-node element mesh is generated according to the set element height of 0.2 mm and circumferential mesh density to ensure that the width of each element is close to the engineering requirements. After the mesh is generated, all node numbers and their three-dimensional coordinates, as well as the correspondence between each element and its eight nodes, are exported and saved to a text file, thereby constructing complete finite element model node information and element topology information.

[0046] Based on the principle of the mesh mapping algorithm and combined with the regional partitioning information obtained in S35, each element in the finite element model of the sample is renumbered to achieve automatic allocation of the element component number. First, each physical height interval obtained in S35 is matched with the height mapping relationship corresponding to each row of elements when meshing. Using the node and element information exported in S206, the center height of each element is compared with the pre-divided regional interval to determine its partition type.

[0047] Specifically, the program iterates through all units in the model sequentially, and for the first... i Each element reads its center height value and searches the interval array for the partition index that satisfies "interval start height ≤ element center height < interval end height". Then, the component number corresponding to that element in LS-PrePost is set to the lithology number corresponding to that partition. The result is as follows: Figure 9 As shown, Figure 9 This is a schematic diagram of the finite element model of the siltstone-mudstone interbedded sample in Example 3.

[0048] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to the embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the claims and their equivalents.

Claims

1. A method for reconstructing a finite element model of layered rock structure based on digital core images, characterized in that, The following steps are involved: S1: Acquire digital images of the rock sample and convert the digital images into grayscale images; S2: Using the center line of the grayscale image as the axis of symmetry, extract a vertical rectangular strip as the analysis area, calculate the horizontal mean of the analysis area, and obtain the vertical grayscale distribution curve; S3: Apply moving average filtering to eliminate high-frequency jitter in the vertical grayscale distribution curve, and then apply smoothing filtering to further denoise, resulting in a smooth grayscale curve. S4: Perform first-order difference on the smooth grayscale curve, take the positions where the absolute difference is greater than the threshold as the initial strip boundary candidates, apply the minimum interval and minimum width filtering rules to filter, and determine the start and end pixels, width and average grayscale of each strip. S5: The smooth gray-scale curve is partitioned based on the K-means clustering algorithm to obtain several clustering label intervals. Each clustering label interval corresponds to a physical region type. Based on the rock sample height and unit size, each clustering label interval is mapped to a physical region based on FEM grid partitioning. S6: Based on the rock sample size, construct a finite element model of the sample. Based on the mesh mapping algorithm, renumber each element in the finite element model of the sample to finally generate a reconstructed finite element model with real partition information.

2. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 1, characterized in that, In S2, the analysis area includes lithological segments that represent the typical bedding orientation and vertical variation of cross-banding. The average horizontal pixel value is calculated along each row of the extracted analysis area to obtain a vertical grayscale profile, including: Let the grayscale matrix of the original image be... I ( i , j After selecting the central vertical band range, the average grayscale value of each row is expressed by the formula: （1）； In the formula, This represents the set of horizontal pixels of the selected vertical band. I Represents the grayscale matrix of the original image. i =1,…, H For vertical pixel rows, j =1,…, W Horizontal pixel column, y i Indicates the first i The average gray value of the row. N j This indicates the number of pixels horizontally.

3. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 1, characterized in that, In step S3, a first-level moving average is performed on the longitudinal grayscale distribution curve. The moving average uses a symmetrical window to obtain a first-level smooth curve. The formula is expressed as: （2）； In the formula, x i Indicates the first i The original grayscale values ​​of each pixel row. x i+j Indicates the current line i Centered on, offset to the left and right k The grayscale value corresponding to the row m This represents the half-length parameter of the sliding window. m The value ranges from 2 to 5, 2 m +1 indicates the length of the symmetrical window used in the moving average. Indicates the first i The output value is the sliding average of the pixel rows.

4. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 3, characterized in that, The first-level smooth curve The input is smoothed using a Savitzky–Golay filter with second-order smoothing, employing a local window of odd length and a polynomial order. p By performing local least-squares polynomial fitting at each center location, a smooth grayscale curve is obtained, expressed by the formula: （3）； In the formula, c k The selected window length is 2 M +1 and polynomial order p The weighted combination weights of the center points are pre-calculated using least squares fitting, representing the weighted combination weights of the center points under this local fitting. y i SG This represents the grayscale value after two levels of smoothing. M Indicates half the window width. k Indicates the relative offset within the window. Indicates the first i + k The grayscale values ​​of each pixel row after first-level smoothing.

5. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 1, characterized in that, In step S4, the smooth grayscale curve is subjected to first-order difference processing, and the difference formula is expressed as: d i = y i － y i-1 , i =2,……, N （4）； In the formula, y i Indicates the first i Smooth grayscale values ​​of each pixel row, N represents Section length, d i Indicates the first i The difference result of each pixel row y i-1 Indicates the first i- Smoothed grayscale values ​​for a single pixel row; The adaptive threshold is calculated by taking a proportion of the maximum absolute value of the difference, as expressed by the formula: （5）； In the formula, T Indicates an adaptive threshold. α This represents the adjustment coefficient. α Typical values ​​range from 0.05 to 0.2, used to distinguish significant mutations from background fluctuations.

6. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 5, characterized in that, All conditions are met | d i The position of │>T is used as an initial candidate point for the strip boundary. Minimum interval and minimum width filtering rules are applied to eliminate erroneous detections where the distance is less than the set minimum separation pixel number or the corresponding strip width is shorter than a set value. Based on the final boundary set, the start and end pixels of each strip are constructed. s k , e k ] Calculate its width w k = e k - s k +1, calculate the average gray level of this interval, expressed by the formula: （6）； In the formula, s k Indicates the first k The starting pixel row number of each strip. e k Indicates the first k The row number of the terminating pixel of each stripe. w k Indicates the first k The width of each strip Indicates the first k The average gray value of each strip.

7. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 6, characterized in that, In S5, the K-means clustering algorithm is used to partition the smooth grayscale curve after smoothing and noise reduction to obtain several clustering label intervals. Each clustering label interval corresponds to a physical region type. Based on the clustering results, the position of clustering label change is detected in sequence to determine the upper and lower boundaries of pixels in each region. Let the total pixel height of the rock sample in the image be... H p The actual physical height is H The conversion factor from pixel to physical height is λ = H / H p The starting pixel of each interval s k With terminating pixel e k Mapped to actual height range λ [( s k -1) λ ,( e k -1)], combined with the preset element height in the finite element model h elem Alignment is performed so that each image partition matches one or more cell layers of the FEM mesh, thereby achieving automatic mapping of grayscale image partitions to the finite element mesh and assignment of physical region attributes.

8. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 1, characterized in that, In step S6, the node coordinates and element connection relationships of the finite element model of the sample are read, and the first... i The coordinates of the eight nodes of each unit in the longitudinal direction are: z i,1 , z i,2 , ..., z i,8 The formula for calculating the geometric center height of this unit is expressed as: （7）； In the formula, Z i,j Indicates the first i Unit 1 j The coordinates of each node in the vertical direction Z i (c) Indicates the first i The geometric center height of each finite element. j Indicates the index of the node within the cell. i Indicates the index of the finite element.

9. The method for reconstructing a finite element model of layered rock structure based on core digital images according to claim 8, characterized in that, Load the starting height of each partition from the external partition information file. a p Termination height b p and the corresponding partition type number T p To form an ordered triple ( a p , b p , T p ),in, p =1, ..., P This represents all partitions, for each unit. i Iterate through all partitions in sequence, for the first partition... p Each partition is checked to determine if the center height of the unit meets the following requirements: （8）； If satisfied, then number the components of that unit. part i The value is assigned to the type number of the partition, which is represented as part i = T p And terminate the partition traversal of that unit; If satisfied Z i (c) = T p Then let part i = T p After all finite element elements have been processed, a reconstructed finite element model with real partition information is generated.

Images