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

Abstract

The application provides a bedding structure rock finite element model reconstruction method based on a core digital image, and belongs to the technical field of digital cores and comprises the following steps: S1: acquiring a gray image of a rock sample; S2: calculating and analyzing the horizontal mean value of a region to obtain a longitudinal gray distribution curve; S3: applying a sliding average filter and a smoothing filter to obtain a smoothed gray curve; S4: performing first-order differentiation on the smoothed gray curve, taking positions with an absolute difference value greater than a threshold value as initial strip boundary candidates, and determining the start and end pixels, the width and the average gray of each strip; S5: mapping each clustering label interval to a physical interval based on FEM grid partition; and S6: constructing a sample finite element model to generate a reconstructed finite element model with real partition information. The application can directly utilize a photo in combination with gray profile analysis, efficiently identify and extract macroscopic structure characteristics of rock, and realize efficient and low-cost reconstruction of a bedding structure rock finite element model.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of digital core and particularly relates to a bedding structure rock finite element model reconstruction method based on a digital core image. BACKGROUND

[0002] For typical bedding structure rocks such as dolomite containing siliceous bands, the mechanical properties thereof are significantly affected by the internal structure heterogeneity and band distribution characteristics of the rocks, and existing researches are mostly focused on homogeneous or mesoscopic uniform models, and the overall and local mechanical response laws of the layered structure rocks such as siliceous bands are still unclear, which further limits the essential safety and efficient development of related underground engineering (such as oil and gas drilling, reservoir fracturing and underground space utilization, etc.), and therefore, it is urgent to develop a mechanical modeling and numerical analysis method capable of reflecting the actual layered structure characteristics.

[0003] Traditional three-dimensional digital core reconstruction technology is mostly based on high-resolution X-ray CT or FIB-SEM experimental methods to obtain two-dimensional or three-dimensional gray-scale images of rock samples, and then a specific algorithm is used to realize three-dimensional reconstruction. Although this method can restore the distribution of pores, fractures and minerals at the microscale, the details of the reconstruction are mainly concentrated on the micron and below scale, and the restoration capability of the layered structure rocks at the macro scale (such as millimeter level samples) is limited. At the same time, the experimental equipment is expensive, the data processing amount is large, and the model reconstruction cycle is long, resulting in high overall application cost, which is difficult to meet the needs of engineering batch modeling and efficient simulation. SUMMARY

[0004] In order to solve the problems existing in the prior art, the application provides a bedding structure rock finite element model reconstruction method based on a digital core image.

[0005] The technical scheme adopted by the application to solve its technical problems is:

[0006] The technical scheme provides a bedding structure rock finite element model reconstruction method based on a digital core image, including the following steps:

[0007] S1: obtaining a digital image of a rock sample and converting the digital image into a gray-scale image;

[0008] S2: taking a vertical rectangular band as an analysis region with the center line of the gray-scale image as the symmetry axis, calculating the transverse mean value of the analysis region, and obtaining a longitudinal gray-scale distribution curve;

[0009] S3: applying a sliding average filter to eliminate high-frequency jitter of the longitudinal gray-scale distribution curve, and then applying a smoothing filter to further denoise, to obtain a smoothed gray-scale curve;

[0010] S4: first-order difference is made to the smoothed gray scale curve, positions with absolute difference value greater than a threshold value are taken as initial strip boundary candidates, minimum interval and minimum width filtering rules are applied to filter, start and end pixels, width and average gray scale of each strip are determined;

[0011] S5: K-means clustering algorithm is used to partition the smoothed gray scale curve, a plurality of clustering label intervals are obtained, each clustering label interval corresponds to a physical region type, according to the height of the rock sample and the unit size, each clustering label interval is mapped to a physical interval based on FEM grid partition;

[0012] S6: based on the size of the rock sample, a sample finite element model is constructed, based on the grid mapping algorithm, each unit in the sample finite element model is renumbered, and finally a reconstructed finite element model with real partition information is generated.

[0013] Preferably, in S2, the analysis region includes a lithology section that can represent the typical bedding trend and the vertical variation rule of the cross-bedding strip, the horizontal pixel mean value of each row of the intercepted analysis region is calculated to obtain a longitudinal gray scale profile, including:

[0014] Let the original image gray scale matrix be I ( i , j ), after selecting the center vertical strip range, the gray scale mean value of each row is formula (1):

[0015] (1);

[0016] In the formula, represents the horizontal pixel set of the selected vertical strip, I represents the original image gray scale matrix, i =1,…, H is the vertical pixel row, j =1,…, W is the horizontal pixel column, y i represents the gray scale mean value of the i th row, N j represents the number of horizontal pixels.

[0017] Preferably, in S3, the longitudinal gray scale distribution curve is first-order sliding average, the symmetric window is used for sliding average, and the first-order smoothed curve is obtained, and the formula is formula (2):

[0018] (2);

[0019] In the formula, x i represents the gray scale mean value of the iThe 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 first i The output value is the row moving average of the pixels.

[0020] 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:

[0021] (3);

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

[0023] Preferably, in step S4, the smooth grayscale curve undergoes first-order difference processing, and the difference formula is expressed as:

[0024] d i = y i - y i-1 , i =2,……, N (4);

[0025] In the formula, y i Indicates the first i Smooth grayscale values ​​of each pixel row, N denotes Section length, di represents the difference result of the first i pixel row, y i-1 represents the smoothed gray value of the first i- pixel row;

[0026] The maximum value of the difference absolute value is taken as the adaptive threshold, and the formula is as follows:

[0027] (5);

[0028] In the formula, T represents the adaptive threshold, α represents the adjustment coefficient, α Typically, the value is between 0.05 and 0.2, which is used to distinguish significant mutations from background fluctuations.

[0029] Preferably, all positions satisfying |D|>T are taken as preliminary candidate points of the strip boundary, and the minimum interval and minimum width filtering rules are applied to eliminate false detections with a distance less than a set minimum number of separation pixels or a corresponding strip width shorter than a set value. Based on the final boundary set, the starting pixel and the ending pixel of each strip are constructed d i s k , e k ] and the width of the strip is calculated w k = e k - s k +1, and the average gray value of the interval is calculated, and the formula is as follows:

[0030] (6);

[0031] In the formula, s k represents the starting pixel row number of the first k strip, e k represents the ending pixel row number of the first k strip, w k represents the width of the first k strip, represents the average gray value of the first k strip.

[0032] ​Preferably, in S5, the smoothed and denoised smoothed gray scale curve is partitioned using a K-means clustering algorithm to obtain a plurality of clustering label intervals, each clustering label interval corresponding to a physical region type, and a clustering label change position is detected in order based on the clustering result to determine the upper and lower boundaries of pixels of each region.

[0033] Let the total pixel height of the rock sample in the image be H p , the actual physical height be H , and the conversion coefficient of pixel to physical height be λ = H / H p The starting pixel s k and the ending pixel e k of each interval are mapped to the actual height interval λ [ s k -1) λ ,( e k -1)], and are aligned with the preset element height in the finite element model h elem to match each image partition to one or more element layers of the FEM grid, thereby realizing automatic mapping of the gray scale image partition to the finite element grid and attribute assignment of the physical region.

[0034] Preferably, in S6, the node coordinates and element connection relationship of the sample finite element model are read, and the coordinates of the eight nodes of the first i element in the longitudinal direction are z i,1 , z i,2 , …, z i,8 , and the calculation formula of the geometric center height of the element is:

[0035] (7);

[0036] In the formula, Z i,j represents the coordinates of the first i node of the first j element in the longitudinal direction, Z i (c) represents the geometric center height of the first i finite element, j represents the index of the node in the element, i represents the index of the finite element.

[0037] Preferably, the starting height of each partition is loaded from the external partition information file a p , the ending height b p and the corresponding partition type number T p , constitute an ordered triplet a p , b p , T p , wherein p =1, …, P represents all partitions, for each unit i , all partitions are traversed in turn, for the p th partition, it is judged whether the unit center height satisfies:

[0038] (8) ;

[0039] If yes, the part number part i of the unit is assigned to the type number of the partition, represented as part i = T p , and the partition traversal of the unit is terminated;

[0040] If Z i (c) = T p , let part i = T p , after all finite element units are processed, a reconstructed finite element model with real partition information is generated.

