Material pore multi-scale characterization based on electron microscope images and regional b-spline curve compression penetration path prediction method

Abstract

A method for multi-scale characterization of material pores based on electron microscopy images and prediction of compressive propagation paths using B-spline curves is characterized by the following steps: First, the scanning electron microscope images are discretized and binarized using MATLAB. Then, pore identification and multi-parameter extraction are performed based on the binary matrix. Next, undetected and newly formed pores are generated based on fractal dimension. Then, the pore regions are divided using K-means clustering. Finally, the crack propagation path is predicted based on the mechanical principle that large pores preferentially propagate under pressure. This invention abandons simple linear fitting methods and uses weighted fitting of B-spline curves based on pore type and size weights. This method not only has clear physical meaning, reflecting the principle of macropore preferential propagation, but also generates a smooth curve as the predicted path, which can more realistically simulate the tortuous propagation behavior of internal cracks in materials, significantly improving prediction accuracy and reliability.

Description

Technical Field

[0001] This invention relates to a method for predicting the propagation of micro-regional cracks, and more particularly to a method for multi-scale characterization of material porosity based on electron microscopy images and prediction of the compressive penetration path of regional B-spline curves, belonging to the technical field of material microstructure analysis and performance prediction. Background Technology

[0002] The pore structure of porous materials (such as pore size, distribution, morphology, and connectivity) directly determines their key indicators such as mechanical properties and permeability. Accurately characterizing the pore structure and predicting fracture pathways under pressure is of great significance for material selection and structural optimization in engineering applications. Currently, pore analysis methods based on microscopic observation techniques such as scanning electron microscopy (SEM) are widely used, but existing technologies still have many shortcomings:

[0003] In terms of pore characterization, traditional binarization processing is prone to missed or misjudged pores due to unreasonable mesh division. Moreover, the classification of pores often relies on a single size threshold, lacking comprehensive consideration of multi-dimensional parameters such as the fractal characteristics of pore distribution and directional consistency, making it difficult to fully reflect the multi-scale structural characteristics of pores. At the same time, existing methods often ignore unscanned micropores and newly formed pores that may be generated during the compression process, resulting in incomplete pore count and affecting the accuracy of subsequent analysis.

[0004] In predicting penetration paths, existing technologies mostly rely on simple fitting based on pore center coordinates, failing to fully consider the impact of differences in pore size and type on penetration priority. Furthermore, the fitted curves are prone to problems such as unevenness and unclear physical meaning, making it difficult to accurately reflect the true pattern of pore penetration under pressure. In addition, existing methods do not divide pores into regions and calculate independent fractal dimensions, failing to accurately capture the differentiated influence of pore structures in different regions on the penetration path, resulting in insufficient reliability of prediction results. The inability to accurately predict the damage evolution process of materials under pressure conditions introduces significant uncertainty into the safety assessment of engineering structures, and this technical challenge has remained unresolved. Therefore, a penetration path prediction method is needed that can achieve accurate multi-scale characterization of pores, supplement information on newly formed pores, and perform scientifically weighted fitting based on pore characteristics to overcome the shortcomings of existing technologies. Summary of the Invention

[0005] The purpose of this invention is to address the lack of effective quantitative means for the complex characteristics of pores in existing methods. These methods rely solely on a single size threshold for pore classification, failing to form a comprehensive quantitative system that covers multiple dimensions such as the fractal characteristics of pore distribution and directional consistency. Furthermore, they have not established a correlation prediction model from microscopic pore structure characteristics to macroscopic mechanical behavior. In particular, they cannot accurately capture the differentiated influence of different types of pores on crack propagation in specific microscopic regions, making it difficult to scientifically predict regional crack propagation paths. This severely restricts the depth of material performance research and the reliability of engineering applications.

[0006] To address this, this invention provides a method for characterizing the microscopic pore structure of defective materials and predicting regional crack propagation that is algorithmically sound, comprehensively representative, and accurately predictive. This method employs multi-scale image processing and machine learning algorithms. First, through refined binarization and pore selection of electron microscopy images, a classification system for macropores, mesopores, and micropores based on area thresholds is established. Then, a fractal feature extraction system based on pore classification is constructed, calculating the distribution fractal dimension and structural fractal dimension for the three types of pores and newly formed pores under pressure, achieving accurate quantification of complex pore features. Based on this, multi-dimensional information such as the fractal features, center coordinates, and size weights of multi-level pores is integrated to construct a crack propagation prediction model based on a single electron microscopy image. Through key steps such as K-means clustering region division and weighted B-spline curve fitting, the prediction accuracy of crack propagation paths within specific microscopic regions is significantly improved, providing a new technical solution for solving the problem of material damage prediction.

[0007] To achieve the above-mentioned objectives, the technical solution of this invention is: a method for multi-scale characterization of material porosity based on electron microscopy images and prediction of the compressive penetration path of B-spline curves, characterized by comprising the following steps:

[0008] Its characteristics include the following steps:

[0009] Step 1: Discretize and binarize the scanning electron microscope image using MATLAB: Set the lower left corner of the image as the origin of the coordinate system, divide it into several square grids with side lengths smaller than the minimum pore size, construct the SEM image discrete matrix, perform binary assignment based on the grid gray values, and identify potential pore regions through the matrix row and column distribution characteristics;

[0010] Step 2: Pore identification and multi-parameter extraction based on binary matrix: Calculate the coordinates of the pore center point, classify pores into three levels—macropores, mesopores, and micropores—accumulate the number of pores in the three categories, and then systematically calculate geometric parameters such as porosity, void ratio, pore perimeter, and equivalent diameter, as well as three types of fractal dimensions—pore structure, equivalent diameter, and distribution. At the same time, analyze the pore orientation characteristics through second-order central moments to obtain weighted average orientation angles, orientation consistency measures, and regional pore density, thereby systematically obtaining the geometric morphology and spatial distribution parameters of the pore system.

