Porous medium sub-resolution porosity estimation method based on single-energy CT (Computed Tomography)

Abstract

The invention relates to a porous medium sub-resolution porosity estimation method based on single-energy CT, which comprises the following steps: based on a known database, acquiring a mass attenuation coefficient and density of a sample mineral, and calculating a linear attenuation coefficient, a relative linear attenuation coefficient upper limit and a relative linear attenuation coefficient lower limit of a sample mineral component; calculating initial porosity distribution according to the linear attenuation coefficient obtained by CT and the relative linear attenuation coefficient upper limit and the relative linear attenuation coefficient lower limit obtained by calculation; and generating porosity distribution data, and outputting the three-dimensional space model to improve the porosity estimation precision of the CT image. According to the method, the porosity is calculated by using the linear attenuation coefficient of the mineral component, the assumption of uniformity of the mass attenuation coefficient is broken through, the porosity estimation precision is improved, the application range is widened, and the stability of subsequent work is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of porous medium characterization, and particularly to a method for estimating the sub-resolution porosity of a porous medium based on monoenergetic CT, which is particularly suitable for the analysis of porous media with complex mineral compositions and cross-scale pore structures. Background Art

[0002] Porosity is a core parameter for evaluating the pore characteristics of rock materials. It not only affects the accuracy of pore flow analysis and mechanical simulation, but also the description of its distribution heterogeneity is crucial for predicting microscale flow, heat conduction, and chemical erosion processes. Traditional methods usually use macroscopic porosity or the total porosity at the sample scale as evaluation indicators, but it is difficult to reflect the complexity of the pore distribution at the mesoscale. X-ray CT technology can achieve three-dimensional visualization of pore structures through micron-scale resolution imaging, becoming the mainstream method for mesoscale porosity analysis. In CT data processing, traditional methods quantify porosity through binary segmentation (i.e., distinguishing solids and pores). However, the existence of sub-resolution pores (i.e., pores smaller than the CT resolution) leads to significant errors. If they are classified as solids, the pore connectivity is underestimated; if they are classified as pores, the total porosity is overestimated. Although such pores are small, they play a key role in processes such as fluid transport and chemical erosion, and there is an urgent need for a cross-scale characterization method.

[0003] To address this challenge, according to the representative volume element (REV) concept, the description of the micro-structure is simplified by defining the average properties at a specific scale. Based on the gray value or linear attenuation coefficient (μl) of CT images, existing studies have proposed estimating sub-resolution porosity by hysteresis segmentation of components. However, such methods have the following key defects: (1) The assumption of the uniformity of the mass attenuation coefficient fails: Existing methods default that the mass attenuation coefficient (μm) of rock minerals is uniform, ignoring the differences between minerals. With the development of high-resolution CT technology, as the X-ray energy increases or heavy metal elements are included, the μm differences between mineral grains and the matrix are significant, resulting in an exacerbation of the porosity calibration error based on a single μm. (2) Calibration benchmark deviation: Assuming that the linear attenuation coefficient is only related to density and not considering the differences in μm in a multi-mineral system, the upper limit calibration of the relative linear attenuation coefficient of the matrix is incorrect, ultimately affecting the accuracy of porosity estimation.

[0004] In view of this, the present invention is specifically proposed. Summary of the Invention

[0005] In view of the deficiencies of the prior art, the present invention provides a method for estimating the sub-resolution porosity of a porous medium based on monoenergetic CT. By calculating the linear attenuation coefficient of the sample minerals, and then calculating the porosity of the sample minerals to generate porosity distribution data and output a three-dimensional space model, it solves the problem in the prior art that the pores of rock materials cannot be accurately analyzed, thereby making it difficult to carry out subsequent work.

[0006] The technical solution of the present invention is as follows:

[0007] The present invention provides a method for estimating the sub-resolution porosity of a porous medium based on monoenergetic CT, which includes the following steps:

[0008] Based on a known database, obtain the mass attenuation coefficient and density of the sample mineral, and calculate the linear attenuation coefficient, upper limit of relative linear attenuation coefficient, and lower limit of relative linear attenuation coefficient of the sample mineral components;

[0009] Calculate the initial porosity distribution through the upper limit of relative linear attenuation coefficient and the lower limit of relative linear attenuation coefficient obtained from the above calculations;

[0010] Generate porosity distribution data and output a three-dimensional space model to improve the porosity estimation accuracy of the CT image.

[0011] The present invention calculates the upper and lower limits of the relative linear attenuation coefficient of the sample mineral components through the mass attenuation coefficient and density of the sample mineral, then calculates the porosity according to the calculation results, and finally outputs a three-dimensional space model according to the porosity, which can detect the sub-resolution pores of the matrix, greatly improving the calculation accuracy of the sample mineral porosity and providing more accurate data support for subsequent work.

[0012] Preferably, before obtaining the mass attenuation coefficient and density of the sample mineral, the following steps are further included:

[0013] Search for representative slices of the sample mineral, and use the k-means algorithm to perform clustering analysis on the attenuation of the representative slices to obtain K ranges of linear attenuation coefficients and cluster centers.

[0014] Preferably, calculating the linear attenuation coefficient of the sample mineral components includes the following steps:

[0015] Divide the sample mineral into three groups, namely macropores, large particles, and matrix. Then, according to the corresponding X-ray energy in the database, select the corresponding mass attenuation coefficient and density, and substitute the obtained mass attenuation coefficient and density into formula (1) to obtain the linear attenuation coefficient of each component of the sample mineral in the crystal state (theoretical pore-free state):

[0016] μ l =μ m ·ρ (1);

[0017] Where, μl is the linear attenuation coefficient; μm is the mass attenuation coefficient; ρ is the density, with the unit g / cm 3 .

[0018] Preferably, the calculation of the lower limit of the relative linear attenuation coefficient includes the following steps:

[0019] Scanning the air near the sample mineral, the obtained result is the lower limit of the relative linear attenuation coefficient;

[0020] Preferably, the calculation of the upper limit of the relative linear attenuation coefficient includes the following steps:

[0021] Under the same X-ray energy, the ratio of the linear attenuation coefficients of the large particles and the matrix in the crystalline state is constant. Combining formula (1) and formula (2), the upper limit of the relative linear attenuation coefficient of the matrix is obtained:

[0022]

[0023] where, μ l_max is the upper limit of the relative linear attenuation coefficient of the matrix, with the unit of cm 2 / g; μ m_m is the mass attenuation coefficient of the matrix; μ m_p is the mass attenuation coefficient of the large particles; ρ m is the density of the matrix mineral in the crystalline state, with the unit of g / cm 3 ; ρ p is the density of the large particle mineral in the crystalline state, with the unit of g / cm 3 ; μ l_p is the clustering center of the linear attenuation coefficient of the large particles, with the unit of cm -1 .

[0024] Preferably, the calculation of the initial porosity distribution includes the following steps:

[0025] Calculating the porosity through formula (3):

[0026]

[0027] where, is the porosity; μ l_j is the linear attenuation coefficient at position j, with the unit of cm -1 ; μ l_max is the linear attenuation coefficient corresponding to the matrix in the crystalline state, with the unit of cm -1 ; μ l_a is the maximum value of the linear attenuation coefficient corresponding to the macropores, with the unit of cm -1 .

[0028] Preferably, before outputting the three-dimensional space model, the following steps are further included:

[0029] Analyzing the pore size distribution information of the sample mineral and correcting the macropore ratio;

[0030] Correcting the total porosity.

[0031] Preferably, the correction of the macropore ratio includes the following steps:

[0032] Define the pore size range of the macropores;

[0033] Determine the relative volume content of the macropores, and calculate the pore content of the macropore spatial resolution, i.e., the macropore porosity;

[0034] Mark the positions of the macropores.

[0035] Preferably, the correction of the total porosity includes the following steps:

[0036] Deduct the macropore porosity from the total porosity of the sample mineral to determine the measured matrix porosity;

[0037] Calculate the measured matrix pore volume by multiplying the pixel volume of the sample mineral by the measured matrix porosity;

[0038] Calculate the uncorrected matrix pore volume according to formula (4);

[0039] Calculate the matrix porosity correction coefficient according to formula (5);

[0040] Apply the matrix porosity correction coefficient to the uncorrected matrix porosity according to formula (6) for correction;

[0041] Where:

[0042]

[0043]

[0044] Where, V mtx is the uncorrected matrix pore volume; is the sub-resolution pore rate spatial distribution matrix of the matrix; V mtx_c is the measured matrix pore volume; C is the matrix porosity correction coefficient; is the sub-resolution pore rate spatial distribution matrix of the corrected matrix.

[0045] Preferably, when the matrix contains multiple mineral components, the weighted average method is used to calculate the equivalent mass attenuation coefficient of the matrix sub-class according to the volume fraction of different components, and the calculation method is shown in formula (7):

[0046]

[0047] Where, μ m_c is the equivalent mass attenuation coefficient of the matrix sub-class, with the unit cm 2 / g; μ m_k is the mass attenuation coefficient of the kth mineral under this matrix sub-class, with the unit cm 2 / g; r k is the proportion of the kth mineral.

[0048] Preferably, when the matrix contains multiple mineral components, the crystal state density of the same element ratio of the matrix subclass is calculated according to the element ratio in the chemical formula of the component, and the calculation method is as shown in formula (8):

[0049]

[0050] where ρ c is the crystal state density of the same element ratio of the matrix subclass, with the unit g / cm 3 ; ρ k is the density of the k-th mineral under this matrix subclass, with the unit g / cm 3 ; S T refers to the total area of all minerals in the slice, with the unit cm 2 ; S K refers to the area of the k-th mineral, with the unit cm 2 .

[0051] The beneficial effects of the present invention compared with the prior art are as follows:

[0052] 1. By setting the upper and lower limits of the relative linear attenuation coefficient of the matrix and calculating the porosity of the matrix in combination with the CT scan results, the present invention breaks through the assumption of the uniformity of the mass attenuation coefficient, improving the estimation accuracy and applicable range of the porosity.

[0053] 2. The present invention integrates the elimination of closed pores and cross-scale characterization, comprehensively considering pores below the micron level. The present invention proposes a multi-source data correction method to ensure the stability of the method. The present invention is applicable to high-energy CT and cases containing heavy metal elements, solving the bottleneck of emerging technologies.

[0054] It should be understood that the implementation of any embodiment of the present invention does not mean that multiple or all of the above beneficial effects need to be simultaneously achieved or reached. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only exemplary, and for those of ordinary skill in the art, without creative efforts, other implementation drawings can be obtained according to the provided drawings.

[0056] The structures, ratios, sizes, etc. shown in this specification are only used to cooperate with the content disclosed in the specification for those familiar with this technology to understand and read, and are not used to limit the limiting conditions for the implementation of the present invention. Therefore, they do not have technical substance. Any modification of the structure, change of the proportional relationship, or adjustment of the size, without affecting the effects that the present invention can produce and the purposes that can be achieved, should still fall within the scope covered by the technical content disclosed in the present invention.

[0057] Figure 1 Process flow diagram of the sub-resolution porosity estimation method for porous media based on monoenergetic CT provided by the embodiments of the present invention;

[0058] Figure 2 Schematic diagrams of three connected structures, where (a) is the D3Q7 connected structure, (b) is the D3Q19 connected structure, and (c) is the D3Q27 connected structure;

[0059] Figure 3 Schematic diagram of the three-dimensional voxel model of the sample mineral, where (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2;

[0060] Figure 4 Schematic diagram of the representative slices of the sample mineral, where (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2;

[0061] Figure 5 Schematic diagram of the clustering results of the sample mineral, where (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2;

[0062] Figure 6 Schematic diagram of the variation trend of the mass attenuation coefficient of the mineral with the X-ray energy;

[0063] Figure 7 Schematic diagram of the porosity distribution of the representative slices, where (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2;

[0064] Figure 8 Schematic diagram of the three-dimensional space and porosity distribution of the sample mineral, where (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2. Detailed implementation manners

[0065] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer and more understandable, the embodiments of the present invention will be further described in detail below in conjunction with the embodiments and the drawings. Herein, the illustrative embodiments of the present invention and their descriptions are used to explain the present invention, but not to limit the present invention.

[0066] In the present invention, unless otherwise clearly specified or defined, terms such as "installed", "connected", "linked", "fixed", etc. shall be understood in a broad sense. For example, it can be a fixed connection, a detachable connection, or integrated; it can be a mechanical connection or an electrical connection; it can be directly connected or indirectly connected through an intermediate medium, and it can be the communication inside two components or the interaction relationship between two components. For those of ordinary skill in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0067] It should be understood that the terms "comprising / including", "consisting of", or any other variant is intended to cover non-exclusive inclusion, so that a product, device, process, or method including a series of elements not only includes those elements, but also includes other elements that are not explicitly listed when needed, or further includes elements inherent to such product, device, process, or method. Without further limitation, the elements defined by the statement "comprising / including..." or "consisting of..." do not exclude the existence of additional identical elements in the product, device, process, or method including the said elements.

[0068] It is also necessary to understand that the terms "upper", "lower", "front", "rear", "left", "right", "top", "bottom", "inner", "outer", etc. indicating the orientation or position relationship are based on the orientation or position relationship shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device, component, or structure referred to must have a specific orientation, be constructed or operated in a specific orientation, and cannot be understood as a limitation to the present invention.

[0069] In addition, the terms "first" and "second" are only used for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, the features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the present invention, the meaning of "a plurality" is two or more, unless otherwise clearly and specifically defined.

[0070] The present invention relates to a method for estimating the sub-resolution porosity of porous media based on monoenergetic CT, which calculates its linear attenuation coefficient through the mass attenuation technique and density of the sample minerals, and then calculates the porosity through the linear attenuation coefficient, and outputs a three-dimensional space model, solving the problem that the measurement accuracy of porosity in the prior art is inaccurate and affects the implementation of subsequent work.

[0071] The following describes the implementation of the present invention in detail in combination with preferred embodiments.

[0072] The present invention provides a method for estimating the sub-resolution porosity of porous media based on monoenergetic CT, as Figure 1As shown, the method includes the following steps:

[0073] (1) Perform component segmentation - divide the components of the sample into three types: macropores, large particles, and matrix based on the linear attenuation coefficient of X-ray CT:

[0074] First, search for representative slices for analysis. This is because the CT data volume is large. If clustering and partitioning are performed on all slices, the computational workload is too large. By analyzing a small number of representative slices and applying it to the component partitioning of the entire 3D model, the workload can be reduced; statistically analyze the range of linear attenuation coefficients of each slice, and use the slice with the widest coverage range and obvious characteristics as the representative slice (if there are multiple slices meeting the conditions, all are included in the analysis scope; considering the influence of extreme linear attenuation coefficients, consider including the surrounding area of the representative slice in the consideration option) to ensure that the representative slice is representative in terms of component segmentation;

[0075] Next, perform clustering analysis on the linear attenuation coefficients of the representative slices, including the following steps:

[0076] a. Expand the linear attenuation coefficient matrix of the representative slice into multiple groups of arrays;

[0077] b. Cluster the data through the k-means algorithm, perform clustering at least twice in the order of natural numbers, and obtain the linear attenuation coefficient range and clustering center of each category;

[0078] c. Select the number of clusters according to the maximum clustering principle, and identify the regions where large particles and their subclasses are located (in actual operation, dividing the sample mineral components into three parts: macropores, matrix, and large particles is sufficient for analyzing the porosity distribution, but in order to calibrate the upper limit of the relative linear attenuation coefficient of the matrix, the large particle subclasses are distinguished in the present invention). It should be noted that since in actual operation, large particles may contain multiple types, it is difficult to determine the number of clusters without performing mineral analysis. Therefore, when selecting the clustering result, as much as possible select a high number of clusters while also ensuring that a single mineral is not divided into multiple classes. Determine the clustering number of large particles.

[0079] d. After determining the clustering number of large particles, the large particles and the matrix can be distinguished. Select the large particles of the known sample mineral type and the corresponding linear attenuation coefficients for matrix boundary calibration. At the same time, in order to ensure that the matrix with a lower porosity can be distinguished, after selecting the clustering number of large particles, it is necessary to distinguish the particle gaps through the clustering merging results considering a larger K value (number of clusters).

[0080] Since the particle gaps are too narrow, their linear attenuation coefficients are affected by the surrounding solid phase and increase. The number of clusters selected according to the maximum clustering principle may cause adjacent particles to bond. Therefore, selecting the clustering result with a larger K value can make the particle gaps significantly more obvious to a certain extent.

[0081] e. After classifying into three categories of macropores, large particles, and matrix, further detect and label the closed pores in the macropores and matrix. This is because the macropores, large particles, and matrix are classified by the linear attenuation coefficient, and the closed pores inside the large particles are included in the macropores and matrix. Although neither NMR (Nuclear Magnetic Resonance, a spectroscopy technique based on the interaction between the spin characteristics of atomic nuclei and an external magnetic field, widely used in fields such as chemistry, materials science, and medicine) nor mercury intrusion porosimetry for structure detection can include closed pores, especially since these closed pores mainly originate from the interior of the sample minerals. Therefore, in the present invention, all the units connected to the outside are marked in the macropore and matrix arrays, and the remaining unmarked units are the closed pores.

[0082] For the connected units, in the present invention, referring to the definition of the discrete velocity model in the lattice method, three connection modes are provided, as Figure 2 shown. These connection modes respectively correspond to the D3Q7, D3Q19, and D3Q27 discrete velocity models in LBM (Lattice Boltzmann Method), where (a) is D3Q7, (b) is D3Q19, and (c) is D3Q27. In actual operation, the connection mode can be selected according to specific requirements for identifying closed pores.

[0083] (2) Calculate the linear attenuation coefficient of the matrix:

[0084] First, identify the mineral components and proportions contained in the sample minerals. Different minerals contained in the matrix can be distinguished by dual-energy X-ray CT (X-ray Computed Tomography) or XRD (X-ray Diffraction). Use the NIST XCOM database (a standard reference database developed by the National Institute of Standards and Technology in the United States, focusing on the calculation of photon-matter interaction data, covering an energy range from 1 keV to 100 GeV, and supporting photon cross-section analysis of elements, compounds, and mixtures) to query the mass attenuation coefficient of the identified minerals in the crystalline state, and select the corresponding mass attenuation coefficient of the minerals according to the X-ray energy used for detecting the minerals; use the Mindat.org database (with the core advantages of open collaboration, multi-language coverage, and dynamic update, becoming an authoritative data platform in the field of mineralogy) to query the density of the identified minerals. Then, based on the Beer-Lambert law (providing the core theoretical basis for spectrophotometric analysis through the linear relationship between absorbance, concentration, and optical path), since the linear attenuation coefficient of large particles is relatively concentrated and stable, while the linear attenuation coefficient of the matrix is usually more dispersed, due to the porosity difference, substitute the mass attenuation coefficient and density of large particles and the matrix in the crystalline state into formula (1) to calculate the linear attenuation coefficient in the theoretically pore-free state:

[0085] μ l= μ m ·ρ (1);

[0086] Wherein, μ l is the linear attenuation coefficient; μ m is the mass attenuation coefficient; ρ is the density, with the unit g / cm 3 .

[0087] (3) Calculate the upper and lower limits of the matrix linear attenuation coefficient:

[0088] In CT image analysis, the lower limit μ l_min of the matrix linear attenuation coefficient uses the result in air during the current scan as a reference.

[0089] The calculation method of the upper limit μ l_max of the matrix linear attenuation coefficient includes the following steps:

[0090] If the mineral component in the matrix is a single component, the calculation method includes the following steps:

[0091] Since under the same X-ray energy, the ratio of the linear attenuation coefficients of large particles and the matrix in the crystal state is constant, combining formula (1) and formula (2) in step (2) can calculate the upper limit of the matrix linear attenuation coefficient:

[0092]

[0093] Wherein, μ l_max is the upper limit of the relative linear attenuation coefficient of the matrix, with the unit cm 2 / g; μ m_m is the mass attenuation coefficient of the matrix; μ m_p is the mass attenuation coefficient of the large particles; ρ m is the density of the matrix mineral in the crystal state, with the unit g / cm 3 ; ρ p is the density of the large particle mineral in the crystal state, with the unit g / cm 3 ; μl _p is the clustering center of the large particle linear attenuation coefficient, with the unit cm -1 .

[0094] When the mineral components in the matrix consist of multiple components, the backscattered color technique can be combined to more accurately identify the mineral components and proportions of the matrix. According to the volume fractions of different components, the weighted average method is used to calculate the equivalent mass attenuation coefficient of the matrix subclass (Formula (7)). Similarly, the crystal state density of the same element ratio in the matrix subclass is also calculated using the weighted average method (Formula (8)). It should be noted that if slices are used for mineral composition identification, the volume fraction should be changed to the area fraction accordingly. After calculating the equivalent mass attenuation coefficient of the matrix subclass and the crystal state density of the same element ratio in the matrix subclass, substitute them into Formulas (1) and (2) for calculation:

[0095]

[0096] where, μ m_c is the equivalent mass attenuation coefficient of the matrix subclass, with the unit cm 2 / g; μ m_k is the mass attenuation coefficient of the k-th mineral under this matrix subclass, with the unit cm 2 / g; r k is the proportion of the k-th mineral.;

[0097]

[0098] where, ρ c is the crystal state density of the same element ratio in the matrix subclass, with the unit g / cm 3 ; ρ k is the density of the k-th mineral under this matrix subclass, with the unit g / cm 3 ; S T refers to the total area of all minerals in the slice, with the unit cm 2 ; S K refers to the area of the k-th mineral, with the unit cm 2 .

[0099] (4) Calculate the porosity:

[0100] By determining the interval range of the matrix linear attenuation coefficient and its corresponding porosity, linear interpolation calculation can be performed on the porosity within this interval, thereby realizing the estimation of sub-resolution.

[0101] Based on the Beer-Lambert law, the linear attenuation coefficient shows a linear relationship with the density of the substance. When the X-ray energy is fixed, there is also a linear correlation between the linear attenuation coefficient (or the gray value of the CT image) and the porosity: when the matrix mineral is in the crystal state, its porosity can be regarded as 0%, and the corresponding linear attenuation coefficient is: When the porosity reaches 100%, the matrix mineral is completely occupied by macropores, and the corresponding linear attenuation coefficient is: Based on this linear relationship, by defining the mapping relationship between the upper and lower limits of the relative linear attenuation coefficient of the matrix and the porosity, the sub-resolution porosity of the matrix can be calculated as shown in formula (3):

[0102]

[0103] Wherein, is the porosity; μ l_j is the linear attenuation coefficient at position j, with the unit of cm -1 ; μ l_max is the linear attenuation coefficient corresponding to the matrix in the crystal state, with the unit of cm -1 ; μ l_a is the maximum value of the linear attenuation coefficient corresponding to the macropores, with the unit of cm -1 .

[0104] Wherein, in this embodiment, μ l_max is the linear attenuation coefficient corresponding to the matrix in the crystal state, where the crystal state is equivalent to the pore-free state; μ l_a is the upper limit of the relative linear attenuation coefficient corresponding to the macropores. The macropores correspond to air, and its value is the same as the lower limit of the linear attenuation coefficient of the matrix in step (2) of this embodiment.

[0105] (5) Correct the porosity:

[0106] 1) Correct the macropore ratio through pore size distribution information:

[0107] By using techniques such as NMR (nuclear magnetic resonance) and MIP (maximum density projection / mercury intrusion porosimetry) to obtain the pore size distribution information of the sample mineral, the limitations of the segmentation algorithm can be effectively compensated, thereby improving the accuracy of macropore and matrix segmentation. Segmenting the pores of the macropores can reduce the complexity of subsequent component division, improve the accuracy of division, and at the same time avoid including the macropore pore region in the calculation range of the matrix, thereby improving the feasibility of sub-resolution porosity estimation. The specific steps are as follows:

[0108] a. Define the pore size range of the macropores: Macropores are pores with a diameter higher than the image resolution. Affected by the resolution, CT images cannot identify all the pores of the sample mineral. Determine the spatial resolution of the CT device used to scan the sample mineral, and divide the true length of the sample mineral by the number of pixels included in this length to obtain the true scale value represented by each pixel;

[0109] b. Determine the relative volume content of the macropores: Detect the pore size distribution of the sample mineral through NMR or MIP, and calculate the pore content greater than the spatial resolution;

[0110] c. Mark the positions of macropores: Sort the linear attenuation coefficients of all slices according to their magnitudes. Take the linear attenuation coefficient corresponding to the relative volume content of macropores at the lower value end as the threshold. Those below this threshold are macropores.

[0111] 2) Correct the total porosity:

[0112] After obtaining the estimation result of the sub-resolution porosity of the matrix according to step (4), compare the estimation result with the actually measured total porosity:

[0113] If there is no difference in the comparison result, generate the final porosity distribution data and output a three-dimensional spatial model to improve the porosity estimation accuracy of the high-resolution CT image;

[0114] If there is a deviation in the comparison result, correct the estimation result based on the actual value to improve the accuracy. During the correction process, the segmentation of macropores and the matrix remains unchanged, where the porosity of macropores is always regarded as 100%. The correction is only for adjusting the porosity distribution of the matrix part to ensure the accuracy of the correction result. That is to say, during the correction process, the relative distribution characteristics of the sub-resolution porosity of the matrix part remain unchanged, that is, the correction is only for the overall numerical value of the matrix porosity and does not change the relative size relationship of the internal porosity. The specific steps are as follows:

[0115] a. Deduct the porosity of macropores from the total porosity of the sample mineral to determine the actually measured matrix porosity;

[0116] b. Use the product of the pixel volume of the sample mineral and the actually measured matrix porosity to calculate the actually measured matrix pore volume, and calculate the uncorrected matrix pore volume. The calculation method is shown in formula (4);

[0117] c. Calculate the proportionality coefficient between the actually measured matrix void volume and the uncorrected matrix pore volume, that is, the matrix porosity correction coefficient. The calculation method is shown in formula (5);

[0118] d. Apply the matrix porosity correction coefficient calculated in step c to the uncorrected matrix porosity to obtain the corrected matrix porosity distribution. The calculation method is shown in formula (6);

[0119]

[0120] where, V mtx is the uncorrected matrix pore volume; is the sub-resolution porosity spatial distribution matrix of the matrix; V mtx_c is the actually measured matrix pore volume; C is the matrix porosity correction coefficient; is the sub-resolution porosity spatial distribution matrix of the corrected matrix.

[0121] In order to further deepen the understanding of the present invention, the following specific implementation cases are given. The case data used in this implementation case are all from the open platform Digital Rocks Portal, which is an online open source platform focusing on porous media research and provides rich CT scanning data and related research materials. This implementation case mainly includes the following steps:

[0122] This implementation case selected CT data of one Belgian Fieldstone (naturally formed loose stone or large pebble) and two Bentheimer sandstone samples (Bentheim sandstone) for verification. The detailed information of the sample minerals is shown in Table 1. Figure 1 and Figure 2 In all the figures of , (a) is Belgian Fieldstone; (b) is Bentheimer_1; (c) is Bentheimer_2. The samples used in these datasets are widely studied, with complete mineral composition and pore information, and have verification value. The CT data format is RAW (.raw) or NetCDF (.nc), and the 3D voxel model is as follows Figure 3 shown.

[0123] Table 1 Sample information used in the present invention

[0124]

[0125] Component segmentation results: By counting the distribution range of the linear attenuation coefficients of all slices in the CT data, the slice with the widest distribution is selected as the representative slice. The representative slice numbers of the Belgian Fieldstone and two Bentheimer samples used in the present invention are 8,400 and 442, respectively. Figure 4 shown.

[0126] The linear attenuation coefficients of representative slices were subjected to k-means clustering analysis, and the clustering number K was tried as [2, 7] respectively. The component clustering results of the Belgian Fieldstone and two Bentheimer samples used are shown in the figure. Figure 5 According to the principle of selecting the highest number of clusters possible, while at the same time ensuring that a single mineral is not divided into multiple categories, the clustering results with larger k values are usually combined to distinguish the large grain gaps, and finally the number of clusters, the number of large grains, the number of matrix classes, and the number of macropores are determined. The results are shown in Table 2.

[0127] Table 2 Statistics of component segmentation results for each case

[0128]

[0129] The results of cluster analysis were applied to the entire 3D model, successfully segmenting out large particles, matrix, and macropores. After the initial segmentation of large particles, matrix, and macropores, the closed pore regions were detected and removed through connectivity analysis. In this implementation case, the D3Q19 model commonly used in the lattice Boltzmann method was selected to define connectivity, and the units connected to the boundary in the macropores and matrix were marked. The unmarked units were regarded as closed pores, mainly distributed inside the intact particles.

[0130] The process for determining the lower limit of the relative linear attenuation coefficient of the matrix is as follows: According to the clustering results with K = 3, the linear attenuation coefficient ranges and clustering center results of each component are shown in Table 3.

[0131] Table 3 Linear attenuation coefficient ranges and clustering centers of each component

[0132]

[0133] Among them, the upper limit of the relative linear attenuation coefficient of the macropores is μl _a , and the calculation process of the linear attenuation coefficient of the matrix part is as follows:

[0134] It is known that the mineral components of Belgian Fieldstone are quartz, glauconite, and illite clay minerals, and the mineral components of the Bentheimer sample are quartz, feldspar, and kaolinite clay minerals. Through the Mindat.org database, the chemical formulas, densities, etc. of quartz, feldspar, illite, and kaolinite are shown in Table 4.

[0135] Table 4 Densities and X-ray mass attenuation coefficients of minerals in the samples

[0136]

[0137] Through the NIST XCOM database, the variation trend of the mass attenuation coefficient of minerals with X-ray energy was queried, as Figure 6 shown (mass attenuation coefficients of common minerals in the X-ray energy range of 100 - 400 keV). In these samples, the content of quartz has an absolute advantage and the linear attenuation coefficient (or gray value) is stable. Therefore, quartz was selected as a reference to calculate the relative relationship of the linear attenuation coefficient with the matrix extraction. For example, when the X-ray energy of Bentheimer_1 is 150 keV, according to formula (2), the linear attenuation coefficient of the matrix crystal (kaolinite) is calculated to be 1.23 times that of the quartz crystal. According to the method of the present invention, taking the clustering center of the linear attenuation coefficient of the particles (μl _p ) as a reference, μ l_max = 1.23μ l_p . Similarly, μ of the matrix in other cases can be obtained, l_max specifically shown in Table 5.

[0138] Table 5 Upper limits of the relative linear attenuation coefficients of the matrix in each sample (μ l_max )

[0139]

[0140] Combining the upper and lower limits of the linear attenuation coefficient of the matrix, the porosity of the matrix part is calculated according to formula (3), and the following results are obtained: On the representative slice, the spatial distribution of the sub-resolution porosity is as Figure 7 shown, and the results show that there are obvious differences in porosity between the macropore and matrix regions.

[0141] The distribution of the porosity of the three-dimensional whole explains that the distribution of the porosity within the sample shows significant heterogeneity. The porosity in the macropore region is significantly higher than that in the matrix region and shows a certain fluctuation along the slice sequence, as Figure 8 shown, and the fluctuation of the average porosity may be related to the particle distribution of the sample minerals, the characteristics of the particle contact interface, and the local distribution of the closed pores.

[0142] The method for estimating the sub-resolution porosity of porous media based on monoenergetic CT provided by the present invention has the following advantages:

[0143] 1. Break through the assumption of the uniformity of the mass attenuation coefficient, and improve the accuracy and application range of porosity estimation;

[0144] 2. Integrate the elimination of closed pores and cross-scale characterization, and comprehensively consider the pores below the micron scale;

[0145] 3. Propose a multi-source data correction method to ensure the stability of the method.

[0146] 4. Applicable to high-energy CT and the case containing heavy metal elements, and solve the bottleneck of emerging technologies.

[0147] The above has described the embodiments of the present invention. The above description is exemplary, not exhaustive, and is not limited to the disclosed embodiments. Many modifications and variations are obvious to those of ordinary skill in the art in the technical field without departing from the scope and spirit of the described embodiments.

Claims

1. A method for estimating the sub-resolution porosity of a porous medium based on monoenergetic CT, characterized in that, It includes the following steps: Based on a known database, obtain the mass attenuation coefficient and density of the sample mineral, and calculate the linear attenuation coefficient, upper limit of relative linear attenuation coefficient, and lower limit of relative linear attenuation coefficient of the sample mineral components; Calculate the initial porosity distribution through the upper limit of relative linear attenuation coefficient and the lower limit of relative linear attenuation coefficient obtained from the above calculations; Generate porosity distribution data and output a three-dimensional space model to improve the porosity estimation accuracy of CT images.

2. The sub-resolution porosity estimation method according to claim 1, wherein Before obtaining the mass attenuation coefficient and density of the sample mineral, the following steps are also included: Search for representative slices of the sample mineral, and use the k-means algorithm to perform cluster analysis on the attenuation of the representative slices to obtain K ranges of linear attenuation coefficients and cluster centers.

3. The sub-resolution porosity estimation method according to claim 2, characterized in that The calculation of the linear attenuation coefficient of the sample mineral components includes the following steps: Divide the sample mineral into three groups, namely macropores, large particles, and matrix. Then, select the corresponding mass attenuation coefficient and density according to the X-ray energy in the database, and substitute the obtained mass attenuation coefficient and density into formula (1) to obtain the linear attenuation coefficient of each component of the sample mineral in the crystal state: μ l = μ m · ρ (1); where, μ l is the linear attenuation coefficient; μ m is the mass attenuation coefficient; ρ is the density, with the unit g / cm 3 .

4. The sub-resolution porosity estimation method according to claim 3, characterized in that The calculation of the lower limit of relative linear attenuation coefficient includes the following steps: Scan the air near the sample mineral, and the result obtained is the lower limit of relative linear attenuation coefficient; Preferably, the calculation of the upper limit of relative linear attenuation coefficient includes the following steps: Under the same X-ray energy, the ratio of the linear attenuation coefficients of large particles and matrix in the crystal state is constant. Combine formula (1) and formula (2) to obtain the upper limit of relative linear attenuation coefficient of the matrix: Among them, μ l_max is the upper limit of the relative linear attenuation coefficient of the matrix, with the unit of cm -1 ; μ m_m is the mass attenuation coefficient of the matrix, with the unit of cm 2 / g; μ m_p is the mass attenuation coefficient of the large particles, with the unit of cm 2 / g; ρ m is the density of the matrix mineral in the crystal state, with the unit of g / cm 3 ; ρ p is the density of the large particle mineral in the crystal state, with the unit of g / cm 3 ; μ l_p is the clustering center of the linear attenuation coefficient of the large particles, with the unit of cm -1 .

5. The sub-resolution porosity estimation method according to claim 4, wherein The calculation of the initial porosity distribution includes the following steps: Calculate the porosity through formula (3): Among them, is the porosity; μ l_j is the linear attenuation coefficient at position j, with the unit of cm -1 ; μ l_max is the linear attenuation coefficient corresponding to the matrix in the crystal state, with the unit of cm -1 ; μ l_a is the maximum value of the linear attenuation coefficient corresponding to the macropores, with the unit of cm -1 .

6. The sub-resolution porosity estimation method according to claim 1, characterized in that, Before outputting the three-dimensional space model, the following steps are also included: Analyze the pore size distribution information of the sample mineral and correct the macropore ratio; Correct the total porosity.

7. The sub-resolution porosity estimation method according to claim 6, wherein The correction of the macropore ratio includes the following steps: Define the pore size range of macropores; Determine the relative volume content of macropores, and calculate the pore content of the macropore spatial resolution, that is, the macropore porosity; Mark the positions of macropores.

8. The sub-resolution porosity estimation method according to claim 7, characterized in that The correction of the total porosity includes the following steps: Deduct the macropore porosity from the total porosity of the sample mineral to determine the measured matrix porosity; Use the product of the pixel volume of the sample mineral and the measured matrix porosity to calculate the measured matrix pore volume; Calculate the uncorrected matrix pore volume according to formula (4); Calculate the matrix porosity correction coefficient according to formula (5); Apply the matrix porosity correction coefficient to the uncorrected matrix porosity according to formula (6) for correction; Where: Among them, V mtx is the uncorrected matrix pore volume; is the sub-resolution porosity spatial distribution matrix of the matrix; C is the matrix porosity correction coefficient; V mtx_c is the measured matrix pore volume; is the sub-resolution porosity spatial distribution matrix of the corrected matrix.

9. The sub-resolution porosity estimation method according to claim 3, wherein When there are multiple mineral components in the matrix, calculate the equivalent mass attenuation coefficient of the matrix subclass by weighted average according to the volume fractions of different components, and the calculation method is shown in formula (7): Among them, μ m_c is the equivalent mass attenuation coefficient of the matrix subclass, with the unit of cm 2 / g; μ m_k is the mass attenuation coefficient of the k-th mineral under this matrix subclass, with the unit of cm 2 / g; r k is the proportion of the k-th mineral.

10. The sub-resolution porosity estimation method according to claim 3, characterized in that, When there are multiple mineral components in the matrix, calculate the crystal state density of the same element ratio of the matrix subclass according to the element ratio in the chemical formula of the components, and the calculation method is shown in formula (8): Among them, ρ c is the crystal state density of the same element ratio of the matrix subclass, with the unit of g / cm 3 ; ρ k is the density of the k-th mineral under this matrix subclass, with the unit of g / cm 3 ; S T refers to the total area of all minerals in the slice, with the unit of cm 2 ; S K refers to the area of the k-th mineral, with the unit of cm 2 .

Images