CN102867302A - Core fracture identification method based on three-dimensional image information processing - Google Patents

Abstract

The invention relates to a core fracture identification method based on three-dimensional image information processing, belonging to the technical field of image segmentation. The method comprises the following steps of: scanning a core to be processed by using a CT (computerized tomography) machine, storing a scanned core CT sequence image in a computer, reading the well scanned core CT sequence image by using core three-dimensional image fracture identification software, performing three-dimensional image reconstruction on the core CT sequence image by utilizing a marching cubes algorithm, and performing a series of calculation, getting and screening processed by utilizing the differences among pores, random noise and fracture in morphological features in a three-dimensional data field to achieve the purpose of identifying a core fracture. According to the method provided by the invention, the core fracture can be fast and accurately identified; and the problems of causing error and even mistake due to negligence of communication information of the core fracture between layers in the existing isolate analysis of each core CT sequence image can be solved.

Description

Rock core crack identification method based on the three-dimensional image information processing

Technical field

The present invention relates to a kind of image processing techniques, particularly a kind of rock core crack identification method of processing based on three-dimensional image information belongs to the image Segmentation Technology field.

Background technology

In the oil geology field, crack identification is the important content of oil reservoir and reservoir research, and is significant to the raising of oil and gas development and oil recovery factor.Rock produces mechanical destruction under effect of stress, the rift structure of the obvious displacement of nothing is in the crack.Because crannied existence, and fracture development are connective good, just can make mud shale and carbonate reservoir become valuable reservoir.

China carries out very early to the research work in rock core crack, the later stage fifties in last century, the geologist of China's oil gas field just begins to identify the crack that exists in the rock core by experience according to features such as the morphological feature of rock core and tomography positions, and this is a kind of recognition methods based on qualitative analysis.Since the eighties in last century, mainly by appear, the data such as rock core, well logging, earthquake, well testing, formation testing pilot production, drilling well, well logging, identify crack in the rock core by means of the analysis-by-synthesis such as determination method, well test analysis method, crack Statistics Method and well-log information method in geology way of qualitative analysis, the core chamber.By this method, can analyze qualitatively the rock core crack; Simultaneously obtain fracture number in the rock core, the sxemiquantitative information such as fracture length according to statistics.But said method all is crack existence and the distribution situation of analyzing whole rock core inside by the feature on some surfaces of observing rock core, does not obtain the real structure of rock core inside.From the end of the nineties in last century so far, increasingly mature along with the computer image treatment means moves towards the pattern recognition stage to the identification in rock core crack from the numerical simulation calculation stage.Oil geology research department is widely used the computerized tomography technology now, namely claim again computed tomography, CT, scan containing crannied rock core by the CT machine, obtain rock core CT sequence image, just can be in the situation that harmless rock core utilizes rock core CT sequence image described above to observe, analyze the real structure of this rock core inside.The features of different material on rock core CT sequence image such as crack, hole and rock there are differences, can extract crack and hole according to image segmentation algorithm, identify crack in the rock core according to morphological feature again, not only can analyze the crack on the core surface, and can go deep into inner its crack of observing of rock core in the situation of change of different aspects.

Wang Renyi, Ao pioneer (extract [J] based on the rock core crack pattern picture of image processing techniques. Xinjiang Geology, 2006, (04) .) etc. utilize the Binarization methods of maximum variance between clusters to extract the crack from core profile; Superb, Wang Zhengyong (based on the rock core crack extract algorithm [J] of Beamlet. information and electronic engineering, 2010, (02) .) a kind of rock core crack pattern based on multiple dimensioned Beamlet conversion is proposed as extraction algorithm; Analogy continue industry (some Study on Problems [D] that rock core CT sequence image is processed. Sichuan University's Dissertation Database, 2011) propose to extract the crack based on the phase place unification algorism.But, said method has plenty of the core surface of analysis image, have plenty of and analyze each Zhang Yanxin CT sequence image isolatedly and identify its crack, only utilized the rock core two-dimensional image information, and all CT sequence chart are not combined analysis, ignored the communication information of rock core crack between the different rock core CT sequence images more than two or two, so be easy to and become long and narrow linear in the two dimensional image, just isolated and disconnected hole or random noise are mistaken for the crack mutually in 3 d data field, perhaps subcircular will be in the two dimensional image, should be that the target of the part in crack is mistaken for and is not the crack in 3 d data field, so often cause the recognition result in rock core crack deviation to occur.Therefore purpose of the present invention just provides a kind of method, can identify accurately and rapidly the rock core crack.

Summary of the invention

Purpose of the present invention is to overcome existing defective and deficiency in the prior art just, and a kind of rock core crack identification method based on the three-dimensional image information processing that can identify quickly and accurately the rock core crack is provided.The method is to use first the pending rock core of CT machine scanning, obtains rock core CT sequence image; According to mobile cube (Marching Cubes, MC) algorithm the rock core CT sequence image that obtains is carried out three-dimensional image reconstruction again; Then the difference of utilizing the feature of the different materials such as crack, hole and rock on rock core CT sequence image to exist by a series of calculating, obtain and screening process, thereby reaches the purpose in identification rock core crack.

For achieving the above object, the present invention adopts following technical scheme to realize.

The present invention is based on the rock core crack identification method that three-dimensional image information is processed, its invention main points mainly are that the rock core CT sequence image that scanning obtains is carried out three-dimensional image reconstruction, and carry out crack identification and calculating in the 3-D view of rebuilding, and its ultimate principle is as follows:

In CT rock core sequence image, except containing rock, crack and hole information, also comprise inevitably a large amount of random noises.Wherein, the feature of rock is that gray level is lower; The feature of hole is that gray level is very high, generally shows as irregular polygon, ellipse, sub-circular; Although the gray level in crack is also very high, it is that rock texture loses a successional plane, so different from hole, it is long and narrow linear in two dimensional image or dendritic distributes.Random noise does not then have fixing gray level, generally has an area little in two dimensional image, the characteristics such as discrete distribution.In order to obtain quantitative information from rock core CT sequence image, an important step will be cut apart rock core CT sequence image exactly, to isolate hole target and crack target from the rock core CT sequence image of reality.But the partitioning algorithm of intensity-based information unavoidably can extract picture noise simultaneously.Present method is usually according to hole and random noise and crack modal different differentiation in two dimensional image.Yet, being the target of long and narrow linear in the two dimensional image, probably just mutually isolated and disconnected hole or random noise in the 3 d data field of reality be not so it is the crack.The target of described sub-circular also is likely the part in crack in 3 d data field.Therefore, analyze each Zhang Yanxin CT sequence image comes crack identification deviation even mistake probably to occur isolatedly.

Described crack and hole morphological feature in 3 d data field differs larger, simultaneously, learns by observing rock core CT sequence image, and for real noise, then it does not have continuity between the different rock core CT sequence images more than two or two; For the crack, although be the impact point that isolates from two dimensional image, between the different rock core CT sequence images more than two or two, there is very large probability to interconnect.Therefore, the present invention carries out the rock core crack identification according to the difference of crack, hole and random noise morphological feature in 3 d data field.Namely use the pending rock core of CT machine scanning, obtain rock core CT sequence image and be stored in the computing machine, carry out three-dimensional image reconstruction by mobile cube (Marching Cubes) algorithm.In the 3 d data field that three-dimensional image reconstruction consists of, obtain one or more than one target.Each target in the described 3 d data field has comprised the information of crack, hole and the random noise of rock core inside.

The present invention is based on the rock core crack identification method that three-dimensional image information is processed, it is characterized in that comprising successively following concrete operation step:

Step 1: use the pending rock core of CT machine scanning, obtain rock core CT sequence image and be stored in the computing machine;

Step 2: use rock core 3-D view crack identification software to read in the good rock core CT sequence image of step 1 scanning, and utilize marching cubes algorithm that rock core CT sequence image is carried out three-dimensional image reconstruction, in the 3 d data field that three-dimensional image reconstruction consists of, obtain one or more than one target;

Step 3: each target that three-dimensional image reconstruction obtains in the scanning step 2, obtain the grown form parameter of each target, comprise target volume, target surface area; Calculate the value of the 3D shape factor F of each target, the relational expression of F and described target volume and target surface area is following formula:

$F=\frac{36\mathrm{\π}{V}_p}{S_p}---\left(1\right)$

Wherein, V _pBe target volume, S _pBe target surface area, calculate the value of the F of each target according to formula (1); Again the F value is screened, if satisfy F＜0.05, illustrate that this target is planar distribution, keep this area target, otherwise the non-area target of this target is described, then this non-area target of deletion;

Step 4: each area target that keeps in the scanning step 3, obtain the minimal circumscribed sphere radius R of each area target _MinWith equivalent spheres radius R _Eq, the ratio of calculating minimal circumscribed sphere radius and equivalent spheres radius

According to Ratio screen again, if satisfy

Illustrate that this area target possesses ductility simultaneously, then keep the area target that this possesses ductility, otherwise illustrate that this area target does not possess ductility, then deletes the area target that this does not possess ductility;

Step 5: each that keeps in the calculation procedure 4 possesses the length of the external rectangular parallelepiped of approximate minimum of the area target of ductility, and length and width and height are compared between any two, obtains the longest edge D among them _MaxMinor face D _Min, according to

Ratio screen, if satisfy Illustrate that this area target that is ductile is the crack, it is kept as the crack of identifying, if this area target that is ductile does not satisfy described condition, illustrate that this area target that is ductile is not the crack, then deletes this area target.

Mobile cube described in the step 2 of the present invention (Marching Cubes) algorithm is the classic algorithm in the iso-surface patch algorithm, and it is also referred to as contour surface and extracts (Isosurface Extraction).Be applied among the present invention, be actually each width of cloth figure in the rock core CT sequence image is merged together, regard a 3 d data field as, the voxel that therefrom will have the same grayscale value extracts, and couples together with topological form, forms contour surface.The method of extracting contour surface is each voxel of processing one by one in the above-mentioned 3 d data field, and deciding the structural form of the inner contour surface of this voxel according to the gray-scale value of each each angle point of voxel, the contour surface of all voxel inside is combined and gets final product the most at last.

Described voxel generally has two kinds of definition: a kind of be with common two dimensional image in the definition of pixel similar, directly the sampled point in the 3 d data field as voxel; Another kind of then be eight adjacent sampled points in the 3 d data field, i.e. summit on the square of unit length size, the zone definitions that comprises is voxel, what the present invention adopted is the second definition.

Described mobile cube (Marching Cubes) algorithm, its reference is as follows:

Qian Feng, horse is beautiful, Yang Shengqi, Wan Wanggen. the research of marching cubes algorithm and improvement [J] computer engineering and application, 2010, (34).

Neilson G M.On marching cubes[J].IEEE Transactions on Visualization and Computer Graphics，2003，9（3）：283-297.

Miu Bin and, Deng Yuanmu, Huang Feizeng. based on corresponding point matching faultage image three-dimensional interpolation method [J] Chinese medicine physics magazine, 2000,17(1): 14-16.

The ginger Yun. based on the image three-dimensional visualization algorithm research [D] of iso-surface patch. Chengdu: University of Electronic Science and Technology, 2007.

The minimal circumscribed sphere of the described target of step 4 of the present invention is to make all voxels of this target all be positioned at the ball of the radius minimum of ball, and the radius of this ball is exactly the minimal circumscribed sphere radius R of described target _MinDescribed equivalent spheres radius R _EqIt is the radius of the ball that equates with described target volume.

The method of length that the described calculating of step 5 of the present invention possesses the approximate minimum external rectangular parallelepiped of area target of ductility may further comprise the steps:

(1) calculates the described height that possesses the minimum external rectangular parallelepiped of area target of ductility

Adopt least square method to seek the fit Plane of described target, its ultimate principle is to make the quadratic sum of deviation minimum, specifically is exactly to find a plane, so that each the distance from the point to the plane sum is minimum on the described target; The three dimensions midplane formula that adopts cartesian coordinate system to determine is as follows:

Ax+By+Cz+D=0 （2）

Wherein (x, y, z) is a bit on the plane, A, and B, C are x, y, the coefficient of z, D are constant factor.

Any point (t in the three dimensions ₁, t ₂, t ₃) as follows to the range formula of plane Ax+By+Cz+D=0:

$d=\frac{|{\mathrm{At}}_1+{\mathrm{Bt}}_2+{\mathrm{Ct}}_3+D|}{\sqrt{A^2+B^2+C^2}}---\left(3\right)$

The error formula on data point and plane is as follows in the Calculation of Three Dimensional space:

e=Ax+By+Cz+D （4）

According to the least square ratio juris, the formula of plane least square fitting is as follows:

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{e_i}^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D)}^2---\left(5\right)$

(x wherein _i, y _i, z _i) for needing a bit in the fit Plane, e _iThe error that represents this data point and plane.

For formula (2), establish D=1, namely get plane Ax+By+Cz=-1, according to formula (5), if point (x _i, y _i, z _i) in the plane, Ax then _i+ By _i+ Cz _i=-1, e _i=Ax _i+ By _i+ Cz _i+ 1, the plane least square fitting, namely

Minimum, A, B, C are got respectively partial differential:

$\frac{\∂}{\∂A}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)}^2=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)x_i=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\mathrm{Ax}}_i^2+{\mathrm{By}}_i{x}_i+{\mathrm{Cz}}_i{x}_i+x_i)=0$

$\frac{\∂}{\∂B}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)}^2=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)y_i=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(A{x}_i{y}_i+{\mathrm{By}}_i^2+{\mathrm{Cz}}_i{y}_i+y_i)=0$

$\frac{\∂}{\∂B}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)}^2=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+1)Z_i=2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}({\mathrm{Ax}}_i{z}_i+{\mathrm{By}}_i{z}_i+{\mathrm{Cz}}_i^1+z_i)=0---\left(6\right)$

Get system of equations after the arrangement:

$\left\{\begin{array}{c}A\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^2+B\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{x}_i+C\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{z}_i{x}_i+\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{x}_i=0\\ A\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i+B\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i^2+C\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{z}_i{y}_i+\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i=0\\ A\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{z}_i+B\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{y}_i{z}_i+C\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{z}_i^2+\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{z}_i=0\end{array}\right.---\left(7\right)$

Solve an equation:

$\left[\begin{array}{c}A\\ B\\ C\end{array}\right]={\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{y}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i^2& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i{z}_i\\ \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i{z}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i{z}_i& \underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i^2\end{array}\right]}^{-1}\left[\begin{array}{c}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i\\ -\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\\ -\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{z}_i\end{array}\right]---\left(8\right)$

According to (7) formula as can be known, by the computing to the point set coordinate, can access the parameters of fit Plane, and described target is that three-dimensional point set consists of by a series of voxels, and be planar distribution, therefore, select to utilize the method for least square plane match to ask its fit Plane.

According to formula (2), the equation of establishing fit Plane is Ax+By+Cz+D=0, some P _i(x _i, y _i, z _i) be any point on the described target, its subpoint P on fit Plane then _jCoordinate be P _j(x _i-Ak, y _i-Bk, z _i-Ck), wherein $k=\frac{D+{\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i}{A^2+B^2+C^2},$ According to formula (3) but invocation point P _iDistance to the plane is $\frac{|{\mathrm{Ax}}_i+{\mathrm{By}}_i+{\mathrm{Cz}}_i+D|}{\sqrt{A^2+B^2+C^2}},$ Travel through each point on the described target, can obtain each some projection point set PointArray in the plane, and each the distance from the point to the plane set M.Relatively gather each numerical value among the M, judge k＜0 or k 〉=0, obtain respectively being positioned at the point of both sides, plane to the ultimate range d on plane in above-mentioned target _Max1, d _Max2, the formula of height H that calculates the approximate minimum external rectangular parallelepiped of above-mentioned target is as follows:

H=d _max1+d _max2 （9）

(2) calculate the described length that possesses the minimum external rectangular parallelepiped of area target of ductility

Ask any two points p in the known three dimensions ₁(a ₁, b ₁, c ₁), p ₂(a ₂, b ₂, c ₂) between the formula of distance as follows:

$f=\sqrt{{(a_1-a_2)}^2+{(b_1-b_2)}^2+{(c_1-c_2)}^2}---\left(10\right)$

Projection point set PointArray as a convex closure, is used and rotates 2 P that the algorithm that gets stuck is obtained distance maximum on the convex closure ₁(x ₁, y ₁, z ₁) and P ₂(x ₂, y ₂, z ₂), the length of the approximate minimum external rectangular parallelepiped of described target is P ₁To P ₂Distance, the formula of length L that can get the approximate minimum external rectangular parallelepiped of described target according to (10) formula is as follows:

$L=\sqrt{{(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2}---\left(11\right)$

(3) calculate the described width that possesses the minimum external rectangular parallelepiped of area target of ductility

Ask vector

To vector

The formula of projection as follows:

$n=\frac{|\stackrel{\→}{C}\×\stackrel{\→}{D}|}{\left|\stackrel{\→}{D}\right|}---\left(12\right)$

P sets up an office _k(x _k, y _k, z _k) be any point on the convex closure, vector Vector

By formula (12) but invocation point P _kTo P ₁, P ₂The straight line of determining apart from d _kFormula as follows:

$d_k=\frac{|\stackrel{\→}{A}\×\stackrel{\→}{B}|}{\left|\stackrel{\→}{B}\right|}---\left(13\right)$

Traversal each point on the convex closure can obtain on the convex closure each point to P ₁And P ₂The distance set of the straight line that determines.Relatively each numerical value in the distance set is judged

Value be less than 0 or more than or equal to 0, obtain respectively being positioned at by P on the convex closure ₁And P ₂The point of the straight line both sides that determine is to the ultimate range d of this straight line _Max3And d _Max4, the width of the approximate minimum external rectangular parallelepiped of above-mentioned target is two ultimate range sums, the formula of width W that calculates the approximate minimum external rectangular parallelepiped of above-mentioned target is as follows:

W=d _max3+d _max4 （14）

Through said process, finally calculated height, the length and wide of the approximate minimal circumscribed sphere of described target.

The present invention has following characteristics and useful technique effect:

Rock core crack identification method of the present invention is used the pending rock core of CT machine scanning, obtains rock core CT sequence image.Carry out three-dimensional image reconstruction by the rock core CT sequence chart that scanning is obtained, and utilize hole, random noise and the crack difference that morphological feature exists in 3 d data field, can identify accurately and rapidly the rock core crack; Solve each the Zhang Yanxin CT sequence image of analysis that in the past isolated, but ignored error even wrong problem that rock core crack communication information between the layers causes.

Description of drawings

The rock core CT sequence image that Fig. 1 the present invention is pending; Wherein, (a) being 4 pictures among the rock core CT sequence image A, (b) is 4 pictures among the rock core CT sequence image B;

Fig. 2 is the result behind the three-dimensional image reconstruction under the rock core CT sequence image different points of view among Fig. 1; Wherein, (a) being the three-dimensional image reconstruction result of rock core CT sequence image A, (b) is the three-dimensional image reconstruction result of rock core CT sequence image B;

Fig. 3 is the result in the three-dimensional rock core crack that extracts in the rock core CT sequence image among Fig. 2; Wherein, (a) being the extraction result of rock core CT sequence image A, (b) is the extraction result of rock core CT sequence image B; From then on figure can intuitively clearly see the crack that exists in the rock core CT sequence image, illustrates that namely rock core crack identification effect of the present invention is accurately.

Embodiment

Below in conjunction with embodiment and as a result figure the present invention is described in further detail, but protection scope of the present invention is not limited to the following stated; Below all " rock core 3-D view crack identification softwares " of mentioning all by seminar of the present invention independent research, it uses operation steps to be operation steps for the rock core crack identification method of processing based on three-dimensional image information.

Embodiment 1

Step 1: use the pending rock core of CT machine scanning, obtain 400 rock core CT sequence images and be stored in the computing machine, such as (a) figure among Fig. 1, observe picture and find, the position that the figure center takes over has the target of long and narrow linear;

Step 2: the picture function that reads that uses rock core 3-D view crack identification software, read in the good rock core CT sequence image of step 1 scanning, this software carries out three-dimensional image reconstruction according to mobile cube (Marching Cubes) algorithm to this CT sequence chart, obtains the 3 d data field that is made of above-mentioned rock core CT sequence image.

Step 3: use the targeted scans function of rock core 3-D view crack identification software to carry out stroke scanning to rebuilding the 3 d data field that obtains in the step 2, whenever obtain a target just with its preservation.Scanning obtains 53 mutual disjunct targets after finishing altogether, such as (a) figure among Fig. 2, from then on can be observed long and narrow planar target among the figure, and granular little target;

Step: 4: use the calculating form factor function of rock core 3-D view crack identification software, the grown form parameter of 53 targets that targeted scans obtains in the calculation procedure 3 comprises target volume, target surface area; Calculate the value of the 3D shape factor F of above-mentioned each target by formula (1); After calculating is finished, find to only have goal satisfaction F＜0.05, be area target, this software keeps this area target, and all the other targets are all deleted;

Step 5: use the calculating radius of a ball function of rock core 3-D view crack identification software, the minimal circumscribed sphere radius R of the area target that keeps in the calculation procedure 4 _MinWith equivalent spheres radius R _Eq, and the ratio of calculating minimal circumscribed sphere radius and equivalent spheres radius

Find this goal satisfaction

For possessing malleable area target, this software keeps the area target that this possesses ductility;

Step 6: use the digital simulation plane function of rock core 3-D view crack identification software, by the fit Plane of the area target that possesses ductility of reservation in formula (2) and formula (8) calculation procedure 5.

Step 7: use the minimum external rectangular parallelepiped height function of calculating of rock core 3-D view crack identification software, by the external rectangular parallelepiped height of approximate minimum of the area target that possesses ductility of reservation in formula (1) and formula (9) calculation procedure 5.

Step 8: use the minimum external rectangular parallelepiped length function of calculating of rock core 3-D view crack identification software, by the external rectangular parallelepiped length of approximate minimum of the area target that possesses ductility of reservation in formula (11) calculation procedure 5.

Step 9: use the minimum external rectangular parallelepiped width function of calculating of rock core 3-D view crack identification software, by the external wide side's height of approximate minimum of the area target that possesses ductility of reservation in formula (13) and formula (14) calculation procedure 5.

Step 10: the minimum external rectangular parallelepiped length of side ratio function of calculating of using rock core 3-D view crack identification software, to step 7, obtain respectively length, the width of approximate minimum external rectangular parallelepiped in the step 8, step 9 and compare highly between any two, obtain the longest edge D among them _MaxMinor face D _Min, calculate

Ratio, find that the area target that possesses ductility that keeps in the step 5 satisfies

This software keeps it as the crack of identifying, shown in (a) figure among Fig. 3, co-exist in a crack through this rock core of identification.

Recorded crack 1 target that extracts through above-mentioned steps and three hole targets---hole 1 wherein in the following table 1, the three-dimensional configuration parameter comparison of hole 2 and hole 3, by formula (1),

With

The calculating acquired results as shown in table 1 below.

Table 1 crack 1 and hole 1, the three-dimensional configuration parameter comparison of hole 2 and hole 3

Observe table 1 and find, step 3 is satisfied in crack 1 wherein, and three screening conditions of step 4 and step 5 are successfully identified; And hole 1,3 in hole 2 and hole do not satisfy step 3, and three screening conditions of step 4 and step 5 final are deletedly fallen.

Embodiment 2

Step 1: use the pending rock core of CT machine scanning, obtain 205 rock core CT sequence images and be stored in the computing machine, such as (b) figure among Fig. 1, observe picture and find, also have sub-circular and granular target among the figure except the target of long and narrow linear, their gray level is very approaching;

Step 2: the picture function that reads that uses rock core 3-D view crack identification software, read in the good rock core CT sequence image of step 1 scanning, this software carries out three-dimensional image reconstruction according to mobile cube (Marching Cubes) algorithm to this CT sequence chart, obtains the 3 d data field that is made of above-mentioned rock core CT sequence image.

Step 3: use the targeted scans function of rock core 3-D view crack identification software to carry out stroke scanning to rebuilding the 3 d data field that obtains in the step 2, whenever obtain a target just with its preservation.Scanning obtains 340 mutual disjunct targets after finishing altogether, such as (b) figure among Fig. 2, learns that from figure long and narrow planar target is arranged, and graininess, approximate spherical little target;

Step: 4: use the calculating form factor function of rock core 3-D view crack identification software, the grown form parameter of 340 targets that targeted scans obtains in the calculation procedure 3 comprises target volume, target surface area; Calculate the value of the 3D shape factor F of above-mentioned each target by formula (1); After calculating was finished, finding had 3 goal satisfaction F＜0.05, is area target, and this software keeps these 3 area targets, and all the other targets are all deleted;

Step 5: use the calculating radius of a ball function of rock core 3-D view crack identification software, the minimal circumscribed sphere radius R of 3 area targets that keep in the calculation procedure 4 _MinWith equivalent spheres radius R _Eq, and the ratio of calculating minimal circumscribed sphere radius and equivalent spheres radius

Discovery has 2 goal satisfactions

For possessing malleable area target, this software keeps these 2 area targets that possess ductility, deletion another one target;

Step 6: use the digital simulation plane function of rock core 3-D view crack identification software, by 2 fit Plane that possess the area target of ductility that keep in formula (2) and formula (8) calculation procedure 5.

Step 7: use the minimum external rectangular parallelepiped height function of calculating of rock core 3-D view crack identification software, by 2 external rectangular parallelepiped height of approximate minimum that possess the area target of ductility that keep in formula (1) and formula (9) calculation procedure 5.

Step 8: use the minimum external rectangular parallelepiped length function of calculating of rock core 3-D view crack identification software, by 2 external rectangular parallelepiped length of approximate minimum that possess the area target of ductility that keep in formula (11) calculation procedure 5.

Step 9: use the minimum external rectangular parallelepiped width function of calculating of rock core 3-D view crack identification software, by 2 external wide side's heights of approximate minimum that possess the area target of ductility that keep in formula (13) and formula (14) calculation procedure 5.

Step 10: the minimum external rectangular parallelepiped length of side ratio function of calculating of using rock core 3-D view crack identification software, to step 7, obtain respectively length, the width of approximate minimum external rectangular parallelepiped in the step 8, step 9 and compare highly between any two, obtain the longest edge D among them _MaxMinor face D _Min, calculate

Ratio, find that 2 area targets that possess ductility that keep in the step 5 only have 1 to satisfy

This software keeps it as the crack of identifying, deletion another one target shown in (b) figure among Fig. 3, co-exists in a crack through this rock core of identification.Recorded crack 2 targets that extract through above-mentioned steps and three hole targets---hole 4 wherein in the following table 2, the three-dimensional configuration parameter comparison of hole 5 and hole 6, by formula (1), With

The calculating acquired results as shown in table 2 below.

The three-dimensional configuration parameter comparison of table 2 crack 2 and hole 4, hole 5 and hole 6

	Crack 2	Hole 4	Hole 5	Hole 6
					Form factor	0.00156837	0.205297	0.252467	0.894843
The ratio of two radius	5.528	1.3812	2.6641	1.3476
					The ratio on both sides	4.303	2.18996	1.16752	2.47214

Observe table 2 and find, step 3 is satisfied in crack 2 wherein, and three screening conditions of step 4 and step 5 are successfully identified; And hole 4,6 in hole 5 and hole do not satisfy step 3, and three screening conditions of step 4 and step 5 final are deletedly fallen.

Claims (1)

1. recognition methods based on the rock core crack of three-dimensional image information processing is characterized in that comprising successively following concrete operation step:

Step 1: use the pending rock core of CT machine scanning, obtain rock core CT sequence image and be stored in the computing machine;

Step 2: use rock core 3-D view crack identification software to read in the good rock core CT sequence image of step 1 scanning, and utilize marching cubes algorithm that rock core CT sequence image is carried out three-dimensional image reconstruction, in the 3 d data field that three-dimensional image reconstruction consists of, obtain one or more than one target;

Step 3: each target that three-dimensional image reconstruction obtains in the scanning step 2, obtain the grown form parameter of each target, comprise target volume, target surface area; Calculate the value of the 3D shape factor F of each target, the relational expression of F and described target volume and target surface area is following formula:

$F=\frac{36\mathrm{\π}{V}_p}{S_p}---\left(1\right)$

Wherein, V _pBe target volume, S _pBe target surface area, calculate the value of the F of each target according to formula (1); Again the F value is screened, if satisfy F＜0.05, illustrate that this target is planar distribution, keep this area target, otherwise the non-area target of this target is described, then this non-area target of deletion;

Step 4: each area target that keeps in the scanning step 3, obtain the minimal circumscribed sphere radius R of each area target _MinWith equivalent spheres radius R _Eq, the ratio of calculating minimal circumscribed sphere radius and equivalent spheres radius According to

Ratio screen again, if satisfy

Illustrate that this area target possesses ductility simultaneously, then keep the area target that this possesses ductility, otherwise illustrate that this area target does not possess ductility, then deletes the area target that this does not possess ductility;

Step 5: each that keeps in the calculation procedure 4 possesses the length of the external rectangular parallelepiped of approximate minimum of the area target of ductility, and length and width and height are compared between any two, obtains the longest edge D among them _MaxMinor face D _Min, according to Ratio screen, if satisfy

Illustrate that this area target that is ductile is the crack, it is kept as the crack of identifying, if this area target that is ductile does not satisfy described condition, illustrate that this area target that is ductile is not the crack, then deletes this area target.

Images