[0011] Step 3: Generate undetected and newly formed pores based on fractal dimension: Calculate the Euclidean distance between all pairs of pore center points and obtain the global average distance. Determine the number of newly formed pores according to the scaling factor. Generate micropores and newly formed pores that conform to the fractal law based on the existing pore distribution, and ensure that the minimum distance between the new pores and the original pores is not less than the grid side length to avoid overlap.

[0012] Step 4: Use the K-means clustering algorithm to divide the pores into regions: Based on the weighted Euclidean distance formula, initialize the cluster centers of macropores, mesopores, and micropores according to the pore area and spatial distribution characteristics. Iteratively assign pores to the nearest cluster center and update the center coordinates until convergence, aggregating the pores into three feature regions (dominant region (macropores), secondary dominant region (mesopores), and auxiliary region (micropores)). Repeat the fractal dimension calculation process for each region to obtain the pore distribution fractal dimension, structural fractal dimension, and equivalent diameter fractal dimension of each region.

[0013] Step 5: Based on the mechanical principle of preferential penetration of large pores under pressure, predict the crack propagation path: Assign weights to the four regions according to their priority (macropores > mesopores > micropores > newly formed pores), sort the pore center point set using the nearest neighbor method and parameterize the chord length, solve for the control points by weighted least squares fitting of B-spline curves, uniformly sample in the parameter domain [0,1] and use the de Boer algorithm to calculate the coordinates of each point on the curve, and fuse them into a smooth crack propagation path prediction curve, realizing multi-scale prediction from micropore structure to macro crack behavior.

[0014] Furthermore, the definition and weighted fitting method of the B-spline curve in step five are as follows:

[0015] (1) Define the curve structure:

[0016] Curve order: p=3 (cubic B-spline)

[0017] Number of control points: n (must satisfy p+1) <n<N)

[0018] Node vector U (U has n+p+1 elements): Using a uniform type, in the form of:

[0019] (40)

[0020] In the formula: internal node {u p+1 ,…,u m-p-1} are uniformly distributed in the interval [0,1].

[0021] (Example: p=3, n=5, then U=(0,0,0,0,0.5,1,1,1,1))

[0022] (2) Cox-de Boor recursive formula:

[0023] This formula should be defined recursively in layers:

[0024] 1) Basis functions of degree 0 (p=0):

[0025] (41)

[0026] 2) Basis functions of degree p (p>0):

[0027] It is obtained by recursive calculation using the lower-order basis function;

[0028] (42)

[0029] In the formula: if the denominator is 0, the result of the entire fraction is defined as 0;

[0030] (3) Construction of the fundamental function matrix:

[0031] (43)

[0032] (4) Normalization verification:

[0033] 1) Verification conditions:

[0034] (44)

[0035] 2) Numerical tolerance check:

[0036] (45)

[0037] 3) Handling verification failures:

[0038] If the sum of any row is not equal to 1 (out of tolerance), then:

[0039] a. Re-examine the definition of node vector U;

[0040] b. Verify the implementation of the Cox-de-Boor recursive algorithm;

[0041] c. Check parameter t i Whether the range of values ​​is within the domain of the node vector;

[0042] (5) Weighted least squares fitting:

[0043] To find the control point P (of size n×2), we need to solve a system of linear equations:

[0044] (46)

[0045] In the formula: N is the basis function matrix; W is an N×N diagonal weight matrix, where the diagonal elements are the weights of each pore corresponding to the ordered point set, and Q... s It is an ordered set of points.

[0046] Furthermore, in step one, the lower left corner of the image is set as the origin of the coordinate system, and the image is divided into a square grid with a side length of a, where the side length of the grid a is less than the minimum pore size. The grid gray value g is then binarized and assigned a value: when g ≤ the threshold, a value of 1 is assigned to represent a pore, and when g > the threshold, a value of 0 is assigned to represent the matrix.

[0047] Furthermore, in step two, the coordinates of the pore center (X...) r Y r The pore area S is obtained by calculating the average of the coordinates of all grids with a value of 1 within the pore; the pore area S is obtained by multiplying the number of grids with a value of 1 N within the pore by the grid area a².

