CN104794709B - A kind of dividing method of three-dimensional core image hole and venturi - Google Patents
A kind of dividing method of three-dimensional core image hole and venturi Download PDFInfo
- Publication number
- CN104794709B CN104794709B CN201510170889.2A CN201510170889A CN104794709B CN 104794709 B CN104794709 B CN 104794709B CN 201510170889 A CN201510170889 A CN 201510170889A CN 104794709 B CN104794709 B CN 104794709B
- Authority
- CN
- China
- Prior art keywords
- point
- region
- hole
- skeleton
- dimensional
- 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
Abstract
The invention discloses a kind of three-dimensional core image hole and the dividing method of venturi, belong to technical field of image segmentation.This method is that connection mark is first done to three-dimensional core model；Three-dimensional distance conversion is done to each connected region again；Result after conversion is sought into gradient modulus value；Then look for key point；Skeleton generation is carried out to the key point found；The skeleton that each connected region generates is added to obtain whole three-dimensional core model skeleton C again；Seek the aperture cross section area of each skeletal point in skeleton C；Hole seed point is found again；Region growing is carried out to hole seed point, so as to mark off each hole of three-dimensional core model, the position that contacted with each other between hole is venturi；So just complete the segmentation of three-dimensional core model hole and venturi.The method of the present invention can reduce conventional segmentation methods some errors are big and phenomena such as over-segmentation and less divided, make three-dimensional core image pore throat segmentation result more accurate and effective.
Description
Technical field
The present invention relates to a kind of image Segmentation Technology, more particularly to a kind of segmentation of the hole and venturi of three-dimensional core image
Method, belong to technical field of image segmentation.
Background technology
Reservoir typically refer to can Reservoir Fluid, and can the rock stratum that is percolated of convection body.It is micro- inside reservoir
View hole gap structure directly affects the ability and seepage characteristic of reservoir storage oil gas, and this micropore structure simultaneously finally influences stone
Oil, natural gas etc. preserve content.So-called pore structure refer to rock core mesopore and venturi (abbreviation pore throat) geometric properties,
Distribution situation and the interconnection between them；Its mesopore refers to the space not being filled in rock, and it reflects rock
The ability of the storage oil gas of stone.And venturi is exactly to connect the long and narrow position between rock adjacent pores.And form, the size of venturi
The seepage characteristic of hole is reflected again.The characteristic of research rock core mesopore and venturi is predicting oil/gas formation, the regularity of distribution with carrying
The main method of the high rate of oil and gas recovery.And split hole and venturi, then it is the geometry for studying rock core mesopore and venturi
The premise of the characteristics such as feature, distribution situation.
It is complicated due to hole, it is traditional larger error to be present to three-dimensional core image pore throat dividing method.Such as
Over-segmentation phenomenon be present in existing watershed algorithm；Minimax ball algorithm can not accurately determine the position of venturi and can not be accurate
Really reflect the true radius of long and narrow hole；And there is less divided phenomenon in the algorithm based on skeleton segment.In recent years, researcher
In constantly research and discovery to three-dimensional core image pore throat segmentation problem, but current research conditions can also without proposition one kind
To carry out efficient and accurate dividing method to three-dimensional core image.Based on this, it is desirable to be able to which there have to be a kind of efficient, accurately right
The dividing method of three-dimensional core image hole and venturi, this is exactly the task place of the present invention.
The content of the invention
The purpose of the present invention is exactly the defects of being to overcome in the presence of prior art and deficiency, there is provided a kind of three-dimensional rock core
The dividing method of image hole and venturi, the dividing method are first to build three-dimensional core model, connection mark are done, to each connection
The steps such as three-dimensional distance conversion are done in region, are eventually found the seed point of hole, and the seed point can be understood as the centre bit of hole
Put, region growing is carried out according to the growing strategy of hole seed point region growing, so as to mark off each hole；Phase between hole
The position mutually contacted i.e. venturi.The method of the present invention can reduce conventional segmentation methods some errors are big and over-segmentation and owe
Phenomena such as segmentation, make three-dimensional core image pore throat segmentation result more accurate and effective.
The dividing method of a kind of three-dimensional core image hole and venturi provided by the invention, including following operating procedure：
Step 1：Prepare one group of rock core binary sequence figure to be split, three-dimensional core model is thus built, to three-dimensional rock core
Model does connection mark；To each connected region A of connection marki2 operated according to the following steps to step 5, wherein i=1,
2nd, 3 ... n, n are the number of connected region；
Step 2：To each connected region A in step 1iDo three-dimensional distance conversion；
Step 3：The result after three-dimensional distance conversion is done to step 2 and seeks gradient modulus value；
Step 4：The three-dimensional distance transformed value and step 3 gained gradient modulus value obtained according to step 2 finds key point；
Step 5：The key point found according to step 4 carries out skeleton generation；
Step 6：Each connected region A that step 5 is generatediSkeleton be added to obtain the skeleton of whole three-dimensional core model
C；
Step 7：All points in obtained skeleton C are added to step 6 and carry out aperture cross section area calculating；
Step 8：The aperture cross section area calculated according to step 7 carries out hole seed point searching；
Step 9：Region is carried out according to the growing strategy of hole seed point region growing with the hole seed point that step 8 is found
Growth, so as to mark off each hole of three-dimensional core model, the position to be contacted with each other between hole is venturi；Complete three
Tie up the segmentation of core model hole and venturi.
In above-mentioned technical proposal, gradient modulus value is asked to the result of three-dimensional distance conversion described in step 3, using following methods
Calculate：
To three-dimensional core model x, tri- directions of y, z calculate dx respectively using SOBEL operators, tri- components of dy, dz,
SOBEL operators are shown in that Rafael C.Gonzalez that Electronic Industry Press publishes etc. write, what Ruan Qiuqi etc. was translated《At digital picture
Reason》The third edition；
If current point is p (x, y, z), the distance value of current point is d (x, y, z)；Its three components represent as follows respectively：
Dx=2 × d (x-1, y, z)+d (x-1, y-1, z)+d (x-1, y+1, z) -2 × d (x+1, y, z)
- d (x+1, y-1, z)-d (x+1, y+1, z) (1)
Dy=2 × d (x, y-1, z)+d (x-1, y-1, z)+d (x+1, y-1, z) -2 × d (x, y+1, z)
- d (x-1, y+1, z)-d (x+1, y+1, z) (2)
Dz=2 × d (x, y, z-1)+d (x, y-1, z-1)+d (x, y+1, z-1) -2 × d (x, y, z+1)
- d (x, y-1, z+1)-d (x, y+1, z+1) (3)
Required gradient modulus value is obtained by below equation (4),
Dx, dy, dz are respectively the partial derivative in three directions of x, y, z in formula, and D (x, y, z) is current point p's (x, y, z)
The gradient modulus value of distance value.
In above-mentioned technical proposal, key point is found according to three-dimensional distance transformed value and gradient modulus value described in step 4, used
Following methods are found：
(1) local maximum point for the three-dimensional distance transformed value that step 2 obtains is found, the local maximum point refers to：This point
Range conversion value be more than in local neighborhood centered on the point three-dimensional distance transformed value a little；
(2) connection mark, the local maximum point in same connected region are done to three-dimensional distance transformed value local maximum point
Middle reservation gradient modulus value minimum point is as key point.
In above-mentioned technical proposal, the key point found described in step 5 carries out skeleton generation, is generated using following methods：
Each connected region AiAll key points in three-dimensional distance transformed value maximum point as source point, other key points
As point of destination, skeleton is generated using dijkstra critical path method (CPM)s；Converted in dijkstra critical path method (CPM)s with three-dimensional distance
The gradient modulus value of value cube is used as weights.
It is described to having a calculated hole diameters cross section area in obtained skeleton C described in step 7 in above-mentioned technical proposal
Aperture cross section area is the aperture section of skeletal point and the intersecting area of hole in skeleton C, is calculated using following methods：
(1) to each skeletal point on skeleton C, the method for aperture section is used as using tangential direction of the current point on skeleton
To, on known normal and plane a little in the case of, obtain aperture tangental equation；
(2) the aperture cross section area each put in skeleton C is calculated.
It is described to seek aperture tangental equation in above-mentioned technical proposal, seek method with normal：Current point is p on known skeletal point
(x0,y0,z0), former point is p1 (x1,y1,z1), latter point is p2 (x2,y2,z2), then normal direction is：
N (A, B, C)=normalize (p2-p1) (5)
(5) normalize is normalization operator in formula；Then aperture tangental equation is：
Ax+Ey+Cz+D=0 (6)
(6) A, B, C are three components of the normal vector tried to achieve above in formula, D=-A × x0-B×y0-C×z0。
Described to calculate the aperture cross section area each put in skeleton C in above-mentioned technical proposal, step is as follows：
1. with current skeletal point p (x0,y0,z0) it is that seed point does region growing,
I, by current skeletal point p (x0,y0,z0) add the growth point set of the skeletal point；
All 26 neighborhood points of the institute of the last growth of II, traversal a little, if satisfaction：(a) neighborhood point is porosity points, i.e., non-
Rock point.(b) neighborhood point is less than 1.0 to the distance of aperture section, then neighborhood point can be grown, and neighborhood point is included in into the skeletal point
Growth point set；Otherwise neighborhood point is abandoned；During the first secondary growth, it is described " last time growth institute a little it is all 26 neighbours
Domain point " is current skeletal point p (x0,y0,z0) 26 neighborhood points；
III, repeat step II, until 26 neighborhood points of all growing points of the skeletal point are all unsatisfactory in step II
(a) when with (b) condition, the region growing of the skeletal point terminates；
After 2. region growing terminates, statistics, which has grown, counts out, as current skeletal point p (x0,y0,z0) aperture section
Area Si；By identical computational methods obtain in skeleton C aperture cross section area a little, wherein i=1,2,3 ... M, M are
The number put on skeleton C.
In above-mentioned technical proposal, to the searching of hole seed point described in step 8, found using following methods：
(1) the aperture cross section area of all skeletal points in skeleton C is counted, calculates average equivalent radius of circle R；
(2) all skeletal points in skeleton C are traveled through, to each skeletal point, calculate the radius of equivalent circle r of the skeletal pointi；Such as
Fruit ri<λ R then skip the point, otherwise using the point as the search center of circle, with riFor search radius；If the aperture section of current skeletal point
Area is more than the aperture cross section area of all skeletal points in current search scope, then current skeletal point is added hole seed point collection
Close；Untill having traveled through all skeletal points, so as to find all hole seed points；Wherein 0<λ<1, general λ takes 1/2,1/2
It is an empirical value；The M of i=1,2,3 ..., M are the number put on skeleton C.
In above-mentioned technical proposal, region growing is carried out to hole seed point described in step 9, carried out using following methods：
(1) hole seed point is as a region BiInitial point, i.e. each hole seed point forms a region
Bi, hole seed point is added into region BiGrowth point set, the L of wherein i=1,2,3 ..., L be hole seed point number；
The region BiFor region AiSubset, i.e., each region AiOne or more hole seed point may be found out；
(2) to each region BiThe point of last time growth, statistical average three-dimensional distance transformed value, initial average distance become
Change the distance value that value is each hole seed point；
(3) with each growth district BiAverage three-dimensional distance transformed value as sort by, sort from big to small；
(4) region growing velocity stages is carried out according to ranking results in step (3)；
(5) D takes a region B in orderi, region growing is carried out according to region growing speed, order D is any, but with
All grown afterwards by this order D, i.e., after regional has sequentially grown first time respectively, first region can just start second
Growth；
(6) region BiAfter the completion of growth, judge whether that all porosity points have all been divided into regional, if institute
There are porosity points to be all divided into regional, then segmentation terminates；Otherwise, judge this region whether be order D last,
If it is step (2) is carried out；Otherwise, step (5) is carried out.
Advantage possessed compared with prior art of the invention and beneficial technique effect are as follows：
The dividing method of three-dimensional core image hole and venturi of the present invention, with the rock core binary sequence figure structure of reality
Three-dimensional core model is built, its hole and venturi are split using method of the present invention, and segmentation result figure is carried out
Observation, so as to demonstrate the reliability of method therefor of the present invention and practicality.The hole and venturi of the three-dimensional core image of the present invention
Segmentation, the error for reducing conventional segmentation methods is larger and the phenomenon of over-segmentation and less divided, its segmentation result are very accurate
Really.
Brief description of the drawings
Fig. 1 is one of rock core binary sequence figure in the embodiment of the present invention；
Fig. 2 is three-dimensional core model in the embodiment of the present invention；
Fig. 3-1 is a rock core bianry image in the embodiment of the present invention；
Fig. 3-2 is the result for carrying out range conversion in the embodiment of the present invention to Fig. 3-1；
Fig. 3-3 is the result that gradient modulus value is sought Fig. 3-2 in the embodiment of the present invention；
Fig. 3-4 is that the result that Fig. 3-2 and Fig. 3-3 finds key point is combined in the embodiment of the present invention；
Fig. 4 is the design sketch of the key point that three-dimensional core model is found in the embodiment of the present invention；
Fig. 5 is the design sketch of the skeleton C of three-dimensional core model generation in the embodiment of the present invention；
Fig. 6 is the schematic diagram of the aperture section of three-dimensional core model skeletal point in the embodiment of the present invention；Wherein, in cylinder
Have a guilty conscience line --- it is skeleton line, cylinder is hole, and the area of shadow surface is the aperture section of skeletal point and intersecting for hole
Area, the cylinder for not drawing skeleton line are adjacent hole.
Fig. 7 is the schematic diagram of three-dimensional core model hole seed point in the embodiment of the present invention；
Fig. 8 is the three-dimensional final segmentation result of core model in the embodiment of the present invention；
Fig. 9-1 is a side of the three-dimensional final segmentation result of core model in the embodiment of the present invention；
Fig. 9-2 is the picture of Fig. 9-1 partial enlargements in the embodiment of the present invention；
Fig. 9-3 is another side of the hole that Fig. 9-1 is indicated in the embodiment of the present invention；
Fig. 9-4 is the partial enlargement picture of one of three venturis of hole that Fig. 9-1 is indicated in the embodiment of the present invention；
Fig. 9-5 is the partial enlargement picture of one of three venturis of hole that Fig. 9-1 is indicated in the embodiment of the present invention；
Fig. 9-6 is the partial enlargement picture of one of three venturis of hole that Fig. 9-1 is indicated in the embodiment of the present invention；
Figure 10 is method flow diagram of the present invention.
Embodiment
The present invention is described in further detail below by specific embodiment and with reference to accompanying drawing, it is necessary to it is pointed out here that
It is that the embodiment is only intended to further describe the present invention, but is not meant to be to any of the scope of the present invention
Limit.
Embodiment：In order that dividing method of the present invention easily facilitates understanding and close to true application, below from rock
Heart bianry image starts to three-dimensional core model to segment as only, whole operation flow is illustrated, concrete operation step
It is as follows：
(1) one group of rock core binary sequence figure to be split of preparation, picture number 128, form BMP, size 128 ×
128, as shown in Figure 1；By the three-dimensional core model of rock core binary sequence figure structure as shown in Fig. 2 grey parts are hole in figure,
Hyalomere in white cube frame in addition to hole is divided into rock, i.e. rock is not shown；
(2) connection mark is done to three-dimensional core model, to each connected region A of connection marki(n of i=1,2,3 ...,
N is the number of connected region) three-dimensional distance conversion is done respectively；Because three-dimensional distance transformation results display effect is not in three dimensions
It is too obvious, the effect that key point is found by range conversion can be preferably showed using two dimensional image, therefore take a rock core two
It is worth image, as shown in figure 3-1, two-dimensional distance is done to it and converts obtained result as shown in figure 3-2；Brighter local distance in figure
Transformed value is bigger；All operations done to Fig. 3-1 to Fig. 3-4 are intended merely to show to find key point by range conversion
Effect, it is unrelated with the practical operation of the embodiment of the present invention；
(3) to each connected region AiThe result for doing three-dimensional distance conversion seeks gradient modulus value respectively；It is logical preferably to show
The effect that key point is found in range conversion is crossed, gradient modulus value is sought Fig. 3-2, obtains result as shown in Fig. 3-3；
The described pair of result done after three-dimensional distance conversion seeks gradient modulus value, is calculated using following methods：
To three-dimensional core model x, tri- directions of y, z calculate dx respectively using SOBEL operators, tri- components of dy, dz,
SOBEL operators are shown in that Rafael C.Gonzalez that Electronic Industry Press publishes etc. write, what Ruan Qiuqi etc. was translated《At digital picture
Reason》The third edition；If current point is p (x, y, z), the distance value of current point is d (x, y, z)；Three component distances represent such as respectively
Under：
Dx=2 × d (x-1, y, z)+d (x-1, y-1, z)+d (x-1, y+1, z) -2 × d (x+1, y, z)
- d (x+1, y-1, z)-d (x+1, y+1, z) (1)
Dy=2 × d (x, y-1, z)+d (x-1, y-1, z)+d (x+1, y-1, z) -2 × d (x, y+1, z)
- d (x-1, y+1, z)-d (x+1, y+1, z) (2)
Dz=2 × d (x, y, z-1)+d (x, y-1, z-1)+d (x, y+1, z-1) -2 × d (x, y, z+1)
- d (x, y-1, z+1)-d (x, y+1, z+1) (3)
Required gradient modulus value is obtained by below equation (4),
Dx, dy, dz are respectively the partial derivative in three directions of x, y, z in formula, and D (x, y, z) is current point p's (x, y, z)
The gradient modulus value of distance value.
(4) each connected region A is combinediThree-dimensional distance transformation results and gradient modulus value find key point respectively, obtain
As a result as shown in figure 4, wherein translucent grey parts are whole rock core hole, bead is key point, is shown in Fig. 4 whole
The key point of three-dimensional core model；In order to will not can not visually be seen because key point is too small, by key point using radius as 2
The form of the bead of individual pixel, true key point are the centre of sphere of bead in figure；Found preferably to show by range conversion
The effect of key point, key point is found with reference to Fig. 3-2 and Fig. 3-3, obtains result such as Fig. 3-4；
It is described that key point is found according to three-dimensional distance transformed value and gradient modulus value, found using following methods：
1. finding the local maximum point for the three-dimensional distance transformed value that step (2) obtains, the local maximum point refers to：This point
Range conversion value be more than the three-dimensional distance transformed value of institute a little in local neighborhood centered on the point, local neighborhood is using standing
Cube frame, for example length is used all for the cube frame of 5 pixels；When three-dimensional core model outer surface edge very not
Smoothly, such as there is situations such as zigzag, the size of cube frame can be adjusted, and such as use length all for 7 pixels
The cube of point；The outer surface of three-dimensional core model is more coarse, and cube frame is bigger；
2. connection mark is done to three-dimensional distance transformed value local maximum point, in the local maximum point in same connected region
Retain gradient modulus value minimum point as key point；
(5) each connected region AiAll key points in three-dimensional distance transformed value maximum point as source point, other passes
Key point generates skeleton as point of destination using dijkstra critical path method (CPM)s；Three-dimensional distance is used in dijkstra critical path method (CPM)s
The gradient modulus value of transformed value cube is used as weights.The skeleton that each connected region obtains is added to obtain three-dimensional core model
Skeleton C；As a result as shown in figure 5, curved lines are skeleton in Fig. 5, bead is key point；
(6) all points obtained in (5) on skeleton C are carried out calculating its aperture cross section area, a section schematic diagram is such as
Shown in Fig. 6, it is emphasized that the aperture section of any on skeleton C is determined by the tangent line of point and this on skeleton
's；
On the skeleton C aperture cross section area a little be calculated as follows：
(I) to each point on skeleton C, using tangential direction of the current point on skeleton as the normal direction of aperture section,
On known normal and plane a little in the case of, obtain aperture tangental equation；
Normal seeks method：Current point is p (x on known skeletal point0,y0,z0), former point is p1 (x1,y1,z1), latter point is
p2(x2,y2,z2), then normal direction is：
N (A, B, C)=normalize (p2-p1) (5)
(5) normalize is normalization operator in formula；Then aperture tangental equation is：
Ax+By+Cz+D=0 (6)
(6) A, B, C are three components of normal vector tried to achieve above in formula, D=-A × x0-B×y0-C×z0；
(II) each aperture section of skeletal point and the intersecting area of hole in skeleton C are calculated：
1. with current skeletal point p (x0,y0,z0) it is that seed point does region growing, step is as follows：
I, by current skeletal point p (x0,y0,z0) add the growth point set of the skeletal point；
All 26 neighborhood points of the institute of the last growth of II, traversal a little, if satisfaction：(a) neighborhood point is porosity points, i.e., non-
Rock point (b) neighborhood point is less than 1.0 to the distance of aperture section, then neighborhood point can be grown, and neighborhood point is included in into the skeletal point
Growth point set；Otherwise neighborhood point is abandoned；During the first secondary growth, it is described " last time growth institute a little it is all 26 neighbours
Domain point " is current skeletal point p (x0,y0,z0) 26 neighborhood points；
III, repeat step II, until 26 neighborhood points of all growing points of the skeletal point are all unsatisfactory in step II
(a) when with (b) condition, the region growing of the point terminates；
After 2. region growing terminates, statistics, which has grown, counts out, as current skeletal point p (x0,y0,z0) aperture section
Area Si；By identical computational methods obtain in skeleton C aperture cross section area a little；M, M are for wherein i=1,2,3 ...
The number put on skeleton C；
(7) the hole seed point of three-dimensional core model is found；
The searching to hole seed point, found using following methods：
1. counting the cross section area of all skeletal points in skeleton C, average equivalent radius of circle R is calculated；
2. traveling through all skeletal points in skeleton C, to each skeletal point, the radius of equivalent circle r of the skeletal point is calculatedi；Such as
Fruit ri<λ R (wherein 0<λ<It is an empirical value that 1, general λ, which take 1/2,1/2) point is then skipped, otherwise justify by search of the point
The heart, with riFor search radius；If the cross section area of current skeletal point is more than the section of all skeletal points in current search scope
Area, then current skeletal point is added hole seed point set；Untill having traveled through all skeletal points, so as to find all holes
Gap seed point.
Described the step of finding method is carried out to hole seed point 1. middle average equivalent radius of circle R calculating using following
Method：
First calculate average pore size cross section area S：
Wherein S is average pore size cross section area, and M is the number put on skeleton C, SiCut in aperture each to be put on skeleton C
Face area；
Average equivalent radius of circle R is calculated by S again：
Wherein R is average equivalent radius of circle, and S is average pore size cross section area, and π is pi.
In the radius of equivalent circle r of the step of carrying out finding method to hole seed point 2. middle skeleton pointiCalculating make
With following methods：
Wherein riFor the radius of equivalent circle of skeletal point, SiFor the aperture cross section area of the skeletal point, π is pi.Tied
For fruit as shown in fig. 7, translucent grey parts are whole rock core hole in figure, bead is hole seed point；In order to visually will not
It can not see, bead of the hole seed point using radius as 5 pixels be shown, true hole because hole seed point is too small
Seed point is the centre of sphere of bead in figure.
(8) region growing is carried out to obtained hole seed point；
Described pair of obtained hole seed point carries out region growing, is carried out using following methods：
1. a hole seed point is as a region BiInitial point, i.e. each hole seed point forms a region Bi,
Seed point is added into region BiGrowth point set, the L of wherein i=1,2,3 ..., L be hole seed point number；The area
Domain BiFor region AiSubset, i.e., each region AiOne or more hole seed point may be found out；
2. to each region BiThe point of last time growth, statistical average three-dimensional distance transformed value, it is initial it is average it is three-dimensional away from
From the distance value that transformed value is each hole seed point；
3. with each growth district BiAverage three-dimensional distance transformed value as sort by, sort from big to small；
4. according to step, 3. middle ranking results carry out region growing velocity stages；
Region growing velocity stages of the sorting position in preceding 1/3 region is FAST, if a total of 30 region BiRow
Sequence, i.e. i=1,2,3,, 30；Preceding 1/3 region of that sorting position refers to region of the sequence preceding 10；Cited below 1/3
To 4/5, last 1/5 by that analogy；If region BiIt is not integer that total number, which is multiplied by 1/3 or is multiplied by 4/5, is taken immediate whole
Number replaces.Sorting position is MEDIUM in 1/3 to 4/5 region growing velocity stages；Sorting position is given birth in last 1/5 region
Long velocity stages is SLOW.Speed is FAST's, and 26 neighborhoods are used when neighborhood point is determined whether in growing strategy；Speed is
MEDIUM's, when neighborhood point is determined whether in growing strategy, using 18 neighborhoods；Speed is SLOW's, using 6 neighborhoods.Growth
Velocity sorting section can be adjusted according to actual conditions, but ensure that the more forward speed of growth of sequence is faster.
5. D takes a region B in orderi, region growing is carried out according to region growing speed, order D is any, but after
All grown by this order D, i.e., after regional has sequentially grown first time respectively, it is secondary that first region can just start second
It is long；
A 6. region BiAfter the completion of growth, judge whether that all porosity points have all been divided into regional, if institute
There are porosity points to be all divided into regional, then segmentation terminates；Otherwise, judge this region whether be order D last,
If it is step is carried out 2.；Otherwise, step is carried out 5..
Described pair of obtained hole seed point carries out region growing the step of 5. in region growing is carried out to regional
Growing strategy it is as follows：
All neighborhood points of the institute of the last growth of traversal a little, according to the difference of the speed of growth, neighborhood may be 26 neighbours
Domain, 18 neighborhoods, 6 neighborhoods；If meet：(a) neighborhood point is porosity points, i.e., non-rock point (b) neighborhood point not by any region
The set of growing point includes；The neighborhood point is then included in the set of region growing point.During the first secondary growth, described " upper one
Secondary growth all neighborhood points a little " be seed neighborhood of a point point.
Segmentation each region B after terminatingiA hole is represented, the position to be contacted with each other between hole is venturi.
Obtain without display venturi in result such as Fig. 8, Fig. 8, that is, the venturi for connecting hole is set to background colour.Split in figure
Different zones B outiIndicated with different gray scales, each one hole of region representation, the position to contact with each other between hole
It is set to venturi.
Fig. 9-1 show a side of the final segmentation result of three-dimensional core model, and wherein white ovals sign is point
Cut a hole in result；It is the picture of Fig. 9-1 partial enlargements shown in Fig. 9-2, wherein oval sign is the three of the hole
Individual venturi；Another side for the hole that Fig. 9-1 is indicated shown in Fig. 9-3, wherein oval sign be the hole three larynxs
Road；And the hole is illustrate only in figure, remaining hole is not shown；It is Fig. 9-1 signs shown in Fig. 9-4, Fig. 9-5 and Fig. 9-6
Hole three venturis partial enlargement picture, white arrow pointed location is venturi position.From segmentation result
As can be seen that the dividing method of the present invention is more accurate to the ration of division of hole and venturi.All with unified ash in Fig. 9-1 to Fig. 9-6
Degree level shows venturi, and for the position of clearer performance venturi, the position cut-off that venturi is connected with hole is set into the back of the body
Scenery；Therefore see that hole and venturi are to disconnect in Fig. 9-1 to Fig. 9-6.
In the present embodiment, it is entered with the rock core binary sequence figure structure threedimensional model of reality and using the method for the present invention
Row segmentation.Segmentation result is observed, so as to demonstrate the reliability of dividing method of the present invention and practicality.
The above embodiment of the present invention is the preferred embodiments of the present invention, is not that technical scheme of the present invention is appointed
What is limited, as long as the technical scheme realized without creative work on the basis of above-described embodiment, is regarded as
Enter the present invention protect in the range of.
Claims (2)
1. the dividing method of a kind of three-dimensional core image hole and venturi, it is characterised in that comprise the following steps：
Step 1：Prepare one group of rock core binary sequence figure to be split, three-dimensional core model is thus built, to three-dimensional core model
Connection mark is done, to each connected region A of connection marki2 operated by the following step to step 5；Wherein, i=1,2,
3 ... n, n are the number of connected region；
Step 2：To each connected region A in step 1iDo three-dimensional distance conversion；
Step 3：The result after three-dimensional distance conversion is done to step 2 and seeks gradient modulus value；
Step 4：The three-dimensional distance transformed value and step 3 gained gradient modulus value obtained according to step 2 finds key point；
Step 5：The key point found according to step 4 carries out skeleton generation；
Step 6：The skeleton that step 5 is generated is added to obtain the skeleton C of whole three-dimensional core model；
Step 7：All points carry out aperture cross section area calculating in the skeleton C obtained to step 6；
Step 8：The aperture cross section area calculated according to step 7 carries out hole seed point searching；
Step 9：Region life is carried out according to the growing strategy of hole seed point region growing with the hole seed point that step 8 is found
It is long, so as to mark off each hole and venturi of three-dimensional core model；Complete the segmentation of three-dimensional core model hole and venturi；
The aperture cross section area calculated described in step 8 according to step 7 carries out hole seed point searching, uses following methods：
(1) the aperture cross section area of all skeletal points in skeleton C is counted, calculates average equivalent radius of circle R；
(2) all skeletal points in skeleton C are traveled through, to each skeletal point, calculate the radius of equivalent circle r of the skeletal pointi；If ri<
λ R then skip the point, otherwise using the point as the search center of circle, with riFor search radius；If the aperture cross section area of current skeletal point
More than the aperture cross section area of all skeletal points in current search scope, then current skeletal point is added hole seed point set；
Untill having traveled through all skeletal points, so as to find all hole seed points；Wherein 0<λ<1, general λ, which take, 1/2,1/2 is
One empirical value；The M of i=1,2,3 ..., M are the number put on skeleton C.
2. according to the method for claim 1, it is characterised in that the hole seed point root found described in step 9 with step 8
Region growing is carried out according to the growing strategy of hole seed point region growing, uses following methods：
(1) hole seed point is as a region BiInitial point, by hole seed point add region BiGrowth point set
Close, the L of wherein i=1,2,3 ..., L are the number of hole seed point；The region BiFor region AiSubset, i.e., each region Ai
One or more hole seed point may be found out；
(2) to each region BiThe point of last time growth, statistical average three-dimensional distance transformed value, initial average distance transformed value
The distance value of as each hole seed point；
(3) with each growth district BiAverage three-dimensional distance transformed value as sort by, sort from big to small；
(4) region growing velocity stages is carried out according to ranking results in step (3)；
(5) D takes a region B in orderi, region growing is carried out according to region growing speed step, order D is any, but after
All grown by this order D, i.e., after regional has sequentially grown first time respectively, it is secondary that first region can just start second
It is long；
(6) region BiAfter the completion of growth, judge whether that all porosity points have all been divided into regional, if all holes
Gap point has all been divided into regional, then segmentation terminates；Otherwise, judge this region whether be order D last, if
It is then to carry out step (2)；Otherwise, step (5) is carried out.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510170889.2A CN104794709B (en) | 2015-04-10 | 2015-04-10 | A kind of dividing method of three-dimensional core image hole and venturi |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510170889.2A CN104794709B (en) | 2015-04-10 | 2015-04-10 | A kind of dividing method of three-dimensional core image hole and venturi |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104794709A CN104794709A (en) | 2015-07-22 |
| CN104794709B true CN104794709B (en) | 2018-03-09 |
Family
ID=53559489
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510170889.2A Active CN104794709B (en) | 2015-04-10 | 2015-04-10 | A kind of dividing method of three-dimensional core image hole and venturi |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104794709B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105427383B (en) * | 2015-11-23 | 2017-04-05 | 中国石油大学(华东) | A kind of pore throat cross-sectional configuration method of the blowhole network model for considering concavity and convexity |
| CN105551018B (en) * | 2015-12-03 | 2018-04-20 | 长江大学 | Towards each component extracting method of carbonate rock complicated reservoirs digital cores |
| CN106383078B (en) * | 2016-09-20 | 2019-05-07 | 中国石油天然气股份有限公司 | The determination method and device of rock waterflood efficiency |
| CN106803253B (en) * | 2017-01-17 | 2019-06-04 | 四川大学 | A kind of three-dimensional rock image crack identification method |
| CN107993261B (en) * | 2017-11-02 | 2019-02-12 | 中国科学院地质与地球物理研究所 | A kind of hole based on three-dimensional Core Scanning Image and pore throat recognition methods |
| CN110246138B (en) * | 2018-03-09 | 2021-06-15 | 中国石油化工股份有限公司 | Method for segmenting pore throat of digital core image |
| CN109087301B (en) * | 2018-06-29 | 2021-09-21 | 中国石油大学(华东) | Core throat segmentation method based on centering axis and surface model |
| CN108491677A (en) * | 2018-07-04 | 2018-09-04 | 河海大学 | Pore character statistical method based on the micro pore model for improving maximum ball |
| CN109242985B (en) * | 2018-10-29 | 2020-06-05 | 中国科学院力学研究所 | Method for determining key parameters of pore structure from three-dimensional image |
| CN110018108B (en) * | 2019-05-14 | 2020-04-28 | 中国科学院地质与地球物理研究所 | Method and equipment for measuring pore diameter of rock pore |
| CN113139302B (en) * | 2021-05-20 | 2023-03-10 | 电子科技大学 | Area growth-based solution breaking identification method |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101556703A (en) * | 2009-05-16 | 2009-10-14 | 中国石油大学(华东) | Method for establishing network model based on serial section image |
| CN102222359A (en) * | 2011-05-24 | 2011-10-19 | 中国石油天然气股份有限公司 | Method for remodeling three-dimensional pore structure of core |
| CN104376556A (en) * | 2014-10-31 | 2015-02-25 | 四川大学 | Rock CT image target segmentation method |
-
2015
- 2015-04-10 CN CN201510170889.2A patent/CN104794709B/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101556703A (en) * | 2009-05-16 | 2009-10-14 | 中国石油大学(华东) | Method for establishing network model based on serial section image |
| CN102222359A (en) * | 2011-05-24 | 2011-10-19 | 中国石油天然气股份有限公司 | Method for remodeling three-dimensional pore structure of core |
| CN104376556A (en) * | 2014-10-31 | 2015-02-25 | 四川大学 | Rock CT image target segmentation method |
Non-Patent Citations (4)
| Title |
|---|
| Geometric and Topological Analysis of Three-Dimensional Porous Media:Pore Space Partitioning Based on Morphological Skeletonization;Z. Liang等;《Journal of Colloid and Interface Science》;20001231;第13-24页 * |
| 低渗砂岩储层孔隙结构表征及应用研究;朱洪林;《中国博士学位论文全文数据库工程科技I辑》;20150315(第3期);参见第2章 * |
| 储层微观孔喉网络图形识别方法;张婷等;《吉 林 大 学 学 报 （地 球 科 学 版）》;20110930;第４１卷(第５期);参见第1646-1650页 * |
| 应用VTK 技术进行三维图像显示及参数计算;杨晏春等;《计算机与数字工程》;20091130;第37卷(第11期);参见第4-5节 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104794709A (en) | 2015-07-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104794709B (en) | A kind of dividing method of three-dimensional core image hole and venturi | |
| CN107025685B (en) | Airborne building roof point cloud modeling method under topology perception | |
| CN103279983B (en) | The modeling method of China Tang dynasty style ancient building | |
| CN101650890A (en) | Method for expressing processing of roads on map | |
| CN102567492A (en) | Method for sea-land vector map data integration and fusion | |
| CN105303612B (en) | A kind of extract digital network method based on Triangulated irregular network model | |
| CN111858810B (en) | Modeling elevation point screening method for road DEM construction | |
| CN106097450A (en) | A kind of contour lines creation method and device | |
| CN102903139B (en) | Accelerated rendering method for contours | |
| CN109118588B (en) | Automatic color LOD model generation method based on block decomposition | |
| CN107240073A (en) | A kind of 3 d video images restorative procedure merged based on gradient with clustering | |
| CN106409129A (en) | Road condition drawing method and road condition drawing device | |
| CN105005580B (en) | A kind of method for showing reservoir landform and device thereof | |
| CN104318605A (en) | Parallel lamination rendering method of vector solid line and three-dimensional terrain | |
| CN106373175A (en) | Terrain height graph data loading method | |
| CN104794308B (en) | Domain image based on Image Edge-Detection is converted to CIF document methods | |
| CN101763356B (en) | Acquisition and organization method of topographic data in virtual map | |
| CN109213763A (en) | The organization and management method and system of Vehicle-borne Laser Scanning point cloud | |
| CN105468375A (en) | Surface structure light point cloud data oriented corresponding point search structure construction method | |
| CN104143214A (en) | Electronic map polygon triangulation method and device | |
| CN104036096B (en) | Method for mapping bump features on inclined face to manufacturing feature bodies | |
| CN107169080A (en) | A kind of geospatial analysis system being combined based on GIS and spatial database | |
| CN110260876A (en) | A kind of road model generation method and system based on oblique photograph and GIS technology | |
| CN103968781B (en) | High accuracy fast phase unwrapping method based on structure limit | |
| CN1273940C (en) | Fast drawing forest method of graded hierarchical assembling depth paste-up atlas |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |