CN117274533A - Three-dimensional geometric model construction method and performance prediction method for porous rock - Google Patents
Three-dimensional geometric model construction method and performance prediction method for porous rock Download PDFInfo
- Publication number
- CN117274533A CN117274533A CN202311314325.2A CN202311314325A CN117274533A CN 117274533 A CN117274533 A CN 117274533A CN 202311314325 A CN202311314325 A CN 202311314325A CN 117274533 A CN117274533 A CN 117274533A
- Authority
- CN
- China
- Prior art keywords
- rock
- model
- image
- tetrahedral
- cluster
- 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.)
- Pending
Links
- 239000011435 rock Substances 0.000 title claims abstract description 115
- 238000000034 method Methods 0.000 title claims description 38
- 238000010276 construction Methods 0.000 title abstract description 12
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 33
- 238000002591 computed tomography Methods 0.000 claims abstract description 20
- 239000011148 porous material Substances 0.000 claims abstract description 13
- 230000011218 segmentation Effects 0.000 claims abstract description 13
- 238000003708 edge detection Methods 0.000 claims abstract description 9
- 238000007405 data analysis Methods 0.000 claims abstract description 6
- 230000007246 mechanism Effects 0.000 claims abstract description 5
- 238000012360 testing method Methods 0.000 claims description 16
- 238000004088 simulation Methods 0.000 claims description 15
- 230000003044 adaptive effect Effects 0.000 claims description 9
- 238000002372 labelling Methods 0.000 claims description 9
- 239000011159 matrix material Substances 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 238000005260 corrosion Methods 0.000 claims description 3
- 230000007797 corrosion Effects 0.000 claims description 3
- 230000006872 improvement Effects 0.000 claims description 3
- 238000012669 compression test Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 238000007906 compression Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 238000011068 loading method Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 201000004569 Blindness Diseases 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000000919 ceramic Substances 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000013210 evaluation model Methods 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 230000003746 surface roughness Effects 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
- G06T17/205—Re-meshing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/13—Edge detection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/181—Segmentation; Edge detection involving edge growing; involving edge linking
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- 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/10—Image acquisition modality
- G06T2207/10004—Still image; Photographic image
- G06T2207/10012—Stereo images
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
- G06T2207/30184—Infrastructure
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
The invention provides a three-dimensional geometric model construction method of porous rock, which comprises the following steps: selecting a rock sample, and acquiring an internal scanning image of the rock sample by using an X-ray CT scanner; based on a rock CT image, based on a total variation image denoising model, acquiring self-adaptive parameters by utilizing an edge detection operator, so that CT scanning image noise is removed, and edge contour and texture information of the CT scanning image are reserved; performing threshold segmentation on the CT scanning image by adopting an iterative self-organizing data analysis threshold segmentation algorithm to obtain a binary image which accords with the interior of a real rock and comprises a rock mechanism and pores; and reconstructing a tetrahedral model of the rock, namely a rock finite element reconstruction grid model by taking the established surface model as an input domain and adopting a Delong tetrahedral algorithm according to boundary limiting conditions in the surface model. The geometric model construction method can accurately reflect the real characteristics of the rock and provides the related model to predict the mechanical behavior of the rock.
Description
Technical Field
The invention relates to the technical field of petroleum and natural gas exploration and development, in particular to a three-dimensional geometric model construction method and a performance prediction method of porous rock.
Background
The mechanical properties of the rock are necessary foundations for engineering design and implementation such as safe drilling, high-efficiency fracturing and the like; the structural characteristics of rock are important factors affecting the mechanical properties of rock. Therefore, the scientific construction of the rock mechanical parameter evaluation model with the complex structure is a precondition and foundation for rock mechanical property evaluation, and has important significance for the safe and efficient exploration and development of oil and gas reservoirs.
Natural rock contains a large number of discrete, multi-scale and irregular pores, and the microstructure of the rock has a direct relationship with not only the porosity and permeability of the rock, but also an indirect relationship with young's modulus, poisson's ratio, uniaxial compressive strength, etc. However, in conventional numerical rock modeling, the microstructure in the rock is often eliminated or simplified, and thus calculation errors may be caused when pores or surface roughness of the microstructure exist inside the rock. The simplified measuring method has the advantages of shallow research on the mechanical properties of rock force and large error.
Therefore, there is a need for a more accurate method of rock geometric model construction that more accurately reflects the true characteristics of the rock and more accurately predicts its mechanical behavior on the rock.
Disclosure of Invention
In view of the above, the present invention aims to provide a three-dimensional geometric model construction method and a performance prediction method for porous rock. The invention aims to solve the problems of shallow research on the rock mechanical property and large error of the existing method.
In order to achieve the above purpose, the invention provides a three-dimensional geometric model construction method of porous rock, comprising the following steps:
s1, selecting a rock sample, and acquiring an internal scanning image of the rock sample by using an X-ray CT scanner;
s2, based on the CT image obtained in the step S1, based on a total variation image denoising model, obtaining self-adaptive parameters by utilizing an edge detection operator, so that the noise of the CT scanning image is removed, and the edge contour and texture information of the CT scanning image are reserved;
s3, based on the CT image obtained in the step S1, performing threshold segmentation on the CT scanning image by adopting an iterative self-organizing data analysis threshold segmentation algorithm so as to obtain a binary image which is consistent with the interior of the real rock and contains a rock mechanism and pores;
s4, reconstructing a tetrahedral model of the rock, namely a rock finite element reconstruction grid model, by taking the established surface model as an input domain and adopting a Delong tetrahedral algorithm according to boundary limiting conditions in the surface model.
Further, the expression of the adaptive parameter ζf obtained by using the edge detection operator in the step S2 is as follows:
wherein:represents a convolution operator, d θ Represents an edge detection operator, d θ Comprises d 0 、d 45 、d 90 And d 135 ；
Wherein:the rest of the edge operators can be obtained by rotating θ∈θ (θ= {45 °,90 °,135 ° }).
Further, the expression of the CT scan image noise in step S2 is as follows:
wherein: u is the original image; f is a total variation image denoising model; x is e.OMEGA, where OMEGA represents image space; μ >0 is a regularization parameter.
Further, the specific algorithm flow in step S4 is as follows:
step 1: the following initial control parameters are determined: similarity function D ij (X i ,X j ) Number of cluster centers k, number of cluster centers N c Minimum number of samples θ N Standard deviation theta of sample distance distribution S Minimum distance theta of clustering center C Iterative operation times I, clustering center logarithm L;
step 2: the following loop is performed
For i＝1:n
For j＝i:n
M(i,j)＝||1i-1j||
end
end
Wherein 1 is a feature vector of pixels, n represents the number of pixels, and M is a matrix for storing the distance between each pixel point;
using average value of square distance between every two pixel points in matrix, i.e. average clustering distance ave, to obtain R 1 And R is R 2 Wherein R is 1 =ave×x (x ranges from 0.5 to 2), R 2 ＝2*R 1 The method comprises the steps of carrying out a first treatment on the surface of the Taking each pixel point as a circle center, taking R as a circle center 1 Grouping pixels for radius, counting the number of pixels in each group, sorting each group from large to small, and taking the largest number of groups as an initial clustering center C 1 Judging the pixel point group arranged at the second positionIf the distance between the first center and the second center is larger than R 2 Then take it as the second initial cluster center C 2 Performing similar judgment on the rest pixel point groups until no new cluster center appears, and then outputting an initial cluster center and cluster number;
step 3: sample NAssigning to nearest clusters S j The classification criteria are as follows:
if it isWhen->Is the smallest, then->
The clustering center value is modified according to the following formula:
the cluster domains S are calculated according to the following formula j Average distance of each sample pixel to the cluster center:
the total average distance of all sample pixels to their respective cluster centers is calculated according to the following formula:
step 4: if N C K/2, performing a splitting step, and calculating the standard of the sample clusters in each cluster according to the formula shown belowDifference vector sigma j ：
σ j ＝(σ 1j ,σ 2j ,...σ nj ) T
Finding the maximum value sigma in the standard deviation vector jmax If sigma is satisfied jmax >θ S And meet the followingOr N > 2 (theta) N +1), then its corresponding cluster center is +.>Split into two new cluster centers +.>And->
If the distance between the two cluster centers is smaller than theta C I.e.
Will D ij The values of (2) are arranged in descending order and are combined according to the following formula:
and from N C Is subtracted byAnd->
Step 5: if it isIf t is less than l, returning to the step 3, if +.>The algorithm is finished, and a cluster center { C } is output 1 ,C 2 ,...,C k }
Step 6: preliminary labeling is performed in the image according to the following formula
The labeling is performed again after the labeling is completed, and in addition, the expansion and corrosion methods in the image morphology are needed to process each segmented region.
Further, the specific flow of the delouli tetrahedral algorithm in step S4 is as follows:
step 1: defining a curved surface model as a convex domain containing all points to be inserted;
step 2: extracting constraint information of a boundary edge set and a triangle set in the curved surface model;
step 3: inserting one vertex in an initial tetrahedron vertex set (ITV) into the convex domain, deleting elements violating the empty boundary sphere criterion to form a cavity, and generating tetrahedron elements by using the inserted vertex of the cavity and the surface triangle;
step 4: storing the generated tetrahedral element, repeating the step 3, and inserting all different vertexes into the convex domain to obtain an Initial Tetrahedral Set (ITS);
step 5: generating an Interior Point Set (IPS) by adopting an adaptive algorithm, and inserting the IPS into the ITS to form a Delould tetrahedral set;
step 6: restoring boundary edges and boundary triangles;
step 7: tetrahedral elements outside the domain are eliminated;
step 8: adopting a quality improvement algorithm to improve the quality of the Deluxe Tetrahedral Set (DTS);
step 9: and removing the slider elements in the deluo inner tetrahedral set (DTS) to finally obtain a finite element tetrahedral model, namely a rock finite element reconstruction grid model.
The invention also provides a method for predicting the performance of the porous rock, which utilizes the finite element reconstruction grid model of the rock, obtains the functional relation between the mechanical property of the rock and the random pore structure through a numerical simulation test, and combines the obtained functional relation with the reconstruction model to numerically simulate and predict the performance of the rock in engineering.
The invention has the beneficial effects that:
1. the invention provides a three-dimensional geometric model construction method of porous rock, which adopts a denoising algorithm based on total variation, and compared with other denoising algorithms, the method has the advantages that the details of images are reserved, and meanwhile, the denoising effect is better. Meanwhile, in the invention, the impact filter and the anisotropic diffusion filter are combined into the self-adaptive total variation denoising algorithm, and the edge operator detection method is used for distinguishing edge details and the region to be smoothed, so that the self-adaptive parameters for controlling the diffusion degree of the variation model can be obtained. The sensitivity of the adaptive parameters to noise is reduced, thereby enabling application of morphological principles of image processing to the denoising model. More importantly, the method avoids the easy-to-occur piecewise constant effect while realizing smooth denoising and detail protection of the image.
2. The invention provides a three-dimensional geometric model construction method of porous rock, which adopts an improved iterative self-organizing data analysis threshold segmentation algorithm, can automatically determine the optimal clustering center and the optimal clustering number, and avoids randomness and blindness in the traditional algorithm so as to obtain a binary image which accords with the interior of real rock and comprises a rock mechanism and pores.
3. The invention provides a performance prediction method of porous rock, which researches the relation between the porosity and the mechanical property of the rock in detail, mainly the relation between the porosity and the compressive strength of the rock; according to the method, the rock finite element reconstruction grid model is utilized, the functional relation between the rock mechanical property and the random pore structure is obtained through a numerical simulation test, and the performance of the rock in engineering can be simulated and predicted through the obtained functional relation and the reconstruction model numerical value.
Additional advantages, objects, and features of the invention will be set forth in part in the description which follows and in part will become apparent to those having ordinary skill in the art upon examination of the following or may be learned from practice of the invention. The objects and other advantages of the invention may be realized and obtained by means of the instrumentalities and combinations particularly pointed out in the specification.
Drawings
FIG. 1 is a flow chart of a method for constructing a three-dimensional geometric model of porous rock according to the present invention;
FIG. 2 is a raw CT image of a Chongqing sandstone sample of example 1;
FIG. 3 is a schematic diagram showing the results of the mechanical compression test of Chongqing sandstone sample in example 1;
FIG. 4 is a graph comparing rock CT scan images before and after denoising;
FIG. 5 is a graph comparing rock CT scan images before and after threshold segmentation;
FIG. 6 is a finite element model of rock reconstruction;
FIG. 7 is a graph comparing stress-strain curves obtained by data simulation and test;
FIG. 8 is a partial three-dimensional reconstruction model of different porosities of a rock;
FIG. 9 is a theoretical model and numerical simulation results;
FIG. 10 shows the results of the model and numerical simulation according to the present invention.
Detailed Description
In order to make the technical scheme, advantages and objects of the present invention more clear, the technical scheme of the embodiment of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiment of the present invention. It will be apparent that the described embodiments are some, but not all, embodiments of the invention. All other embodiments, which can be obtained by a person skilled in the art without creative efforts, based on the described embodiments of the present invention belong to the protection scope of the present application.
As shown in fig. 1, the invention provides a three-dimensional geometric model construction method of porous rock, which comprises the following steps:
s1, selecting a rock sample, and acquiring an internal scanning image of the rock sample by using an X-ray CT scanner;
s2, based on the CT image obtained in the step S1, based on a total variation image denoising model, obtaining self-adaptive parameters by utilizing an edge detection operator, so that the noise of the CT scanning image is removed, and the edge contour and texture information of the CT scanning image are reserved;
s3, based on the CT image obtained in the step S1, performing threshold segmentation on the CT scanning image by adopting an iterative self-organizing data analysis threshold segmentation algorithm so as to obtain a binary image which is consistent with the interior of the real rock and contains a rock mechanism and pores;
s4, reconstructing a tetrahedral model of the rock, namely a rock finite element reconstruction grid model, by taking the established surface model as an input domain and adopting a Delong tetrahedral algorithm according to boundary limiting conditions in the surface model.
The invention also provides a performance prediction method of the porous rock, which utilizes the constructed rock finite element reconstruction grid model to obtain the functional relation between the rock mechanical property and the random pore structure through a numerical simulation test, and combines the obtained functional relation with the reconstruction model to numerically simulate and predict the performance of the rock in engineering.
Example 1
The rock sample selected in this embodiment is Chongqing sandstone, and the specific steps are as follows:
first, selecting a sandstone sample, and acquiring an internal scanning image of the sandstone sample by using an X-ray CT scanner.
The X-ray CT scanner of this embodiment is a somatim oscilloscope manufactured by SIEMENS corporation. Basic working principle of SOMATOM oscilloscope: in the working process of the CT device, the X-ray source continuously emits X-ray beams which can penetrate through a sample on the sample table, the emitted X-rays attenuate to a certain extent according to the components of the sample, and the attenuated X-rays are received by the receiver to serve as original information of imaging and reconstruction. In this process, the detector and receiver are rotated to acquire raw data or image data from various directions. The whole process is monitored by a Syngo scope one console.
The X-ray tube focal spot size of the somom oscilloscope is 0.8X 0.5mm. At 2% Modulation Transfer Function (MTF), the high contrast resolution is 15.6lp/cm. The X-ray beam penetrating the rock sample is measured by 24 line detectors. The SOMATOM oscilloscope provides excellent detailed images using a Siemens-specific ultra-fast ceramic detector. The electron current was 81mA, the acceleration voltage was 130kV, and the scan time was 1500ms. The reconstruction matrix of the Syngo scope one console is 512×512. The rock sample was scanned by a somom oscilloscope to obtain CT images as shown in fig. 2.
And secondly, carrying out mechanical test on the rock sample, thereby obtaining the physical and mechanical parameters of the rock sample.
First, the mechanical properties of Chongqing sandstone were determined experimentally. The mechanical properties of the reconstructed or reference model are then determined by a numerical program. Finally, the mechanical properties obtained by numerical simulation are compared with those obtained by experiments.
The mechanical test is classified into a uniaxial compression test and a triaxial compression test.
Uniaxial compression test: the Chongqing sandstone sample C1 was uniaxially compressed using a hydraulically controlled uniaxial loading system to obtain its Uniaxial Compressive Strength (UCS), and the test results are shown in FIG. 3 (a). The diameter and height of the cylindrical rock sample were 50mm and 100mm, respectively. The UCS of Chongqing sandstone sample C1 is 103.684MPa.
Triaxial compression test: by using the multifunctional triaxial rock test system 600-50HT Plus (compared with other triaxial rock test systems, the multifunctional triaxial rock test system adopts the self-adaptive confining pressure compensation technology developed by university e Lille 1, so that the influence of piston motion on confining pressure can be eliminated, the confining pressure acting on a rock sample is ensured to be stable, in an experiment, the maximum axial stress can reach 1000KN, the maximum confining pressure can be loaded to 60 MPa), and triaxial compression is carried out on the counter-weight rock samples C2, C3 and C4, and the experimental result is shown in figures 3 (b), (C) and (d). In order to obtain the mechanical parameters of Chongqing sandstone, namely cohesion c and internal friction angle, at least three different confining pressures are required to obtain the corresponding representative molar coulomb circles and breaking envelope. The uniaxial compressive strength of the rock can then be calculated using a linear molar coulomb relationship. The whole process of triaxial compression of the rock accords with the loading method recommended by the international rock mechanics society. In an example, fifteen different confining pressures were applied to the rock triaxial experiment from 2MPa to 30MPa in turn, with each different rock sample confining pressure differing by 2MPa. Young's modulus at different confining pressures may be calculated from the slope of the stress-strain curve at different confining pressures, while Poisson's ratio may be calculated from the ratio of axial strain and hoop strain. The cohesion and internal friction angle can be calculated from the slope and intercept of the tangent line according to the molar coulomb circle.
And thirdly, based on the CT image obtained in the first step, based on a total variation image denoising model, utilizing an edge detection operator to obtain self-adaptive parameters, thereby removing the noise of the CT scanning image, retaining the edge contour and texture information of the CT scanning image, and enabling the processed CT image to reflect the spatial distribution rule of pores in the rock as accurately as possible, as shown in fig. 4.
f＝u+n f (1)
u is the original image, f is the image degraded by noise, and n is the gaussian noise white with zero mean and standard deviation δ.
Four edge operators are introduced to calculate the adaptive parameter ζ f When θ=0, respectively:
the rest of the edge operators can be obtained by rotating θ∈θ (θ= {45 °,90 °,135 ° }). Adaptive parameters ζ according to the definition of the edge operator f Can be written as follows:
the method noise can be written as:
wherein: u is the original image; f is a total variation image denoising model; x is e.OMEGA, where OMEGA represents image space; μ >0 is a regularization parameter.
In the edge region of the image,
at mu>In the case of 0, the number of the cells,therefore, the traditional full-variance model can remove noise and blur edge detail information, but the self-adaptive full-variance model can remove noise of an image and retain the edge detail information. The bragg gate algorithm is used to solve quickly while also avoiding piecewise constant effects.
For the case of k >0 and c e IR, the operator cut can be defined as follows:
let u be 1 ＝f,And c=ζ f Let λ >0, respectively, be the auxiliary variables, then equation (5) can be equated to a fast iterative algorithm expressed as follows:
jia has demonstrated that for k=0, 1..There is->Because ofTherefore->
For any C E (0, 1)]，Can be determined, i.e. an optimal solution can be determined, from which an optimal denoising scheme can be obtained.
And fourthly, based on the CT image obtained in the step S1, performing threshold segmentation on the rock CT scanning image by adopting an iterative self-organizing data analysis threshold segmentation algorithm, as shown in fig. 5.
Step 1: the following initial control parameters are determined: e.g. similarity function D ij (X i ,X j ) Number of cluster centers k, number of cluster centers N c Minimum number of samples θ N Standard deviation theta of sample distance distribution S Minimum distance theta of clustering center C Iterative operation times I, clustering center logarithm L;
step 2: the following loop is performed
For i＝1:n
For j＝i:n
M(i,j)＝||1i-1j||
end
end
Wherein 1 is a feature vector of pixels, n represents the number of pixels, and M is a matrix for storing the distance between each pixel point;
using average value of square distance between every two pixel points in matrix, i.e. average clustering distance ave, to obtain R 1 And R is R 2 Wherein R is 1 =ave×x (x ranges from 0.5 to 2), R 2 ＝2*R 1 The method comprises the steps of carrying out a first treatment on the surface of the Taking each pixel point as a circle center, taking R as a circle center 1 As a radius pairGrouping the pixels, counting the number of the pixels in each group, sequencing each group from large to small, and taking the group with the largest number as an initial clustering center C 1 Judging the pixel point group arranged at the second position, if the distance between the pixel point group and the first center is larger than R 2 Then take it as the second initial cluster center C 2 Performing similar judgment on the rest pixel point groups until no new cluster center appears, and then outputting an initial cluster center and cluster number;
step 3: sample NAssigning to nearest clusters S j The classification criteria are as follows:
if it isWhen->Is the smallest, then->
The clustering center value is modified according to the following formula:
the cluster domains S are calculated according to the following formula j Average distance of each sample pixel to the cluster center:
the total average distance of all sample pixels to their respective cluster centers is calculated according to the following formula:
step 4: if N C If k/2 is less than or equal to, performing a splitting step, and calculating a standard deviation vector sigma of sample clusters in each cluster according to the following formula j ：
σ j ＝(σ 1j ,σ 2j ,...σ nj ) T
Finding the maximum value sigma in the standard deviation vector jmax If sigma is satisfied jmax >θ S And meet the followingOr N > 2 (theta) N +1), then its corresponding cluster center is +.>Split into two new cluster centers +.>And->
If the distance between the two cluster centers is smaller than theta C I.e.
Will D ij The values of (2) are arranged in descending order and are combined according to the following formula:
and from N C Is subtracted byAnd->
Step 5: if it isIf t is less than l, returning to the step 3, if +.>The algorithm is finished, and a cluster center { C } is output 1 ,C 2 ,...,C k }
Step 6: preliminary labeling is performed in the image according to the following formula
The labeling is performed again after the labeling is completed, and in addition, the expansion and corrosion methods in the image morphology are needed to process each segmented region.
And fifthly, reconstructing a tetrahedral model of the rock, namely reconstructing a grid model of the finite element of the rock by taking the established surface model as an input domain and adopting a delould tetrahedral algorithm according to boundary limiting conditions in the surface model.
The specific flow of the deluo inner tetrahedral algorithm in the step S4 is as follows:
step 1: defining a curved surface model as a convex domain containing all points to be inserted;
step 2: extracting constraint information of a boundary edge set and a triangle set in the curved surface model;
step 3: inserting one vertex in an initial tetrahedron vertex set (ITV) into the convex domain, deleting elements violating the empty boundary sphere criterion to form a cavity, and generating tetrahedron elements by using the inserted vertex of the cavity and the surface triangle;
step 4: storing the generated tetrahedral element, repeating the step 3, and inserting all different vertexes into the convex domain to obtain an Initial Tetrahedral Set (ITS);
step 5: generating an Interior Point Set (IPS) by adopting an adaptive algorithm, and inserting the IPS into the ITS to form a Delould tetrahedral set;
in step 5, adding some self-adaptive points into the finite element tetrahedral model to generate internal tetrahedral elements. At the same time, the added points should avoid co-sphericity of more than four points. The adaptive algorithm is summarized as follows: 1. calculating diameter spheres of triangle elements in the curved surface model to define a Protected Area Set (PAS); 2. generating an Initial Interior Point Set (IIPS) using a grid-based method; 3. carrying out a random arrangement algorithm on the IIPS to avoid more than 4 point co-sphericity; 4. checking the IIPS to eliminate points located in the PAS to obtain an Interior Point Set (IPS); 5. IPS is inserted into the Initial Tetrahedral Set (ITS), obtaining the Delustered Tetrahedral Set (DTS).
Step 6: restoring boundary edges and boundary triangles;
step 7: tetrahedral elements outside the domain are eliminated;
step 8: adopting a quality improvement algorithm to improve the quality of the DTS;
step 9: and removing the slider elements in the DTS to finally obtain a finite element tetrahedral model, namely a rock finite element reconstruction grid model.
And sixthly, utilizing the constructed rock finite element reconstruction grid model, obtaining a functional relation between the rock mechanical property and the random pore structure through a numerical simulation test, and predicting the performance of the rock in engineering through the numerical simulation of the obtained functional relation and the reconstruction model.
Step 1, model verification
The test sample was three-dimensionally reconstructed by the above-described reconstruction method, as shown in fig. 6, fig. 6 (a) is a schematic structural view of the reconstructed finite element model, and fig. 6 (b) is a perspective view of the reconstructed finite element model. And then mechanical parameters obtained through a mechanical test are substituted into model attributes as initial calibration parameters, and numerical simulation calculation is carried out through ABAQUS. The stress-strain curve obtained by numerical simulation can be highly matched with the stress-strain curve obtained by test, as shown in fig. 7, so that the reconstructed model can safely predict the mechanical properties of the porous sandstone.
Step 2, numerical simulation test
After denoising and threshold segmentation are carried out on images obtained by CT scanning of 19 rock samples, respectively reconstructing finite element models with the porosities of 3.027%,3.979%,5.945%,7.396%,8.692%,9.385%,10.021%,11.312%,11.934%,12.582%,13.243%,13.972%,14.747%,15.538%,16.341%,17.153%,18.107%,19.030% and 20.261%, and establishing a finite element model with the same size with the porosity of 0, as shown in fig. 8, carrying out numerical simulation uniaxial compression test to obtain the respective compressive strengths, wherein fig. 8 (a) is a finite element model with the porosity of 0, fig. 8 (b) is a three-dimensional reconstruction model with the porosity of 3.979%, fig. 8 (c) is a three-dimensional reconstruction model with the porosity of 15.538%, and fig. 8 (d) is a three-dimensional reconstruction model with the porosity of 20.261%.
The numerical results were fitted with the proposed power model and exponential model, respectively, as shown in fig. 9, with correlation coefficients 0.98852 and 0.97995, respectively.
In view of this, a combined model between compressive strength and porosity for porous rock was proposed in combination with the two models described above:
σ＝σ 0 γ p -3pe -2.5p
wherein σ represents the compressive strength of the porous rock; sigma (sigma) 0 Represents compressive strength at a rock porosity of 0; gamma is a control factor, and the value range is 0.98322 +/-9.8639 multiplied by 10 -5 The method comprises the steps of carrying out a first treatment on the surface of the p represents the porosity of the porous rock.
The correlation coefficient of the compressive strength and porosity combined model described above reached 0.99453, as shown in fig. 10, which suggests that the model can accurately describe the correlation between compressive strength and porosity of an acceptable porous rock. I.e. the model proposed by the present invention is reliable for predicting the compressive strength of porous rock.
Finally, it is noted that the above embodiments are only for illustrating the technical solution of the present invention and not for limiting the same, and although the present invention has been described in detail with reference to the preferred embodiments, it should be understood by those skilled in the art that the technical solution of the present invention may be modified or substituted without departing from the spirit and scope of the technical solution, and the present invention is intended to be covered in the scope of the present invention.
Claims (6)
1. The method for constructing the three-dimensional geometric model of the porous rock is characterized by comprising the following steps of:
s1, selecting a rock sample, and acquiring an internal scanning image of the rock sample by using an X-ray CT scanner;
s2, based on the CT image obtained in the step S1, based on a total variation image denoising model, obtaining self-adaptive parameters by utilizing an edge detection operator, so that the noise of the CT scanning image is removed, and the edge contour and texture information of the CT scanning image are reserved;
s3, based on the CT image obtained in the step S1, performing threshold segmentation on the CT scanning image by adopting an iterative self-organizing data analysis threshold segmentation algorithm so as to obtain a binary image which is consistent with the interior of the real rock and contains a rock mechanism and pores;
s4, reconstructing a tetrahedral model of the rock, namely a rock finite element reconstruction grid model, by taking the established surface model as an input domain and adopting a Delong tetrahedral algorithm according to boundary limiting conditions in the surface model.
2. The method for constructing a three-dimensional geometric model of porous rock according to claim 1, wherein in step S2, the adaptive parameter ζ is obtained by using an edge detection operator f The expression of (2) is as follows:
wherein:represents a convolution operator, d θ Represents an edge detection operator, d θ Comprises d 0 、d 45 、d 90 And d 135 ；
The rest of the edge operators can be obtained by rotating θ∈θ (θ= {45 °,90 °,135 ° }).
3. The method for constructing a three-dimensional geometric model of porous rock according to claim 2, wherein the expression of the CT scan image noise in step S2 is as follows:
wherein: u is the original image; f is a total variation image denoising model; x is e.OMEGA, where OMEGA represents image space; μ >0 is a regularization parameter.
4. The method for constructing a three-dimensional geometric model of porous rock according to claim 1, wherein the specific algorithm flow in step S4 is as follows:
step 1: the following initial control parameters are determined: similarity function D ij (X i ,X j ) Number of cluster centers k, number of cluster centers N c Minimum number of samples θ N Standard deviation theta of sample distance distribution S Minimum distance theta of clustering center C Iterative operation times I, clustering center logarithm L;
step 2: the following loop is performed
For i＝1:n
For j＝i:n
M(i,j)＝||1i-1j||
end
end
Wherein 1 is a feature vector of pixels, n represents the number of pixels, and M is a matrix for storing the distance between each pixel point;
using the average of the square distances between all pixels in the matrix, i.e. average clusteringDistance ave, find R 1 And R is R 2 Wherein R is 1 =ave×x (x ranges from 0.5 to 2), R 2 ＝2*R 1 The method comprises the steps of carrying out a first treatment on the surface of the Taking each pixel point as a circle center, taking R as a circle center 1 Grouping pixels for radius, counting the number of pixels in each group, sorting each group from large to small, and taking the largest number of groups as an initial clustering center C 1 Judging the pixel point group arranged at the second position, if the distance between the pixel point group and the first center is larger than R 2 Then take it as the second initial cluster center C 2 Performing similar judgment on the rest pixel point groups until no new cluster center appears, and then outputting an initial cluster center and cluster number;
step 3: sample NAssigning to nearest clusters S j The classification criteria are as follows:
if it isWhen->Is the smallest, then->
The clustering center value is modified according to the following formula:
the cluster domains S are calculated according to the following formula j Average distance of each sample pixel to the cluster center:
the total average distance of all sample pixels to their respective cluster centers is calculated according to the following formula:
step 4: if N C If k/2 is less than or equal to, performing a splitting step, and calculating a standard deviation vector sigma of sample clusters in each cluster according to the following formula j ：
σ j ＝(σ 1j ,σ 2j ,...σ nj ) T
Finding the maximum value sigma in the standard deviation vector jmax If sigma is satisfied jmax >θ S And meet the followingOr N > 2 (theta) N +1), then its corresponding cluster center is +.>Split into two new cluster centers +.>And->
If the distance between the two cluster centers is smaller than theta C I.e.
Will D ij The values of (2) are arranged in descending order, according toThe following formula is combined:
and from N C Is subtracted byAnd->
Step 5: if it isReturn to step 3 if->The algorithm is finished, and a cluster center { C } is output 1 ,C 2 ,...,C k }
Step 6: preliminary labeling is performed in the image according to the following formula
The labeling is performed again after the labeling is completed, and in addition, the expansion and corrosion methods in the image morphology are needed to process each segmented region.
5. The method for constructing a three-dimensional geometric model of porous rock according to claim 1, wherein the specific flow of the deluster tetrahedral algorithm in step S4 is as follows:
step 1: defining a curved surface model as a convex domain containing all points to be inserted;
step 2: extracting constraint information of a boundary edge set and a triangle set in the curved surface model;
step 3: inserting one vertex in an initial tetrahedron vertex set (ITV) into the convex domain, deleting elements violating the empty boundary sphere criterion to form a cavity, and generating tetrahedron elements by using the inserted vertex of the cavity and the surface triangle;
step 4: storing the generated tetrahedral element, repeating the step 3, and inserting all different vertexes into the convex domain to obtain an Initial Tetrahedral Set (ITS);
step 5: generating an Interior Point Set (IPS) by adopting an adaptive algorithm, and inserting the IPS into the ITS to form a Delould tetrahedral set;
step 6: restoring boundary edges and boundary triangles;
step 7: tetrahedral elements outside the domain are eliminated;
step 8: adopting a quality improvement algorithm to improve the quality of the Deluxe Tetrahedral Set (DTS);
step 9: and removing the slider elements in the deluo inner tetrahedral set (DTS) to finally obtain a finite element tetrahedral model, namely a rock finite element reconstruction grid model.
6. A performance prediction method of porous rock is characterized in that: the rock finite element reconstruction grid model constructed by the method of claims 1-5 is utilized to obtain a functional relation between the rock mechanical property and the random pore structure through a numerical simulation test, and the performance of the rock in engineering is predicted through the obtained functional relation and the numerical simulation of the reconstruction model.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311314325.2A CN117274533A (en) | 2023-10-11 | 2023-10-11 | Three-dimensional geometric model construction method and performance prediction method for porous rock |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311314325.2A CN117274533A (en) | 2023-10-11 | 2023-10-11 | Three-dimensional geometric model construction method and performance prediction method for porous rock |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117274533A true CN117274533A (en) | 2023-12-22 |
Family
ID=89221325
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311314325.2A Pending CN117274533A (en) | 2023-10-11 | 2023-10-11 | Three-dimensional geometric model construction method and performance prediction method for porous rock |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117274533A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117557742A (en) * | 2024-01-12 | 2024-02-13 | 武汉大学 | 3D rock reservoir modeling method based on digital image and machine learning |
-
2023
- 2023-10-11 CN CN202311314325.2A patent/CN117274533A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117557742A (en) * | 2024-01-12 | 2024-02-13 | 武汉大学 | 3D rock reservoir modeling method based on digital image and machine learning |
| CN117557742B (en) * | 2024-01-12 | 2024-03-22 | 武汉大学 | 3D rock reservoir modeling method based on digital image and machine learning |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107449707B (en) | Three-dimensional characterization determination method and device for quantification of pores with different scales in shale reservoir | |
| US10223782B2 (en) | Digital rock physics-based trend determination and usage for upscaling | |
| JP6884104B2 (en) | A learning contour identification system that uses portable contour metrics derived from contour mapping | |
| US8170799B2 (en) | Method for determining in-situ relationships between physical properties of a porous medium from a sample thereof | |
| CN117274533A (en) | Three-dimensional geometric model construction method and performance prediction method for porous rock | |
| Zhao et al. | Digital measurement of 2D and 3D cracks in sandstones through improved pseudo color image enhancement and 3D reconstruction method | |
| Zhang et al. | GPU-accelerated 3D reconstruction of porous media using multiple-point statistics | |
| US20140052420A1 (en) | Digital Rock Analysis Systems and Methods that Estimate a Maturity Level | |
| Zhou et al. | Analysis of fracture properties of three-dimensional reconstructed rock model using hierarchical-fractal annealing algorithm | |
| US11530997B2 (en) | Material properties from two-dimensional image | |
| Li et al. | Reconstructing the 3D digital core with a fully convolutional neural network | |
| Jouini et al. | Texture segmentation of 3D x-ray micro-computed tomography images using U-NET | |
| Kahl et al. | Microfabric and anisotropy of elastic waves in sandstone–An observation using high-resolution X-ray microtomography | |
| Ziabari et al. | Simurgh: A Framework for CAD-Driven Deep Learning Based X-Ray CT Reconstruction | |
| Vengrinovich | Bayesian image and pattern reconstruction from incomplete and noisy data | |
| Suuronen et al. | Enhancing industrial X-ray tomography by data-centric statistical methods | |
| Wu et al. | Damage evolution characteristics of 3D reconstructed bedding-containing shale based on CT technology and digital image processing | |
| Al-Mestarehi et al. | Creating a Complete Model of the Wooden Pattern from Laser Scanner Point Clouds Using Alpha Shapes. | |
| CN115808352B (en) | Rock fracture extraction and complexity characterization method and device based on digital rock core | |
| Zhang et al. | 3D pore space reconstruction using deep residual deconvolution networks | |
| RU2774959C1 (en) | Method for determining filtration properties of non-homogeneous porous samples | |
| Lavrukhin et al. | The Influence of Image Morphology on Neural Network-Based Segmentation Results. | |
| He | Reconstruction of 3D microstructure of the rock sample basing on the CT images | |
| CN116930219A (en) | Rock pore throat parameter analysis method, throat length prediction method, device and equipment | |
| CN118150441A (en) | Method for evaluating rock micro-pore structure based on fluid flow characteristics |
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 |