Rock CT image mineral identification method based on Latin hypercube random search

Abstract

The invention relates to a rock CT image mineral identification method based on Latin hypercube random search, and the method comprises the steps: carrying out the CT scanning of a to-be-identified rock, and obtaining a slice gray-scale map of the to-be-identified rock; according to the slice grey-scale map of the to-be-identified rock, obtaining a two-dimensional distribution map of each mineral in a target area of the to-be-identified rock; according to the two-dimensional distribution diagram of each mineral, obtaining the area proportion and the area proportion range of each mineral, and the dispersion index and the dispersion index range of each mineral; determining constraint conditions according to an area proportion and an area proportion range, and constructing a mathematical model based on a loss function of a dispersion index; the mathematical model is solved based on a Latin hypercube random search method, and the gray threshold segmentation interval of each mineral is obtained; and recognizing the mineral components of the rock to be recognized according to the gray threshold segmentation interval of each mineral. According to the invention, the mineral particles in the rock can be accurately identified.

Description

Technical Field

[0001] The invention relates to the technical field of CT image processing, and in particular to a rock CT image mineral recognition method based on Latin hypercube random search. Background Art

[0002] CT imaging technology has become a key tool in geological exploration, mineral resource development, geotechnical engineering, oil and gas exploration and development, and other fields. By generating three-dimensional digital core models, researchers can visually observe the internal structure of rocks. However, existing methods are still insufficient in mineral composition identification, especially in complex multi-mineral rock systems. Accurate mineral composition identification is crucial for mineral resource development, rock stability assessment, and rock mechanics analysis. It not only helps to accurately estimate mineral reserves and optimize mining strategies, but also provides basic data for rock mechanics models, revealing the role of different minerals in rock deformation, strength, and failure mechanisms, thereby effectively predicting geological disasters and guiding safety design.

[0003] Current CT image mineral recognition methods mostly rely on grayscale-based algorithms, such as K-means clustering and generative adversarial networks (GANs). Although these methods have certain advantages in processing speed, they have obvious limitations in accuracy and model authenticity. The main problems include:

[0004] (1) Randomness of cluster centers leads to inaccurate recognition: The K-means clustering algorithm is very sensitive to the selection of initial cluster centers during CT image segmentation. The randomness of the initial centers often leads to inaccurate clustering results, especially when the distribution of mineral particles is complex, and it is difficult to effectively distinguish the boundaries of different minerals. In addition, although the clustering algorithm can identify different types of minerals, its processing will modify the original image data, resulting in the loss of some image details, further affecting the accuracy of recognition.

[0005] (2) Lack of optimization algorithms and prominent local optimal problems: Although the K-means algorithm is easy to calculate, it is easy to fall into local optimality and cannot guarantee global optimal results. Existing CT image processing methods have not effectively introduced advanced optimization algorithms and cannot fully improve the recognition accuracy of complex mineral structures. At the same time, the simple processing of data structure by clustering methods limits its accurate recognition of complex mineral morphology and makes it difficult to show sufficient flexibility in heterogeneous rock systems.

[0006] (3) Falsification and complexity of high-resolution image processing: Although GAN-based technology can improve the resolution of low-resolution images and has made some progress in the characterization of mineral structure, clay morphology and pore structure, its model training process is complex, computationally expensive and time-consuming, making it difficult to adapt to application scenarios that require rapid processing and accurate identification of mineral components. In addition, the high-resolution images generated by GAN are not entirely derived from the original data, but are generated through model speculation, which may cause the generated detailed information to be inconsistent with the actual mineral structure, making it impossible to accurately identify the mineral category, especially in complex rock systems.

[0007] (4) Impact on image authenticity and insufficient generalization ability: These methods are highly dependent on training data, and the generalization ability of the models in different data sets and application scenarios is insufficient, making it difficult to cope with the recognition of a variety of complex mineral forms. In addition, although GAN and clustering algorithms have improved the ability to characterize mineral structures to a certain extent, their image processing often modifies the characteristics of the original data and cannot maintain a complete reflection of the mineral structure, thus limiting the practicality of these methods in actual mineral identification and classification. High computational costs and opaque model mechanisms further limit the performance of these methods in real-time applications.

[0008] In addition, the current dual-energy CT-based core mineral identification and porosity distribution joint characterization method also has technical defects:

[0009] (1) Limited mineral resolution and identification accuracy: The imaging quality of dual-energy CT is limited by the scanning resolution. Usually, tiny mineral particles in the core (especially submicron-sized particles) are difficult to accurately distinguish, which affects the precise identification of minerals.

[0010] (2) Difficulty in distinguishing minerals with similar density and atomic number: Dual-energy CT mainly relies on the difference in density and atomic number of the material to distinguish different minerals. However, for minerals with similar density and atomic number (such as some iron minerals and silicate minerals), it is difficult to effectively distinguish them through dual-energy CT, resulting in inaccurate results.

[0011] (3) Radiation dose issue: Dual-energy CT scanning involves a higher radiation dose, which may affect the health and safety of operators when used under laboratory conditions and also increase the cost and complexity of equipment use.

[0012] (4) Complexity of mineral composition reconstruction algorithms: Converting dual-energy CT image data into mineral composition information usually requires accurate physical models and algorithm support. Some existing algorithms are highly dependent on parameters such as mineral density and thickness, and are easily affected by errors, which limits the recognition accuracy.