[0041] Compared with the prior art, the present application has the following advantages:

[0042] 1. Based on the macroscopic rock surface image shot by a digital camera, the present application can accurately identify various types of strip and matrix intervals through image gray feature analysis and partition algorithm, and can realize efficient FEM grid unit placement by combining automatic data processing and physical parameter mapping, which can not only significantly reduce the reconstruction cost, but also fully retain the macroscopic partition characteristics and engineering related structural details of the bedded structure rock, thereby providing a low-cost, efficient and strong structure authenticity model basis for the numerical study of the mechanical behavior of bedded structure rock.

[0043] 2. The present application effectively overcomes the dependence on expensive CT / SEM equipment and micro-scale reconstruction in the prior art, and can directly utilize ordinary photos combined with gray profile analysis and automatic partitioning algorithm to efficiently identify and extract the macroscopic structural characteristics of bedded structure rocks such as siliceous band dolomite, greatly improve the accuracy of macroscopic structure partitioning such as band and matrix through noise reduction processing and clustering partitioning, and can directly output physical intervals and types suitable for FEM modeling, realize automatic mapping from pixels to finite element model, greatly reduce the modeling cost and technical threshold, and fill the gap in the macroscopic structure FEM reconstruction link of the prior art. BRIEF DESCRIPTION OF DRAWINGS

[0044] The above and / or additional aspects and advantages of the present application will become apparent and more readily appreciated from the following description of the embodiments, taken in conjunction with the accompanying drawings in which:

[0045] Figure 1 is a whole flowchart of the present application;

[0046] Figure 2 is a schematic diagram of the central vertical analysis area of the siliceous band dolomite sample in Example Two;

[0047] Figure 3 is a gray original profile partitioning diagram of the siliceous band dolomite sample in Example Two;

[0048] Figure 4 is a gray profile partitioning diagram after noise reduction of the siliceous band dolomite sample in Example Two;

[0049] Figure 5 is a schematic diagram of the finite element model of the siliceous band dolomite sample in Example Two;

[0050] Figure 6 is a schematic diagram of the central vertical analysis area of the siltstone-mudstone interbedded sample in Example Three;

[0051] Figure 7 is a gray original profile partitioning diagram of the siltstone-mudstone interbedded sample in Example Three;

[0052] Figure 8 is a gray profile partitioning diagram after noise reduction of the siltstone-mudstone interbedded sample in Example Three;

[0053] Figure 9 is a schematic diagram of the finite element model of the siltstone-mudstone interbedded sample in Example Three. DETAILED DESCRIPTION

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

[0055] Example 1

[0056] 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:

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

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

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

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

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

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

[0063] In S2, the analysis area includes a lithology section that can represent the vertical variation law of typical layering trend and cross-banded, to further ensure the representativeness of the central strip, in the selection and verification process of the candidate strip, the adaptive check is introduced, the transverse pixel mean value of the intercepted analysis area is calculated along each row, the longitudinal gray profile is obtained, including:

[0064] Let the original image gray matrix be I i j , after selecting the central vertical strip range, the gray mean value of each row is

[0065] (1);

[0066] In the formula, represents the transverse pixel set of the selected vertical strip, I represents the original image gray matrix, i =1,…, H is the longitudinal pixel row, j =1,…, W is the transverse pixel column, y i represents the gray mean value of the i th row, N j represents the number of transverse pixels.

[0067] In S3, the first-order sliding average is performed on the longitudinal gray distribution curve to suppress high-frequency fluctuations, and the symmetric window is used for sliding average to obtain the first-order smoothed curve , which is expressed by the formula:

[0068] (2);

[0069] In the formula, x i represents the original gray value of the i th pixel row, x i+j represents the gray value corresponding to the left and right shift of i rows with the current row k as the center, m represents the half-length parameter of the sliding window, m the value range is 2 to 5, 2 m +1 represents the length of the symmetric window used for sliding average, and the corresponding window length 2m+1 is between 5 and 11, to balance between denoising and trend keeping, represents the sliding average output value of the i th pixel row.

[0070] The first-order smoothed curve ​​The input is smoothed by a Savitzky-Golay filter, and the two-stage noise reduction process gradually weakens high-frequency and medium-frequency disturbances while retaining the main trend of the gray scale to the maximum extent, providing a stable and high signal-to-noise ratio input signal for subsequent strip partition boundary positioning and lithology differentiation. An odd-length local window is used, and the polynomial order p A local least squares polynomial fitting is performed at each center position to obtain a smoothed gray scale curve, which is represented by the formula:

[0071] (3) ;

[0072] In the formula, c k The selected window length 2 M +1 and the polynomial order p (2 to 4 are selected) are pre-calculated by least squares fitting, representing the weighted combination weight of the center point under this local fitting, y i SG represents the gray scale value after two-stage smoothing, M represents the half-window width, M The value range is 7 to 25, corresponding to a window of 15 to 51, k represents the relative offset within the window, represents the i + k pixel row after one-stage smoothing.

[0073] In S4, the smoothed gray scale curve is subjected to first-order difference processing, and the difference formula is represented as:

[0074] d i = y i - y i-1 , i =2,……, N (4) ;

[0075] In the formula, y i represents the smoothed gray scale value of the i pixel row, N denotes profile length, d i represents the difference result of the i pixel row, y i-1 represents the smoothed gray scale value of the i- 1 pixel row;

[0076] The maximum value of the absolute value of the difference is taken as a proportion of the adaptive threshold, and the formula is represented as:

[0077] (5);

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

[0079] 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:

[0080] (6);

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

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

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

[0084] 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:

[0085] (7);

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

[0087] 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:

[0088] (8);

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

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

[0091] Example 2

[0092] 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:

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

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

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

[0096] S23: Apply sliding average filter to eliminate high-frequency jitter of longitudinal gray distribution curve, and then apply smoothing filter to further denoise, to obtain a smoothed gray curve that retains the main trend and suppresses medium-high frequency interference, providing high-quality signals for partition boundary identification.

[0097] As shown in Figure 3 , Figure 3 is the gray original profile partitioning diagram of the siliceous banded dolomite sample in Example Two. First, the original gray profile (horizontal average value sequence in pixel units) is processed by two smoothing methods: one is short window sliding average, and the other is Savitzky-Golay filter with a larger window.

[0098] The sliding average method can remove high-frequency “spur” noise, while the Savitzky-Golay filter retains a more complete overall trend. The combination of the two can, on the premise of ensuring the fidelity of the main trend of the curve, maximally suppress local random fluctuations, providing high-quality signals for subsequent boundary detection and partitioning. The processing results are shown in Figure 4 , Figure 4 is the gray profile partitioning diagram after noise reduction of the siliceous banded dolomite sample in Example Two.

[0099] S24: Based on the average gray of each band, automatically set the threshold value to preliminarily distinguish all bands into siliceous bands (low gray), siliceous dolomite (medium gray), and dolomite matrix (high gray), to achieve initial division of macroscopic lithology types;

[0100] First, calculate the average gray of each curve according to the upper and lower boundaries of each band interval, and take the overall average of the average gray of all bands as the division scale, to obtain two automatic threshold lines corresponding to the gray levels of 105.7 and 162.5, respectively. Then, color the area below 105.7 as blue, representing siliceous bands; color the area with gray values between 105.7 and 162.5 as green, representing siliceous dolomite; and color the area with gray values above 162.5 as red, representing dolomite matrix, thus completing the preliminary division of macroscopic lithology types.

[0101] S25: In the multi-interval automatic segmentation link, K-means clustering is used to divide the gray curve into three zones, and combined with the physical height and unit size parameters, the image pixel segments are mapped to the physical intervals required for finite element modeling, to realize automatic partitioning from image to model;

[0102] Combined with the total height of the sample 5 cm and the pixel-centimeter ratio of the image row number, each pixel row is mapped to the physical height, and the clustering partition number is assigned to the corresponding unit of the finite element model, to automatically complete the “part number” reconstruction from two-dimensional image to three-dimensional FEM unit.

[0103] S26: Constructing the finite element model of the sample in accordance with the actual sample size, saving the information of each unit and node of the sample, and using LS-PrePost to newly create a cylindrical geometry with a diameter of 50 mm and a height of 50 mm in the "Geometry" module thereof, which is completely consistent with the actual core sample size.