[0048] Furthermore, in step two, the box size sequence ε is defined. k = a 2 k (k=0,1,2,…,K), where K satisfies εk≤ min(image length L, width W); for each εk, count the number N(εk) boxes covering the apertures. k ); for the dataset (log ε) k ,log N(ε k A linear fit is performed, and the negative slope of the fitted line is the fractal dimension D.

[0049] Step two also includes fractal property verification: the correlation coefficient R² of the linear fit is required to be ≥ 0.95, and the fluctuation of the fractal dimension D in the scale-free interval is less than 5%.

[0050] Furthermore, in step three, the spacing E is generated. 新 = E(1 + β(D - D 基准 ), where E is the global average pore distance, D is the fractal dimension of the system, and D 基准 The baseline fractal dimension is β, where β is the scaling factor; the number of newly formed pores, N 生成 = γN, where N is the number of identified pores and γ is the proportionality coefficient.

[0051] Furthermore, in step four, initializing cluster centers includes: selecting the macropore center with the largest equivalent diameter, the medium pore center adjacent to the macropore, and the densely distributed micropore center; and using the weighted Euclidean distance formula to perform pore allocation and cluster center update.

[0052] Furthermore, in step five, the weight allocation method is as follows: for known pore sizes in the macropore, mesopore, and micropore regions, the weight ω is... i= ω 区域基础 (S i / S max ), where ω 宏基础 > ω 中基础 > ω 微基础 S i S is the area of ​​the pore. max The maximum area among all known pores is used; newly generated pores are assigned a uniformly lower weight.

[0053] Furthermore, in step five, the nearest neighbor method is used to sort the set of pore center points, and the chord length is parameterized to obtain the parameters. The B-spline curve is a cubic B-spline curve; by solving the linear equation system (N T WN)P =N T WQ s Obtain the control points P of the B-spline curve, where N is the basis function matrix, W is the diagonal weight matrix, and Q is the control point P. s The system consists of an ordered set of pore center points; a predicted curve is generated by uniformly sampling in the parameter domain and calculating the coordinates of the curve points using the de Boer algorithm.

[0054] The beneficial effects of this invention are:

[0055] 1. This invention achieves a comprehensive quantitative characterization of pore structures by refining the mesh and binarizing the data, combined with multi-scale parameter calculations (size, orientation, fractal dimension, etc.), thus avoiding the problems of missed detection and misjudgment in traditional methods.

[0056] 2. This invention innovatively introduces a new pore generation mechanism, and combines fractal dimension to adjust the generation spacing, supplementing the information of unscanned micropores and pressure-induced new pores, improving the completeness of pore statistics, and providing a more realistic microstructure basis for subsequent prediction.

[0057] 3. This invention uses K-means clustering to divide pore regions and calculates the independent fractal dimension of each region, accurately capturing the differentiated characteristics of pore structures in different regions.

[0058] 4. This invention abandons the simple linear fitting method and adopts a weighted fitting of B-spline curves based on pore type and size weights. This method not only has clear physical meaning and reflects the principle of prioritizing the penetration of macropores, but also generates a smooth curve as the predicted path, which can more realistically simulate the tortuous propagation behavior of internal cracks in materials, significantly improving prediction accuracy and reliability.

[0059] 5. The algorithm of this invention is reasonable, the method is standardized and highly repeatable, and it is applicable to a variety of porous materials. It provides reliable technical support for the evaluation of material mechanical properties and the safety design of engineering structures. It is applicable to the quantitative analysis of pore structure of porous materials such as rocks and composite materials, and can accurately predict the potential crack penetration path of materials under pressure conditions. Attached Figure Description

[0060] Figure 1 This is a binarized image from an embodiment of the present invention.

[0061] Figure 2 This is a diagram showing the pore shape and spatial distribution according to an embodiment of the present invention.

[0062] Figure 3 This is a fractal dimension diagram of the pore distribution according to an embodiment of the present invention.

[0063] Figure 4 This is a fractal dimension diagram of the equivalent diameter of the pores in an embodiment of the present invention.

[0064] Figure 5 This is a diagram of the pore fractal dimension parameters according to an embodiment of the present invention.

[0065] Figure 6 This is a pore distribution diagram of pores of different sizes relative to the center of the pore group in an embodiment of the present invention.

[0066] Figure 7 This is a pore statistics diagram of the size of the pores relative to the center of the pore group in an embodiment of the present invention.

[0067] Figure 8 This is a graph showing the pore characterization parameters of an embodiment of the present invention.

[0068] Figure 9 This is a graph showing the B-spline prediction results of an embodiment of the present invention. Detailed Implementation

[0069] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0070] See Figures 1 to 9 The present invention provides a method for multi-scale characterization of material porosity and prediction of the compressive penetration path of B-spline curves based on electron microscopy images, characterized by comprising the following steps:

[0071] Step 1: Discretize and binarize the scanning electron microscope image using MATLAB: Set the lower left corner of the image as the origin of the coordinate system, divide it into several square grids with side lengths smaller than the minimum pore size, construct the SEM image discrete matrix, perform binary assignment based on the grid gray values, and identify potential pore regions through the matrix row and column distribution characteristics;

[0072] Step 2: Pore identification and multi-parameter extraction based on binary matrix: Calculate the coordinates of the pore center point, classify pores into three levels—macropores, mesopores, and micropores—accumulate the number of pores in the three categories, and then systematically calculate geometric parameters such as porosity, void ratio, pore perimeter, and equivalent diameter, as well as three types of fractal dimensions—pore structure, equivalent diameter, and distribution. At the same time, analyze the pore orientation characteristics through second-order central moments to obtain weighted average orientation angles, orientation consistency measures, and regional pore density, thereby systematically obtaining the geometric morphology and spatial distribution parameters of the pore system.

[0073] Step 3: Generate undetected and newly formed pores based on fractal dimension: Calculate the Euclidean distance between all pairs of pore center points and obtain the global average distance. Determine the number of newly formed pores according to the scaling factor. Generate micropores and newly formed pores that conform to the fractal law based on the existing pore distribution, and ensure that the minimum distance between the new pores and the original pores is not less than the grid side length to avoid overlap.

[0074] Step 4: Use the K-means clustering algorithm to divide the pores into regions: Based on the weighted Euclidean distance formula, initialize the cluster centers of macropores, mesopores, and micropores according to the pore area and spatial distribution characteristics. Iteratively assign pores to the nearest cluster center and update the center coordinates until convergence, aggregating the pores into three feature regions (dominant region (macropores), secondary dominant region (mesopores), and auxiliary region (micropores)). Repeat the fractal dimension calculation process for each region to obtain the pore distribution fractal dimension, structural fractal dimension, and equivalent diameter fractal dimension of each region.

[0075] Step 5: Based on the mechanical principle of preferential penetration of large pores under pressure, predict the crack propagation path: Assign weights to the four regions according to their priority (macropores > mesopores > micropores > newly formed pores), sort the pore center point set using the nearest neighbor method and parameterize the chord length, solve for the control points by weighted least squares fitting of B-spline curves, uniformly sample in the parameter domain [0,1] and use the de Boer algorithm to calculate the coordinates of each point on the curve, and fuse them into a smooth crack propagation path prediction curve, realizing multi-scale prediction from micropore structure to macro crack behavior.

[0076] The definition and weighted fitting method of the B-spline curve in step five are as follows:

[0077] (1) Define the curve structure:

[0078] Curve order: p=3 (cubic B-spline)

[0079] Number of control points: n (must satisfy p+1) <n<N)

[0080] Node vector U (U has n+p+1 elements): Using a uniform type, in the form of:

[0081] (40)

[0082] In the formula: internal node {u p+1 ,…,u m-p-1} are uniformly distributed in the interval [0,1].

[0083] (Example: p=3, n=5, then U=(0,0,0,0,0.5,1,1,1,1))

[0084] (2) Cox-de Boor recursive formula:

[0085] This formula should be defined recursively in layers:

[0086] 1) Basis functions of degree 0 (p=0):

[0087] (41)

[0088] 2) Basis functions of degree p (p>0):

[0089] It is obtained by recursive calculation using the lower-order basis function;

[0090] (42)

[0091] In the formula: if the denominator is 0, the result of the entire fraction is defined as 0;

[0092] (3) Construction of the fundamental function matrix:

[0093] (43)

[0094] (4) Normalization verification:

[0095] 1) Verification conditions:

[0096] (44)

[0097] 2) Numerical tolerance check:

[0098] (45)

[0099] 3) Handling verification failures:

[0100] If the sum of any row is not equal to 1 (out of tolerance), then:

[0101] a. Re-examine the definition of node vector U;

[0102] b. Verify the implementation of the Cox-de-Boor recursive algorithm;

[0103] c. Check parameter t iWhether the range of values ​​is within the domain of the node vector;

[0104] (5) Weighted least squares fitting:

[0105] To find the control point P (of size n×2), we need to solve a system of linear equations:

[0106] (46)

[0107] In the formula: N is the basis function matrix; W is an N×N diagonal weight matrix, where the diagonal elements are the weights of each pore corresponding to the ordered point set, and Q... s It is an ordered set of points.

[0108] In step one, the lower left corner of the image is set as the origin of the coordinate system and divided into a square grid with a side length of a. The side length of the grid, a, is less than the minimum pore size. The grid gray value g is then binarized and assigned a value: when g ≤ the threshold, a value of 1 is assigned to represent a pore, and when g > the threshold, a value of 0 is assigned to represent the matrix.

[0109] In step two, the pore center coordinates (X...) r , Y r The pore area S is obtained by calculating the average of the coordinates of all grids with a value of 1 within the pore; the pore area S is obtained by multiplying the number of grids with a value of 1 N within the pore by the grid area a².

[0110] In step two, the box size sequence ε is defined. k = a 2 k (k=0,1,2,…,K), where K satisfies εk ≤ min(image length L, width W); for each εk, count the number N(εk) boxes covering the apertures. k ); for the dataset (log ε) k , log N(ε k A linear fit is performed, and the negative slope of the fitted line is the fractal dimension D.

[0111] Step two also includes fractal property verification: the correlation coefficient R² of the linear fit is required to be ≥ 0.95, and the fluctuation of the fractal dimension D in the scale-free interval is less than 5%.

[0112] In step three, the spacing E is generated. 新 = E(1 + β(D - D 基准 ), where E is the global average pore distance, D is the fractal dimension of the system, and D 基准 The baseline fractal dimension is β, where β is the scaling factor; the number of newly formed pores, N 生成 = γN, where N is the number of identified pores and γ is the proportionality coefficient.

[0113] In step four, initializing cluster centers includes: selecting the macropore center with the largest equivalent diameter, the medium pore center adjacent to the macropore, and the densely distributed micropore center; and using the weighted Euclidean distance formula to perform pore allocation and update cluster centers.

[0114] In step five, the weight allocation method is as follows: for the known pore sizes in the macropore, mesopore, and micropore regions, the weight ω is... i = ω 区域基础 (S i / S max ), where ω 宏基础 > ω 中基础 > ω 微基础 S i S is the area of ​​the pore. max The maximum area among all known pores is used; newly generated pores are assigned a uniformly lower weight.

[0115] In step five, the nearest neighbor method is used to sort the set of pore center points, and the chord length is parameterized to obtain the parameters. The B-spline curve is a cubic B-spline curve; by solving the linear equation system (N T WN)P = N T WQ s Obtain the control points P of the B-spline curve, where N is the basis function matrix, W is the diagonal weight matrix, and Q is the control point P. s The system consists of an ordered set of pore center points; a predicted curve is generated by uniformly sampling in the parameter domain and calculating the coordinates of the curve points using the de Boer algorithm.

[0116] The following is a detailed description: The binarization method in step one is as follows:

[0117] The lower left corner of the quantized electron microscope image is set as the origin of the rectangular coordinate system, and the image is divided into three parts. For a square grid with side length 'a', it is necessary to ensure that the grid side length 'a' is less than the minimum pore size to avoid missed detections.

[0118] set up( If the origin of the rectangular coordinate system is (0, 0), then the coordinates of the center point of each small grid are ( ),and =(N-0.5)a, =(M-0.5)a.

[0119] Equation (1) is the discrete matrix of the SEM image, where the grayscale value at the center point of each discrete grid is [value missing]. .

[0120] (1)

[0122] Equation (2) is the binary assignment, with the black area assigned a value of 1 and the white area assigned a value of 0:

[0123] (2)

[0124] In the formula: g represents the grayscale value;

[0125] When the top row and the left and right columns of the matrix are all 0, it represents a gap, as shown in equation (3):

[0126] (3)

[0127] The binarized image is obtained by using MATLAB software in the following steps, and the binarized image is shown below. Figure 1 As shown.

[0128] The method for aperture identification and parameter calculation in step two is as follows:

[0129] 2.1) Pore screening conditions: In formula (3), except for the top row and the left and right columns which are all 0, the other parts are randomly distributed between 0 and 1, and there is no other part where a row or column is all 0;

[0130] 2.2) Extraction of pore center coordinates:

[0131] Assuming that the number of 1s in equation (3) is N, the coordinates of the pore center point are as shown in equation (4);

[0132] (4)

[0133] In the formula: and The coordinates of each pore that is 1;

[0134] 2.3) Pore Classification and Area Statistics:

[0135] Assuming that the number of 1s in the matrix of equation (3) is N, then the pore area is:

[0136] S=Na 2 (5)

[0137] In the formula: a is the grid side length;

[0138] 2.4) Pore perimeter

[0139] If the grid above a grid of 1 is 0 or does not exist (at the image edge), then N(i,j) is incremented by 1; if the grid below a grid of 1 is 0 or does not exist (at the image edge), then N(i,j) is incremented by 1; if the grid to the left of a grid of 1 is 0 or does not exist (at the image edge), then N(i,j) is incremented by 1; if the grid to the right of a grid of 1 is 0 or does not exist (at the image edge), then N(i,j) is incremented by 1.

[0140] (6)

[0141] In the formula: a is the grid side length; N(i,j) is the total number of calculations mentioned above.

[0142] 2.5) Pore equivalent diameter

[0143] From equation (5), the equivalent diameter d can be obtained. e for:

[0144] (7)

[0145] In the formula: S is the pore area.

[0146] From equation (7), we can obtain .

[0147] 2.6) Number of pores

[0148] To set multiple threshold levels:

[0149] Macrohole: d e ≥d3

[0150] Central hole: d2≤d e <d3

[0151] Micropores: d e <d1

[0152] d3, d2, and d1 are set based on experience.

[0153] 2.7) Fractal dimension of pore structure

[0154] 1) Extract the area S of each pore. i and its perimeter C i .

[0155] 2) Calculate for each pore and construct a dataset. .

[0156] 3) Calculate the fractal dimension D:

[0157] (8)

[0158] 2.8) Pore equivalent diameter and fractal dimension

[0159] 1) Data sorting:

[0160] The equivalent diameter set {d} derived from equation (7) i Arranged in descending order, the result is {d1 ≥ d2 ≥ d3 ≥ … ≥ d N (If pores with the same equivalent diameter appear, they should also be placed in the above sorting. For example, if {3, 2, 2, 1} is {d3, d2, d2, d1}.)

[0161] 2) Calculate the cumulative quantity:

[0162] For the j-th diameter d in the sorted sequence j, Define aperture threshold r j =d j Cumulative quantity N j For an equivalent diameter greater than or equal to r j The number of pores is used to obtain the data sequence {(r1,N1), (r2,N2),…,(r Ntotal N total 3) Logarithmic transformation:

[0163] Taking the natural logarithm of each data pair in the above data sequence yields a new sequence {(lnr1,lnN1),(lnr2, lnN2),…,(ln r... Ntotal , lnN total 4) Univariate linear regression: For data points {(lnr1,lnN1),(lnr2, lnN2),…,(ln r...}. Ntotal , lnN total If a least-squares linear fit is performed, the fractal dimension D is calculated using the following formula:

[0164] (9)

[0165] 2.9) Calculate the fractal dimension of pore distribution using the box counting method.

[0166] (1) Define the box size sequence

[0167] 1) Select box size

[0168] Define a series of box sizes It is usually a multiple of the initial grid side length 'a':

[0169] (10)

[0170] In the formula: K satisfies L and W are the length and width of the image.

[0171] 2) Coverage verification

[0172] Ensure that the box size covers a multi-scale range, from microscopic (individual pores) to macroscopic (overall image).

[0173] (2) Count the number of boxes covering the pores.

[0174] 1) Divide the boxes

[0175] For each box size Divide the image into A box.

[0176] 2) Determine if the box is covered

[0177] Define indicator functions :

[0178] (11)

[0179] 3) Calculate the coverage quantity

[0180] Statistics for each Number of boxes covering the pores:

[0181] (12)

[0182] (3) Calculation of fractal dimension D

[0183] 1) Logarithmic transformation:

[0184] Regarding box size and coverage Taking the natural logarithm (or the common logarithm) yields the dataset:

[0185] (13)

[0186] 2) Calculation:

[0187] The equation of the straight line is fitted using the least squares method, and the fractal dimension D is the negative slope of the fitted line:

[0188] (14)

[0189] In the formula: .

[0190] (4) Verification of fractal properties

[0191] 1) Linear correlation:

[0192] Requirement R 2 ≥0.95, otherwise it is necessary to check whether the box size selection or pore distribution has fractal characteristics.

[0193] 2) Scale-free intervals:

[0194] Within a certain box size range (e.g.) arrive The fluctuation of the fractal dimension D should be less than 5%.

[0195] 3) Local calculation of fractal dimension:

[0196] For each box size (k=2, 3, ..., K-2), using 3 adjacent points ( Fitting local slope D k .

[0197] (15)

[0198] 4) Calculate the mean and standard deviation:

[0199] Mean: (16)

[0200] Standard deviation: (17)

[0201] Fluctuation verification: (18)

[0202] 2.10) Angle of the pore's own orientation

[0203] The set of coordinates of the center of each grid with a pore size of 1:

[0204] (1) Calculate the second-order center distance:

[0205] (19)

[0206] (2) Calculate the pore direction angle:

[0207] The direction of the major axis of the best-fit ellipse for the porosity is calculated using the second-order central moment:

[0208] (20)

[0209] Note: The four-quadrant arctangent function should be used in actual calculations.

[0210] (twenty one)

[0211] (3) Directional statistical analysis:

[0212] 1) Directional distribution rose diagram:

[0213] a. Divide the range of -180° to 180° into equal parts.

[0214] b. Count the number of pores in each interval.

[0215] c. Draw a before-and-after histogram.

[0216] 2) Calculation of average direction:

[0217] a. Convert each angle to a unit vector:

[0218] (twenty two)

[0219] b. Calculate the average vector:

[0220] (twenty three)

[0221] c. Calculate the average direction:

[0222] (twenty four)

[0223] 3) Measurement of directional consistency:

[0224] a. Calculate the average vector length:

[0225] (25)

[0226] b. Calculate the circular variance:

[0227] (26)

[0228] In the formula: the closer R is to 1, the better the consistency of direction; the closer R is to 0, the more random the direction.

