CN109242985B - Method for determining key parameters of pore structure from three-dimensional image - Google Patents
Method for determining key parameters of pore structure from three-dimensional image Download PDFInfo
- Publication number
- CN109242985B CN109242985B CN201811265575.0A CN201811265575A CN109242985B CN 109242985 B CN109242985 B CN 109242985B CN 201811265575 A CN201811265575 A CN 201811265575A CN 109242985 B CN109242985 B CN 109242985B
- Authority
- CN
- China
- Prior art keywords
- solid
- void
- radius
- voxels
- space
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T19/00—Manipulating 3D models or images for computer graphics
- G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Architecture (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Image Generation (AREA)
Abstract
The invention provides a method for determining key parameters of a pore structure from a three-dimensional image, which comprises the following steps: firstly, converting a test rock core to obtain a three-dimensional binary image formed by solid voxels and gap voxels, and respectively establishing a solid sphere structure element and a gap sphere structure element with corresponding radiuses; and converting the corresponding void space to obtain each void space and a communicated space block thereof, and identifying the pore and the pore throat to obtain the pore and pore throat characteristics of the void space in the current test core. The pore throat boundary in the whole identification process is clear, no human factors exist, on one hand, the microstructure characteristics of the porous medium pore structure can be effectively represented, on the other hand, the obtained key parameters of the pore structure are used as constraint conditions to establish a pore network model, the flow properties such as absolute permeability, relative permeability and the like can be effectively predicted under the condition of no coefficient adjustment, and the method has important significance for revealing pore size flow mechanism.
Description
Technical Field
The invention relates to the field of hydrogeology and petroleum geology, in particular to a method for determining key parameters of a pore structure from a three-dimensional image reflecting the pore structure of a reservoir.
Background
In the fields of petroleum engineering, hydrology, environmental engineering and the like, the representation of the micro-pore structure of the porous medium has important significance for establishing a mechanical model to reveal a flow mechanism of the pore size, and can guide the practical engineering application. The pore radius distribution of the porous medium can be obtained through high-pressure mercury injection, nitrogen adsorption and other core experiments, but the core experiments are high in time and economic cost, and the obtained pore structure parameters are limited. With the development of 3D imaging techniques and image reconstruction techniques, porous media can be rapidly scanned into three-dimensional digital images, the images being represented by a set of binary voxels, each voxel having a value of 1 (void voxel) or 0 (solid voxel), as shown in fig. 2. Digital images of suitable resolution may embody detailed microscopic information of the porous medium.
Based on three-dimensional digital images, the complex void space is reduced to a system consisting of pores and pore throats, and key parameters of the pore structure can be determined, including: pore number, pore throat number, coordination number distribution, pore radius distribution, pore throat shape factor distribution, tortuosity, and the like. The primary task in achieving this is how to identify complex void spaces as pores and pore throats.
There are three main methods of pore and pore throat identification that are widely used at present: the central axis Method (MA), the maximum sphere Method (MB), and the central axis-large sphere method (AB).
The MA method marks void voxels which do not affect the connectivity and topological properties of the void space, and the void voxels which are not marked form the central axis of the void space. The central axis is made up of a series of intersecting links, each representing a throat, and the intersections of the links representing pores. The MA method is very sensitive to irregular geometric details of the pore surface, and redundant links are formed, but it is difficult to define a definite standard to determine which links are redundant, excessive deletion of redundant links results in a small number of pores and pore throats, and partial deletion of redundant links results in a large number of pores and pore throats. In addition, after the pore and the pore throat are positioned through the central axis, the volume of the pore and the pore throat is required to be distributed by dividing the pore space, the pore space is guaranteed to occupy most of the pore space, the clear boundary of the pore and the pore throat is not clear in the process, and the dividing process is somewhat random and empirical.
The MB method comprises the steps of approaching a void space through a series of maximum balls, establishing a family tree relationship according to the size and the relative position of each maximum ball, marking the positions of ancestor balls corresponding to a void position and children balls corresponding to a pore throat position, marking all the maximum balls as void balls or pore throat balls, and dividing the space occupied by void blocks and pore throat blocks. The radius threshold for determining whether a maximum sphere is an ancestor sphere when the MB method locates pores and pore throats through family tree relationships is an empirical value. The segmentation of the space occupied by the pore block and the pore throat block has certain experience and roughness without clear pore and pore throat boundaries in the segmentation process.
The AB method absorbs the advantages of the MA method and the MB method, the maximum sphere is established under the constraint of the central axis, the family tree relationship is established, the position of the ancestor sphere corresponds to the position of the pore, the position of the child sphere corresponds to the position of the pore throat, and finally the space occupied by the pore block and the pore throat block is divided by adopting a double-speed expansion algorithm. AB method the threshold radius to determine whether a maximum sphere is an ancestor sphere in locating pores and pore throats is also an empirical value. In the process of dividing the space occupied by the pores and the pore throats, the clear definition of the pore and pore throat boundaries is not provided, and the division of the pore blocks and the pore throat blocks has certain experience.
Disclosure of Invention
Based on the problems in the prior art, the method for determining the key parameters of the pore structure from the three-dimensional image of the porous medium is established, the pore space can be identified as the pore and the pore throat according to the morphological characteristics of the pore and the pore throat, and the boundary of the pore and the pore throat is clear and clear without human factors in the identification process.
In particular, the present invention provides a method for determining key parameters of pore structure from a three-dimensional image, comprising the steps of:
and 500, synthesizing all the three-dimensional binary images after identifying the pores and the pore throats to obtain the characteristics of the pores and the pore throats in the pore space in the current test core.
In one embodiment of the present invention, 0 in the three-dimensional binary image represents a solid voxel of a solid space, 1 represents a void voxel of a void space, the unit of radius is one voxel, and one solid voxel or one void voxel corresponds to one voxel.
In one embodiment of the invention, the solid spherical structural element takes a solid voxel as a center, all solid voxels are in the radius range of the solid spherical structural element, and the rest solid voxels are void voxels; the gap spherical structural elements take a gap voxel as a center, all the gap voxels are in the radius range, and the rest are solid voxels.
In one embodiment of the present invention, when converting solid voxels or void voxels located at a boundary, a layer of solid voxels of a corresponding radius size needs to be added outside the corresponding boundary.
In one embodiment of the invention, when converting solid voxels, the conversion is performed only for solid voxels at the inner edge of the void space.
In one embodiment of the invention, the morphology of the pores is such that the radii of their connected void spaces are not all greater than their own radii, while the morphology of the pore throats is such that the radii of their connected void spaces are all greater than their own radii.
In one embodiment of the present invention, the process of determining the connected space block of the current radius void space in the three-dimensional binary image is as follows:
step 411, converting the original three-dimensional binary image by using the solid sphere structural element and the void sphere structural element corresponding to the current radius to obtain a three-dimensional binary image showing all unfilled void spaces larger than the radius and filled void spaces smaller than or equal to the radius by the solid voxel;
step 412, selecting a solid sphere structural element and a gap sphere structural element corresponding to a radius smaller than a voxel unit of the current radius to convert on the basis of the original three-dimensional binary image, so as to obtain a three-dimensional binary image formed by all unfilled gap spaces larger than or equal to the current radius and the filled gap spaces smaller than the current radius and filled by solid voxels;
step 413, performing boolean operation on the three-dimensional binary image in the steps 411 and 412, namely completely separating the void space of the current radius in the step 411, and performing connectivity analysis on the void space in a six-connectivity manner by using a bwleaeln function of MATLAB to obtain a connected space block of the void space of the radius; and repeating the steps to obtain the communicated space blocks with all the radius void spaces.
In one embodiment of the present invention, the process of identifying all connected spatial blocks of the current radius void space as pores and pore throats based on their morphological characteristics is as follows:
step 421, projecting each surface of each connected space block to obtain a projection plane which is away from the surface by one pixel from the outer normal direction of the corresponding surface;
step 422, if there are only solid voxels in the projection plane, this surface is not connected to larger void spaces; if there are only void voxels in the projection plane, this surface is connected to a larger void space; if the surface is not projected, the surface is a boundary;
if the projection plane contains both solid voxels and void voxels, the surface is extended by one or more voxel distances in the normal direction to obtain further projection planes, and then the following three conditions can be obtained according to the voxel conditions in the projection planes:
1) if a plurality of further projection surfaces represent all void voxels, judging that the surfaces are connected with a large void space;
2) if further several projection planes are represented as both void voxels and solid voxels, performing voxel analysis on the projection planes of other five surfaces of the connected block, and if the remaining five projection planes are all judged to have both solid voxels and void voxels, directly judging that the connected block is a pore throat;
3) if the two conditions are not met, judging that the surface is not connected with a large void space;
step 423, according to the above judgment result, if a connected block has two surfaces connected with a larger void space, or the projection planes corresponding to the six surfaces simultaneously contain solid voxels and void voxels, the void space block in the connected block is identified as the throat; for a boundary connected block, if a surface connects larger void spaces, or each existing projection plane contains both void and solid voxels, then the connected block is identified as a pore throat; the rest of the communicating blocks are the pores.
The invention provides a Pore-throat morphology recognition method PTM (Pore-ThroatMorphology) based on a porous medium three-dimensional image, and pores and Pore throats are recognized according to morphological characteristics of the pores and the Pore throats to calculate key parameters of a Pore structure. The pore throat boundary of the pore is clear in the whole identification process, and no human factor exists. The method has the advantages that more comprehensive pore structure parameters can be obtained, on one hand, the micro characteristics of the pore structure of the porous medium can be effectively represented, on the other hand, the obtained key parameters of the pore structure are used as constraint conditions to establish a pore network model, the flow properties such as absolute permeability and relative permeability can be effectively predicted under the condition that the coefficient is not adjusted, and the method has important significance for revealing pore size flow mechanisms. Can guide oil gas exploration and development in the field of petroleum engineering.
Drawings
FIG. 1 is a flow chart of a method for determining key parameters of pore structure from three-dimensional images, according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of a three-dimensional binary image according to an embodiment of the invention;
FIG. 3 is a schematic representation of solid voxels in a two-dimensional map of solid spherical structural elements in accordance with one embodiment of the present invention;
FIG. 4 is a schematic representation of void voxels in a transformed two-dimensional map of void spherical structural elements in accordance with an embodiment of the present invention;
FIG. 5 is a schematic diagram of a transition with a radius of 2 for a solid spherical structural element in accordance with one embodiment of the present invention;
FIG. 6 is a schematic diagram of a process for solid voxel and void voxel conversion at a current radius according to one embodiment of the present invention;
FIG. 7 is a schematic diagram of projection analysis of a connected block in accordance with an embodiment of the present invention;
FIG. 8 is a schematic illustration of a process for identifying pores and pore throats in accordance with an embodiment of the present invention;
FIG. 9 is a schematic illustration of a simplified representation of pores and pore throats in accordance with an embodiment of the present invention;
FIG. 10 is a schematic representation of the pores and pore throats after total analysis in one embodiment of the invention;
FIG. 11 is a schematic view of a pore throat section taken in accordance with an embodiment of the present invention;
FIG. 12 is a schematic view of the communication paths of the pores in different orientations in one embodiment of the invention.
Detailed Description
As shown in FIG. 1, the novel method for determining key parameters of pore structure from three-dimensional images of one embodiment of the present invention generally comprises the steps of:
for a rock sample, the rock can be scanned by imaging technical means such as CT or FIB and the like to obtain a three-dimensional image corresponding to the rock, the three-dimensional image is a gray image and consists of numbers in the range of 0-255, and the numbers in different ranges correspond to different gray scales to represent different components, such as: the components of the gray scale representation corresponding to the numbers ranging from 190 to 255 are gaps, and the components of the gray scale representation corresponding to the numbers ranging from 0 to 189 are rock skeletons. And then carrying out threshold division on the three-dimensional gray-scale image, and setting the value of any voxel in the image to be 1 if the value of the voxel is in the range of 190-255 so as to represent a gap, and setting the value of the voxel to be 0-189 if the value of the voxel is in the range of 0-189. Setting the voxel value to be 0 to represent a rock skeleton, and finally obtaining a three-dimensional binary image as shown in fig. 2, wherein the value of each voxel in the three-dimensional binary image is no longer a number within the range of 0-255 but 0 or 1, if 0, the rock skeleton is represented, if 1, a gap is represented, and the finally formed state is the three-dimensional binary image, which is the basis of the following analysis.
In the following description, 0 in the three-dimensional binary image represents a solid voxel of a solid space, 1 represents a void voxel of a void space, and the unit of radius is one voxel, for example, radius 10 is ten voxels; one solid voxel or one void voxel corresponds to each voxel.
the solid spherical structural element takes a solid voxel as a center, all solid voxels are arranged in the radius range of the solid spherical structural element, the rest solid voxels are gap voxels, and the size of the solid spherical structural element is 2 x the radius + 1. The gap spherical structural element takes a gap voxel as a circle center, all gap voxels are in the radius range, the rest solid voxels are, and the size of the gap spherical structural element is 2 x the radius + 1.
For example, the central voxel of the solid spherical structural element is 0, the voxel values within the voxel range with the distance from the central voxel being the radius are also 0, and the rest of the voxel values are 1. If the diameter is 10, the voxel values within a distance of 10 voxels from the central voxel are all 0. The interstitial spherical structural elements are the same, only 0 and 1 are interchanged.
as shown in fig. 3, a represents a two-dimensional image represented by solid voxel 0 and void voxel 1, b is a two-dimensional representation of a solid spherical ball structural element, and c is a two-dimensional image converted by the solid voxel. The conversion process is illustrated as follows, taking as an example the solid voxel conversion: respectively corresponding the center of a solid sphere structural element with a corresponding radius to each solid voxel in a three-dimensional binary image, then corresponding other voxels in the solid sphere structural element to voxels beside the corresponding solid voxel one by one, and then performing logic and operation, wherein if the value of a certain gap voxel in the three-dimensional binary image within the range of the solid sphere structural element is 1 and the value of a corresponding voxel in the solid sphere structural element is 0, 1&0 is 0, and the gap voxel in the solid sphere structural element in the three-dimensional binary image is changed into a solid voxel; if a certain solid voxel value in the solid spherical ball structural element in the three-dimensional binary image is 0, and the corresponding voxel value in the solid spherical ball structural element is 0, 0&0 is 0, and the solid voxel in the solid spherical ball structural element in the three-dimensional binary image is still a solid voxel.
As shown in fig. 4, a represents a two-dimensional image represented by a solid voxel 0 and a pore voxel 1, b is a two-dimensional representation of a pore spherical structural element, and c is a two-dimensional image after transformation of a pore voxel. The conversion process of the interstitial spherical structural elements is the same as that of the solid spherical structural elements and is not repeated here.
As shown in fig. 5 and 6, in the above conversion, the void space less than or equal to the radius of the current solid sphere structural element is completely filled with the solid voxel; since the void space larger than the radius has the void voxels, after the transformation of the void sphere structural element with the same radius, the replaced void voxels are restored to be original, so that the current three-dimensional binary image accurately takes the space with the current radius as a boundary line, and the existing void spaces are all clearly represented to be the void spaces larger than the current radius. In the iteration step with the minimum radius, performing Boolean operation on the void space in the original three-dimensional binary image and the void space currently converted into the solid space to obtain the void space with the current iteration radius, and in the iteration steps with other radii, performing Boolean operation on the void space obtained in the iteration step with the previous radius and the void space currently converted into the solid space to obtain the void space with the current iteration radius;
after the void space corresponding to each radius is converted and subjected to Boolean operation, the converted three-dimensional binary images are respectively and independently stored, the void space in the three-dimensional binary images only has the void space corresponding to the radius, and the void spaces with the rest radii are converted into solid spaces
Fig. 5 (a) illustrates a state of a solid voxel corresponding to a solid spherical structural element, in which three void voxels enclosed inside the solid spherical structural element are converted into solid voxels for a solid voxel i to which the solid spherical structural element having a radius of 2 is attached; and adding solid voxels outside the image boundary to prevent overflow of the solid spherical structural elements, wherein the number of layers of the added solid voxels is equal to the radius of the solid spherical structural elements. (b) The norm of each voxel in the solid spherical structural element is the distance from the voxel to the center of the solid spherical structural element, and the maximum norm is the radius of the solid spherical structural element.
In fig. 6 a is a schematic diagram of solid voxel conversion; b is a schematic of void voxel transition; c is a schematic diagram of the converted three-dimensional binary image; d is a schematic diagram of the result of Boolean operation with c.
And finally, combining all independently stored three-dimensional binary images to obtain a void space map under each radius, and laying a foundation for later analysis.
During specific conversion, a radius range can be visually observed according to the size of the actual void space, then the three-dimensional binary image is converted from the radius 1, after the current radius conversion is finished, the radius of one unit (one voxel) is increased each time, and the conversion is continued until the radius reaches the maximum radius.
If the maximum radius of the actual rock core is 40 voxels, the maximum radius is visually selected to be 45 voxels, solid spherical structural elements and pore spherical structural elements with the radius corresponding to the pore space of the 45 voxels are established, the three-dimensional binary image is converted, corresponding space gaps do not exist in the state, therefore, no effect is generated, and the conversion effect can be realized only when the radius is within 40.
Furthermore, to improve conversion efficiency, solid spherical ball structural elements may be converted only for solid voxels at the inner edge of each void space, with solid voxels in other locations not being converted.
At the time of conversion, if a solid voxel and a pore voxel at the boundary are encountered, solid voxels of the corresponding radius layer are added outside the boundary on the three-dimensional binary image. When converting the void space of radius 10 at the boundary, a solid voxel layer of 10 pixels is added outside the boundary to avoid overflow.
as shown in fig. 7 and 8, the morphology of the pores is such that the radii of the void spaces to which they are attached are not all greater than their own radii, while the morphology of the pore throats is such that the radii of the void spaces to which they are attached are all greater than their own radii. For the pore throat at the boundary, since one section is the boundary of the core, the radius of the void space connected at the other end is larger than the radius of the void space.
Since 1 in the three-dimensional value image represents a pore space, 0 represents a rock skeleton. Therefore, the void space is determined, and in the subsequent pore throat identification, only the voxel with the value of 1 needs to be searched in the three-dimensional binary image.
The void spaces to be analyzed can be classified by radius, since the void spaces of the corresponding radius can be located by interconversion and boolean operations of the void voxels and the solid voxels at each radius; then, connectivity analysis is carried out on the void spaces to obtain a series of connected space blocks, and the connected space blocks under the same radius can be numbered according to the obtained sequence.
The process of determining each connected space block according to the three-dimensional binary image is as follows:
step 411, converting the original three-dimensional binary image by using the solid sphere structural element and the void sphere structural element corresponding to the current radius to obtain a three-dimensional binary image showing all unfilled void spaces larger than the radius and filled void spaces smaller than or equal to the radius by the solid voxel;
the current void space radius is the void space to be processed which is selected by increasing the radius; as can be seen from the foregoing description, after the solid spherical structural elements and the void spherical structural elements corresponding to the solid spherical structural elements are converted, the void spaces with the radius smaller than or equal to the radius are all filled with the solid voxels. While void spaces larger than this radius remain intact.
Step 412, selecting a solid sphere structural element and a gap sphere structural element corresponding to a radius smaller than a voxel unit of the current radius to convert on the basis of the original three-dimensional binary image, so as to obtain a three-dimensional binary image formed by all unfilled gap spaces larger than or equal to the current radius and the filled gap spaces smaller than the current radius and filled by solid voxels;
step 413, performing boolean operation on the three-dimensional binary image in the steps 411 and 412, namely completely separating the void space of the current radius in the step 411, and performing connectivity analysis on the void space in a six-connectivity manner by using a bwleaeln function of MATLAB to obtain a connected space block of the void space of the radius; and repeating the steps to obtain the communicated space blocks with all the radius void spaces.
Wherein a in FIG. 7 is a schematic diagram of a linear pore throat and a pore connected at both ends, and b is a schematic diagram of a diagonal pore throat and a pore connected at both ends; the above analysis process all needs to be based on the three-dimensional binary image after the previous corresponding radius down-conversion; through the logical calculation process of Boolean operation, the content of the corresponding radius void space can be isolated, and connectivity analysis is carried out, so that all connected space block distribution maps of the current radius void space in the three-dimensional binary image are established.
The process of identifying all connected spatial blocks of the current radius void space as pores and pore throats from their morphological features is as follows:
step 421, projecting each surface of each connected space block to obtain an image slice which is one pixel away from the surface from the outer normal direction of the corresponding surface;
here it is necessary to process all connected space blocks at one radius at a time and then process all connected space blocks at the next radius. Each connected block, having been determined by the bwleaeln function of MATLAB, can have its corresponding position determined directly by the regionprops function of MATLAB.
As shown in fig. 9, in the analyzed three-dimensional binary image, only the connected space blocks corresponding to the radii of the pore space to be analyzed are obtained, and connectivity analysis is performed on the pore spaces of the radii in a 6-way connection manner, so that the corresponding connected space blocks are obtained. When positioning the connected space block, the connected space block can be regarded as a regular cuboid, i.e. can be composed of connected void voxels and some solid voxels around them, regardless of whether there are projections or pits. Each pore space block is positioned by a vertex and side length of the corresponding communicating block. And selecting a solid sphere structure element and a gap sphere structure operation corresponding to a voxel radius smaller than the current radius from the original three-dimensional binary image to convert the original three-dimensional binary image, filling all the gap spaces smaller than the current radius by the solid voxels, and only the gap spaces larger than or equal to the current radius in the image. And then determining all connected space blocks corresponding to the current radius according to the determined positions of the connected blocks in the processed original three-dimensional binary image. During projection, taking the front surface of the connected block as an example, a projection surface of which the front surface advances by one unit in the normal direction is obtained, and theoretically, the projection surface is beyond the range of the current connected block but can just reflect the connection state of the connected block.
The specific state of each surface connection includes the following cases:
step 422, if there are only solid voxels in the projection plane, this plane is not connected with larger void space; if only a void voxel exists in the projection plane, the plane is connected with a larger void space; if the surface is not projected, the surface is a boundary;
if the projection plane contains both solid voxels and void voxels, then the surface is extended by one or more voxel distances in the normal direction to obtain further projection planes, and then the following three conditions can be obtained according to the content of the voxel conditions in the projection planes:
1) if a plurality of further projection surfaces represent all void voxels, judging that the surfaces are connected with a large void space;
2) if further several projection surfaces are represented by existing gap voxels and solid voxels, acquiring projection surfaces of other five surfaces of the communicating block for judgment, if the remaining five projection surfaces are the existing solid voxels and the gap voxels, the communicating block is in an inclined state, two ends of the communicating block are respectively connected with large gap spaces, and the communicating block is directly judged as a pore throat; .
3) If the two conditions are not met, the surface is judged not to be connected with large void space.
Step 423, according to the above judgment result, if a connected block has two surfaces connected with a larger void space, or the projection planes corresponding to the six surfaces simultaneously contain solid voxels and void voxels, the void space block in the connected block is identified as the throat; for a boundary connected block, if a surface connects larger void spaces, or each existing projection plane contains both void and solid voxels, then the connected block is identified as a pore throat; the rest of the communicating blocks are the pores.
The above process is repeated, and the analysis results of the pores and pore throats of the void spaces with all radii in the three-dimensional binary image can be obtained.
Specific examples are as follows: let six surfaces of the connected block be marked as N, S, W, E, U, D respectively. The plane of projection of each surface is derived from an image slice one pixel from the surface in the direction normal to the corresponding surface. The voxel condition in each projection plane is as follows:
1) if there are only solid voxels in the projection plane, the corresponding surfaces of the void space blocks in the connected block are either pore walls or interfaces with void spaces that are transformed into solid spaces with radii smaller than the iteration radius. The void space blocks are not connected to larger void spaces in the surface direction, see fig. 3 a; the projection plane corresponding to the S, N, U, D surface of the connecting block C1 and C2).
2) If there are only void voxels in the projection plane, the corresponding surface of the void space block in the connected block is the interface with the larger void space. The reason is that the void spaces with a radius smaller than the iteration radius have been converted into solid spaces, and after connectivity analysis these void space blocks corresponding to the iteration radius cannot be connected to void spaces of the same radius, and therefore the void space blocks can only be connected to larger void spaces in the surface direction.
3) There are three cases if the projection plane contains both solid voxels and void voxels. First, the surface of the connected block corresponding to the pore space block may be the interface with the larger pore space, where another interface radius in the projection plane belonging to the larger pore space is smaller than the iteration radius. To further determine this possibility, we extended the acquired image slices in the surface outside normal direction to determine the corresponding voxel cases, and if some of the projection planes contain only void voxels, the void space block is connected to a larger void space in the surface direction (FIG. 3 a; projection plane corresponding to surface W of connected block C2). Otherwise, we analyze the projection planes corresponding to the 6 surfaces of the connected block simultaneously, and if each projection plane contains a void voxel and a solid voxel, the block of void space connects more than one larger void space in different surface directions, i.e., the block of void space is the throat (FIG. 3; six corresponding projection planes of connected block C3). Furthermore, if the voxel condition does not satisfy both of the above conditions, the pore space block does not have a pore space with a larger boundary in the surface direction.
And 500, synthesizing all the three-dimensional binary images after identifying the pores and the pore throats to obtain the characteristics of the pores and the pore throats of all the pore spaces in the current test rock core, and finally obtaining the pore number and the pore throat number of the whole rock core.
As shown in FIG. 10, the invention proposes a Pore-Throat Morphology recognition method PTM (Pore-through Morphology) based on a three-dimensional image of a porous medium, and the pores and the Pore throats are recognized according to the morphological characteristics of the pores and the Pore throats to calculate key parameters of the Pore structure. The pore throat boundary of the pore is clear in the whole identification process, and no human factor exists. The method has the advantages that more comprehensive pore structure parameters can be obtained, on one hand, the micro characteristics of the pore structure of the porous medium can be effectively represented, on the other hand, the obtained key parameters of the pore structure are used as constraint conditions to establish a pore network model, the flow properties such as absolute permeability and relative permeability can be effectively predicted under the condition that the coefficient is not adjusted, and the method has important significance for revealing pore size flow mechanisms. Can guide oil gas exploration and development in the field of petroleum engineering.
After the number of pores and the number of pore throats of the whole core are obtained, the relation parameters related to the pore structure can be accurately determined, and the specific key parameter contents are as follows:
1. pore volume-radius distribution and pore throat volume-radius distribution
In each iteration, the radius of the identified pores and pore throats is equal to the iteration radius, and the volume of the identified pores and pore throats can be determined by counting the number of voxels in the corresponding pore space block. The final series of pore volume-radii and pore throat volume-radii form a pore volume-radius distribution and a pore throat volume-radius distribution, respectively.
2. Distribution of coordination number
The coordination number of a pore means the number of pore throats attached to the pore. After the identification process of the pores and pore throats is finished, the coordination number of each pore can be determined by calculating the number of the pore throats connected to each pore, and finally, the coordination number distribution is generated according to the statistical result of the pore number-pore coordination number.
3. Pore throat shape factor distribution
In fig. 11, a is a grayscale image of the core, b is a perspective view of the pore throat extracted from a, and c is an enlarged schematic cross-sectional view of b. It is shown that the key to determining the pore throat shape factor is to find the true pore throat cross-section perpendicular to the flow path, and in order to be able to obtain this as much as possible, we use the interface between the pore and the pore throat to achieve this. The interface to the pore in each pore throat is first determined, then all pore throat cross-sections parallel to the interface are obtained, and finally the shape factor of each cross-section is calculated by the following formula:
G＝A/P2
wherein G is a shape factor;
a is the area of the cross-section;
p is the perimeter of the cross-section;
and taking the average value of all the section shape factors as the shape factor of the pore throat, and finally obtaining the distribution of the pore throat shape factors according to the statistical result of the pore throat shape factors and the pore throat number.
In the identified pore and pore throat system, each pore throat is numbered from 1 in sequence, and then each voxel value in each pore throat is changed from 1 to the number value of the pore throat; and determining the voxel condition on the projection plane of the six surfaces of each pore according to the method, wherein if the projection plane of a certain direction of a pore has some voxels with the value of n, the plane of the certain direction is the interface of the pore and the pore throat with the number of n, and finally, the interfaces of all pores and pore throats can be obtained.
4. Tortuosity
In fig. 12, a is an X-direction communication path, b is a Y-direction communication path, and c is a Z-direction communication path. Tortuosity is the average extension of the flow path. The preferential flow path in the pore space can be approximated by the central axis, so we use the central axis to calculate the tortuosity of the pore space; for one direction, firstly, a central axis in the middle of a pore in each image slice perpendicular to the direction is determined, then, a starting point of the communication path is obtained through the intersection point of the central axes in boundary slices of the image slices, and finally, the communication path is obtained by conducting guidance search from each starting point through the central axis. Once all the communication paths in the direction are determined, the tortuosity in the direction is calculated by the following equation:
in the formula taudIs the tortuosity of that direction; n is the number of communication paths in the direction; liIs the actual length of the communication path i; lsIs a straight length corresponding to the communication path i; finally, the average tortuosity across the sample pore space was calculated as:
τ＝(τx+τy+τz)/3
wherein τ is the average tortuosity of the sample pore space; tau isxIs tortuosity in the X direction; tau isyIs tortuosity in the Y direction; tau iszIs the tortuosity in the Z direction.
Thus, it should be appreciated by those skilled in the art that while a number of exemplary embodiments of the invention have been illustrated and described in detail herein, many other variations or modifications consistent with the principles of the invention may be directly determined or derived from the disclosure of the present invention without departing from the spirit and scope of the invention. Accordingly, the scope of the invention should be understood and interpreted to cover all such other variations or modifications.
Claims (6)
1. A method for determining key parameters of pore structure from a three-dimensional image, comprising the steps of:
step 100, selecting a test rock core, scanning by using a three-dimensional imaging technology to obtain a three-dimensional gray image, and performing threshold division conversion on the three-dimensional gray image to obtain a three-dimensional binary image formed by solid voxels and void voxels;
step 200, respectively establishing a solid spherical ball structure element of a conversion solid voxel and a gap spherical ball structure element of a conversion gap voxel with corresponding radiuses according to the radiuses from the minimum gap to the maximum gap in the three-dimensional binary image;
step 300, starting iteration from the minimum gap radius, respectively corresponding each solid voxel by using the circle center of the solid spherical structural element with the corresponding radius, converting the gap voxels in the range of the solid spherical structural element into solid voxels, respectively corresponding each remaining gap voxel by using the circle center of the gap spherical structural element with the same radius, converting the solid voxels in the range of the gap spherical structural element into gap voxels, converting the gap space smaller than or equal to the current radius into a solid space, and keeping the gap space larger than the current radius unchanged; in the iteration step with the minimum radius, performing Boolean operation on the void space in the original three-dimensional binary image and the void space currently converted into the solid space to obtain the void space with the current iteration radius, and in the iteration steps with other radii, performing Boolean operation on the void space obtained in the iteration step with the previous radius and the void space currently converted into the solid space to obtain the void space with the current iteration radius; then storing the currently converted three-dimensional binary image, and then performing next radius conversion by taking the original three-dimensional binary image as a basis;
step 400, after all radii are processed, performing connectivity analysis on the void space in the three-dimensional binary image of each radius, determining all connected space blocks of the void space of the current radius in the three-dimensional binary image, then identifying all connected space blocks of the void space of the current radius as pores and pore throats according to morphological characteristics of the pores and the pore throats, and storing the identified three-dimensional binary image;
500, synthesizing all three-dimensional binary images after identifying the pores and the pore throats to obtain the characteristics of the pores and the pore throats of the pore space in the current test core;
the process of determining the connected space block of the current radius void space in the three-dimensional binary image is as follows:
step 411, converting the original three-dimensional binary image by using the solid sphere structural element and the void sphere structural element corresponding to the current radius to obtain a three-dimensional binary image showing all unfilled void spaces larger than the radius and filled void spaces smaller than or equal to the radius by the solid voxel;
step 412, selecting a solid sphere structural element and a gap sphere structural element corresponding to a radius smaller than a voxel unit of the current radius to convert on the basis of the original three-dimensional binary image, so as to obtain a three-dimensional binary image formed by all unfilled gap spaces larger than or equal to the current radius and the filled gap spaces smaller than the current radius and filled by solid voxels;
step 413, performing boolean operation on the three-dimensional binary image in the steps 411 and 412, namely completely separating the void space of the current radius in the step 411, and performing connectivity analysis on the void space in a six-connectivity manner by using a bwleaeln function of MATLAB to obtain a connected space block of the void space of the radius; repeating the steps to obtain communicated space blocks of all the radius void spaces;
the process of identifying all connected spatial blocks of the current radius void space as pores and pore throats from their morphological features is as follows:
step 421, projecting each surface of each connected space block to obtain a projection plane which is away from the surface by one pixel from the outer normal direction of the corresponding surface;
step 422, if there are only solid voxels in the projection plane, this surface is not connected to larger void spaces; if there are only void voxels in the projection plane, this surface is connected to a larger void space; if the surface is not projected, the surface is a boundary;
if the projection plane contains both solid voxels and void voxels, the surface is extended by one or more voxel distances in the normal direction to obtain further projection planes, and then the following three conditions can be obtained according to the voxel conditions in the projection planes:
1) if a plurality of further projection surfaces represent all void voxels, judging that the surfaces are connected with a large void space;
2) if further several projection planes are represented as both void voxels and solid voxels, performing voxel analysis on the projection planes of other five surfaces of the connected block, and if the remaining five projection planes are all judged to have both solid voxels and void voxels, directly judging that the connected block is a pore throat;
3) if the two conditions are not met, judging that the surface is not connected with a large void space;
step 423, according to the above judgment result, if a connected block has two surfaces connected with a larger void space, or the projection planes corresponding to the six surfaces simultaneously contain solid voxels and void voxels, the void space block in the connected block is identified as the throat; for a boundary connected block, if a surface connects larger void spaces, or each existing projection plane contains both void and solid voxels, then the connected block is identified as a pore throat; the rest of the communicating blocks are the pores.
2. The method of claim 1,
0 in the three-dimensional binary image represents a solid voxel of a solid space, 1 represents a void voxel of a void space, the unit of radius is one voxel, and one solid voxel or one void voxel corresponds to one voxel.
3. The method of claim 1,
the solid spherical structural element takes a solid voxel as a center, all solid voxels are in the radius range of the solid spherical structural element, and the rest solid voxels are void voxels; the gap spherical structural elements take a gap voxel as a center, all the gap voxels are in the radius range, and the rest are solid voxels.
4. The method of claim 1,
when solid voxels or void voxels located at the boundary are converted, a layer of solid voxels of a corresponding radius size needs to be added outside the corresponding boundary.
5. The method of claim 1,
when solid voxels are transformed, only solid voxels at the inner edge of the void space are transformed.
6. The method of claim 1,
the morphology of the pores is characterized in that the radii of their associated void spaces are not all greater than their own radii, while the morphology of the pore throats is characterized in that the radii of their associated void spaces are all greater than their own radii.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811265575.0A CN109242985B (en) | 2018-10-29 | 2018-10-29 | Method for determining key parameters of pore structure from three-dimensional image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811265575.0A CN109242985B (en) | 2018-10-29 | 2018-10-29 | Method for determining key parameters of pore structure from three-dimensional image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109242985A CN109242985A (en) | 2019-01-18 |
| CN109242985B true CN109242985B (en) | 2020-06-05 |
Family
ID=65078673
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811265575.0A Active CN109242985B (en) | 2018-10-29 | 2018-10-29 | Method for determining key parameters of pore structure from three-dimensional image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109242985B (en) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109993786B (en) * | 2019-03-08 | 2021-05-18 | 中国石油大学(北京) | Tortuosity acquisition method, device, equipment and storage medium |
| CN110018108B (en) * | 2019-05-14 | 2020-04-28 | 中国科学院地质与地球物理研究所 | Method and equipment for measuring pore diameter of rock pore |
| CN110992483B (en) * | 2019-11-19 | 2024-04-09 | 中国石油大学(华东) | Method for printing real three-dimensional fracture-cavity type oil reservoir physical model based on reverse modeling |
| CN112082917B (en) * | 2020-08-03 | 2021-04-30 | 西南石油大学 | Gas-water unsteady two-phase seepage simulation method based on dynamic network simulation |
| CN112967238A (en) * | 2021-02-23 | 2021-06-15 | 广东工业大学 | Porous medium permeability prediction method, electronic device and storage medium |
| CN113405965B (en) * | 2021-06-08 | 2022-06-28 | 浙江广天构件集团股份有限公司 | Method for analyzing pore connectivity of cement-based material particle stacking system |
| CN116337708B (en) * | 2022-12-06 | 2023-11-17 | 武汉理工大学 | Method and system for extracting porous material pore characteristics and fluid distribution |
Family Cites Families (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102087676B (en) * | 2010-12-13 | 2012-07-04 | 上海大学 | Pore network model (PNM)-based bionic bone scaffold designing method |
| US9070049B2 (en) * | 2013-03-15 | 2015-06-30 | Bp Corporation North America Inc. | Systems and methods for improving direct numerical simulation of material properties from rock samples and determining uncertainty in the material properties |
| CN103822865B (en) * | 2014-03-20 | 2016-05-04 | 中国石油大学(华东) | A kind of high-resolution three-dimension digital cores modeling method |
| CN104794709B (en) * | 2015-04-10 | 2018-03-09 | 四川大学 | A kind of dividing method of three-dimensional core image hole and venturi |
| CN106127816A (en) * | 2016-03-08 | 2016-11-16 | 中国石油大学(华东) | A kind of shale matrix reservoirs interstitial space characterizing method |
| CN106127777B (en) * | 2016-06-27 | 2017-08-29 | 中山大学 | A kind of three dimensions crack separation identification and characterizing method |
| CN107817199B (en) * | 2016-09-14 | 2019-09-03 | 中国石油化工股份有限公司 | A kind of construction method of tight sand multi-scale porosity model and application |
| CN106546521B (en) * | 2016-10-12 | 2019-02-19 | 北京师范大学 | A method of soil macropore spacial framework is quantified based on CT scan technology |
| CN108389193B (en) * | 2018-02-12 | 2021-03-16 | 中国人民解放军第三军医大学 | Method and system for analyzing three-dimensional voxel model behind carotid artery stent |
| CN108491677A (en) * | 2018-07-04 | 2018-09-04 | 河海大学 | Pore character statistical method based on the micro pore model for improving maximum ball |
-
2018
- 2018-10-29 CN CN201811265575.0A patent/CN109242985B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109242985A (en) | 2019-01-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109242985B (en) | Method for determining key parameters of pore structure from three-dimensional image | |
| Delerue et al. | New algorithms in 3D image analysis and their application to the measurement of a spatialized pore size distribution in soils | |
| CN110084304B (en) | Target detection method based on synthetic data set | |
| CN113850825A (en) | Remote sensing image road segmentation method based on context information and multi-scale feature fusion | |
| CN109447994A (en) | In conjunction with the remote sensing image segmentation method of complete residual error and Fusion Features | |
| CN105261068A (en) | Micro-CT technology-based reservoir core three-dimensional entity model reconstruction method | |
| CN111899353A (en) | Three-dimensional scanning point cloud hole filling method based on generation countermeasure network | |
| CN107977992A (en) | A kind of building change detecting method and device based on unmanned plane laser radar | |
| CN105654486A (en) | Complex reservoir rock pore structure parameter extraction method | |
| CN109345625B (en) | Rock core image self-adaptive partition three-dimensional reconstruction method | |
| CN107709699A (en) | Generate the three-dimensional micro model of porous rock sample | |
| CN108122221A (en) | The dividing method and device of diffusion-weighted imaging image midbrain ischemic area | |
| JP2004334640A (en) | Method for identifying multimedia data and its program | |
| CN112634429B (en) | Rock core three-dimensional image reconstruction method based on mixed depth generation model | |
| CN111028335B (en) | Point cloud data block surface patch reconstruction method based on deep learning | |
| CN112199886A (en) | Processing method of PRB data deep learning geological map prediction model | |
| CN112818777B (en) | Remote sensing image target detection method based on dense connection and feature enhancement | |
| CN111833432B (en) | Three-dimensional reconstruction method based on core two-dimensional gray scale image | |
| CN115546649B (en) | Single-view remote sensing image height estimation and semantic segmentation multi-task prediction method | |
| CN106226831A (en) | A kind of porosity inverse matrix method of rock core | |
| CN114494637A (en) | Sandstone three-dimensional real model reconstruction method based on structural body matrix | |
| CN114612315A (en) | High-resolution image missing region reconstruction method based on multi-task learning | |
| CN111986078A (en) | Multi-scale core CT image fusion reconstruction method based on guide data | |
| CN113447419B (en) | Dividing system for porous medium communicated pore structure unit | |
| CN111612869A (en) | Method for analyzing geological mapping based on grid data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |