CN106803253B - A kind of three-dimensional rock image crack identification method - Google Patents

A kind of three-dimensional rock image crack identification method Download PDF

Info

Publication number
CN106803253B
CN106803253B CN201710030343.6A CN201710030343A CN106803253B CN 106803253 B CN106803253 B CN 106803253B CN 201710030343 A CN201710030343 A CN 201710030343A CN 106803253 B CN106803253 B CN 106803253B
Authority
CN
China
Prior art keywords
triangle gridding
triangle
class
unit
crack
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
Application number
CN201710030343.6A
Other languages
Chinese (zh)
Other versions
CN106803253A (en
Inventor
滕奇志
夏晨木
何小海
张耀
吴晓红
卿粼波
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan University
Original Assignee
Sichuan University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Sichuan University filed Critical Sichuan University
Priority to CN201710030343.6A priority Critical patent/CN106803253B/en
Publication of CN106803253A publication Critical patent/CN106803253A/en
Application granted granted Critical
Publication of CN106803253B publication Critical patent/CN106803253B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0004—Industrial image inspection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G06F18/232—Non-hierarchical techniques
    • G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • G06F18/23213—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions with fixed number of clusters, e.g. K-means clustering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/00—3D [Three Dimensional] image rendering
    • G06T15/005—General purpose rendering architectures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20036—Morphological image processing

Abstract

The invention discloses a kind of three-dimensional rock image crack identification methods.This method takes turns doing resurfacing, Laplce's Mesh Smoothing and lattice simplified operation to the three-dimensional each connected component of blowhole model；It is then based on triangle gridding unit triangle area and unit normal vector feature, calculates cluster centre using K- mean cluster；Triangle gridding unit is classified again according to cluster centre；Hole repair is executed to the three-dimensional surface structure that the triangle gridding unit in each triangle gridding class is constituted, to obtain sealing surface structured set；It whether is crack according to the current sealing surface structured set of the parameter decisions such as form factor, area ratio；Morphological dilations processing is executed to sealing surface structured set surface and voxel of object point set and carries out logical AND operation with the voxel point set of initial three-dimensional blowhole model connected component, the voxel point set obtained with operation is crack.It is compared with the traditional method, the present invention can be identified correctly and accurately extract rock fracture.

Description

A kind of three-dimensional rock image crack identification method
Technical field
The present invention relates to a kind of image segmentation and identification technology more particularly to a kind of crack identification sides of three-dimensional rock image Method belongs to image segmentation and identification technology field.
Technical background
Rock generates mechanical destruction under stress, without the rift structure being obviously displaced crack.Crack is a kind of Special hole.In field of petroleum geology, crack identification is oil reservoir and the important content for storing heterogeneous research, to oil and gas development Raising with oil recovery factor is of great significance.
The present invention is directed the crack identification of rock CT sequence image.The present invention needs first to every in rock CT sequence chart Width image carries out binaryzation, recycles the sequence chart of binaryzation to construct three-dimensional blowhole model, and identify in three-dimensional space Rock fracture out.It is one group of continuous two-dimensional sequence figure via the rock image that CT is obtained, is divided using conventional two-dimensional image and calculated Method successively can make crack extract to width CT image each in sequence chart, and then can obtain the 3-D image in crack.But above-mentioned side Method has ignored related information of the rock fracture between multiple continuous rock CT images merely with two-dimensional image information.This will The common hole that CT image is in long and narrow linear and three-dimensional space is in non-fracture pattern is caused easily to be mistaken for crack or CT image In non-long and narrow linear, three-dimensional space is easily mistaken for non-crack in the crack of fracture pattern.
The current existing crack identification method based on three-dimensional space is according to the minimum external of each isolated three-dimensional pore space The radius of a ball and the ratio of equivalent spheres radius, form factor, the ginseng such as ratio of approximate minimum circumscribed rectangular body longest edge and most short side Number is to identify rock fracture.But since crack is not usually self-existent, it is associated between crack and common hole closely, so The misclassification rate of method fracture is higher, and the crack being connected with common hole is such as judged to non-crack, or general by what is be connected with crack Through-hole gap is judged to crack.
In recent years, researcher is in constantly research and discovery three-dimensional rock image crack identification problem, but at present also not It is proposed that one kind can carry out efficiently and accurately to three-dimensional rock image crack and know method for distinguishing.It is existing that the present invention is based on this research Shape, propose it is a kind of can the efficiently and accurately method that identifies three-dimensional rock image crack.
Summary of the invention
It is an object of the invention to overcome defect and deficiency of the existing technology, a kind of three-dimensional rock image crack is provided Recognition methods.
It is a kind of special hole that three-dimensional rock image crack identification method provided by the invention, which is based on the thought crack, Gap is often connected with common hole, therefore cannot directly determine whether some hole is crack.The present invention is to three-dimensional blowhole mould Each connected component of type executes resurfacing, Laplce's Mesh Smoothing, lattice simplified, grid cell unit normal vector and calculates Equal sequence of operations.Since crack is there are planar feature, the probability being simplified positioned at the triangle gridding unit of fracture faces is high.Three Angle grid cell, which will be simplified rear mesh triangles area, to become larger, and grid units normal vector direction tends to be unified.According to grid three Edged surface product and grid units normal vector direction character execute K- mean cluster to triangle gridding unit, according to what is obtained after cluster Cluster centre repartitions triangle gridding class, and the triangle gridding unit positioned at same fracture faces will gather in same class.Into And it is whether the three-D space structure that can determine that each triangle gridding class is constituted has FRACTURE CHARACTERISTICS by parameters such as form factors. In order to improve the accuracy of crack identification, 3-D image morphological dilations behaviour is executed to the three-D space structure with FRACTURE CHARACTERISTICS Make, the voxel point set for the three-D space structure for later obtaining expansion and the voxel of initial three-dimensional blowhole model connected component Point set carries out logical AND operation, and the result with operation is rock fracture.
The recognition methods in three-dimensional rock image crack provided by the invention, including following operating procedure:
Pretreatment: its corresponding three-dimensional blowhole model is constructed to one group of rock binary sequence figure to be identified, and is mentioned All connected components of modulus type；Reject the wherein very few connected component of tissue points；Use marching cube isosurface extraction (Marching Cube) algorithm carries out resurfacing to each connected component, obtains the corresponding triangle gridding of each connected component Model Aq, wherein q=1,2,3 ... k, k are triangle grid model number；Laplce's grid is made to obtained triangle grid model Smoothing processing, and it is lattice simplified to triangle grid model progress using Mesh simplification algorithm；Then to each AqIt independently executes Following step 1 arrives the operation of step 10；
Step 1: calculating the triangle grid model A that pretreatment obtainsqIn each triangle gridding unit unit normal vector with Triangle area, and adjust unit normal vector direction；Followed by the unit normal vector and triangle area feature of triangle gridding unit To AqTriangle gridding unit carry out K- mean cluster；
Step 2: the cluster centre obtained according to step 1 is by AqTriangle gridding unit reclassify, classification obtain it is multiple Triangle gridding class；
Step 3: extracting the triangle gridding connected component of each triangle gridding class in step 2；
Step 4: rejecting the noise contribution in the triangle gridding connected component that step 3 obtains, recalculate triangle gridding class Cluster centre, and the adjacent all triangle gridding classes of the cluster centre recalculated are merged into one kind；Then to every A triangle gridding class independently executes the operation that step 5 arrives step 10；
Step 5: all triangle gridding units in each triangle gridding class that step 4 obtains form a three-dimensional surface knot Structure executes the operation of threedimensional model hole repair to this three-dimensional surface structure, one or more closed three-dimensionals can be obtained after operation Face structure, one sealing surface structured set of these three-dimensional surface structure compositions；
Step 6: calculating the minimum circumscribed rectangular body for the sealing surface structured set that step 5 obtains；
Step 7: calculating the form factor for the sealing surface structured set that step 5 obtains；It is external to calculate the minimum that step 6 obtains The form factor of cuboid；Reference area compares parameter；
Step 8: the parameters such as two form factors and area ratio for being obtained using step 7 determine current sealing surface structure collection Whether close is crack.If it is determined that result is very, then to continue to execute subsequent operation；If it is determined that result is false, then current sealing surface knot Structure set is not crack；
Step 9: determining that result is the tissue points on genuine sealing surface structured set surface and inside in extraction step 8, to mentioning The voxel point set taken makees the processing of 3-D image morphological dilations；
Step 10: by the voxel of the voxel point set expanded in step 9 and initial three-dimensional blowhole model connected component Point set carries out logical AND operation, and the voxel point set obtained after operation is crack.
In above-mentioned technical proposal, the very few connected component of rejecting tissue points described in pretreatment, the tissue points preferentially used Number threshold value is 100.
In above-mentioned technical proposal, Laplce's Mesh Smoothing that triangle grid model executes is handled described in pretreatment, Preferentially use following parameter setting: for Laplce's Mesh Smoothing index for 0.2, smooth number is 100.
In above-mentioned technical proposal, the processing of Triangular Mesh Simplification described in pretreatment preferentially uses mesh extraction (Decimation of Triangle Meshes) algorithm, the control of mesh extraction ratio remove 90% triangulation network 0.9 Lattice unit.
It is preferential to use to the adjustment in triangle gridding unit of cells normal vector direction described in step 1 in above-mentioned technical proposal Following adjustment mode:
If it is (x, y, z) that unit normal vector, which is calculated, if when z < 0, adjusting unit normal vector is (- x ,-y ,-z).
In above-mentioned technical proposal, to the K- mean cluster of triangle gridding unit described in step 1, cluster centre is calculated, it is excellent First calculated using following manner:
K- means clustering algorithm process is as follows:
(1) initial cluster center is determined；
(2) all triangle gridding units are assigned to away from the affiliated triangle gridding class of its nearest cluster centre；
(3) cluster centre is updated；
(4) (2) (3) step is executed repeatedly until meeting termination condition.
Relevant calculation formula is as follows in K- means clustering algorithm:
Initial cluster center V in process (1)iCalculation formula it is as follows, wherein i=1,2,3:
Wherein σ is AqThe unit normal vector of the maximum triangle gridding unit of intermediate cam area；ξ is AqIntermediate cam area is minimum Triangle gridding unit unit normal vector.
Enable ΔjFor AqIn j-th of triangle gridding unit, λjFor ΔjUnit normal vector, wherein [1, m] j ∈, m AqContain Triangle gridding unit number；For AqIn i-th of triangle gridding class, cluster centre Vi, wherein i=1,2,3；It is i-th A triangle gridding classIn j-th of triangle gridding unit,ForUnit normal vector,ForTriangle area, wherein j ∈ [1, n], n are i-th of triangle gridding classIntermediate cam grid cell number.
Process (2) intermediate cam grid cell ΔjWith each cluster centre ViDistance thetaijCalculation formula it is as follows:
Wherein
Cluster centre calculation formula is as follows in process (3):
Wherein
Add is defined as:
Wherein n is i-th of triangle gridding classIntermediate cam grid cell number.
Stopping criterion for iteration in process (4) are as follows: the number of iterations is more than 50 times or cluster centre amplitude of variation less than 0.001； Kth iteration cluster centre variance valuesCalculation formula it is as follows:
WhereinCluster centre at the end of for kth time iterationThree coordinate values.
In above-mentioned technical proposal, triangle gridding unit described in step 2 is reclassified, and is preferentially carried out using following methods Classification:
Triangle gridding unit ΔjWith cluster centre ViDistance be θij, and
θkj=min { θij|i∈[1,2,3]} (10)
If θkj< 0.59359877, then ΔjBelong to kth triangle gridding classOtherwise ΔjIt is not belonging to any kind.
In above-mentioned technical proposal, the triangle gridding connected component of triangle gridding class is extracted described in step 3, under preferential use State method:
If triangle gridding unit ΔiTriangle vertex and triangle gridding unit ΔjTriangle vertex have at least one vertex phase Together, then ΔiWith ΔjSpace is adjacent.Each triangle gridding classIn, spatially neighbouring relations extract triangle gridding connected componentWherein [1, k] t ∈, k are triangle gridding classesTriangle gridding connected component number.Triangle gridding classWith the triangulation network Lattice connected componentRelationship it is as follows:
In above-mentioned technical proposal, to the method for triangle gridding connected component noise eliminating described in step 4, under preferential use State method:
Triangle gridding classThe sum of triangle gridding unit triangle area beTriangle gridding connected componentTriangle The sum of grid cell triangle area isWherein [1, k] t ∈, k are triangle gridding classesTriangle gridding connected component number. Calculation formula is as follows:
WhereinFor triangle gridding unitTriangle area.
It is rightIn all triangle gridding connected componentsIf
Then by triangle gridding connected componentIt is considered as noise jamming, and from triangle gridding classMiddle rejecting.
In above-mentioned technical proposal, the method for the close triangle gridding class of merging described in step 4, preferentially using following methods into Row merges:
If cluster centre ViWith cluster centre VjDistance be χij, calculation formula is as follows:
Wherein
If χij< 0.7853981, then it is adjacent to regard the affiliated triangle gridding class of two cluster centres.It will be in three triangle gridding classes All adjacent classes merge into a class.Triangle gridding class after merging is still denoted asWherein 1≤i≤3.
In above-mentioned technical proposal, the minimum circumscribed rectangular body of sealing surface structured set described in step 6 is calculated, preferential to use Method described in following documents calculates: three-dimensional core image crack automatic recognition [journal article] Deng Zhiqiu, Teng Qizhi-" meter Calculation machine and digital engineering " 1 phase in 2013.
In above-mentioned technical proposal, the calculating of form factor described in step 7, area ratio preferentially uses following methods:
1. the form factor Φ of sealing surface structured set ψψCalculation formula is as follows:
2. the form factor Φ of minimum circumscribed rectangular bodily form ∏∏Calculation formula is as follows:
3. area ratio Λ calculation formula is as follows:
Wherein VψFor the volume of sealing surface structured set ψ；SψFor triangle gridding unit triangle areas all before ψ hole repair The sum of；V∏For the volume of sealing surface structured set minimum circumscribed rectangular body ∏；S∏For the surface area of ∏.
In above-mentioned technical proposal, the criterion in crack described in step 8 preferentially uses following criterion:
When following condition is set up:
Or when the following conditions are set up:
Wherein
Then determine that the current surface sealing surface structured set ψ and the voxel point set δ of inside belong to a part in crack, and executes Subsequent step；Otherwise it is determined as non-crack.
In above-mentioned technical proposal, 3-D image morphological dilation described in step 9 preferentially uses following methods:
Wherein it is located at (0,0,0), the cube having a size of 3*3 centered on B；Z is space arbitrary point (z1,z2,z3)
(B)z=c | c=b+z, b ∈ B } (25)
In above-mentioned technical proposal, the operation of logical AND described in step 10 preferentially uses following methods:
Crack voxel point set ζ calculation formula are as follows:
ζ=p | p ∈ Hq∩δ′} (26)
Wherein HqFor AqThe voxel point set of corresponding initial three-dimensional blowhole model connected component, p are any tissue points. δ ' is the voxel point set obtained after morphological dilations processing in step 9.
Compared with prior art, the present invention possessed advantage and beneficial technical effect are as follows:
Three-dimensional rock image crack identification method of the present invention constructs three-dimensional rock with actual rock binary sequence figure Stone pore model identifies its crack using method of the present invention, and observes recognition result, to verify The reliability and practicability of method therefor of the present invention.The present invention solves that crack identification in conventional method is incorrect and crack The problem of inaccuracy is extracted, the present invention can be identified correctly and accurately extract rock fracture.
The present invention is by state natural sciences fund " rock Micro Heterogeneous structure three-dimensional image reconstruction and increase resolution skill Art research (61372174) " is subsidized.
Detailed description of the invention
Fig. 1-1 is one of rock binary sequence figure in the embodiment of the present invention；
Fig. 1-2 is the three-dimensional blowhole model constructed in the embodiment of the present invention by rock binary sequence figure；
Fig. 2-1 is in the embodiment of the present invention to the three-dimensional blowhole model after model denoising in Fig. 1-2；
Fig. 2-2 is the result for carrying out resurfacing in the embodiment of the present invention to model in Fig. 2-1；
Fig. 3-1 is to execute the smoothed out partial enlargement picture of Laplce to model in Fig. 2-2 in the embodiment of the present invention；
Fig. 3-2 is to execute the smoothed out partial enlargement picture of Laplce to model in Fig. 2-2 in the embodiment of the present invention；
Fig. 4-1 is to execute the partial enlargement picture after mesh extraction in the embodiment of the present invention to model in Fig. 3-1；
Fig. 4-2 is to execute the partial enlargement picture after mesh extraction in the embodiment of the present invention to model in Fig. 3-1；
Fig. 5 is the partial schematic diagram that grid units normal vector is calculated model in Fig. 4-1 in the embodiment of the present invention；
Fig. 6-1 is in the embodiment of the present invention to reclassifying after the triangle gridding unit K- mean cluster of model in Fig. 5 Result figure；
Fig. 6-2 is in the embodiment of the present invention to reclassifying after the triangle gridding unit K- mean cluster of model in Fig. 5 Result figure；
Fig. 7-1 is to execute connected component to triangle gridding unit sorted in Fig. 6-1 in the embodiment of the present invention to mark Result figure；
Fig. 7-2 is to execute connected component to triangle gridding unit sorted in Fig. 6-1 in the embodiment of the present invention to mark Result figure；
Fig. 7-3 is to execute connected component to triangle gridding unit sorted in Fig. 6-1 in the embodiment of the present invention to mark Result figure；
Fig. 8-1 is the effect picture of three-dimensional blowhole model crack identification in the embodiment of the present invention；
Fig. 8-2 is the minimum circumscribed rectangular body schematic diagram of the first crack in Fig. 8-1 in the embodiment of the present invention；
Fig. 8-3 is the minimum circumscribed rectangular body schematic diagram of the second crack in Fig. 8-1 in the embodiment of the present invention；
Fig. 8-4 is the minimum circumscribed rectangular body schematic diagram of third crack in Fig. 8-1 in the embodiment of the present invention；
Fig. 8-5 is the effect picture of three-dimensional blowhole model crack identification in the embodiment of the present invention；
Fig. 8-6 is the effect picture that crack and original hole are shown simultaneously in Fig. 8-5 in the embodiment of the present invention；
Fig. 8-7 is the effect picture of three-dimensional blowhole model crack identification in the embodiment of the present invention；
Fig. 8-8 is the effect picture that crack and original hole are shown simultaneously in Fig. 8-7 in the embodiment of the present invention；
Fig. 8-9 is the effect picture of three-dimensional blowhole model crack identification in the embodiment of the present invention；
Fig. 8-10 is the effect picture that crack and original hole are shown simultaneously in Fig. 8-9 in the embodiment of the present invention；
Fig. 9 is method flow diagram of the present invention.
Specific embodiment
Below by specific embodiment and in conjunction with attached drawing, invention is further described in detail, it is necessary to indicated herein It is that the embodiment is only intended to further describe of the invention, but is not meant to be to any of the scope of the present invention It limits.
Embodiment: in order to make recognition methods of the present invention easily facilitate understanding and close to true application, below from rock Stone binary sequence image starts until three-dimensional blowhole model crack identification is completed, and is illustrated to entire operation process. Specific steps are as follows:
(1) one group of rock binary sequence figure to be identified of preparation, picture number 256, format BMP, size 256 × 256, wherein representing as Figure 1-1；The three-dimensional blowhole model constructed by rock binary sequence figure, as shown in Figs. 1-2, figure Middle grey parts are blowhole.
(2) all connected components for extracting three-dimensional blowhole model, reject wherein connection of the tissue points number less than 100 Component, denoising result is as shown in Fig. 2-1.With marching cube isosurface extraction (Marching Cube) algorithm to each connection Component carries out resurfacing, obtains the corresponding triangle grid model A of each connected componentq, wherein q=1,2,3 ... k, k tri- Angle grid model number.Resurfacing result is as shown in Fig. 2-2.
(3) to each resurfacing result AqExecute Smoothness Index 0.2, Laplce's Mesh Smoothing of smooth number 100 Processing；Mesh Smoothing effect is as shown in model partial enlarged view 3-1 and 3-2.
(4) lattice simplified processing is executed to smoothed out triangle grid model using triangle gridding extraction algorithm, grid is taken out Take ratio control 0.9.Lattice simplified effect is as shown in model partial enlarged view 4-1 and 4-2, and straightway is triangle gridding in figure Boundary.
(5) unit normal vector of all triangle gridding units is calculated, if it is (x, y, z) that unit normal vector, which is calculated, if z When < 0, then adjusting unit normal vector is (- x ,-y ,-z).Unit normal vector calculated result is as shown in figure 5, figure intermediate cam grid list Arrow in member is its unit normal vector.
(6) using the unit normal vector of triangle gridding unit and triangle gridding unit triangle area to AqIn triangle gridding Unit carries out K- mean cluster, obtains three cluster centre V1、V2、V3。
K- means clustering algorithm process is as follows:
1) A is soughtqThe unit normal vector σ and A of the maximum triangle gridding unit of intermediate cam areaqIntermediate cam area is the smallest by three The unit normal vector ξ of angle grid cell.Then enabling σ, σ × ξ, σ × (σ × ξ) is respectively three initial cluster centres；
2) triangle gridding unit Δ is calculatedjWith each cluster centre ViDistance thetaij, calculation formula are as follows:
Wherein
Wherein λjFor triangle gridding ΔjUnit normal vector.
Triangle gridding unit is assigned to away from the nearest affiliated triangle gridding class of cluster centre.
3) Δ is enabledjFor AqIn j-th of triangle gridding unit, λjFor ΔjUnit normal vector, wherein [1, m] j ∈, m AqContain Triangle gridding unit number；For AqIn i-th of triangle gridding class, cluster centre Vi, wherein i=1,2,3；For I-th of triangle gridding classIn j-th of triangle gridding unit,ForUnit normal vector,ForTriangle area, wherein j ∈ [1, n], n are i-th of triangle gridding classIntermediate cam grid cell number；
Recalculate cluster centre Vi, calculation formula is as follows:
Wherein
Add is defined as:
Wherein n is i-th of triangle gridding class intermediate cam grid cell number.
4) at the end of calculating kth time iteration, cluster centre variance values
WhereinCluster centre at the end of for kth time iterationThree coordinate values.
If the number of iterations k is more than 50 or cluster centre amplitude of variationTerminate iteration less than 0.001, otherwise returns to stream Journey 2).
(7) each triangle gridding unit Δ is calculatedjWith cluster centre ViDistance thetaij, and
θkj=min { θij|i∈[1,2,3]} (9)
If θkj< 0.59359877 is then by ΔjIt is classified as kth triangle gridding classOtherwise ΔjIt is not belonging to any kind.
(8) triangle gridding class is extractedTriangle gridding connected component: if triangle gridding unit ΔiTriangle vertex with ΔjAt least one vertex of triangle vertex it is identical, then ΔiWith ΔjSpace is adjacent.Each triangle gridding classIn, spatially Neighbouring relations extract triangle gridding connected componentWherein [1, k] t ∈, k are triangle gridding classesTriangle gridding connection point Measure number.Triangle gridding connected component extracts result as shown in Fig. 7-1,7-2 and 7-3.
(9) each triangle gridding class is calculatedIn all triangle gridding units triangle area andCalculate the triangulation network Lattice connected componentIn all triangle gridding units triangle area andWherein [1, k] t ∈, k are triangle gridding classes's Triangle gridding connected component number.Calculation formula is respectively as follows:
WhereinFor triangle gridding unitTriangle area.
It is rightIn all triangle gridding connected componentsIf
Then by this triangle gridding connected componentIt is considered as noise jamming, and from triangle gridding classMiddle rejecting.
After cancelling noise, the cluster centre of each class is recalculated, and merges similar class.Calculate cluster centre ViWith it is poly- Class center VjDistance χij, calculation formula is as follows:
Wherein
If χij< 0.7853981, then it is adjacent to regard the affiliated triangle gridding class of two cluster centres.It will be in three triangle gridding classes All adjacent classes merge into a class.Triangle gridding class after merging is still denoted asWherein 1≤i≤3.The merging of adjacent class As a result as shown in Fig. 6-1,6-2, different gray scales respectively indicate a class in figure, and white triangles grid is not belonging to any kind.
(10) three-dimensional blowhole connected component AqIn each triangle gridding class in all triangle gridding units composition One three-dimensional surface structure.The operation of threedimensional model hole repair is executed to this three-dimensional surface structure, obtains one or more after operation Closed three-dimensional surface structure, one sealing surface structured set of these three-dimensional surface structure compositions.
(11) the minimum circumscribed rectangular body of sealing surface structured set: 3-D image crack automatic recognition is calculated by as follows offer advice [journal article] Deng Zhiqiu, Teng Qizhi-" computer and digital engineering " 1 phase in 2013.
(12) the form factor Φ of sealing surface structured set ψ is calculatedψ:
Calculate the form factor Φ of minimum circumscribed rectangular bodily form ∏∏:
Reference area ratio Λ:
Wherein VψFor the volume of sealing surface structured set ψ；SψFor all triangle gridding unit triangular facets before ψ hole repair The sum of product；VΠFor the volume of sealing surface structured set minimum circumscribed rectangular body ∏；SΠFor the surface area of Π.
(13) when following condition is set up:
Or when the following conditions are set up:
Wherein
Then determine that voxel point set that the current surface sealing surface structured set ψ and inside include belongs to a part in crack, and It performs the next step rapid；Otherwise it is determined as non-crack.
(14) all tissue points for extracting the surface and inside that are determined as genuine sealing surface structured set in (13), to mentioning The voxel point set obtained carries out the processing of 3-D image morphological dilations.
(15): by the tissue points of expansion obtains in (14) voxel point set and initial three-dimensional blowhole model connected component Collection carries out logical AND operation, and the voxel point set obtained after operation is crack.With operating result as shown in Fig. 8-1 to 8-10:
Fig. 8-1 is the crack schematic diagram extracted, shares 3 cracks in figure；
Fig. 8-2 to 8-4 is the minimum circumscribed rectangular body schematic diagram in each crack；
Fig. 8-5 to 8-10 is crack different direction display figure, and the effect picture that it is shown simultaneously with original hole.
In the present embodiment, with actual rock binary sequence figure construct threedimensional model and using method of the invention to its into Row crack identification.Recognition result is observed, to demonstrate the reliability and accuracy of recognition methods of the present invention.
The above embodiment of the present invention is the preferred embodiment of the present invention, is not to appoint to technical solution of the present invention What is limited, as long as the technical solution realized on the basis of the above embodiments without creative work, is regarded as falling Enter the present invention protect in the range of.

Claims (6)

1. a kind of three-dimensional rock image crack identification method, it is characterised in that after image preprocessing, algorithm includes following step It is rapid:
Step 1: calculating the triangle grid model A that pretreatment obtainsqIn each triangle gridding unit unit normal vector and triangular facet Product, and unit normal vector direction is adjusted, wherein q=1,2,3 ... k, k are triangle grid model number；Followed by triangle gridding The unit normal vector and triangle area feature of unit are to AqTriangle gridding unit carry out K- mean cluster；
Step 2: the cluster centre obtained according to step 1 is by AqTriangle gridding unit reclassify, classification obtain multiple triangulation networks Lattice class；
Step 3: extracting the triangle gridding connected component of each triangle gridding class in step 2；
Step 4: rejecting the noise contribution in the triangle gridding connected component that step 3 obtains, recalculate the poly- of triangle gridding class Class center, and the adjacent all triangle gridding classes of the cluster centre recalculated are merged into one kind；Then to each three Angle grid class independently executes the operation that step 5 arrives step 10；
Step 5: all triangle gridding units in each triangle gridding class that step 4 obtains form a three-dimensional surface structure, right This three-dimensional surface structure executes the operation of threedimensional model hole repair, and one or more closed three-dimensional surface knots can be obtained after operation Structure, one sealing surface structured set of these three-dimensional surface structure compositions；
Step 6: calculating the minimum circumscribed rectangular body for the sealing surface structured set that step 5 obtains；
Step 7: calculating the form factor for the sealing surface structured set that step 5 obtains；Calculate the minimum circumscribed rectangular that step 6 obtains The form factor of body；Reference area compares parameter；
Step 8: whether two form factors and area ratio parameter obtained using step 7 determine current sealing surface structured set For crack；If it is determined that result is very, then to continue to execute subsequent operation；If it is determined that result is false, then current sealing surface structured set It is not crack；
Step 9: determining that result is the tissue points on genuine sealing surface structured set surface and inside in extraction step 8, to extraction Voxel point set makees the processing of 3-D image morphological dilations；
Step 10: by the voxel point set of the voxel point set expanded in step 9 and initial three-dimensional blowhole model connected component Logical AND operation is carried out, the voxel point set obtained after operation is crack.
2. according to the method described in claim 1, it is characterized in that in step 1 using triangle gridding unit unit normal vector and The K- mean cluster of triangle area feature calculates triangle gridding unit cluster centre, using following calculation methods:
K- means clustering algorithm process is as follows:
1. determining initial cluster center；
2. being assigned to all triangle gridding units away from the affiliated triangle gridding class of its nearest cluster centre；
3. updating cluster centre；
4. process is executed 2. repeatedly, 3. until meeting termination condition；
Relevant calculation formula is as follows in K- means clustering algorithm:
Process 1. in initial cluster center ViCalculation formula it is as follows, wherein i=1,2,3:
Wherein σ is AqThe unit normal vector of the maximum triangle gridding unit of intermediate cam area；ξ is AqIntermediate cam area is the smallest by three The unit normal vector of angle grid cell；
Enable ΔjFor AqIn j-th of triangle gridding unit, λjFor ΔjUnit normal vector, wherein [1, m] j ∈, m AqThe triangle contained Grid cell number；For AqIn i-th of triangle gridding class, cluster centre Vi, wherein i=1,2,3；It is i-th three Angle grid classIn j-th of triangle gridding unit,ForUnit normal vector,ForTriangle area, wherein [1, n] j ∈, n For i-th of triangle gridding classIntermediate cam grid cell number；
Process 2. intermediate cam grid cell ΔjWith each cluster centre ViDistance thetaijCalculation formula it is as follows:
Wherein
3. middle cluster centre calculation formula is as follows for process:
Wherein
Add is defined as:
Wherein n is i-th of triangle gridding classIntermediate cam grid cell number；
Process 4. in stopping criterion for iteration are as follows: the number of iterations is more than 50 times or cluster centre amplitude of variation less than 0.001；K Secondary iteration cluster centre variance valuesCalculation formula it is as follows:
WhereinCluster centre at the end of for kth time iterationThree coordinate values.
3. according to the method described in claim 1, using it is characterized in that step 2 intermediate cam grid cell reclassifies method Following calculation methods:
Triangle gridding unit ΔjWith cluster centre ViDistance be θij, and
θkj=min { θij|i∈[1,2,3]} (10)
If θkjThe then Δ of < 0.59359877jBelong to kth triangle gridding classOtherwise ΔjIt is not belonging to any kind.
4. according to the method described in claim 1, it is characterized in that step 4 intermediate cam mesh connectivity component noise ingredient reject and Merge adjacent all triangle gridding classes, using following calculation methods:
Triangle gridding classTriangle gridding unit triangle area and beTriangle gridding connected componentTriangle gridding list First triangle area and it isWherein [1, k] t ∈, k are triangle gridding classesTriangle gridding connected component number；Calculation formula It is as follows:
WhereinFor triangle gridding unitTriangle area；
It is rightIntermediate cam mesh connectivity componentIf
Then by triangle gridding connected componentIt is considered as noise jamming, and from triangle gridding classMiddle rejecting；
If cluster centre ViWith cluster centre VjDistance be χij, calculation formula is as follows:
Wherein
If χij< 0.7853981, then the view affiliated triangle gridding class of two cluster centres is adjacent, will own in three triangle gridding classes Adjacent class merges into a class；Triangle gridding class after merging is still denoted asWherein 1≤i≤3.
5. according to the method described in claim 1, it is characterized in that in step 8 crack judgement, using following criterion: when When following condition is set up:
Or when the following conditions are set up:
Wherein
VψFor the volume of sealing surface structured set ψ, SψFor the sum of triangle gridding unit triangle areas all before ψ hole repair, Φψ For the form factor of ψ；VПFor the volume of sealing surface structured set ψ minimum circumscribed rectangular body Π, SΠFor the surface area of Π, Φ∏For The form factor of ∏；
Then determine that voxel point set δ that the current surface sealing surface structured set ψ and inside include belongs to a part in crack；Otherwise sentence It is set to non-crack.
6. according to the method described in claim 1, it is characterized in that in step 10 crack accurate extraction, use following extraction side Method:
Crack voxel point set ζ calculation formula are as follows:
ζ=p | p ∈ Hq∩δ′} (23)
Wherein HqFor AqThe voxel point set of corresponding initial three-dimensional blowhole model connected component, p are any tissue points；δ ' is The voxel point set obtained after morphological dilations processing in step 9.
CN201710030343.6A 2017-01-17 2017-01-17 A kind of three-dimensional rock image crack identification method Active CN106803253B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710030343.6A CN106803253B (en) 2017-01-17 2017-01-17 A kind of three-dimensional rock image crack identification method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710030343.6A CN106803253B (en) 2017-01-17 2017-01-17 A kind of three-dimensional rock image crack identification method

Publications (2)

Publication Number Publication Date
CN106803253A CN106803253A (en) 2017-06-06
CN106803253B true CN106803253B (en) 2019-06-04

Family

ID=58984691

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710030343.6A Active CN106803253B (en) 2017-01-17 2017-01-17 A kind of three-dimensional rock image crack identification method

Country Status (1)

Country Link
CN (1) CN106803253B (en)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109521180B (en) * 2018-10-09 2020-08-11 水利部交通运输部国家能源局南京水利科学研究院 Analysis method for cause of degradation and damage of expansive soil channel in cold region and application
CN109580649B (en) * 2018-12-18 2020-11-27 清华大学 Engineering structure surface crack identification and projection correction method and system
CN110018108B (en) * 2019-05-14 2020-04-28 中国科学院地质与地球物理研究所 Method and equipment for measuring pore diameter of rock pore
CN112307803A (en) * 2019-07-25 2021-02-02 中国石油天然气股份有限公司 Digital geological outcrop crack extraction method and device
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
CN113935928B (en) * 2020-07-13 2023-04-11 四川大学 Rock core image super-resolution reconstruction based on Raw format
CN112163621A (en) * 2020-09-29 2021-01-01 成都理工大学 Compact sandstone reservoir pore structure classification and characterization method based on micro ct technology
CN112132966A (en) * 2020-09-29 2020-12-25 成都理工大学 Shale fracture network connectivity characterization method based on topological structure
CN113034346B (en) * 2020-10-10 2023-08-11 清能艾科（深圳）能源技术有限公司 Method and device for restoring broken rock core, electronic equipment and storage medium
CN113379783B (en) * 2021-03-29 2022-09-06 南京理工大学 Two-dimensional fracture reproduction method based on cross-correlation function
CN113222923B (en) * 2021-04-30 2021-11-12 广东石油化工学院 Method for identifying crack type based on core photo
CN117152187B (en) * 2023-10-30 2024-01-26 山东中科冶金矿山机械有限公司 Crack contour extraction method in geological mapping

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101303770A (en) * 2008-05-28 2008-11-12 中山大学 Method for constructing triangle gridding as well as application thereof for geotechnical engineer modeling
CN102867302A (en) * 2012-08-30 2013-01-09 四川大学 Core fracture identification method based on three-dimensional image information processing
CN104376556A (en) * 2014-10-31 2015-02-25 四川大学 Rock CT image target segmentation method
CN104794709A (en) * 2015-04-10 2015-07-22 四川大学 Three-dimensional core image pore and throat segmentation method
WO2016006283A1 (en) * 2014-07-10 2016-01-14 三菱電機株式会社 Structure maintenance management system
CN105261068A (en) * 2015-11-16 2016-01-20 中国石油大学(华东) Micro-CT technology-based reservoir core three-dimensional entity model reconstruction method
CN105279794A (en) * 2015-11-25 2016-01-27 中国石油大学(华东) Reservoir stratum rock core multi-organizational model constructing method based on Micro-CT technology

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101303770A (en) * 2008-05-28 2008-11-12 中山大学 Method for constructing triangle gridding as well as application thereof for geotechnical engineer modeling
CN102867302A (en) * 2012-08-30 2013-01-09 四川大学 Core fracture identification method based on three-dimensional image information processing
WO2016006283A1 (en) * 2014-07-10 2016-01-14 三菱電機株式会社 Structure maintenance management system
CN104376556A (en) * 2014-10-31 2015-02-25 四川大学 Rock CT image target segmentation method
CN104794709A (en) * 2015-04-10 2015-07-22 四川大学 Three-dimensional core image pore and throat segmentation method
CN105261068A (en) * 2015-11-16 2016-01-20 中国石油大学(华东) Micro-CT technology-based reservoir core three-dimensional entity model reconstruction method
CN105279794A (en) * 2015-11-25 2016-01-27 中国石油大学(华东) Reservoir stratum rock core multi-organizational model constructing method based on Micro-CT technology

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
基于中轴线的岩心3维图像孔喉分割算法;龚小明等;《四川大学学报（工程科学版）》;20160630;第48卷;全文
基于模糊C均值聚类算法的岩屑;万红吉等;《四川大学学报（自然科学版）》;20100531;第47卷(第3期);全文
岩心3维图像孔喉分割优化算法;刘雨晨等;《四川大学学报（工程科学版）》;20120630;第44卷;全文
改进相位一致的岩心图像裂缝提取方法;李朝辉等;《计算机与数字工程》;20151231;第43卷(第4期);全文

Also Published As

Publication number Publication date
CN106803253A (en) 2017-06-06

Similar Documents

Publication Publication Date Title
CN106803253B (en) A kind of three-dimensional rock image crack identification method
CN105741355B (en) A kind of block dividing method of triangle grid model
CN101540000B (en) Iris classification method based on texture primitive statistical characteristic analysis
CN104008553B (en) Crack detection method with image gradient information and watershed method conflated
CN105975913B (en) Road network extraction method based on adaptive cluster learning
CN104850825A (en) Facial image face score calculating method based on convolutional neural network
CN109165674A (en) A kind of certificate photo classification method based on multi-tag depth convolutional network
CN105139379B (en) Based on the progressive extracting method of classified and layered airborne Lidar points cloud building top surface
CN105321176A (en) Image segmentation method based on hierarchical higher order conditional random field
CN109191471A (en) Based on the pancreatic cell image partition method for improving U-Net network
CN109919159A (en) A kind of semantic segmentation optimization method and device for edge image
CN103996018A (en) Human-face identification method based on 4DLBP
CN106651865B (en) A kind of optimum segmentation scale automatic selecting method of new high-resolution remote sensing image
CN105809113B (en) Three-dimensional face identification method and the data processing equipment for applying it
CN104331716A (en) SVM active learning classification algorithm for large-scale training data
CN111681274A (en) 3D human skeleton recognition and extraction method based on depth camera point cloud data
CN109447100A (en) A kind of three-dimensional point cloud recognition methods based on the detection of B-spline surface similitude
CN110021028A (en) A kind of automatic clothing method based on garment fashion drawing
CN108961429A (en) A kind of cultural relic fragments model divides automatically and joining method
CN102629321A (en) Facial expression recognition method based on evidence theory
CN104036550A (en) Laser radar point-cloud interpretation and reconstruction method for building elevations on basis of shape semantics
CN109993213A (en) A kind of automatic identifying method for garment elements figure
CN110349159A (en) 3D shape dividing method and system based on the distribution of weight energy self-adaptation
Liu et al. A novel rock-mass point cloud registration method based on feature line extraction and feature point matching
CN111886630A (en) Three-dimensional cell and tissue image analysis for cellular and subcellular morphological modeling and classification

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