[0104] Subsequently, in the "Mesh" module, a three-dimensional entity eight-node unit grid is generated according to the set unit height of 0.04 cm and the circumferential grid density, ensuring that the width of each unit is close to the engineering requirement. After the meshing is completed, all node numbers and their three-dimensional coordinates, as well as the correspondence between each unit and its eight nodes, are exported and saved into a text file, thereby constructing the complete finite element model node information and unit topology information.

[0105] According to the principle of the grid mapping algorithm, combined with the regional partition information obtained in S25, each unit in the sample finite element model is renumbered to realize the automatic allocation of the unit component number. First, each physical height interval obtained in step S25 is matched with the height mapping relationship corresponding to each row of units during meshing. By using the exported node and unit information, the center height of each unit is compared with the pre-divided regional interval to determine its type of partition.

[0106] Specifically, the program iterates through all units in the model in sequence, reads the center height value of the first unit, and finds the partition index that satisfies "interval start height < unit center height < interval end height" in the interval array. Then, the component number corresponding to the unit in LS-PrePost is set to the lithology number corresponding to the partition. The result is shown in i Figure 5 Figure 5 is a schematic diagram of the finite element model of the siliceous banded dolomite sample in Example 2.

[0107] Example Three

[0108] As shown in Figures 6-9 , the present embodiment proposes a method for reconstructing the finite element model of a bedded structure rock based on a core digital image. Taking a certain layer series of siltstone-mudstone interbedded sample as an example, the method includes the following steps:

[0109] S31: A digital image is obtained by shooting the surface of the siltstone-mudstone interbedded sample with a digital camera. After being imported, the image is uniformly converted into a gray mode for subsequent processing.

[0110] S32: The vertical strip on the left side of the center of the image is automatically selected as the analysis area, as shown in Figure 6 Figure 6 ​​​is a schematic diagram of the central vertical analysis area of the siltstone-mudstone interbedded sample in Example Three, and the longitudinal gray scale distribution curve is obtained by calculating the transverse mean value, which reflects the gray scale characteristics of different lithology intervals;

[0111] First, the lateral surface of the siltstone-mudstone interbedded sample is vertically photographed in the laboratory environment using a high-resolution digital camera, ensuring uniform light source and no shadow, and fixing the relative position of the camera and the sample to obtain repeatable images. After the shooting is completed, the original color image is imported into MATLAB and converted into a single-channel gray scale image to remove the interference of color information and only keep the gray scale distribution.

[0112] S33: Apply sliding average filtering to eliminate high-frequency jitter of the longitudinal gray scale distribution curve, and then apply smoothing filtering to further denoise, obtaining a smoothed gray scale curve that retains the main trend and suppresses medium and high frequency interference, providing high-quality signals for boundary identification.

[0113] As shown in Figure 7 , the gray scale original profile zoning diagram of the siltstone-mudstone interbedded sample in Example Three is shown in Figure 7 . First, the original gray scale profile (transverse mean value sequence in pixels) is processed by two smoothing methods: short window sliding average and Savitzky-Golay filtering with a larger window.

[0114] The sliding average method can remove high-frequency "spur" noise, while the Savitzky-Golay filtering retains a more complete overall trend. The combination of the two can maximize the suppression of local random fluctuations while ensuring the integrity of the main trend of the curve, providing high-quality signals for subsequent boundary detection and zoning. The processing results are shown in Figure 8 , and Figure 8 is the denoised gray scale profile zoning diagram of the siltstone-mudstone interbedded sample in Example Three.

[0115] S34: Based on the average gray scale of each strip, the threshold is automatically set, and all strips are preliminarily divided into siltstone (low gray scale) and mudstone (high gray scale), realizing the initial division of macroscopic lithology type;

[0116] First, calculate the average gray scale of each curve according to the upper and lower boundaries of each strip interval, and take the overall mean value of the average gray scale of all strips as the division scale, thus obtaining an automatic threshold line with a gray scale level of 71.3. Then, the area below 71.3 is colored blue, representing siltstone, and the area above 71.3 is colored green, representing mudstone, thus completing the preliminary division of macroscopic lithology type.

[0117] S35: Multi-interval automatic segmentation link, using K-means clustering to divide the gray curve into two zones, and combining physical height and unit size parameters, mapping the image pixel section to the physical interval required for finite element modeling, realizing automatic partitioning from image to model;

[0118] Combined with the total height of the sample 4cm and the pixel-centimeter ratio of the image row number, each pixel row is mapped to the physical height, and the clustering partition number is assigned to the corresponding unit of the finite element model, and the "part number" reconstruction from two-dimensional image to three-dimensional FEM unit is automatically completed.

[0119] S36: Constructing a sample finite element model conforming to the actual sample size, saving the sample unit and node information, using LS-PrePost to create a cylindrical geometry with a diameter of 25mm and a height of 40mm in its "Geometry" module, which is completely consistent with the actual core sample size.

[0120] Then in the "Mesh" module, generate a three-dimensional solid eight-node element grid according to the set unit height of 0.2mm and the circumferential grid density, ensure that the width of each unit is close to the engineering requirement, after the meshing is completed, export all the node numbers and their three-dimensional coordinates, as well as the corresponding relationship between each unit and its eight nodes to a text file, thereby constructing the complete finite element model node information and unit topology information.

[0121] According to the principle of grid mapping algorithm, combined with the regional partition information obtained in S35, renumbering each unit in the sample finite element model to realize the automatic allocation of each unit part number, first match each physical height interval obtained in S35 with the height mapping relationship of each row of units during meshing: use the node and unit information exported in S206 to compare the center height of each unit with the pre-divided regional interval to determine its partition type.

[0122] Specifically, the program traverses all units in the model in turn, reads the center height value of the first i unit, and finds the partition index that satisfies "interval start height < unit center height < interval end height" in the interval array, then sets the part number corresponding to the unit in LS-PrePost to the lithology number corresponding to the partition, the result is shown in Figure 9 , Figure 9 is the finite element model diagram of the siltstone and mudstone interbedded sample in example three.

[0123] While embodiments of the application have been shown and described, it is to be understood that the embodiments described are merely divergences of the principles and application of the present application and that numerous modifications, changes, substitutions, and alterations can be made thereto without departing from the spirit and scope of the present application, which is defined by the following claims and their equivalents.

Claims

1. A method for reconstructing a finite element model of a rock with a layered structure based on a digital image of a core, characterized in that, The method comprises the following steps: S1: obtaining a rock surface image taken by a digital camera and converting the image into a gray image; S2: taking a vertical rectangular strip as an analysis region with the center line of the gray image as a symmetric axis, calculating the lateral mean value of the analysis region, and obtaining a longitudinal gray distribution curve; S3: applying a sliding average filter to eliminate high-frequency jitter of the longitudinal gray distribution curve, and then applying a smoothing filter to further denoise, to obtain a smoothed gray curve; S4: performing first-order differentiation on the smoothed gray curve, taking positions with an absolute difference value greater than a threshold value as initial strip boundary candidates, applying minimum interval and minimum width filtering rules to filter, and determining the start and end pixels, width, and average gray value of each strip; S5: partitioning the smoothed gray curve based on a K-means clustering algorithm to obtain a plurality of clustering label intervals, each clustering label interval corresponding to a physical region type, and mapping each clustering label interval to a physical interval based on FEM grid partitioning according to the height of the rock sample and the unit size; S6: constructing a sample finite element model based on the size of the rock sample, renumbering each unit in the sample finite element model based on a grid mapping algorithm, and finally generating a reconstructed finite element model with real partition information; In S5, the K-means clustering algorithm is used to partition the smoothed gray curve after smoothing and denoising, to obtain a plurality of clustering label intervals, each clustering label interval corresponding to a physical region type, and the clustering label change position is detected in order based on the clustering result to determine the upper and lower boundaries of the pixels of each region. Let the total pixel height of the rock sample in the image be H p , the actual physical height be H , and the conversion coefficient of pixel to physical height be λ = H / H p The starting pixel s k and the ending pixel e k of each interval are mapped to the actual height interval λ [ s k -1) λ , e k -1) ], and are aligned with the preset element height in the finite element model h elem to make each image partition match one or more element layers of the FEM grid, thereby realizing automatic mapping of the gray-scale image partition to the finite element grid and attribute assignment of the physical region. In the S6, the node coordinates and the element connection relationship of the sample finite element model are read, and the coordinates of the eight nodes of the first element in the longitudinal direction are i z i,1 , z i,2 ,… z i,8 The calculation formula of the height of the geometric center of the element is represented as:​ （7）； wherein Z i,j denotes the coordinate of the i th node of the j th element in the longitudinal direction, Z i (c) denotes the height of the geometric center of the i th finite element, j denotes the index of the node within the element, i denotes the index of the finite element. loading the start height of each partition from the outer partition information file a p , the end height b p and the corresponding partition type number T p , to form an ordered triplet ( a p , b p , T p ), wherein p =1, …, P denotes all partitions, for each cell i , all partitions are traversed in sequence, for the p th partition, it is determined whether the cell center height satisfies: （8）； If so, the component number of the unit is assigned part i is assigned the type number of the partition, denoted as part i = T p and the partition traversal of the unit is terminated; If the condition is satisfied Z i (c) = T p , let part i = T p , after the unit processing of all finite elements is completed, a reconstructed finite element model with real partition information is generated.