[0013] In addition to the above-mentioned existing technologies, in the study of the mesoscopic mechanism of the variability of granite mechanical properties and the analysis of granite mineral composition and mesoscopic structure based on CT images, Mao Weize's method only deduces the three-dimensional spatial distribution information from the one-dimensional distribution information of the minerals, and does not have direct spatial resolution. In the construction of rock digital models considering the mesoscopic structure of minerals and their application, Zhang Yanbo's method is completely based on the color transformation of CT images, without considering the physical characteristics of each mineral, and completely relies on image processing methods to subjectively define mineral identification. It has the characteristics of strong subjectivity and lack of physical basis, so it is not suitable for the precise processing of CT minerals.

[0014] In order to solve the above problems, the present invention proposes a rock CT image mineral recognition method based on Latin hypercube random search. Summary of the invention

[0015] The purpose of the present invention is to provide a rock CT image mineral recognition method based on Latin hypercube random search. The CT image processing method based on Latin hypercube random search can not only identify mineral particles more accurately, but also simplify the training process and reduce the consumption of computing resources.

[0016] To achieve the above object, the present invention provides the following solutions:

[0017] The rock CT image mineral recognition method based on Latin hypercube random search includes:

[0018] Perform CT scanning on the rock to be identified to obtain a grayscale image of a slice of the rock to be identified;

[0019] According to the grayscale image of the slice of the rock to be identified, a two-dimensional distribution map of each mineral in the target area of ​​the rock to be identified is obtained;

[0020] According to the two-dimensional distribution map of each mineral, the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral are obtained;

[0021] Determining constraint conditions based on the area proportion and the area proportion range, and constructing a mathematical model based on the loss function of the dispersion index;

[0022] Solving the mathematical model based on the Latin hypercube random search method to obtain the gray threshold segmentation interval of each mineral;

[0023] The mineral composition of the rock to be identified is identified according to the grayscale threshold segmentation interval of each mineral.

[0024] Optionally, obtaining a two-dimensional distribution map of each mineral in the target rock area to be identified includes:

[0025] Slicing, grinding and polishing the target rock area to be identified to obtain a sample;

[0026] Performing X-ray diffraction analysis on the sample to obtain the main mineral components of the rock to be identified, and obtaining the elemental composition of the main mineral according to the chemical formula corresponding to the main mineral components;

[0027] Analyzing the sample using a scanning electron microscope to obtain a two-dimensional image of the surface of the rock to be identified, and obtaining a two-dimensional distribution map of each element according to the elemental composition of the main mineral components;

[0028] The two-dimensional distribution maps of each element are processed by set operations to obtain the two-dimensional distribution maps of each mineral.

[0029] Optionally, according to the two-dimensional distribution map of each mineral, obtaining the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral includes:

[0030] According to the ratio of the number of pixels of a single mineral to the number of pixels of all minerals, the area proportion of each mineral is obtained;

[0031] Calculate the number of connected regions of the i-th mineral in the two-dimensional distribution map of each mineral to obtain the number of connected regions;

[0032] Divide the two-dimensional distribution map of each mineral by using a number of grids equal to the number of connected regions, obtain the number of pixels of the i-th mineral in each grid, and calculate the expectation and variance of the number of pixels;

[0033] Obtaining a dispersion index of each mineral according to a ratio of the variance of the number of pixels to an expectation of the number of pixels;

[0034] The area proportion and fractional index of each mineral in the two-dimensional distribution map of each mineral are statistically analyzed to obtain the area proportion range and dispersion index range of each mineral.

[0035] Optionally, constructing a mathematical model includes:

[0036] Calculating the error value of the dispersion index range of each mineral according to the dispersion index of each mineral to obtain a loss function;

[0037] The loss function is:

[0038]

[0039]

[0040] Taking minimizing the loss function as the goal, the area proportion and the area proportion range determine the constraint conditions, and constructing the mathematical model;

[0041] The mathematical model is:

[0042]

[0043] The constraints are:

[0044]

[0045] Among them, L is the loss function, E i is the error value of the i-th mineral, n is the n-th mineral, a i , b i is the area proportion of the i-th mineral (a i , b i ), x i is the area proportion of the i-th mineral, is the dispersion index range of the ith mineral ρ i is the density of the ith mineral, and The density range of the i-th mineral, a n and b n is the area proportion range of n types of minerals.

[0046] Optionally, solving the mathematical model based on a Latin hypercube random search method includes:

[0047] S1. Obtain the slice grayscale image, the area ratio range of each mineral, the dispersion index range and the density range, and obtain a grayscale histogram based on the slice grayscale image;

[0048] S2, calculate the number of iterations;

[0049] S3, performing the kth iteration, sorting the minerals according to the area proportion range, obtaining the Latin square number of the i-th mineral, and performing Latin hypercube sampling on the area proportion interval of the i-th mineral according to the Latin square number;

[0050] S4, determine whether the random area proportions of all minerals are obtained, if so, calculate the area proportion of the nth mineral; if not, return to S3, wherein the area proportion of the nth mineral is 100% minus the sum of the area proportions of the previous n-1 minerals;

[0051] S5, judging whether the nth mineral is within the area ratio range, if not, returning to S3, if yes, judging whether the number of iterations reaches the calculated number of iterations;

[0052] S6. If the number of iterations does not reach the calculated number of iterations, the number of iterations is increased by one and the process returns to S3. If the number of iterations reaches the calculated number of iterations, the gray value range of each mineral and the new dispersion index are calculated, and the loss function is calculated according to the new dispersion index.

[0053] S7. Obtain the gray value range corresponding to each mineral in the optimal iterative process, and use the gray value range as the gray threshold segmentation interval of each mineral, wherein the optimal iterative process is the iterative process corresponding to minimizing the damage function.

[0054] Optionally, the method for calculating the number of iterations is:

[0055]

[0056]

[0057] Among them, T(1, n, a i , b i , m) is the total number of iterations, P i (x i , a i , b i ) is the number of ways to divide the i-th mineral, m is the number of minerals with crossover parts in mineral density, P k (x k , a i , b i ) is the number of divisions obtained by adding the grayscale frequency of the hth grayscale value when calculating the i-th mineral, 2 m Considering the possibility of arrangement brought by the logarithm m of the density overlap interval, P1(x1, a i , b i ) is the number of division methods obtained by adding the grayscale frequency of the first grayscale value of the ith mineral when calculating the ith mineral, and l is the ith mineral.

[0058] Optionally, calculate a new dispersion index for each mineral including:

[0059] S51, obtaining the area proportion, dispersion index range, grayscale histogram, and density range of each mineral in the kth iteration;

[0060] S52, obtaining the grayscale frequency of the j-th grayscale value, if the grayscale frequency of the j-th grayscale value is not less than the product of the total grayscale frequency and the area ratio, obtaining the grayscale value range of the i-th mineral, otherwise, the grayscale value is increased by one, and step S52 is repeated;

[0061] S53, converting the grayscale value range of the i-th mineral into a corresponding slice grayscale image, obtaining a grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining the number of connected areas, using the same grid as the number of connected areas to divide the grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining a new dispersion index;

[0062] S54. Determine whether the number of minerals reaches the value n. If so, output the result. If not, add one to i and one to j and return to S52.

[0063] Optionally, the gray value range of the i-th mineral obtained in step S52 includes:

[0064] The area proportions of each mineral in the kth iteration are sorted according to the density ranges of the minerals, corresponding to the order of grayscale values ​​from small to large in the grayscale histogram, and the grayscale frequency is allocated to the i-th mineral according to the area proportion to obtain the grayscale value range of the i-th mineral.

[0065] Optionally, the method further includes: determining the RGB color value of each mineral, converting the grayscale threshold segmentation interval of each mineral into a corresponding RGB color value, and generating a three-dimensional model.

[0066] The beneficial effects of the present invention are:

[0067] 1. Strong practicability: This method overcomes the reliance on complex models and high computing resources in existing CT image processing technology, adopts a simple and efficient optimization algorithm, reduces the requirements for training data and high-performance computing equipment, and significantly improves operability in practical applications. Whether in a laboratory environment or on a mining site, the present invention can be quickly deployed and achieve accurate mineral identification, and has wide applicability.

[0068] 2. Efficiency: The present invention significantly improves the processing speed through the uniform sampling strategy of Latin hypercube random search, avoiding the complex training process in the traditional generative adversarial network (GAN) method. Compared with the prior art, the present invention reduces the computational cost and time consumption, allowing mineral identification and three-dimensional reconstruction to be completed in a shorter time, meeting the application scenarios that require efficient data processing.

[0069] 3. Authenticity guarantee: When generating high-resolution three-dimensional images of mineral particles, the present invention avoids the problem of model speculation generating fictitious details. Through the optimized sampling strategy, the method can ensure the true reflection of the three-dimensional reconstruction results, retain the actual physical properties of the mineral particles, enhance the consistency between the three-dimensional image and the actual rock sample, and provide a higher authenticity guarantee for mineral composition identification and rock structure analysis.

[0070] 4. Accuracy of identification results: Through Latin hypercube random search technology, the present invention effectively overcomes the limitations of the traditional K-means clustering algorithm in dealing with complex mineral systems, and significantly improves the recognition accuracy of mineral grain boundaries, especially between minerals with similar grayscale values. This method can accurately distinguish different mineral components, ensure the accuracy of mineral identification and three-dimensional reconstruction, and performs well in complex multi-mineral rock systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0072] Figure 1 The XRD analysis results of the granite samples of the embodiments of the present invention are as follows;

[0073] Figure 2 Part of the two-dimensional mineral distribution map obtained by scanning electron microscope (SEM) analysis of the granite sample of the embodiment of the present invention, wherein (a) is the first two-dimensional mineral distribution map, (b) is the second two-dimensional mineral distribution map, (c) is the third two-dimensional mineral distribution map, and (s) is the fourth two-dimensional mineral distribution map;

[0074] Figure 3 A grayscale histogram obtained by CT scanning analysis of a granite sample according to an embodiment of the present invention;

[0075] Figure 4 The grayscale image of the granite sample of the embodiment of the present invention is compared with the RGB image after conversion, wherein (a) is the grayscale image before conversion, and (b) is the RGB image after conversion;

[0076] Figure 5 3D reconstructed model of a granite sample according to an embodiment of the present invention, wherein (a) is all minerals, (b) is quartz, (c) is plagioclase, (d) is orthoclase, (e) is muscovite, and (f) is biotite;

[0077] Figure 6 This is a flow chart of a method for identifying minerals in rock CT images based on Latin hypercube random search according to an embodiment of the present invention;

[0078] Figure 7 The algorithm flow chart of solving the gray threshold segmentation range of each mineral according to the embodiment of the present invention;

[0079] Figure 8 The figure is a flow chart of an algorithm for calculating a dispersion index according to an embodiment of the present invention. DETAILED DESCRIPTION

[0080] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0081] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0082] This embodiment provides a rock CT image mineral identification method based on Latin hypercube random search, including:

[0083] Perform CT scanning on the rock to be identified to obtain a grayscale image of the slice of the rock to be identified;

[0084] According to the grayscale image of the rock slice to be identified, a two-dimensional distribution map of each mineral in the target area of ​​the rock to be identified is obtained;