[0229] 2.11) Porosity and void ratio

[0230] Total black area (S) p The area of ​​the white region (S) is defined as the superposition of all regions assigned a value of 0; m If S is defined as the superposition of all regions assigned a value of 1, then S = ... p and S m for:

[0231] (27)

[0232] (28)

[0233] In the formula: S ij The area of ​​each discrete region.

[0234] Rock porosity and void ratio can be defined as:

[0235] (29)

[0236] (30)

[0237] Step two of the pore structure parameters was performed using MATLAB software, as follows: Figure 2 , Figure 3 , Figure 4 , Figure 5 , Figure 6 , Figure 7 , Figure 8 As shown in Table 1.

[0238] Table 1 shows the overall fractal dimension of the pores.

[0239] Pore ​​type <![CDATA[Equivalent diameter fractal dimension D1]]> <![CDATA[Fractal dimension D2 of pore structure]]> <![CDATA[Fractal dimension D3 of pore distribution]]> overall 2.519 1.380 1

[0240] The specific method for step three is as follows:

[0241] 3.1) Calculate the global average distance E:

[0242] The Euclidean distance between each pair of pore center points is calculated as shown in equation (31):

[0243] (31)

[0244] In the formula: This represents the distance between the center points of the j-th and k-th pores;

[0245] The Euclidean distances between each center point are calculated according to equation (31) to form a distance matrix. The average distance E is calculated by averaging all distance values ​​in the distance matrix.

[0246] (32)

[0247] In the formula: Number of pores;

[0248] 3.2) Fractal dimension adjustment to generate spacing:

[0249] (33)

[0250] In the formula: D 基准 The fractal dimension D is calculated using a small number of samples and then the mean is taken. The scaling factor is 0.1.

[0251] 3.3) Number of center coordinates generated

[0252] Empirical Proportioning Method: Based on the material type and research objectives, a proportionality coefficient γ is set between the number of newly formed pores and the number of identified pores N, i.e.:

[0253] (34)

[0254] In the formula: γ (0.1~0.3) is adjusted according to the material type and the degree of pressure.

[0255] 3.4) Generation method:

[0256] With E 新 Combining existing pore centers and pore distribution patterns with r 生成 , with E 新 The center coordinates of unscanned micropores and newly formed micropores are generated to determine the spacing; if the existing pores are directionally aligned, they are preferentially generated along that direction; the generated pores must maintain a minimum spacing d with the original pores. min =a, to avoid overlap;

[0257] The results of step three are obtained using MATLAB software, such as... Figure 8 As shown.

[0258] The method for dividing the pores into regions in step four is as follows:

[0259] In the cluster analysis method, the K-means algorithm is used to divide the porosity into a dominant region, a secondary dominant region, and an auxiliary region 1. The auxiliary region 2 (newly formed pores) is automatically divided at the same time as the coordinates of the center point of the new pores are generated.

[0260] Weighted Euclidean distance formula:

[0261] (35)

[0262] In the formula: d p d q This is the corresponding equivalent pore diameter.

[0263] 4.1) Gap region division

[0264] Initialize cluster centers:

[0265] Macropore: Select the pore center P1 with the largest equivalent diameter.

[0266] Medium pore: Select the medium pore center P2 adjacent to the macropore.

[0267] Micropores: Select the densely distributed micropore center P3.

[0268] Assigned to the nearest cluster center P K : Coordinates of the center point of each pore ( , Using equation (35), the distance to each cluster center is calculated, and the pores are assigned to the nearest cluster center.

[0269] Update cluster centers: For each cluster P K Calculate the average coordinates of all its points, as shown in equation (36).

[0270] (36)

[0271] In the formula: For clustering Number of pore coordinates.

[0272] Iteration: Repeat the above allocation and update steps until the cluster centers no longer change or the preset number of iterations is reached.

[0273] Region division: Based on the clustering results, the pore center points are divided into the dominant region (macropores), the secondary dominant region (mesopores), and the auxiliary region 1 (micropores).

[0274] 4.2) Calculation of Independent Fractal Dimension

[0275] Repeat steps 7-10 in step two above for each region to obtain the directional angle distribution, pore structure, pore equivalent diameter, and pore distribution fractal dimension of the three regions.

[0276] Step four using MATLAB software yields the fractal dimensions of the three types of pores, as shown in the following results. Figure 3 , Figure 4 , Figure 5 , Figure 6 , Figure 7 As shown in Table 2.

[0277] Table 2 shows the fractal dimension of pore classification.

[0278] Pore ​​type <![CDATA[Equivalent diameter fractal dimension D1]]> <![CDATA[Fractal dimension D2 of pore structure]]> <![CDATA[Fractal dimension D3 of pore distribution]]> micropores 1.809 1.485 1 central hole 2.020 1.409 1 Hong Kong 2.476 1.314 1

[0279] The specific method for step five is as follows:

[0280] To predict the potential penetration path of materials under pressure, this method comprehensively considers multi-scale information such as pore center coordinates, size and type. By introducing a weight matrix to reflect the principle of macropore priority penetration, and using cubic B-spline curves for fitting, a smooth and physically meaningful predicted path is finally obtained.

[0281] In step five, the example uses cubic B-spline curve prediction, and the specific method is as follows:

[0282] 5.1) Data Processing and Weight Allocation

[0283] Suppose there are N pore center points, and their coordinate set is: This collection includes identified macropores, mesopores, micropores, and newly formed pores.

[0284] Assign weights ω to each data pointi Its calculation is based on the area of ​​the pore region and the area of ​​the (known pores):

[0285] (37)

[0286] In the formula: S i S represents the area of ​​the i-th known pore; max ω represents the maximum area among all known pores. 宏 ω 中 ω 微 ω 新 The regional basic weight coefficients satisfy ω 宏 >ω 中 >ω 微 >ω 新 (For example: 1, 0.7, 0.4, 0.2 can be adjusted according to different situations); z is the scaling factor for newly formed pores (generally 0.1, which can be adjusted according to different situations), used to assign a uniform low weight to newly formed pores without area data.

[0287] 5.2) Point Set Sorting and Parameterization

[0288] (1) The pore point set Q is randomly distributed in two-dimensional space and needs to be sorted to establish the mapping relationship between parameter t and spatial location. The nearest neighbor method is used for sorting:

[0289] 1) Initialization: Start with the coordinates of the pore center with the largest weight and set it to φ1.

[0290] 2) Substitution: Starting from φ1, use equation (31) to find the point with the smallest Euclidean distance among the remaining points as φ2, and then find the point with the smallest Euclidean distance from φ2 among the remaining center coordinate points as φ 3， And so on. This leads to the ordered point Q. s ={φ1,φ2,φ3,….,φ N}

[0291] (2) Perform chord length parameterization on the ordered point sequence and assign parameters to each point. :

[0292] 1) Calculate the chord length: Use equation (31) to calculate the Euclidean distance d between adjacent points. i (i=1,2,…,N-1).

[0293] 2) Calculate the cumulative distance H i Total chord length L:

[0294] (38)

[0295] 3) Calculate the parameter value t i :

[0296] (39)

[0297] 5.3) Definition and Weighted Fitting of B-Spline Curves

[0298] (1) Define the curve structure:

[0299] Curve order: p=3 (cubic B-spline)

[0300] Number of control points: n (must satisfy p+1) <n<N)

[0301] Node vector U (U has n+p+1 elements): Using a uniform type, in the form of:

[0302] (40)

[0303] In the formula: internal node {u p+1 ,…,u m-p-1} are uniformly distributed in the interval [0,1].

[0304] (Example: p=3, n=5, then U=(0,0,0,0,0.5,1,1,1,1))

[0305] (2) Cox-de Boor recursive formula:

[0306] This formula should be a hierarchical recursive definition:

[0307] 1) Basis functions of degree 0 (p=0):

[0308] (41)

[0309] 2) Basis functions of degree p (p>0):

[0310] It is obtained by recursive calculation using the lower-order basis functions.

[0311] (42)

[0312] In the formula: if the denominator is 0, the result of the entire fraction is defined as 0.

[0313] (3) Construction of the fundamental function matrix:

[0314] (43)

[0315] (4) Normalization verification:

[0316] 1) Verification conditions:

[0317] (44)

[0318] 2) Numerical tolerance check:

[0319] (45)

[0320] 3) Handling verification failures:

[0321] If the sum of any row is not equal to 1 (out of tolerance), then:

[0322] a. Re-examine the definition of node vector U.

[0323] b. Verify the implementation of the Cox-de-Boor recursive algorithm.

[0324] c. Check parameter t i Whether the range of values ​​is within the domain of the node vector.

[0325] (5) Weighted least squares fitting:

[0326] To find the control point P (of size n×2), we need to solve a system of linear equations:

[0327] (46)

[0328] In the formula: N is the basis function matrix; W is an N×N diagonal weight matrix (the diagonal elements of the matrix are the weights of each pore corresponding to the ordered point set); Q s It is an ordered set of points.

[0329] 5.4) Predicting the generation of the through path curve

[0330] (1) Sampling: Uniformly sample N within the parameter domain t∈[0,1]. s N points s Ideally, it should be 100.

[0331] (47)

[0332] (2) Calculate the curve coordinates: for each sampling parameter t k The coordinates of points on the B-spline curve are calculated using the de Boer algorithm.

[0333] 1) Given conditions:

[0334] a. Set of control points:

[0335] b. Node vectors:

[0336] c. Curve order: p=3 (cubic B-spline)

[0337] d. Parameter values: Total N ssampling points

[0338] 2) De Boer's algorithm recursive formula: For each t k Sample parameters and perform the following steps.

[0339] a. Identify the active range:

[0340] That is, determine t k The node interval that falls within will be used as the starting interval for recursion: .

[0341] b. Initialize control points:

[0342] Extracting local control points: .

[0343] c. Recursive calculation:

[0344] For r=1~p:

[0345] (48)

[0346] In the formula: j = i - p + r, ..., i; if the denominator is 0, then a j,r =0.

[0347] d. Obtain the points on the curve:

[0348] (49)

[0349] In the formula: if t k =0, then d i,p =P0; if t k =, then d i,p =P n-1 .

[0350] From equation (49), we can obtain {C(t0), C(t1), ..., C(t)} Ns-1 )}, link all points in sequence to form a smooth prediction curve.

[0351] Step four of the MATLAB software was used to obtain the B-spline prediction curve, and the result is as follows: Figure 9 As shown.

[0352] The above description is a further detailed explanation of the present invention in conjunction with specific embodiments. It should not be considered that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple substitutions, improvements and changes can be made to the present invention without departing from the concept of the present invention. All such simple substitutions, improvements and changes should be considered to fall within the protection scope of the present invention.

Claims

1. A method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images, characterized in that... Includes the following steps: Step 1: Discretize and binarize the scanning electron microscope image using MATLAB: Set the lower left corner of the image as the origin of the coordinate system, divide it into several square grids with side lengths smaller than the minimum pore size, construct the SEM image discrete matrix, perform binary assignment based on the grid gray values, and identify potential pore regions through the matrix row and column distribution characteristics; Step 2: Pore identification and multi-parameter extraction based on binary matrix: Calculate the coordinates of the pore center point, classify pores into three levels—macropores, mesopores, and micropores—accumulate the number of pores in the three categories, and then systematically calculate geometric parameters such as porosity, void ratio, pore perimeter, and equivalent diameter, as well as three types of fractal dimensions—pore structure, equivalent diameter, and distribution. At the same time, analyze the pore orientation characteristics through second-order central moments to obtain weighted average orientation angles, orientation consistency measures, and regional pore density, thereby systematically obtaining the geometric morphology and spatial distribution parameters of the pore system. Step 3: Generate undetected and newly formed pores based on fractal dimension: Calculate the Euclidean distance between all pairs of pore center points and obtain the global average distance. Determine the number of newly formed pores according to the scaling factor. Generate micropores and newly formed pores that conform to the fractal law based on the existing pore distribution, and ensure that the minimum distance between the new pores and the original pores is not less than the grid side length to avoid overlap. Step 4: Use the K-means clustering algorithm to divide the pores into regions: Based on the weighted Euclidean distance formula, initialize the cluster centers of macropores, mesopores, and micropores according to the pore area and spatial distribution characteristics. Iteratively assign pores to the nearest cluster center and update the center coordinates until convergence, aggregating the pores into three feature regions: the dominant region (macropores), the secondary dominant region (mesopores), and the auxiliary region (micropores). Repeat the fractal dimension calculation process for each region to obtain the pore distribution fractal dimension, structural fractal dimension, and equivalent diameter fractal dimension of each region. Step 5: Based on the mechanical principle of preferential penetration of large pores under pressure, predict the crack propagation path: Assign weights to the four regions according to their priority (macropores > mesopores > micropores > newly formed pores), sort the pore center point set using the nearest neighbor method and parameterize the chord length, solve for the control points by weighted least squares fitting of B-spline curves, uniformly sample in the parameter domain [0,1] and use the de Boer algorithm to calculate the coordinates of each point on the curve, and fuse them into a smooth crack propagation path prediction curve, realizing multi-scale prediction from micropore structure to macro crack behavior.

2. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: The definition and weighted fitting method of the B-spline curve in step five are as follows: (1) Define the curve structure: Curve order: p=3 (cubic B-spline) Number of control points: n (must satisfy p+1) <n<N) Node vector U (U has n+p+1 elements): Using a uniform type, in the form of: （40） In the formula: internal node {u p+1 ,…,u m-p-1 } are uniformly distributed in the interval [0,1]. (Example: p=3, n=5, then U=(0,0,0,0,0.5,1,1,1,1)) (2) Cox-de Boor recursive formula: This formula should be defined recursively in layers: 1) Basis functions of degree 0 (p=0): （41） 2) Basis functions of degree p (p>0): It is obtained by recursive calculation using the lower-order basis function; （42） In the formula: if the denominator is 0, the result of the entire fraction is defined as 0; (3) Construction of the fundamental function matrix: （43） (4) Normalization verification: 1) Verification conditions: （44） 2) Numerical tolerance check: （45） 3) Handling verification failures: If the sum of any row is not equal to 1 (out of tolerance), then: a. Re-examine the definition of node vector U; b. Verify the implementation of the Cox-de-Boor recursive algorithm; c. Check parameter t i Whether the range of values ​​is within the domain of the node vector; (5) Weighted least squares fitting: To find the control point P (of size n×2), we need to solve a system of linear equations: （46） In the formula: N is the basis function matrix; W is an N×N diagonal weight matrix, where the diagonal elements are the weights of each pore corresponding to the ordered point set, and Q... s It is an ordered set of points.

3. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step one, the lower left corner of the image is set as the origin of the coordinate system and divided into a square grid with a side length of a. The side length of the grid, a, is less than the minimum pore size. The grid gray value g is then binarized and assigned a value: when g ≤ the threshold, a value of 1 is assigned to represent a pore, and when g > the threshold, a value of 0 is assigned to represent the matrix.

4. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step two, the coordinates of the pore center (X) r Y r The pore area S is obtained by calculating the average of the coordinates of all grids with a value of 1 within the pore; the pore area S is obtained by multiplying the number of grids with a value of 1 N within the pore by the grid area a².

5. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step two, the box size sequence ε is defined. k = a 2 k (k=0,1,2,…,K), where K satisfies εk ≤ min(image length L, width W); for each εk, count the number N(εk) boxes covering the apertures. k ); for the dataset (log ε) k , log N(ε k A linear fit is performed, and the negative slope of the fitted line is the fractal dimension D.

6. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: Step two also includes fractal property verification: the correlation coefficient R² of the linear fit is required to be ≥ 0.95, and the fluctuation of the fractal dimension D in the scale-free interval is less than 5%.

7. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step three, the spacing E is generated. 新 = E(1 + β(D - D 基准 ), where E is the global average pore distance, D is the fractal dimension of the system, and D 基准 The baseline fractal dimension is β, where β is the scaling factor; the number of newly formed pores, N 生成 = γN, where N is the number of identified pores and γ is the proportionality coefficient.

8. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step four, initializing cluster centers includes: selecting the macropore center with the largest equivalent diameter, the medium pore center adjacent to the macropore, and the densely distributed micropore center; and using the weighted Euclidean distance formula to perform pore allocation and update cluster centers.

9. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step five, the weight allocation method is as follows: for the known pore sizes in the macropore, mesopore, and micropore regions, the weight ω is... i = ω 区域基础 (S i / S max ), where ω 宏基础 > ω 中基础 > ω 微基础 S i S is the area of ​​the pore. max The maximum area among all known pores is used; newly generated pores are assigned a uniformly lower weight.

10. The method for multi-scale characterization of material porosity and prediction of compressive penetration path of regional B-spline curves based on electron microscopy images according to claim 1, characterized in that: In step five, the nearest neighbor method is used to sort the set of pore center points, and the chord length is parameterized to obtain the parameters. The B-spline curve is a cubic B-spline curve; by solving the linear equation system (N T WN)P = N T WQ s Obtain the control points P of the B-spline curve, where N is the basis function matrix, W is the diagonal weight matrix, and Q is the control point P. s The system consists of an ordered set of pore center points; a predicted curve is generated by uniformly sampling in the parameter domain and calculating the coordinates of the curve points using the de Boer algorithm.

Images