2. The method according to claim 1, wherein, In S2, the analysis region includes a lithology section that can represent the typical bedding strike and the vertical variation law of the cross strip, the lateral pixel mean value of the analysis region is calculated along each row to obtain a longitudinal gray profile, including: Let the original image gray matrix be I ( i , j ), select the center vertical band range, and the average gray value of each row is expressed by the formula: （1）； wherein denotes the set of horizontal pixels of the selected vertical band, I denotes the original image gray scale matrix, i = 1,..., H is the vertical pixel row, j = 1,..., W is the horizontal pixel column, y i denotes the average gray scale value of the i row, N j denotes the number of horizontal pixels.

3. The method according to claim 1, wherein, In the S3, a first-order sliding average is performed on the longitudinal gray scale distribution curve, a symmetric window is adopted for the sliding average, and a first-order smooth curve is obtained The formula is represented 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 j 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 first i The output value is the sliding average of the pixel rows.

4. The method according to claim 3, wherein, The first order smoothing curve is obtained by fitting a polynomial of order 3 to the data points in the local window of length 5. The second order smoothing is performed by inputting the Savitzky-Golay filter with an odd length local window and a polynomial order of 3. p The local least square polynomial fitting is performed at each center position to obtain the smoothed gray level curve, which is expressed as: （3）； wherein c k by the selected window length 2 M + 1 and the polynomial order p by a pre-computed least squares fit, representing the weighted combination of the center points under this local fit, y i SG representing the gray value after the second smoothing, M representing the half window width, k representing the relative offset within the window, representing the first i + k pixel row after the first smoothing.

5. The method of claim 1, wherein, In S4, the smoothed gray curve is subjected to first-order differentiation, and the differentiation formula is represented as: d i = y i － y i-1 , i =2,……, N （4）； wherein y i represents the smoothed gray value of the i th pixel row, N represents profile length, d i represents the difference result of the i th pixel row, y i-1 represents the smoothed gray value of the i- 1th pixel row; A proportion of the maximum value of the absolute difference is taken as an adaptive threshold, and the formula is represented as: （5）； wherein T denotes an adaptive threshold, α denotes an adjustment coefficient, α Typically, the value is between 0.05 and 0.2 for distinguishing significant mutations from background fluctuations.

6. The method according to claim 5, wherein, All meet d i The position of T as a strip boundary preliminary candidate point, the minimum interval and minimum width filtering rules are applied to eliminate the error detection of the distance less than the set minimum separation pixel number or the corresponding strip width shorter than the set value, based on the final boundary set, the starting pixel and the terminal pixel of each strip are constructed s k , e k ], the width of which is calculated w k = e k - s k +1, the average gray scale of the interval is counted, and the formula is expressed as: （6）； wherein s k represents the start pixel row number of the k e k represents the end pixel row number of the k w k represents the width of the k represents the average gray value of the k k represents the average gray value of the​​

Images