[0085] According to the two-dimensional distribution map of each mineral, the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral are obtained;

[0086] The constraints are determined by area ratio and area ratio range, and a mathematical model is constructed based on the loss function of the dispersion index;

[0087] The mathematical model is solved based on the Latin hypercube random search method to obtain the gray threshold segmentation interval of each mineral;

[0088] The mineral composition of the rock to be identified is identified by segmenting the intervals according to the grayscale threshold of each mineral.

[0089] Further, obtaining a two-dimensional distribution map of each mineral in the target rock area to be identified includes:

[0090] Slice, grind and polish the target area of ​​the rock to be identified to obtain samples;

[0091] Perform X-ray diffraction analysis on the sample to obtain the main mineral components of the rock to be identified, and obtain the elemental composition of the main minerals based on the chemical formula corresponding to the main mineral components;

[0092] The samples were analyzed using a scanning electron microscope to obtain a two-dimensional image of the rock surface to be identified, and a two-dimensional distribution map of each element was obtained based on the elemental composition of the main mineral components;

[0093] The two-dimensional distribution map of each element is processed by set operation to obtain the two-dimensional distribution map of each mineral.

[0094] Specifically, in this embodiment, the core is scanned using a CT scanning device to obtain characteristics such as rock pore distribution, crack distribution, and large-particle mineral distribution, and then the same specific area of ​​the core, i.e., the target area, is sliced. After grinding and polishing, the slice is subjected to phase analysis using X-ray diffraction (XRD), and EDS Mapping element analysis is performed using a scanning electron microscope (SEM). The grayscale histogram of the CT scan grayscale image of the specific area is determined, the XRD data is processed to obtain the main mineral component types, and the density range of each mineral component is obtained by querying the standard database. The steps include: according to the grayscale image of the corresponding slice position in the CT scan, OpenCV is used to perform grayscale frequency statistics of each grayscale value to obtain a grayscale histogram, the corresponding XRD data is analyzed by least squares fitting through the standard mineral database PDF card to obtain n types of minerals, and the density range corresponding to the i-th mineral is obtained by querying the standard database.

[0095] Furthermore, according to the two-dimensional distribution map of each mineral, the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral are obtained, including:

[0096] According to the ratio of the number of single mineral pixels to the number of all mineral pixels, the area proportion of each mineral is obtained;

[0097] Calculate the number of connected regions of the i-th mineral in each mineral two-dimensional distribution map to obtain the number of connected regions;

[0098] The two-dimensional distribution map of each mineral is divided into grids with the same number as the number of connected regions, the number of pixels of the i-th mineral in each grid is obtained, and the expectation and variance of the number of pixels are calculated;

[0099] The dispersion index of each mineral is obtained according to the ratio of the variance of the number of pixels to the expected number of pixels;

[0100] The area proportion and fractional index of each mineral in the two-dimensional distribution map of each mineral are statistically analyzed to obtain the area proportion range and dispersion index range of each mineral.

[0101] Specifically, in this embodiment, the chemical formula corresponding to the mineral is obtained from the main mineral components of the XRD analysis results, and then its elemental composition is inferred based on the chemical formula. The two-dimensional distribution map of each element obtained by scanning electron microscopy is subjected to set operations to obtain the two-dimensional distribution map of each mineral. Corresponding to the two-dimensional distribution map of the mineral, OpenCV is used to count the number of pixels of all minerals, and then the number of pixels of a single mineral is counted to obtain the area proportion of each mineral. After statistical analysis of multiple mineral two-dimensional distribution maps, the area proportion range of the i-th mineral (a i , b i), in the order of increasing maximum value of area proportion range, define the area proportion of the i-th mineral as x i , the area proportion of the nth mineral Then the dispersion index d of each mineral is calculated. i , the method to calculate the dispersion index is:

[0102] First, the number of connected areas N of the i-th mineral in the two-dimensional distribution map of minerals is counted, and then the two-dimensional distribution map of minerals is evenly divided into N equal square grids. Then, the number of pixels occupied by the mineral in each square grid is counted, and after calculating the expectation and variance of the number of pixels, the dispersion index of the i-th mineral is obtained using the following formula:

[0103]

[0104] Using this method to count multiple mineral two-dimensional distribution maps, the dispersion index range of the i-th mineral is obtained:

[0105] Furthermore, constructing a mathematical model includes:

[0106] The error value of the dispersion index range of each mineral is calculated according to the dispersion index of each mineral to obtain the loss function;

[0107] The loss function is:

[0108]

[0109] Taking minimizing the loss function as the goal, the area ratio and the area ratio range as constraints, and constructing a mathematical model;

[0110] The mathematical model is:

[0111]

[0112] Among them, L is the loss function, E i is the error value of the i-th mineral, n is the n-th mineral, a i , b i is the area proportion of the i-th mineral (a i , b i ), x i is the area proportion of the i-th mineral, is the dispersion index range of the ith mineral ρ i is the density of the ith mineral, and The density range of the i-th mineral, a n and b n It is the area proportion range of the nth mineral.

[0113] Furthermore, solving the mathematical model based on the Latin hypercube random search method includes:

[0114] S1. Obtain a slice grayscale image, an area ratio range of each mineral, a dispersion index range, and a density range, and obtain a grayscale histogram based on the slice grayscale image;

[0115] S2, calculate the number of iterations;

[0116] S3, performing the kth iteration, obtaining the Latin square number of the i-th mineral, and performing Latin hypercube sampling on the area proportion interval of the i-th mineral according to the Latin square number;

