METHOD OF SEPARATING, IDENTIFYING AND CHARACTERIZING CRACKS IN 3D SPACE

Abstract

The present invention discloses a method of separating, identifying and characterizing cracks in 3D space, which processes as follows to a volumetric image, so as to perform the separation, identification and the characterization of the cracks in the 3D space: 1) preprocessing digital image; 2) statistically analyzing basic information of the digital image: the basic information of the image includes porosity, connectivity of each pore, statistics of pore size, and position, size, orientation and anisotropy of each pore-structure; 3) filtration: removing non-crack structure in the image; 4) smoothening: smoothening and mending the image; 5) thinning: thinning the void structure into a thickness d (d can be any value, but more appropriate to be 2 to 3 voxels generally) in a direction with shortest extension in the 3D space; 6) separation: separating intersected cracks in a crack network by breaking the connections; 7) combination: combining those elongated cracks that are disconnected in the last step, merging tiny structures that are formed during the separation to a nearby large cluster, and restoring cracks to the thickness before thinning, and eventually giving out the characterization of the cracks. In the following expression, the wording “void” is used more, emphasizing the “empty” gap in the image rather than the rock solid. In this patent application, it is mainly for the case where the void appears in a state of crack, not excluding the case where the void appears in a state of small pore.

Claims

1. A method of separating, identifying and characterizing cracks in 3D space for a volumetric image, the method comprising: 1) preprocessing digital image; 2) statistically analyzing basic information of the digital image: the basic information of the image includes porosity, connectivity, and position, size, orientation and anisotropy of each void structure; 3) filtration: removing non-crack structure in the digital image; 4) smoothening: smoothening and mending the digital image; 5) thinning: thinning the void structure into a thickness d (2 to 3 voxels) in a most narrow direction in three-dimensions; 6) separation: separating intersected cracks in a crack network by breaking the intersections; 7) combination: combining the void structures that are disconnected in the last step, merging tiny structures that are formed during the separation to nearby ones, and restoring the thickness of the crack; on this basis, analyzing porosity, and position, size, orientation and anisotropy of each void structure, finishing the identification and characterization of the cracks in the 3D space.

2. The method of separating, identifying and characterizing cracks in 3D space according to claim 1, wherein a compute mode of the anisotropy of the structure mentioned in step 2) is as follows: for a void structure containing n voxels, each voxel forms a vector from a center point of a cluster to a current position i presented as a.sub.i=(a.sub.xi, a.sub.yi, a.sub.zi).sup.T, and coordinates of a center of the void structure is calculated by a mean value of all point coordinates in the void structure; then an anisotropy of this void structure can be represented by an orientation matrix R: T=Σ.sub.i=1.sup.n a.sub.ia.sub.i.sup.T= $\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}^2& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{zi}}\\ \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}^2& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}.Math.a_{\mathrm{zi}}\\ \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{zi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}.Math.a_{\mathrm{zi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{zi}}^2\end{array}\right]\hspace{1em}$ the matrix having three eigenvalues τ.sub.1<τ.sub.2<τ.sub.3 and corresponding eigenvectors μ.sub.1, μ.sub.2, μ.sub.3; an isotropic index I=τ.sub.1/τ.sub.3 and an elongation index E=1−τ.sub.1/τ.sub.3 are used to define the anisotropy of the structure; when E.fwdarw.0 and I.fwdarw.1, it indicates that the structure is isotropic; wherein the cluster is an independent structure composed of voxels of same material connected by common surfaces in binary data; if an individual voxel is material 1, and there is no common plane with other voxels around that are material 1, then this voxel is a cluster; a plurality of voxels that are identified the same material can be connected by common planes one by one into a large and very complex cluster, and a plurality of intersected cracks constitute a cluster; and a cluster is a structure that is interconnected inside but not connected with an outside world, in its original definition.

3. The method of separating, identifying and characterizing cracks in 3D space according to claim 2, wherein a process of the filtration mentioned in step 3) is as follows: when the cluster meets a filter criterion, modifying material identities of all the voxels in the cluster, that is, modifying an original void identity 1 to 0, which represents that this small void structure is changed into a part of a rock, and is not exist in a subsequent analysis, the identity 1 representing the void structure, the identify 0 representing a part of the rock.

4. The method of separating, identifying and characterizing cracks in 3D space according to claim 3, wherein a process of the smoothening mentioned in step 4) is as follows: performing smoothening and mending operations to the filtered image data, which specifically use a closed operation in mathematical morphology to realize smoothening and mending of the image, the morphology closed operation including a dilation operation and a following erosion operation; the dilation operation is a basic operation in morphology, an operation to evaluate local maximum; firstly, a previously defined structuring element B is convolved with an image A, that is, a maximum value of a gray scale of the pixel in a region of the image A covered by the structuring element B is calculated first, and the maximum value is assigned to a pixel specified by an origin of the structuring element B; this causes a highlighted area of the grayscale image to expand gradually, which, for a binary image, it is also equivalent to process an OR operation to the structuring element B and the image A, and an result of the operation determines a value of the pixel specified by the origin; wherein the structuring element is a structure in which the origin is centered, and its shape and size can be any; the erosion operation is a basic operation in morphology, and is mutually dual to the dilation, erosion being essentially an operation to evaluate local minimum, having an operating method similar to the dilation, which uses the local minimum for assigning to the pixel specified by the origin of the structuring element B, and this operation causes the highlighted area of the grayscale image to reduce gradually; for the binary image, it is also equivalent to process an AND operation to the structuring element B and the image A, and an result of the operation determines a value of the pixel specified by the origin; since the erosion and the dilation are mutually dual, a complementary set of the image A being eroded by the structuring element B equals to a complementary set of the image A being dilated by the structuring element B; the morphology closed operation processes a dilation operation to the filtered image first, and then processes an erosion operation.

5. The method of separating, identifying and characterizing cracks in 3D space according to claim 4, wherein a process of the thinning mentioned in step 5) is as follows: determining whether a current point is a void point needed to be processed; for the void point, calculating its extension in three directions, that is, checking point by point in both sides (positive and negative of the axis) of the three directions of x, y, z; if encountering a non-void point on one of the directions, then stopping the check on the direction and recording an extension length; after finishing checking the three directions, selecting the direction with the smallest extension length as the direction to be thinned; then thinning this direction from the boundary to the center to a thickness of d (d could be any number and recommended to be 2 or 3) voxels; except the d (2 or 3) voxels, assignments of other voxels on this direction are all changed during this process, that is, identified as 0, making them a part of the rock, and recording relevant information (thinning direction and thinning voxels); after completing the processing, proceeding to a next point.

6. The method of separating, identifying and characterizing cracks in 3D space according to claim 5, wherein the separation mentioned in step 6) is separating the thinned data, so that the originally intersected structures are separated and become independent structures, a specific process of the separation being as follows: (1) using a local grid for point-by-point process from the three-dimensions; generating two 2D square local grids in the three directions, respectively, with a current processing pixel as a center, and two given analysis radiuses (actually half-side-lengths), a large one and a small one, each point-by-point process uses only image information in the grid; taking the current processing pixel as an origin, with two local grids being divided into four quadrants by the origin; for the local grids with large outer radius, if the images in any two quadrants are linear and the angle between them meets a pre-defined criterion (an adjustable input parameter), while a number of void voxels contained in a corresponding quadrant in the grid with smaller inner radius is not less than a preset number (also an adjustable parameter), then they are taken as a projection of two cracks with different intersecting directions on a current section, thus changing a pixel value at the current origin, making it a part of the rock; and (2) for a 3D image, the processing step is to proceed from the three directions in turn until it is determined whether the current pixel should be classified as the rock.

7. The method of separating, identifying and characterizing cracks in 3D space according to claim 6, wherein a process of the combination mentioned in step 7) is as follows: firstly determining whether each void structure should be combined with another void structure, with determining criteria including: i) whether normal directions of the two void structures are approximately parallel; ii) whether they originally belong to the same void structure before the separation step; iii) whether their normal lines are perpendicular to a line of their centers; if the above three conditions are satisfied, then they are determined as should be combined; by this step, a cluster may be comprised of several pieces and represent an identified crack.

8. The method of separating, identifying and characterizing cracks in 3D space according to claim 7, further comprising: based on the information recorded in the thinning step, restoring the shape of each pore to that before thinning; and eventually, analyzing porosity, connectivity, and position, size, orientation and anisotropy of each void structure, finishing the separation, identification and characterization of the cracks in the 3D space.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0042] FIG. 1A and FIG. 1B are schematic views of the 2D crack network that is intersected, wherein FIG. 1A is a schematic view of that before identification, and FIG. 1B is a schematic view of each crack being performed separation and identification as a mixture of individual cracks, where different grey-scales represent different cracks.