[0117] S4, determine whether the random area proportions of all minerals are obtained, if so, calculate the area proportion of the nth mineral; if not, return to S3, wherein the area proportion of the nth mineral is 100% minus the sum of the area proportions of the previous n-1 minerals;

[0118] S5, judging whether the nth mineral is within the area ratio range, if not, returning to S3, if yes, judging whether the number of iterations reaches the calculated number of iterations;

[0119] S6. If the number of iterations does not reach the calculated number of iterations, the number of iterations is increased by one and the process returns to S3. If the number of iterations reaches the calculated number of iterations, the gray value range of each mineral and the new dispersion index are calculated, and the loss function is calculated according to the new dispersion index.

[0120] S7. Obtain the gray value range corresponding to each mineral in the optimal iteration process, and use the gray value range as the gray threshold segmentation interval of each mineral.

[0121] Furthermore, new dispersion indices are calculated for each mineral including:

[0122] S51, obtaining the area proportion, dispersion index range, grayscale histogram, and density range of each mineral in the kth iteration;

[0123] S52, obtaining the grayscale frequency of the j-th grayscale value, if the grayscale frequency of the j-th grayscale value is not less than the product of the total grayscale frequency and the area ratio, obtaining the grayscale value range of the i-th mineral, otherwise, the grayscale value is increased by one, and step S52 is repeated;

[0124] S53, converting the grayscale value range of the i-th mineral into a corresponding slice grayscale image, obtaining a grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining the number of connected areas, using the same grid as the number of connected areas to divide the grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining a new dispersion index;

[0125] S54. Determine whether the number of minerals reaches the value n. If so, output the result. If not, add one to i and one to j and return to S52.

[0126] The gray value range of the i-th mineral obtained in step S52 includes:

[0127] The area proportions of each mineral in the kth iteration are sorted according to the density range of each mineral, corresponding to the order of grayscale values ​​from small to large in the grayscale histogram, and the grayscale frequency is assigned to the i-th mineral according to the area proportion to obtain the grayscale value range of the i-th mineral.

[0128] Furthermore, the method also includes: determining the RGB color value of each mineral, converting the grayscale threshold segmentation interval of each mineral into a corresponding RGB color value, and generating a three-dimensional model.

[0129] Taking the scanning electron microscope results, XRD analysis and computer tomography (CT) results of granite as an example, the XRD analysis results of granite samples are as follows: Figure 1 Specifically, granite is an igneous rock, and due to the high temperature involved in its formation process, there are a variety of mineral crystals. In order to further study the spatial distribution of these minerals, as shown in FIG. Figure 6 As shown, the method proposed by the present invention is implemented by the following steps:

[0130] 1. CT scanning: First, perform CT scanning on the granite sample to obtain its three-dimensional structural information. The multiple slice grayscale images obtained from the CT scanning results can provide the density distribution inside the sample, providing basic data for subsequent mineral differentiation and three-dimensional reconstruction.

[0131] 2. Scanning electron microscopy and XRD: Based on CT scanning, select representative areas of granite samples, and after slicing, grinding, polishing and other processes, prepare thin slice samples for scanning electron microscopy (SEM) and X-ray diffraction analysis (XRD). SEM analysis can obtain high-resolution two-dimensional images of the granite surface, and combined with EDS Mapping to provide two-dimensional spatial distribution information of each chemical element, and then obtain a two-dimensional distribution map of each chemical element; and through XRD analysis, the mineral components contained can be obtained.

[0132] 3. Analysis data: Through XRD analysis, it can be concluded that the main components of granite include quartz and feldspar, followed by mica. By querying the standard database, the corresponding chemical formula and density are shown in Table 1.

[0133] Table 1

[0134] mineral Chemical formula density quartz <![CDATA[SiO2]]> 2.65-2.66 Plagioclase <![CDATA[NaAlSi3O8、CaAl2Si2O8]]> 2.61-2.76 Orthoclase <![CDATA[KAlSi3O8]]> 2.55-2.63 muscovite <![CDATA[KAl2(AlSi3O 10 )(OH)]]> 2.76-3 Biotite <![CDATA[K(Mg,Fe)3AlSi3O 10 (F,OH)2]]> 2.8-3.4

[0135] The two-dimensional distribution map of each element is obtained by scanning electron microscope EDS Mapping analysis. After the two-dimensional distribution map of each element is processed by set operation, the two-dimensional distribution map of each mineral is obtained, such as Figure 2 The corresponding mineral set operations are shown in Table 2:

[0136] Table 2

[0137]

[0138] 3. Calculation of area proportion and dispersion index in the two-dimensional distribution map of minerals: By processing the two-dimensional distribution map of each mineral, the area proportion and dispersion index of different minerals in the selected area are calculated. The area proportion indicates the relative proportion of the mineral in the sample, while the dispersion index is used to describe the distribution uniformity of mineral particles in two-dimensional space. The first four sets of data corresponding to the calculation of the granite sample are shown in Table 3:

[0139] Table 3

[0140]

[0141] 4. Establish a mathematical model: obtain the grayscale image of a specific position of the CT scan sample and generate a grayscale histogram, such as Figure 3 As shown in the figure, the minimization loss function is defined as the optimization algorithm loss function, and its calculation method is as follows: First, the dispersion index range of each mineral to its SEM results is calculated The error value E i , which is calculated as:

[0142]

[0143] Then calculate the loss function:

[0144]

[0145] Corresponding to the objective function and constraints, the following mathematical model is established:

[0146]

[0147]

[0148] 5. Optimization algorithm processing and optimal solution determination: This embodiment uses an optimization algorithm to process the grayscale image and mineral density combination, such as Figure 7In this example, the calculated number of iterations is 11075328. The number of iterations is used as the Latin square number, so that each mineral is sampled in Latin hypercube within the area ratio range. Then, the minerals are sorted in the order of density from small to large, corresponding to the order of grayscale values ​​from small to large in the grayscale histogram. The grayscale frequency is allocated to each mineral according to the area ratio. After the allocation is completed, it is converted into a grayscale value range. The GPU is called to parallelly calculate the dispersion index of each iteration process, such as Figure 8 When calculating, the gray value range is converted into the corresponding gray map, the number of connected areas N in the gray value range is counted, and then the dispersion index is obtained in the same way as the dispersion index of the two-dimensional distribution map of minerals. While obtaining the dispersion index of each iteration process, the GPU parallel calculation is called again to obtain the loss function, and the optimal iteration process is selected, and the gray value range corresponding to each mineral is used as the final gray threshold segmentation result.

[0149] like Figure 7 As shown, the method for solving the mathematical model based on the Latin hypercube random search method includes:

[0150] S1. Read a CT slice grayscale image, the area ratio range and dispersion index range of each mineral calculated from the two-dimensional EDS Mapping image, and the density range of each mineral;

[0151] S2. Obtaining a grayscale histogram by counting the frequency of each grayscale value of the CT slice grayscale image;

[0152] S3. Sort the minerals according to the upper limit of their area proportion, where there are i types of minerals: i=1, 2, 3, ..., n, where the larger the i, the larger the upper limit of the area proportion of the mineral;

[0153] S4. Using the iteration number calculation method to solve the total number of threshold segmentation division methods, the number of division methods is the iteration number, and the number of division methods of the i-th mineral is obtained, that is, the number of Latin squares of the i-th mineral, and Latin hypercube sampling is performed on the area proportion interval of the i-th mineral according to the Latin square number of the i-th mineral to obtain multiple random area proportions of each interval;

[0154] S5. Let the integer k be a counter of the number of iterations, and add one to k for each iteration until it reaches the calculated number of iterations K;

[0155] S6. Start the kth iteration process, let the integer i be the counter, and add one after each calculation of a mineral, and return to S4 until it reaches the value of the number of minerals n.

[0156] S7. Process the ith mineral, select its area share from the values ​​obtained in Latin hypercube sampling according to the combination arrangement method, and obtain multiple random area shares of each mineral. The area share of each mineral can vary within these multiple random area shares. If the ith mineral has mi random area shares, then the total number of combination arrangements of the first n-1 is Πmi.

[0157] S8. The area proportion of the nth mineral is calculated by subtracting the sum of the area proportions of the previous n-1 minerals from 100%. If the area proportion of the nth mineral does not fall within the range, the iteration process is invalid, the k value is not counted, and the process returns to S4 for re-iteration, wherein the sum of the area proportions of the previous n-1 minerals is one of the total combination arrangements;

[0158] When the S9.k value reaches the number of iterations, the GPU is called to parallelly calculate the dispersion index of each iteration process;

[0159] S10. Parallel calculation of the loss function of the dispersion index for the dispersion index range in each iteration process;

[0160] S11. Use the quick sort algorithm to find the optimal iterative process, which is the segmentation result obtained, wherein the optimal iterative process is the iterative process corresponding to the minimized damage function.

[0161] like Figure 8 As shown, the method for calculating the new dispersion index includes:

[0162] S1. First, read a CT grayscale slice image and the following data: the area ratio combination of the kth iteration process, the dispersion index range of each mineral in the two-dimensional distribution map, the grayscale histogram, and the density ratio range of each mineral;

[0163] S2. Sort the minerals according to the lower limit of the density range, and there are i types of corresponding minerals: i = 1, 2, 3, ..., n, where the larger the i, the larger the lower limit of the mineral density range;

[0164] S3. Let integer i be a counter, and after each calculation of a mineral, i is incremented by one until it reaches the value of the number of minerals n;

[0165] S4. Start assigning a grayscale value range to the i-th mineral, let integer j be a grayscale value counter, and add one after each grayscale value is calculated until the ratio of grayscale frequency to total grayscale frequency reaches the area ratio; where grayscale frequency = grayscale frequency + grayscale frequency of j-th grayscale value, where "=" is an assignment operation symbol;

[0166] Among them, assigning a gray value range to each mineral includes: sorting the area proportion of each mineral in the kth iteration according to the density of each mineral from small to large, corresponding to the gray value in the gray histogram from small to large, assigning the gray frequency to each mineral in turn according to the area proportion according to the mineral density order, and converting it into a gray value range after the assignment is completed. The gray frequency assignment is the assignment of the ratio of the sum of the selected gray frequencies to the sum of the gray frequencies of all minerals closest to the area proportion of the mineral;

[0167] S5. after obtaining the gray value range of the i-th mineral, converting the CT slice gray image into a gray image that only retains the gray value range of the i-th mineral;

[0168] S6. Start counting the number N of connected regions on the grayscale image, divide the grayscale image evenly into N square grids of equal size, and count the number of pixels in each square grid;

[0169] S7. Calculate the expectation and variance of the number of pixels in all grids, and obtain the dispersion index of the ith mineral, which is the ratio of the variance to the expectation, that is, calculate the dispersion index of all minerals.

[0170] S8. Determine whether the number of minerals reaches the value n. If so, output the result. If not, add 1 to i and 1 to j and return to S4.

[0171] The number of iterations calculated is used as the Latin square number. The number of iterations is calculated as follows:

[0172] Among them, P n (1, a i , b i ) is to divide the grayscale histogram into n intervals, and the area ratio of each interval satisfies [a i , b i ] is the number of ways to divide . If Number of permutations with alternating order 2 m It takes into account the possible permutations brought about by the logarithm m of the density overlapping interval.

[0173]

[0174]

[0175] Among them, T(1, n, a i , b i , m) is to divide the grayscale histogram into n intervals, and the area ratio of each interval satisfies [a i , b i ] is divided into the number of ways, that is, the number of iterations, P i (x i , a i , b i) is the number of ways to divide the i-th mineral, representing the number of Latin squares of the i-th mineral, m is the number of minerals with crossover parts in the mineral density, P k (x k , a i , b i ) is the number of divisions obtained by adding the grayscale frequency of the hth grayscale value when calculating the i-th mineral, 2 m Considering the possibility of arrangement brought by the logarithm m of the density overlap interval, P1(x1, a i , b i ) is the number of division methods obtained by adding the grayscale frequency of the first grayscale value of the ith mineral when calculating the ith mineral, and l is the ith mineral.

[0176] 6. 3D reconstruction and visualization: Based on the optimal solution, define the RGB color value for each mineral, map the grayscale image of the CT scan with the corresponding RGB color value, and finally generate a 3D digital model, such as Figure 4-5 In this 3D model, different mineral types are intuitively displayed in different colors, which enables accurate distinction and 3D reconstruction of mineral particles in the digital core. Table 4 shows the RGB values ​​corresponding to the solution results.

[0177] Table 4

[0178] mineral RGB Value Threshold segmentation range quartz 0，255，0 100-151 Plagioclase 0，0，255 80-86，93-99 Orthoclase 255，0，255 50-79，87-92 muscovite 255，255，0 152-167 Biotite 255，0，0 168-255

[0179] 7. Result accuracy verification mechanism: First, the volume fraction and three-dimensional dispersion index of each mineral in the converted three-dimensional digital core are counted. When calculating the volume fraction, the number of stereo pixels of each mineral is counted. When calculating the three-dimensional dispersion index, the number of connected three-dimensional connected areas is counted. Then, cubic grids of equal size are divided, and the expectation and variance of the number of stereo pixels in each cubic grid are counted. The rest of the methods are the same as those for two-dimensional.

[0180] The Kolmogorov-Smirnov test was used to judge the accuracy of the transformation results, and the D statistic was defined as follows:

[0181]

[0182] where F 1，n and F 2，n are the empirical distribution functions of the solution result samples and the two-dimensional distribution samples of each mineral, respectively, and sup is the supremum function. Among them, the empirical distribution function depends on the parameters to be tested. When testing the volume fraction of the solution result, the area ratio of the two-dimensional distribution of each mineral is the empirical distribution function, and when testing the dispersion index of the solution result, the dispersion index of the two-dimensional distribution of each mineral is the empirical distribution function.

[0183] The null hypothesis is rejected at α if the D statistic satisfies the following conditions:

[0184]

[0185] Among them, α is the significance level, which is 0.005 in this example, n and m are the sample sizes of the solution result sample and the two-dimensional distribution samples of each mineral, respectively. In this example, n=m=5, and c(α) is the D statistic corresponding to the significance level.

[0186] The calculation for c(α) is:

[0187]

[0188] The conditions for rejecting the hypothesis are:

[0189]

[0190] If the hypothesis is rejected, it means that the result is unreasonable, otherwise it is reasonable. This allows us to scientifically judge the rationality of the result.

[0191] The test results of the corresponding volume fractions are shown in Table 5:

[0192] Table 5

[0193]

[0194] The test results of the corresponding dispersion index are shown in Table 6:

[0195] Table 6

[0196]

[0197]

[0198] From the above test results, it can be seen that the solution is reasonable.

[0199] Through the above implementation steps, the present invention can not only accurately display the mineral distribution in granite in three-dimensional space, but also provide an important analytical tool for geological research and mineral exploration. This method is particularly suitable for core samples with complex components and diverse mineral combinations, and can significantly improve the accuracy and efficiency of research.

[0200] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A method for identifying minerals in rock CT images based on Latin hypercube random search, characterized in that: include: Perform CT scanning on the rock to be identified to obtain a grayscale image of a slice of the rock to be identified; According to the grayscale image of the slice of the rock to be identified, a two-dimensional distribution map of each mineral in the target area of ​​the rock to be identified is obtained; According to the two-dimensional distribution map of each mineral, the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral are obtained; Determining constraint conditions based on the area proportion and the area proportion range, and constructing a mathematical model based on the loss function of the dispersion index; Solving the mathematical model based on the Latin hypercube random search method to obtain the gray threshold segmentation interval of each mineral; The mineral composition of the rock to be identified is identified according to the grayscale threshold segmentation interval of each mineral.

2. The rock CT image mineral identification method based on Latin hypercube random search according to claim 1 is characterized in that: Obtaining the two-dimensional distribution map of each mineral in the target rock area to be identified includes: Slicing, grinding and polishing the target rock area to be identified to obtain a sample; Performing X-ray diffraction analysis on the sample to obtain the main mineral components of the rock to be identified, and obtaining the elemental composition of the main mineral according to the chemical formula corresponding to the main mineral components; Analyzing the sample using a scanning electron microscope to obtain a two-dimensional image of the surface of the rock to be identified, and obtaining a two-dimensional distribution map of each element based on the elemental composition of the main minerals; The two-dimensional distribution maps of each element are processed by set operations to obtain the two-dimensional distribution maps of each mineral.