[0043] FIG. 2 is a schematic view of the 3D crack structure of a rock sample.

[0044] FIG. 3 is the workflow of the invention of separating, identifying and characterizing of cracks.

[0045] FIG. 4 is a schematic view of a 2D slice before and after binarization.

[0046] FIG. 5 is a comparison chart of a rock sample before and after filtration.

[0047] FIG. 6 is a process diagram of implement of dilation algorithm.

[0048] FIG. 7 is a process diagram of implement of erosion algorithm.

[0049] FIG. 8 is a process diagram of implement of closed operation.

[0050] FIG. 9 is a schematic view of the smoothening effect.

[0051] FIG. 10 is a comparison chart before and after thinning.

[0052] FIG. 11A and FIG. 11B are schematic views of the separation process, wherein FIG. 11A shows a state before disconnection and FIG. 11 B shows that the intersection set of the cracks is disconnected.

[0053] FIG. 12A is a comparison chart before separation.

[0054] FIG. 12B is a comparison chart after separation.

[0055] FIG. 13 is individual structure diagram after separation.

[0056] FIG. 14 is a schematic view of T-type intersection.

[0057] FIG. 15 is a schematic view of turn-type intersection.

[0058] FIG. 16 is a schematic view of cross-shaped intersection.

[0059] FIG. 17 is a process diagram of combination.

[0060] FIG. 1.8 is a schematic view of a typical case 1.

[0061] FIG. 19 is a schematic view of a typical case 2.

[0062] FIG. 20A and FIG. 20B are comparison charts before and after combination (shown in a two dimensional plane).

[0063] FIG. 21 is a schematic view in turn of each structure after the combination.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

[0064] The present invention is further described below in combination with the accompanied drawings, but implementations of the present invention are not limited hereto.

[0065] Micro-CT images are commonly asked for routine preprocessing, including smoothening, denoising, edge cutting, image segmentation (identification of crack or pore). Thus the actual work requires 7 steps, and there are 6 main steps after the image segmentation, as shown in FIG. 3, that is:

[0066] 0) preprocessing digital image

[0067] 1) statistically analyzing basic information: basic information of the digital image, such as porosity, connectivity, structure anisotropy, etc.

[0068] 2) filtration: removing non-crack structure in the data, such as very small structures or structures with strong isotropic

[0069] 3) smoothening: smoothening and mending the digital image

[0070] 4) thinning: thinning the void structure into a certain thickness in a most narrow direction in three-dimensions

[0071] 5) separation: separating intersected cracks in a crack network by detecting the intersection and breaking the connection (independent from enclosed crack-surfaces)

[0072] 6) combination: unifying the labels of segments of a crack that are forcibly disconnected in the last step, merging very small structures that are formed during the separation a nearby large structure, restoring the shape of each void to that before thinning, and eventually analyzing porosity, connectivity, and position, size, orientation and anisotropy of each void structure, finishing the separation, identification and characterization of the cracks in the 3D space.

[0073] Data preprocessing is briefly introduced below, and the 6 main implementation steps are then introduced.

[0074] Preprocessing of the micro-CT images can be implemented by a variety of image processing software, where Avizo® software is the most common. First, a set of 3D CT images of rock samples are loaded into the image software, the original images being generally grayscale images. In the case of the image quality being not good enough, it may require smoothening and denoising of the images. Second, it is necessary to crop a desired part from the original data, or, to remove the non-sample part in the images. Image cutting is generally cut into a parallelepiped or cut into a cylinder. Third, it is necessary to transform the cut image into binary image. This step is called image segmentation, also known as binarization, which identifies the structure or material of interest, and interprets the rest as a matrix. The image after the segmentation only has two values, which are replaced with 0 or 1 in a computer, and the image display is converted to a black-and-white image, as shown in FIG. 4. The 3D images after the binarization can be seen in FIG. 2. There are many ways to segment the image, which are not introduced herein.

[0075] Statistically Analyzing of Basic Information

[0076] The purpose of this step is to obtain the basic information of the image, including: porosity, connectivity, and position, size, orientation, and anisotropy of each cluster. This information can be used for determining the values of the various parameters that require to be used in subsequent analysis steps. The implementation program is ctsta10.f90, written by the patent applicant Liu Jie.

[0077] In the basic information analysis process, it relates to a concept “cluster”. A cluster is an independent structure composed of voxels of same material (e.g. void, with a material identity of 1) having a common surface for neighboring two voxels one by one in binary data. If an individual voxel is material 1, and there is not common plane with other voxels around that are material 1, then this voxel is a cluster. A plurality of voxels that are identified the same material can be connected by common planes one by one into a large and very complex cluster, such as a pore network in sandstone. A plurality of intersected cracks also constitute a cluster. In its original definition, a cluster may also be interpreted as a structure that is interconnected inside but not connected with an outside world. All voxels in one cluster have the same unique label different from labels in other clusters.

[0078] Among many parameters, the anisotropy of the pore is a key point, and the anisotropy is calculated by a method as follows:

[0079] for a cluster containing n voxels, each voxel forms a vector from a center point of the cluster to a current position i presented as a.sub.i=(a.sub.xi, a.sub.yi, a.sub.zi).sup.T, and coordinates of a center of the cluster may be calculated by a mean value of all point coordinates in the cluster. Then an anisotropy of the structure can be represented by an orientation matrix R:

[00002]$R=\underset{i=1}{\overset{n}{.Math.}}.Math.a_i.Math.a_i^T=\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}^2& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{zi}}\\ \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}^2& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}.Math.a_{\mathrm{zi}}\\ \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{zi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}.Math.a_{\mathrm{zi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{zi}}^2\end{array}\right]$

[0080] The matrix has three eigenvalues τ.sub.1<τ.sub.2<τ.sub.3, and corresponding eigenvectors μ.sub.1, μ.sub.2, μ.sub.3. An isotropic index I=τ.sub.1/τ.sub.3 and an elongation index E=1−τ.sub.1/τ.sub.3 are used to define the anisotropy of the pore. When E.fwdarw.0 and I.fwdarw.1, it indicates that the cluster is isotropic.

[0081] There are output files from the statistically analyzing of basic information, with each file recording different information of the image.

[0082] 1) Basic information: provides the basic information of the image data, including porosity, specific surface area, number of the clusters (the cluster herein is the void structure within the sample), connectivity and other statistical information.

[0083] 2) Particle/pore throat distribution.

[0084] 3) Result of connectivity: simply lists all the information of percolating clusters in the current image, including their number and their percolating state. On the three directions of x, y, z, 0 represents not percolating, and 1 represents percolating.

[0085] 4) Statistics of cluster size: provides cluster number, number of voxels contained in each cluster, equivalent radius square, total number of voxels contained by the current cluster, and cumulative percentage of total void voxels.

[0086] 5) Shape of each cluster: provides information including cluster size, surface area, dimensions in principle directions, isotropy index, and elongation index.

[0087] 6) Orientation matrix of each cluster: provides the center coordinates and orientation matrix of each cluster.

[0088] 7) Normalized orientation matrix of each cluster: provides the normalized orientation matrix of each cluster and the direction angle of each eigenvector.

[0089] Filtration

[0090] The purpose of this step is to remove unimportant cluster in the model, making a subsequent calculation and analysis more convenient. Filters parameters that can be chosen include size, position, whether percolating, anisotropy, etc of the cluster.

[0091] For example, in general, natural samples contain a lot of tiny structures, representing small pores in the sample and noise generated when images are obtained. These small pores can be ignored in the crack analysis, so filtering these excessively small structures out contributes to the subsequent analysis.

[0092] The structural parameters obtained from the results of analysis process of the basic information provide the basis for determining filtering thresholds.

[0093] When any one cluster (structure) meets the filter condition, material identities of all voxels within the cluster will be modified, that is, modified from the original void identity 1 to 0, indicating that this small void structure transforms into a part of the rock, and is not exist in a subsequent analysis.

[0094] In addition to the case of tiny pores contained in the rock, there is also a case where a solid rock contained in the void structure is identified as 0. With regard to this, a file with opposite identity is generated by swapping the label in the data file (e.g. labmat=1−labmat), and the filtration is run again. After processing, the data needs to be reversed again and stored. FIG. 5 shows that excessively small voids and image noise are significantly suppressed after the filtration.

[0095] Smoothening

[0096] Filtered data serves as the input data of this step. As the complexity of the structure of the natural sample, the digital image that only has been filtered is still difficult to process, so in this step, the image is performed smoothening and mending operation, so that the image shape is as simple as possible and easy to process. The smoothening and mending of the image are realized by using the closing operation of mathematical morphology. The morphological closing operation includes a dilation operation and a following erosion operation.

[0097] Structuring element: can be any shape and size, and has a reference point called origin. In general, the structuring element refers to a square or a disc having an origin centered on a center, but other shapes may also be used.

[0098] Dilation operation: a basic operation in morphology. It is essentially an operation to evaluate local maximum. A previously defined structuring element B is used to be convolved with an image A, that is, a maximum value of a gray scale of the pixel in a region of the image A covered by the structuring element B is calculated first, and the maximum value is assigned to a pixel specified by an origin of the structuring element B. This causes a highlighted area of the grayscale image to expand gradually. For a binary image, it is also equivalent to process an OR operation to the structuring element B and the image A, and a result of the operation determines a value of the pixel specified by the origin (FIG. 6).

[0099] Erosion operation: a basic operation in morphology. Mutually dual to the dilation, erosion is essentially an operation to evaluate local minimum, having an operation method similar to the dilation, which uses the local minimum for assigning to a pixel specified by the origin of the structuring element B. This operation causes the highlighted area of the grayscale image to reduce gradually. For the binary image, it is also equivalent to process an AND operation to the structuring element B and the image A, and a result of the operation determines a value of the pixel specified by the origin (FIG. 7). Since the erosion and the dilation are mutually dual, a complementary set of the image A being eroded by the structuring element B equals to a complementary set of the image A being dilated by the structuring element B.

[0100] Closing operation: a dilation operation is processed to the image first, and then an erosion operation is processed. Generally, the closing operation may fill small pores and close little cracks, and an overall position and shape of the image remain the same (FIG. 8).

[0101] In the processing of 3D data, the structuring element is changed from square and disc to cube and sphere accordingly. Different shapes and sizes of the structuring element produce different processing results, and thus the structuring element should be larger than the size of all noise blocks to be removed. It also can increase times of the dilation and erosion, such as multiple erosions after multiple dilations, to make the smoothening effect better. However, it should be noted that the structuring element that is too large and erosion after multiple dilations may remove some features of the original image, and thus multiple attempts should be performed to select appropriate type, size of the structuring element, and the number of operations based on the processing effect (FIG. 9).

[0102] Thinning

[0103] For crack analysis, the image is required to be processed using thinning algorithm, and a figure that is approximate to the shape of original object and consists of simple curved surfaces is obtained.

[0104] The thinning algorithm proposed in the present embodiment is somewhat similar, but not identical, to the morphological thinning algorithm in 2D image processing. It should be noted that, the present invention processes 3D data. If the 3D image is processed only from a 2D perspective, it is necessary to regard the 3D image as a set of 2D images and to adopt a 2D method to each 2D slice. However, since the performance of the 3D image in 2D slices is not stable, even if they are adjacent slices, there may be a case that the difference is relatively large. If two 2D slices are used thinning algorithm alone, the image in the two slices will be staggered so that an original structure in the three dimensions is misinterpreted into two separate structures. In addition, if the 3D image is processed only from the 2D perspective, it means that the slices on the three directions are required to be processed simultaneously, and results of the slices on the three directions after processing will inevitably exist inconsistencies. How to synthesize the results of the three directions is also a very thorny problem. Therefore, the thinning method commonly used in 2D image is not applicable herein.

[0105] The thinning algorithm proposed in the present embodiment is performed from a 3D perspective. Firstly, a current point is determined whether a void point needed to be processed. For the void point, its extension in the three directions is calculated, that is, checked point by point in the three directions of x, y, z. If a non-void point is encountered on one of the directions, the check on that direction is stopped and an extension length is recorded. After the three directions are finished checking, the direction with smallest extension length is selected as the direction to be thinned. Then this direction is thinned from the boundary to the center to a thickness of d (d could be any number and recommended to be 2 or 3) voxels. Except the d (2 or 3) voxels, assignments of other voxels on this direction are all changed during this process, making them a part of the rock, and recording relevant information (thinning direction and thinning voxels) (FIG. 10). After the processing is completed, calculation of a next point is performed. It should be noted that, thinning is a manner to facilitate the processing of the image, in the subsequent steps, the thinning information recorded here will be used to restore the thinned pore.

[0106] Separation

[0107] This is the key to the present invention. This step takes the data that is thinned in the last step as an input, and perform the separation to the void structure in the image according to a series of pre-defined conditions, so that the originally intersected structure is separated into independent structures.

[0108] In this step, firstly, each cluster is numbered according to how many voxels are contained, and the largest structure is numbered as 1. Then spatial dimensions occupied by each structure in the 3D space, i.e., a minimum coordinates and a maximum coordinates of the structure spanned in three directions of x, y, z, are calculated separately. Next, in order to avoid being affected by other structures, each void structure will be processed separately, and the space to be processed is the rectangular space constituted by the minimum coordinates and the maximum coordinates on the 3D directions. Search space of the process is less than the original data space, which will greatly save the running time of the program.

[0109] The specific processing step is divided into two parts.

[0110] The first part uses local grids for processing point by point from three-dimensions. Two 2D square local grids are generated, with a current processing pixel as a center, and two given analysis radiuses (actually, half-side-lengths of square grids), a large one and a small one, and each point-by-point process uses only image information in the grid. The current processing pixel is taken as an origin, with two local grids being divided into four quadrants by the origin. If for the local grids with large radius, two void structures in any two quadrants are linear and the angle between them meet a pre-defined criterion (an adjustable input parameter), and for the local grid with small radius there are a number of void voxels (also an adjustable input parameter) contained in a corresponding quadrant, then they are taken as a projection of two cracks with different intersecting directions on a current section, thus changing a pixel value at the current origin, making it a part of the rock (FIG. 11A and FIG. 11B). The reason why using the two radiuses, an inner one and an outer one, as the criteria for determining is that, if the small grid contains too few number of void voxels, the calculated orientation of void will be unreliable. Therefore, we use the local grid with the larger radius to determine the orientation of the image, while use the grid with the smaller radius to avoid a phenomenon of too much removal caused by the large radius.

[0111] Since the processed image is 3D, the above processing steps will be performed sequentially from three directions until it is determined whether the current pixel should be classified as the rock. After this step, most of the intersected void structures can be broken into two or more independent void structures. The effect of disconnection depends on the void structure and parameter setting. It is recommended to test multiple times and select appropriately parameters based on the testing results. For example, when the outer radius is too large, it causes the range of disconnected voxels too large and damages the structure. When the radius is too small, the number of the void voxels contained in a very small grid may be too few, and the calculated anisotropy is unreliable. The range of the intersected angle is also necessary to be defined. In the present embodiment, an angle between the two structures can be regarded as approximately parallel when it is 0° to 30° and 150° to 180°, and the structures have an angle between 30° and 150° are considered as different intersected crack structures.

[0112] In the first part, if the number of the void voxels contained in the current local grid is few, it is possible that the void voxels are not forming a linear structure, so that intersection points are not removed in the first part. Thus there is a possibility that some of large intersected cracks are still connected by only a few voxels and keeping an interconnected cluster. This situation was not processed well in the first part, and it needs to be processed in the second part.

[0113] In the second part, it is also processed from three directions, but no longer point-by-point processing, but slice-by-slice processing. The space occupied by the current analyzed cluster is divided into a plurality of non-overlapping cubic sub-volumes, and then each sub-volume is analyzed slice by slice. In a small cube, the anisotropy of the void structure of a first slice of the 3D local grid is first analyzed along a Z-axis direction, and the is determined whether it is linear, and then the anisotropy of the void structure of a next slice is calculated and is determined whether it is linear. If the images of two adjacent slices are linear and their anisotropies have large difference, that is, the adjacent layers are not approximately parallel, then they are considered as different pore structures, and all the points in the second slice belonging to void are identified as a rock type, making the different void structures no longer in contact. The program then continues to use the same method to process the next two adjacent slices of image until the last slice. Then this implementation is repeated from a Y-axis direction and an X-axis direction, respectively. When it is all completed, an analysis of a next 3D sub-volume is performed.

[0114] The results of these two parts of separation operation are shown in FIG. 12A, FIG. 12B and FIG. 13.

[0115] The principle of the algorithm used to calculate the image direction in the local grid in this step is the same as that of the calculation of the 3D direction matrix R (see statistically analyzing basic information), and in case of two-dimensions, there is accordingly:

[00003]$R=\underset{i=1}{\overset{n}{.Math.}}.Math.a_i.Math.a_i^T=\left[\begin{array}{cc}\underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}^2& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}\\ \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{xi}}.Math.a_{\mathrm{yi}}& \underset{i=1}{\overset{n}{.Math.}}.Math.a_{\mathrm{yi}}^2\end{array}\right]$

[0116] The matrix has two eigenvalues τ.sub.1<τ.sub.2, and corresponding eigenvectors μ.sub.1, μ.sub.2. An isotropic index I=τ.sub.1/—T.sub.2.

[0117] The value of I is used to determine whether the projection of the void in the two-dimensions is linear. The two eigenvectors of the matrix can be calculated and their directions represent the angles between principle axes of the structure in this quadrant with the x-axis and y-axis, respectively, i.e., the direction of void structure can be determined.

[0118] Three typical cases often appear in the actual processing of the image are as follows:

[0119] 1) as shown in FIG. 14, when the cracks are two T-type intersecting cracks, using one voxel in the central blue part as an origin, the entire cluster is considered as consisting of three parts: a first quadrant A, a second quadrant B and a fourth quadrant C. The program will calculate the orientation matrix of A, B, and C, respectively. When the number of the void voxels around the origin is enough, the angles between A, B and C will continue to be calculated. If there is a cosine value of the angle having an absolute value less than a given threshold, then it means that there is an intersection situation that angles between A, B, and C is relatively large, and the void voxel at the origin needs to be changed to a rock voxel.

[0120] 2) When the crack is that shown in FIG. 15, the crack appears a large angle transition, which can also be considered as two cracks are connected, and needed to be disconnected. Using a voxel of the central blue part as an origin, the entire cluster is considered as consisting of two parts: a first quadrant A and a fourth quadrant B. The program will calculate the orientation matrix of A and B, respectively. When the number of the void voxels around the origin is enough, the angle between A and B will continue to be calculated. If there is a cosine value of the angle having an absolute value less than a given threshold, then it means that A and B belong to an intersection situation that the angle is relatively large, and the void voxel at the origin needs to be changed to the rock voxel.

[0121] 3) When the cracks are cross-shaped intersected shown as FIG. 16, using a voxel in the central blue part as an origin, the entire cluster is considered as consisting of four parts: a first quadrant A, a second quadrant B, a third quadrant C, and a fourth quadrant D. The program will calculate the orientation matrix of A, B, C and D, respectively. When the number of the void voxels around the origin is enough, the angles between A, B, C and D will continue to be calculated. If there is a cosine value of the angle between any of two of them having an absolute value less than a given threshold, it means that they belong to an intersection situation that the angle is relatively large, and the void voxel at the origin needs to be changed to the rock voxel.

[0122] Combination

[0123] This step takes the result of the last step as an input data, and is a correction of the previous steps. That is, for the actually elongated cracks, they are disconnected in the separation process. Through this step, they are combined into a cluster, this step also includes processing some small clusters generated during the separation, and restoring the crack to a state before being thinned.

[0124] Firstly, the structures in image is labeled. Since after the separation step, the image this time will contain more clusters than before. After completing labeling, the determination of whether each void structure should be combined with another void structure begins, with determining conditions including: i) whether normal directions of the two void structures are approximately parallel; ii) whether they originally belong to the same void before the separation step; iii) whether their normal lines are perpendicular to a line of their centers. If the above three conditions are satisfied, then they are determined as should be combined.

[0125] The combination operation includes re-labeling a supplier cluster number as a receiver cluster number, re-counting the total number of voxels contained in the two clusters, and recording the receiver cluster number, so that the supplier cluster is incorporated into the receiver cluster. During the program is running, the following situation may be encountered: the cluster used as the receiver is also used as the supplier and has been incorporated into other clusters. At this point, the receiver number recorded in the combination before is required to find the receiver cluster, and this pore is received and combined as a new receiver (FIG. 17).

[0126] As shown in FIG. 18, due to the separation algorithm, the cracks originally belong to the same cluster are now divided into three, and the central part has been changed to the rock, that is, black voxel. For these three independent cracks, their normal lines are L1, L2 and L3, respectively, and L4 is a line of centers of two of the cracks. Since the normal line L1 is not approximately parallel to L2 and L3, the crack of it is firstly determined as being not combined with the other two cracks. L2 and L3 are approximately parallel, and then continue to determine whether the line of the centers of the two cracks is perpendicular to their normal lines. It can be seen from the figure that, L4 is approximately perpendicular to L1 and L2, the two cracks originally belong to the same cluster, and thus they are determined as can be combined.

[0127] In the case shown as FIG. 19, the normal lines L1 and L2 of the two cracks are parallel, but the line L3 of their centers does not form as approximate perpendicular to L1 and L2, so the two cracks in this figure cannot be combined.

[0128] Since the separation process may produce some smaller fragmentary voids (such as the last picture of FIG. 15), these voids often do not have a clear orientation because of their small size, and thus they cannot be calculated their normal directions, thereby the above combination operation being difficult to be properly handled. Therefore, they are incorporated into the larger void structure by using a method of calculating between-class distance. Respective center points of these small voids are taken as representing their coordinates, minimum distances between them to the voxels of other pores are calculated in turn, and a larger void that is nearest is selected for incorporating the small void.

[0129] Finally, the shape of each pore is restored to that before being thinning using the information recorded in the thinning step. A final result can be obtained (FIGS. 20A, 20B and 21).

[0130] Characterization of Cracks

[0131] After completing all the processing, characterization of cracks is performed. Herein, only testing data shown in FIG. 13 is taken as an example, and major results are given in tables. Table 1 shows that the number of voxels contained in the processed data is reduced, which is due to change of identity values of some voxels during the separation of cracks, and the increase in surface area is also due to the disconnection of the intersecting cracks. There are only 3 clusters before the separation, and there are 8 clusters after the separation and identification. The change of size and morphological parameters of the cluster is very obvious. Data in Table 3 provide statistical characteristics of the length and width of each crack in the model. More information can be extracted from the 7 output files of ctsta10.f90.

TABLE-US-00001 TABLE 1 comparison of main parameters before and after separation and identification of cracks Parameter Before separation After separation porosity 3.18 2.847% surface area 177856 177346

TABLE-US-00002 TABLE 2 comparison of cluster size before and after separation and identification of cracks NO. SIZE RS{circumflex over ( )}2 ACCUPORE PERCENT 3 before 1 105891 4599.26 105891 0.98479 A total of 3 separation 2 927 32.56 106818 0.99341 clusters, one of .clus 3 709 112.88 107527 1.00000 which accounted for 98.5% of porosity. 8 after 1 38309 4107.55 38309 0.39865 Containing 8 separation 2 30234 3524.90 68543 0.71327 clusters .clus 3 13378 1490.33 81921 0.85248 (structure). 4 10017 683.37 91938 0.95672 5 2946 170.90 94884 0.98738 6 447 111.62 95331 0.99203 7 701 33.84 96032 0.99932 8 65 24.20 96097 1.00000

TABLE-US-00003 TABLE 3 comparison of cluster form before and after separation and identification of cracks Cluster Voxels Surface Dim-1 Dim-2 Dim-3 Iso-index Elong-index before 1 105891 175392 133.00 146.35 185.69 0.716 0.212 The first cluster separation 2 927 922 5.49 15.59 15.75 0.349 0.010 has similar .out2 3 709 1542 1.00 30.05 30.05 0.033 0.000 dimensions in three directions, and the isotropic index is > 0.7 after 1 30234 70439 2.05 172.83 189.33 0.011 0.087 All clusters have a separation 2 38309 60716 2.62 165.17 170.62 0.015 0.032 low isotropic .out2 3 13378 28887 1.96 92.51 123.63 0.016 0.252 index, indicating 4 10017 11414 3.07 69.21 78.32 0.039 0.116 that the structure is 5 2946 2698 3.28 31.58 41.55 0.079 0.240 in crack 6 701 1526 1.00 29.88 29.88 0.033 0.000 morphology. 7 447 1200 5.42 14.43 17.43 0.311 0.172 8 65 211 2.10 6.37 18.50 0.114 0.655

[0132] The 3D binary image with voxel as the basic unit can be performed the separation, identification and characterization of the cracks by six steps, basic characteristic analysis, filtering, smoothening, thinning, separation and combination. This technique provides important supports for researches such as crack morphology, permeability and mechanical response, and scaling up.

[0133] The above-described implementations of the present invention do not constitute a limitation to the scope of protection of the present invention. Any modifications, equivalent substitutions and improvements within the spirit and principle of the present invention shall be all included in the scope of protection of the claims of the present invention.