3. The rock CT image mineral identification method based on Latin hypercube random search according to claim 2 is characterized in that: According to the two-dimensional distribution map of each mineral, the area proportion and area proportion range of each mineral, the dispersion index and dispersion index range of each mineral are obtained, including: According to the ratio of the number of pixels of a single mineral to the number of pixels of all minerals, the area proportion of each mineral is obtained; Calculate the number of connected regions of the i-th mineral in the two-dimensional distribution map of each mineral to obtain the number of connected regions; Divide the two-dimensional distribution map of each mineral by using a number of grids equal to the number of connected regions, obtain the number of pixels of the i-th mineral in each grid, and calculate the expectation and variance of the number of pixels; Obtaining a dispersion index of each mineral according to a ratio of the variance of the number of pixels to an expectation of the number of pixels; The area proportion and fractional index of each mineral in the two-dimensional distribution map of each mineral are statistically analyzed to obtain the area proportion range and dispersion index range of each mineral.

4. The rock CT image mineral identification method based on Latin hypercube random search according to claim 1 is characterized in that: Constructing mathematical models includes: Calculating the error value of the dispersion index range of each mineral according to the dispersion index of each mineral to obtain a loss function; The loss function is: Taking minimizing the loss function as the goal, the area proportion and the area proportion range determine the constraint conditions, and constructing the mathematical model; The mathematical model is: The constraints are: Among them, L is the loss function, E i is the error value of the i-th mineral, n is the n-th mineral, a i , b i is the area proportion of the i-th mineral (a i , b i ), x i is the area proportion of the i-th mineral, is the dispersion index range of the ith mineral ρ i is the density of the ith mineral, and The density range of the i-th mineral, a n and b n is the area proportion of n minerals, d i Dispersion index of the ith mineral.

5. The rock CT image mineral identification method based on Latin hypercube random search according to claim 4 is characterized in that: Solving the mathematical model based on the Latin hypercube random search method includes: S1. Obtain the slice grayscale image, the area ratio range of each mineral, the dispersion index range and the density range, and obtain a grayscale histogram based on the slice grayscale image; S2, calculate the number of iterations; S3, performing a kth iteration, sorting the minerals according to the area proportion range, obtaining the Latin square number of the i-th mineral, performing Latin hypercube sampling on the area proportion interval of the i-th mineral according to the Latin square number, and obtaining the random area proportion of the i-th mineral; S4, determine whether the random area proportions of all minerals are obtained, if so, calculate the area proportion of the nth mineral; if not, return to S3, wherein the area proportion of the nth mineral is 100% minus the sum of the area proportions of the previous n-1 minerals; S5, judging whether the nth mineral is within the area ratio range, if not, returning to S3, if yes, judging whether the number of iterations reaches the calculated number of iterations; S6. If the number of iterations does not reach the calculated number of iterations, the number of iterations is increased by one and the process returns to S3. If the number of iterations reaches the calculated number of iterations, the gray value range of each mineral and the new dispersion index are calculated, and the loss function is calculated according to the new dispersion index. S7. Obtain the gray value range corresponding to each mineral in the optimal iterative process, and use the gray value range as the gray threshold segmentation interval of each mineral, wherein the optimal iterative process is the iterative process corresponding to minimizing the damage function.

6. The rock CT image mineral identification method based on Latin hypercube random search according to claim 5 is characterized in that: The method for calculating the number of iterations is: Among them, T(1, n, a i , b i , m) is the number of iterations, P i (x i , a i , b i ) is the number of ways to divide the i-th mineral, m is the number of minerals with crossover parts in mineral density, P k (x k , a i , b i ) is the number of divisions obtained by adding the grayscale frequency of the hth grayscale value when calculating the i-th mineral, 2 m Considering the possibility of arrangement brought by the logarithm m of the density overlap interval, P1(x1, a i , b i ) is the number of division methods obtained by adding the grayscale frequency of the first grayscale value of the ith mineral when calculating the ith mineral, and l is the lth mineral.

7. The rock CT image mineral identification method based on Latin hypercube random search according to claim 5 is characterized in that: Calculation of new dispersion indices for each mineral includes: S51, obtaining the area proportion, dispersion index range, grayscale histogram, and density range of each mineral in the kth iteration; S52, obtaining the grayscale frequency of the j-th grayscale value, if the grayscale frequency of the j-th grayscale value is not less than the product of the total grayscale frequency and the area ratio, obtaining the grayscale value range of the i-th mineral, otherwise, the grayscale value is increased by one, and step S52 is repeated; S53, converting the grayscale value range of the i-th mineral into a corresponding slice grayscale image, obtaining a grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining the number of connected areas, using the same grid as the number of connected areas to divide the grayscale image that only retains the grayscale value range of the i-th mineral, and obtaining a new dispersion index; S54. Determine whether the number of minerals reaches the value n. If so, output the result. If not, add one to i and one to j and return to S52.

8. The method for identifying minerals in rock CT images based on Latin hypercube random search according to claim 7, characterized in that: The gray value range of the i-th mineral obtained in step S52 includes: The area proportions of each mineral in the kth iteration are sorted according to the density ranges of the minerals, corresponding to the order of grayscale values ​​from small to large in the grayscale histogram, and the grayscale frequency is allocated to the i-th mineral according to the area proportion to obtain the grayscale value range of the i-th mineral.

9. The method for identifying minerals in rock CT images based on Latin hypercube random search according to claim 1, characterized in that: The method further includes: determining the RGB color value of each mineral, converting the grayscale threshold segmentation interval of each mineral into a corresponding RGB color value, and generating a three-dimensional model.

Images