CN114299098A - Rapid rock permeability analysis method based on throat characteristics of scanned image - Google Patents

Abstract

The invention discloses a rock permeability rapid analysis method based on scanned image throat characteristics, which comprises the steps of carrying out image scanning on a rock sample to obtain a rock pore phase and a solid skeleton phase, and identifying pore and throat distribution in the pore phase by combining a central axis of the pore phase and a Euclidean distance matrix of the pore phase; calculating the maximum and minimum throat diameters based on the throat distribution in the void phase, and calculating the fractal dimension of the throat by a box counting method; and finally, establishing the minimum throat diameter, the maximum throat diameter and the fractal dimension of a throat system, establishing a rock sample permeability analysis expression and accurately predicting the rock sample permeability. The method solves the problems of large calculation amount and long time consumption of the existing rock digital imaging permeability analysis, and provides a basis for accurate characterization of the flow capacity of oil gas in the rock.

Description

Rapid rock permeability analysis method based on throat characteristics of scanned image

Technical Field

The invention belongs to the technical field of oil and gas development, and particularly relates to a rapid rock permeability analysis method for throat characteristics of a scanned image.

Background

Rock permeability is an important parameter for accurately predicting the flow capacity of oil gas in underground rocks and guiding efficient development of the oil gas.

At present, rock permeability is generally measured through a laboratory physical experiment, but the laboratory physical experiment method has long test time and low accuracy on a low-porosity rock sample, and the popularization and application of the method for testing the rock sample permeability through the physical experiment method are limited to a great extent.

The rock scanning imaging method provides a platform for digitally analyzing rock flow properties. At present, the rock permeability testing method mainly focuses on analyzing the rock permeability through rock scanning imaging, and the method mainly comprises the steps of directly solving a Navisotke equation, a lattice Boltzmann method and a pore network model method. The method for directly solving the NaviStokes equation and the lattice Boltzmann method needs to divide high-precision grids for calculation of a rock digital imaging gap phase region, needs to develop flow simulation by means of high-performance super-calculation equipment, and is huge in calculation amount and long in time consumption. Compared with a method for directly solving the NavStokes equation and a method for lattice Boltzmann, the pore network model method has higher calculation speed, but because the pore structure regularization representation is carried out on the rock pore phase, the rock permeability analysis precision of the method is low.

Since the rock void space is comprised of open pores and a narrow throat, rock permeability is primarily controlled by the narrow throat structural characteristics. Based on the method, a rock permeability rapid analysis method is established from the characteristics of the scanned image throat, and the method is convenient to popularize and apply.

Disclosure of Invention

Based on the method, the invention provides a rapid rock permeability analysis method based on throat characteristics of a scanned image, which comprises the steps of scanning a rock sample, carrying out binarization processing, and identifying rock gaps and solid skeleton phases; obtaining a void phase centering axis by continuously removing a void phase boundary voxel, and marking the coordinate position of a connection point on the void phase centering axis representing the center of the pore; calculating the distance between each point in the void phase and the wall surface of the framework to obtain a void phase Euclidean distance matrix; dividing the space phase Euclidean distance matrix by using a watershed algorithm to obtain space phase distribution after a local area is divided; finding a connecting point position on a central axis representing the center of a pore on a void phase after a local region is divided by a watershed algorithm, and identifying the distribution of the pore and the throat in the void phase; calculating the maximum and minimum throat diameters based on the throat distribution in the void phase, and calculating the fractal dimension of the throat by adopting a box counting method; and establishing a rock sample permeability analysis expression based on the minimum throat diameter, the maximum throat diameter and the fractal dimension of the throat system, and accurately predicting the rock sample permeability. The method solves the problems of large calculation amount and long time consumption of permeability analysis based on rock digital imaging at present, accurately and quickly analyzes the permeability of the rock based on the throat characteristics of the scanned image, and provides a basis for accurately characterizing the flowing capacity of oil gas in underground rocks.

Specifically, the rapid rock permeability analysis method provided by the invention comprises the following steps:

s1, carrying out image scanning on the rock sample, carrying out binarization processing on the three-dimensional image of the rock sample, and identifying a rock gap phase and a solid skeleton phase;

s2, obtaining a void phase centering axis by continuously removing rock void phase boundary voxels, and marking the coordinate position of a connection point on the void phase centering axis representing the center of the pore;

s3, calculating the distance between each point in the rock void phase and the wall surface of the solid skeleton phase to obtain a void phase Euclidean distance matrix;

s4, segmenting the space phase Euclidean distance matrix by adopting a watershed algorithm to obtain space phase distribution after a local area is divided; finding a connecting point position on a central axis representing the center of a pore on a void phase after a local region is divided by a watershed algorithm, and identifying the distribution of the pore and the throat in the void phase;

s5, calculating the maximum and minimum throat diameters based on the throat distribution in the void phase, and calculating the fractal dimension of the throat by a box counting method;

s6, establishing the minimum throat diameter, the maximum throat diameter and the fractal dimension of a throat system, establishing a rock sample permeability analysis expression, and accurately predicting the rock sample permeability.

Preferably, the voxel size of the three-dimensional image after rock sample scanning in step S1 is L_x×L_y×L_zFor each voxel point on the rock sample scan, the x coordinate range is 1-L_xY coordinate range of 1-L_yZ coordinate range of 1-L_z(ii) a The rock sample three-dimensional image binarization method comprises the following steps:

(1) initial value of threshold T gives original scanning rock sample V_sThe minimum value of the (x, y, z) gray value is subjected to binarization processing by using the formula (1),

in the formula, V_s(x, y, z) is a gray value corresponding to the coordinate point (x, y, z) position of the scanned rock sample image, G_cFor a given threshold value for partitioning the void phase, V_seg(x, y, z) is a numerical value of 01 corresponding to the position of the coordinate point (x, y, z) of the binarized rock sample image, the solid skeleton phase region is 1, the void phase region is 0, and the coordinate point corresponding to the void phase region is marked as x_b,y_b,z_b；

(2) And (3) calculating the porosity phi after the rock sample is divided by adopting the formula (2):

in the above formula V_pRepresenting the voxel value V of each coordinate point of the binarized rock sample image_segThe sum of (x, y, z);

(3) and (3) comparing the porosity calculated by the formula (2) with the porosity measured by the rock sample laboratory, if the porosity is not matched with the porosity measured by the rock sample laboratory, gradually increasing the threshold value T, and repeating the processes (1) and (2) until the porosity calculated by the formula (2) is matched with the porosity measured by the rock sample laboratory.

Preferably, the method for extracting the central axis of the void phase in step S2 is as follows:

(1) for each coordinate point on the boundary of the pore phase and the solid skeleton phase region, searching a neighbor coordinate point connected with the coordinate point, and setting the neighbor coordinate points as front, back, upper, lower, left and right 6 types according to the spatial position dependency relationship of the neighbor coordinate points.

(2) Searching for a point satisfying the following conditions in each direction according to the sequence of the front direction, the back direction, the up direction, the down direction, the left direction and the right direction, and converting the neighbor coordinate point into a solid skeleton region in sequence, namely converting the coordinate point from a value of a void phase region 0 into a value of a solid skeleton phase region 1, wherein the conditions are as follows:

neighbor points belonging to a void phase region in all the 6 direction neighbor points are mutually communicated;

a void phase region after the neighbor coordinate points are converted into the solid skeleton region at least comprises 2 neighbor points in 6 direction neighbor points;

if there is no point in all six directions that satisfies the condition then the central axis of the void phase has been generated;

(3) repeatedly marking each coordinate point on the boundary of the void phase and the solid skeleton phase region based on the void phase after the boundary point is deleted in the step (2), classifying the neighbor coordinate points according to the step (1), and deleting the boundary point according to the step (2) until the requirement of the centered axis output in the step (2) is met;

(4) marking a generated set of mid-axis connection point spatial coordinate positions, marked as U_p(x_p,y_p,z_p) The set of connection point locations characterizes each aperture center location.

Preferably, the calculation method of the space phase euclidean distance matrix in step S3 is as follows:

x_c,y_c,z_cfor each voxel coordinate point x of the distance void phase in the boundary of the void phase region and the skeleton phase region_b,y_b,z_bThe nearest solid wall coordinate point;

according to the formula (3) and (3), the shortest distance between each point in the void phase region and the solid wall surface, namely a void phase Euclidean distance matrix E can be obtained, wherein the size of the distance matrix E is Lx multiplied by Ly multiplied by Lz, and the Euclidean distance value of the skeleton phase region is 0.

Preferably, the identification method of the pore and throat distribution in the void phase in the step S4 is as follows:

(1) the gap phase Euclidean distance matrix is segmented by adopting a watershed algorithm, the image segmentation principle of the watershed algorithm can refer to [ IEEE Transactions on Pattern Analysis & Machine Intelligence,13(06) ], 583-:

U_total＝{U₁(x,y,z)，U₂(x,y,z)，L，U_n(x,y,z)} (19)

in the formula of U_i(i is 1,2, …, n) represents a subset of x, y, z coordinate positions within each divided region, n represents the number of divided void phase local regions, U_tRepresents the bulk void phase;

(2) the segmented pore phase has both pore area and throat area, and the central axis connection point space coordinate position set U is compared_p(x_p,y_p,z_p) With each of the division areas U_iIf centered on the axis connection point, a set of spatial coordinate positions U_p(x_p,y_p,z_p) In which there is a coordinate point in the divided area U_iIf the region is a pore phase, otherwise, the region is a throat phase;

according to a set U of spatial coordinate positions of connection points with a central axis_p(x_p,y_p,z_p) After comparison, the formula (4) can be rearranged as follows:

in the formula n₁Indicates the number of pores after division, n₂Indicates the number of divided throats, U_pore(i)(x,y,z)(i＝1,2,…,n₁) Representing a set of spatial coordinates, U, of elements within each pore phase_throat(i)(x,y,z)(i＝1,2,…,n₂) The set of voxel space coordinate positions within each throat phase is represented, and the set of voxel space coordinate positions within the entire throat phase can be represented as:

U_{throat_total}＝{U_throat1(x,y,z),U_throat2(x,y,z),L,U_thoratn2(x,y,z)} (21)

the binarized image of the rock sample when considering only the throat distribution can be expressed as:

preferably, in step S5, the maximum and minimum throat diameters and the throat fractal dimension calculation method are as follows:

(1) calculating the equivalent throat diameter by adopting equivalent sphere volume for each region of the throat phase,

in the formula (d)_t(i) Is the ith throat diameter, R_esScanning image resolution, N, for a rock sample_s(i) Set U of spatial coordinate positions of voxels of the facies of the laryngeal tract_throat(i)The number of the inner space coordinates is calculated to obtain the diameter d of the throat, namely the minimum diameter d of the throat_tminAnd maximum throat diameter d_tmax；

(2) Rock sample binarization image V based on consideration of throat distribution_throatThe method adopts a box counting method to calculate the fractal dimension of the throat, and comprises the following specific steps:

filling the throat space with a series of cubes of different side lengths, i_bThe size is given by the following equidistant array:

in the above formula, n_stepRepresenting the cube side length spacing, at a given cube side length l_b(i)(i＝1,2,…,n_step) In the case of side length l_b(i) The number of cubes required to completely fill the throat space (region of value 1) is N_s(i) By plotting N on a log-log plot_sAnd l_bData points are subjected to regression fitting by adopting a formula (10) to obtain a throat fractal dimension D_{f_throat}，

ln(N_s)＝-D_{f_throat}lnl_b+C (25)

In the above formula, C is an intercept, constant.

Preferably, in step S6, based on the minimum throat diameter, the maximum throat diameter and the fractal dimension of the throat system, the rock sample permeability analysis method is as follows:

the flow rate of the fluid in a single throat is as follows:

wherein Δ p is the pressure difference across the throat, μ is the viscosity, L (d)_t) The bending length of the throat can be expressed as:

in the formula, L₀Is the length of the direction of flow of the rock sample, D_torIs a fractal dimension of tortuosity of the throat and is based on a fractal dimension D of the rock sample throat_{f_throat}The throat size distribution can be expressed as:

N(d_t) Indicating a throat diameter greater than d_tThe number of pores (2) is obtained by differentiating equation (13) on both sides:

wherein, -dN represents a throat size of d_tTo d_t+ ddt pore size, the flow of fluid in all throats being able to pass through at the smallest throat diameter d_tminAnd maximum throat diameter d_tmaxThe integration between them yields:

according to the darcy formula:

where A is the cross-sectional area in the flow direction of the rock sample, and the core permeability can be expressed as follows according to equations (14) (15):

compared with the prior art, the invention has the following beneficial effects:

(1) the rapid rock permeability analysis method based on the throat characteristics of the scanned image can accurately consider the structural characteristics of the throat, represents the scale invariance of the throat system structure of a rock sample through the fractal dimension of the throat, overcomes the limitation that the conventional digital imaging rock permeability analysis method is limited by single scale resolution, and improves the accuracy of permeability analysis results.

(2) Compared with the existing permeability analysis method for carrying out numerical calculation based on the rock scanning image, the rapid rock permeability analysis method based on the throat characteristics of the scanning image has the advantages that the analysis precision is guaranteed, the calculation speed is greatly improved, and the rapid rock permeability analysis method based on the throat characteristics of the scanning image is convenient to popularize and apply.

Drawings

FIG. 1 is a schematic flow chart of a rapid rock permeability analysis method according to the present invention;

fig. 2 is a sandstone rock sample CT scanning image, and color bars represent gray values, and the overall voxel size: 400 × 400 × 400, physical size: 10^-3m×10^-3m×10^-3m, resolution 2.5 μm;

fig. 3 is a binarized image of a sandstone rock sample CT scanning image, and the overall voxel size is: 550 × 550 × 500, physical size: 10^-3m×10^-3m×10^-3m, resolution 12.559 μm, black indicating a void phase and white indicating a solid skeletal phase;

FIG. 4 is a central axis space distribution of sandstone rock sample CT scanning image binarization image void phases;

FIG. 5 is a space phase Euclidean distance matrix in a binary image of a CT scanning image of a sandstone rock sample, the size of a voxel is 400 multiplied by 400, a color bar value represents the distance from the wall surface of a solid skeleton, and the unit is represented by the number of the voxel;

FIG. 6 shows a sandstone three-dimensional image void phase segmented by a watershed algorithm, wherein color bar values represent marking values of local areas of the sandstone three-dimensional image void phase segmented, and the range is 1-1081;

FIG. 7 is a visualization of a void phase of a binarized image of a CT scan image of a sandstone rock sample, including pores and throats;

FIG. 8 is a visualization of throat distribution in a void phase in a binarized image of a sandstone rock sample CT scan image;

FIG. 9 shows a calculation result of fitting the fractal dimension of the sandstone rock sample throat.

Detailed Description

The present invention will be further described with reference to the following specific examples.

Referring to fig. 1, in this embodiment, a rock permeability rapid analysis method based on a throat characteristic of a scanned image is provided by taking analysis of a rock permeability rate of a sandstone oil, gas and reservoir block in China as an example, and includes the following steps:

s1, carrying out image scanning on the rock sample, carrying out binarization processing, and identifying rock gaps and solid skeleton phases;

in this embodiment, the voxel size of the selected three-dimensional image (fig. 2) after CT scanning of the sandstone rock sample is 400 × 400 × 400, and the physical size is 10^-3m×10^-3m×10^-3m。

For each voxel point on the rock sample scan, the x coordinate range is 1-400, the y coordinate range is 1-400, and the z coordinate range is 1-400. The binaryzation method of the three-dimensional image of the rock sample comprises the following steps:

(1) initial value of threshold T gives original scanning rock sample V_sThe minimum value 0 of the (x, y, z) gradation value is binarized by the formula (1),

in the formula, V_s(x, y, z) is a gray value corresponding to the coordinate point (x, y, z) position of the scanned rock sample image, G_cFor a given threshold value for partitioning the void phase, V_segAnd (x, y, z) is a numerical value of 01 corresponding to the position of the coordinate point (x, y, z) of the binarized rock sample image, the solid skeleton phase region is 1, and the void phase region is 0. The corresponding coordinate point of the void phase region is marked as x_b,y_b,z_b。

(2) And (3) calculating the porosity phi after the rock sample is divided by adopting the formula (2):

in the above formula V_pRepresenting the voxel value V of each coordinate point of the binarized rock sample image_segThe sum of (x, y, z).

(3) And (3) comparing the porosity calculated by the formula (2) with the porosity measured by the rock sample laboratory to be 0.26, if the porosity is not matched with the porosity measured by the rock sample laboratory, gradually increasing the threshold value T, and repeating the processes (1) and (2) until the porosity calculated by the formula (2) is matched with the porosity measured by the rock sample laboratory. The resulting threshold T is 66. The three-dimensional scanning image of the sandstone rock sample after binarization is shown in figure 3.

S2, continuously removing the boundary voxels of the void phase to obtain a void phase centering axis, and marking the coordinate position of a connection point on the void phase centering axis representing the center of the pore;

in this embodiment, for extracting the central axis of the void phase from the sandstone three-dimensional image binarized by S1, the specific operations are as follows:

(1) and searching a neighbor coordinate point connected with each coordinate point on the boundary of the pore phase and the solid skeleton phase region. And (3) according to the spatial position dependency relationship of the neighbor coordinate points, classifying the neighbor coordinate points into front, back, upper, lower, left and right classes 6.

(2) Searching points meeting the following conditions in each direction according to the sequence of front, back, up, down, left and right directions, and converting the neighbor coordinate points into solid skeleton regions in sequence, namely converting the coordinate points from a value of a void phase region 0 into a value of a solid skeleton phase region 1. If there is no point in all six directions that satisfies the condition then it is an indication that a centered axis of the void phase has been generated.

The conditions were as follows:

neighbor points belonging to a void phase region in all the 6 direction neighbor points are mutually communicated;

and the interstitial phase region after the neighbor coordinate points are converted into the solid skeleton region at least comprises 2 neighbor points in 6 direction neighbor points.

(3) And (3) repeatedly marking each coordinate point on the boundary of the void phase and the solid skeleton phase region based on the void phase after the boundary point is deleted in the step (2), classifying the neighbor coordinate points according to the step (1), and deleting the boundary point according to the step (2) until the requirement of the centered axis output in the step (2) is met.

(4) Marking a generated set of mid-axis connection point spatial coordinate positions, marked as U_p(x_p,y_p,z_p). The set of connection point locations characterizing each wellThe center position of the gap.

The resulting void phase centered axis is shown in FIG. 4, according to the above-described procedure.

S3, calculating the distance between each point in the void phase and the wall surface of the framework to obtain a void phase Euclidean distance matrix;

in this embodiment, the calculation method of the gap phase euclidean distance matrix of the binarized sandstone three-dimensional image is as follows:

x_c,y_c,z_cfor each voxel coordinate point x of the distance void phase in the boundary of the void phase region and the skeleton phase region_b,y_b,z_bNearest solid wall coordinate points. And (4) obtaining the nearest distance between each point in the void phase region and the solid wall surface, namely a void phase Euclidean distance matrix E. The size of the distance matrix E is 400 × 400 × 400, and the euclidean distance value of the skeleton phase region is 0, as shown in fig. 5.

S4, segmenting the space phase Euclidean distance matrix by adopting a watershed algorithm to obtain space phase distribution after a local area is divided; and finding the position of a connecting point on a central axis representing the center of the pore on the void phase after the local area is divided by the watershed algorithm, and identifying the distribution of the pore and the throat in the void phase.

In this embodiment, the method for identifying the distribution of the pore and the throat in the void phase of the binarized sandstone three-dimensional image is as follows:

(1) the method comprises the steps of adopting a watershed algorithm to segment a gap phase Euclidean distance matrix of a sandstone three-dimensional image, wherein the image segmentation principle of the watershed algorithm can refer to [ IEEE Transactions on Pattern Analysis & Machine insight, 13 (06);, 583 + 598]. by finding the minimum value 1 in the gap phase Euclidean distance matrix as an initial threshold value, marking the initial threshold value as an initial point, gradually increasing the initial threshold value, and gradually increasing the threshold value in the threshold value increasing process to be larger than the surrounding Euclidean distance values. As the threshold value is gradually increased, the regions with Euclidean distance values smaller than the threshold value are gradually contacted with each other, and the contact line of the boundaries of different adjacent regions is recorded as the watershed. Through the steps, the void phase after the void phase Euclidean distance matrix is divided into different local areas can be expressed as follows:

U_total＝{U₁(x,y,z)，U₂(x,y,z)，L，U_n(x,y,z)} (34)

in the formula of U_i(i is 1,2, …, n) represents a subset of x, y, z coordinate positions in each divided region, n represents the number of the divided sandstone three-dimensional image void phase local regions, and has a value of 1081, U_tIndicating the overall void phase. The sandstone three-dimensional image void phase after the watershed algorithm segmentation is shown in figure 6.

(2) Both pore and throat regions are present in the segmented sandstone rock-like void phase, as shown in fig. 7. By comparing the set of spatial coordinate positions U of the connection points of the central axes_p(x_p,y_p,z_p) With each of the division areas U_iIf centered on the axis connection point, a set of spatial coordinate positions U_p(x_p,y_p,z_p) In which there is a coordinate point in the divided area U_iThen the region is the pore phase, otherwise the region is the throat phase. Thus, the set U of spatial coordinate positions is based on the point of connection with the central axis_p(x_p,y_p,z_p) After comparison, the formula (4) can be rearranged as follows:

in the formula n₁The number of pores after the division is indicated as 833, n₂The number of throats after division is indicated and is 248. U shape_pore(i)(x,y,z)(i＝1,2,…,n₁) Representing a set of spatial coordinates, U, of elements within each pore phase_throat(i)(x,y,z)(i＝1,2,…,n₂) Representing a set of voxel spatial coordinate locations within each laryngeal channel. The set of voxel spatial coordinate positions within the overall laryngeal tract can be expressed as:

U_{throat_total}＝{U_throat1(x,y,z),U_throat2(x,y,z),L,U_thoratn2(x,y,z)} (36)

therefore, the binarized image of the sandstone rock sample only in the case of considering the throat distribution can be expressed as formula (7), and the result is shown in fig. 8.

S5, calculating the maximum and minimum throat diameters based on the throat distribution in the void phase, and calculating the fractal dimension of the throat by a box counting method;

(1) in the embodiment, the equivalent volume of each area of the sandstone rock sample throat phase is equivalent by adopting the sphere volume, the equivalent throat diameter is calculated,

in the formula (d)_t(i) Is the ith throat diameter, R_es2.5 x 10 resolution for rock sample scanning image^-6m，N_s(i) Set U of spatial coordinate positions of voxels of the facies of the laryngeal tract_throat(i)Number of inner space coordinates. Based on the calculated diameter of the throat, the minimum diameter d of the throat can be obtained_tminIs 2.5 multiplied by 10^-6m and maximum throat diameter d_tmaxIs 112X 10^-6m。

(2) Rock sample binarization image V based on consideration of throat distribution_throatThe method adopts a box counting method to calculate the fractal dimension of the throat, and comprises the following specific steps:

filling a throat space by adopting a series of cubes with different side lengths, wherein the side length l of each cube_bThe size is given by the following equidistant array:

in the above formula, n_stepRepresenting the cube side length spacing, with a value of 10. Given cube side length l_b(i)(i＝1,2,…,n_step) In the case of side length l_b(i) The number of cubes required to completely fill the throat space (region of value 1) is N_s(i) By plotting N on a log-log plot_sAnd l_bData points are subjected to regression fitting by adopting a formula (10) (figure 9), and the fractal dimension D of the throat can be obtained_{f_throat}。

In(N_s)＝-D_{f_throat}Inl_b+C (40)

In the above formula, C is an intercept, constant. Obtaining a throat fractal dimension D according to the formula (40)_{f_throat}It was 1.5759, the intercept C was 10.723, and the fitness R2 value was 0.97.

S6, establishing minimum throat diameter, maximum throat diameter and throat system fractal dimension, establishing rock sample permeability analysis expression, and accurately predicting rock sample permeability

In this embodiment, the flow rate of the fluid in a single throat is:

wherein, Deltap is the pressure difference between two ends of the throat of 0.1Pa, and mu is the viscosity of 10^-3Pa·S，L(d_t) The bending length of the throat can be expressed as:

in the formula, L₀Selecting the physical size of the three-dimensional sandstone sample image with the value of 10 in the x direction for the length of the flowing direction of the rock sample^-3m。D_torFractal dimension of tortuosity of throat, 1.1. Fractal dimension D based on rock sample throat_{f_throat}1.5759, the throat size distribution can be expressed as:

N(d_t) Indicating a throat diameter greater than d_tThe number of pores in the substrate. Two-sided calculation for formula (13)The differential can be found as:

wherein, -dN represents a throat size of d_tTo d_t+ ddt pore number. The flow of fluid in all throats can pass through at the minimum throat diameter d_tmin 2.5×10^-6m and maximum throat diameter d_tmax 112×10^-6Integration between m yields:

according to the darcy formula:

wherein A is the cross-sectional area of the rock sample in the flow direction, and A is L_yL_zValue of 10^-6m². According to equation (14) (15), the core permeability may be expressed as:

the permeability of the scanning sandstone sample is calculated according to the formula (16) and is 1.937 multiplied by 10^-12m²And the permeability value measured by a core physical experiment is 1.97 multiplied by 10^-12m²The coincidence is good, and the rock permeability rapid analysis method based on the throat characteristics of the scanned image can accurately predict the rock sample permeability.

It should be noted that the above-mentioned embodiments are merely examples of the present invention, and it is obvious that the present invention is not limited to the above-mentioned embodiments, and other modifications are possible. All modifications directly or indirectly derivable by a person skilled in the art from the present disclosure are to be considered within the scope of the present invention.

Claims (7)

1. A rock permeability rapid analysis method based on throat characteristics of a scanned image is characterized by comprising the following steps:

s1, carrying out image scanning on the rock sample, carrying out binarization processing on the three-dimensional image of the rock sample, and identifying a rock gap phase and a solid skeleton phase;

s2, obtaining a void phase centering axis by continuously removing rock void phase boundary voxels, and marking the coordinate position of a connection point on the void phase centering axis representing the center of the pore;

s3, calculating the distance between each point in the rock void phase and the wall surface of the solid skeleton phase to obtain a void phase Euclidean distance matrix;

s4, segmenting the space phase Euclidean distance matrix by adopting a watershed algorithm to obtain space phase distribution after a local area is divided; finding a connecting point position on a central axis representing the center of a pore on a void phase after a local region is divided by a watershed algorithm, and identifying the distribution of the pore and the throat in the void phase;

s5, calculating the maximum and minimum throat diameters based on the throat distribution in the void phase, and calculating the fractal dimension of the throat by a box counting method;

s6, establishing the minimum throat diameter, the maximum throat diameter and the fractal dimension of a throat system, establishing a rock sample permeability analysis expression, and accurately predicting the rock sample permeability.

2. The method for rapidly analyzing rock permeability based on the throat characteristic of the scanned image according to claim 1, wherein the voxel size of the three-dimensional image after the rock sample is scanned in the step S1 is L_x×L_y×L_zFor each voxel point on the rock sample scan, the x coordinate range is 1-L_xY coordinate range of 1-L_yZ coordinate range of 1-L_z(ii) a The rock sample three-dimensional image binarization method comprises the following steps:

(1) initial value of threshold T gives original scanning rock sample V_sThe minimum value of the (x, y, z) gray value is subjected to binarization processing by using the formula (1),

in the formula, V_s(x, y, z) is a gray value corresponding to the coordinate point (x, y, z) position of the scanned rock sample image, G_cFor a given threshold value for partitioning the void phase, V_seg(x, y, z) is a numerical value of 01 corresponding to the position of the coordinate point (x, y, z) of the binarized rock sample image, the solid skeleton phase region is 1, the void phase region is 0, and the coordinate point corresponding to the void phase region is marked as x_b,y_b,z_b；

(2) And (3) calculating the porosity phi after the rock sample is divided by adopting the formula (2):

in the above formula V_pRepresenting the voxel value V of each coordinate point of the binarized rock sample image_segThe sum of (x, y, z);

(3) and (3) comparing the porosity calculated by the formula (2) with the porosity measured by the rock sample laboratory, if the porosity is not matched with the porosity measured by the rock sample laboratory, gradually increasing the threshold value T, and repeating the processes (1) and (2) until the porosity calculated by the formula (2) is matched with the porosity measured by the rock sample laboratory.

3. The method for rapidly analyzing the rock permeability based on the throat characteristic of the scanned image according to claim 1, wherein the method for extracting the central axis of the void phase in the step S2 is as follows:

(1) for each coordinate point on the boundary of the pore phase and the solid skeleton phase region, searching a neighbor coordinate point connected with the coordinate point, and setting the neighbor coordinate points as front, back, upper, lower, left and right 6 types according to the spatial position dependency relationship of the neighbor coordinate points.

(2) Searching for a point satisfying the following conditions in each direction according to the sequence of the front direction, the back direction, the up direction, the down direction, the left direction and the right direction, and converting the neighbor coordinate point into a solid skeleton region in sequence, namely converting the coordinate point from a value of a void phase region 0 into a value of a solid skeleton phase region 1, wherein the conditions are as follows:

neighbor points belonging to a void phase region in all the 6 direction neighbor points are mutually communicated;

a void phase region after the neighbor coordinate points are converted into the solid skeleton region at least comprises 2 neighbor points in 6 direction neighbor points;

if there is no point in all six directions that satisfies the condition then the central axis of the void phase has been generated;

(3) repeatedly marking each coordinate point on the boundary of the void phase and the solid skeleton phase region based on the void phase after the boundary point is deleted in the step (2), classifying the neighbor coordinate points according to the step (1), and deleting the boundary point according to the step (2) until the requirement of the centered axis output in the step (2) is met;

(4) marking a generated set of mid-axis connection point spatial coordinate positions, marked as U_p(x_p,y_p,z_p) The set of connection point locations characterizes each aperture center location.

4. The method for rapidly analyzing rock permeability based on the throat characteristic of the scanned image according to claim 1, wherein the calculation method of the Euclidean distance matrix of the void phase in the step S3 is as follows:

x_c,y_c,z_cfor each voxel coordinate point x of the distance void phase in the boundary of the void phase region and the skeleton phase region_b,y_b,z_bThe nearest solid wall coordinate point;

the shortest distance from each point in the void phase region to the solid wall surface, namely the Euclidean pore phase distance matrix E, can be obtained according to the formulas (3) and (3).

5. The method for rapidly analyzing the rock permeability based on the throat characteristics of the scanned image according to claim 1, wherein the method for identifying the pore and throat distribution in the void phase in the step S4 is as follows:

(1) the method comprises the following steps of adopting a watershed algorithm to segment a void phase Euclidean distance matrix, finding a minimum value in the void phase Euclidean distance matrix as an initial threshold value, marking the minimum value as an initial point, gradually increasing the initial threshold value, gradually increasing the threshold value in the threshold value increasing process to be larger than surrounding Euclidean distance values, gradually contacting each region of the Euclidean distance value smaller than the threshold value along with the gradual increase of the threshold value, recording a contact line of boundaries of different adjacent regions as a watershed, and dividing the void phase Euclidean distance matrix into void phases after different local regions through the steps, wherein the void phases can be expressed as follows:

U_total＝{U₁(x,y,z)，U₂(x,y,z)，L，U_n(x,y,z)} (4)

in the formula of U_i(i is 1,2, …, n) represents a subset of x, y, z coordinate positions within each divided region, n represents the number of divided void phase local regions, U_tRepresents the bulk void phase;

(2) the segmented pore phase has both pore area and throat area, and the central axis connection point space coordinate position set U is compared_p(x_p,y_p,z_p) With each of the division areas U_iIf centered on the axis connection point, a set of spatial coordinate positions U_p(x_p,y_p,z_p) In which there is a coordinate point in the divided area U_iIf the region is a pore phase, otherwise, the region is a throat phase;

according to a set U of spatial coordinate positions of connection points with a central axis_p(x_p,y_p,z_p) After comparison, the formula (4) can be rearranged as follows:

in the formula n₁Indicates the number of pores after division, n₂Indicates the number of divided throats, U_pore(i)(x,y,z)(i＝1,2,…,n₁) Representing the spatial coordinates of the elements within each pore phaseSet of positions, U_throat(i)(x,y,z)(i＝1,2,…,n₂) The set of voxel space coordinate positions within each throat phase is represented, and the set of voxel space coordinate positions within the entire throat phase can be represented as:

U_{throat_total}＝{U_throat1(x,y,z),U_throat2(x,y,z),L,U_thoratn2(x,y,z)} (6)

the binarized image of the rock sample when considering only the throat distribution can be expressed as:

6. the method for rapidly analyzing rock permeability based on the throat characteristics of the scanned image according to claim 1, wherein in the step S5, the maximum and minimum throat diameters and the calculation method of the fractal dimension of the throat are as follows:

(1) calculating the equivalent throat diameter by adopting equivalent sphere volume for each region of the throat phase,

in the formula (d)_t(i) Is the ith throat diameter, R_esScanning image resolution, N, for a rock sample_s(i) Set U of spatial coordinate positions of voxels of the facies of the laryngeal tract_throat(i)The number of the inner space coordinates is calculated to obtain the diameter d of the throat, namely the minimum diameter d of the throat_tminAnd maximum throat diameter d_tmax；

(2) Rock sample binarization image V based on consideration of throat distribution_throatThe method adopts a box counting method to calculate the fractal dimension of the throat, and comprises the following specific steps:

filling the throat space with a series of cubes of different side lengths, i_bThe size is given by the following equidistant array:

in the above formula, n_stepRepresenting the cube side length spacing, at a given cube side length l_b(i)(i＝1,2,…,n_step) In the case of side length l_b(i) The number of cubes required for completely filling the throat space is N_s(i) By plotting N on a log-log plot_sAnd l_bData points are subjected to regression fitting by adopting a formula (10) to obtain a throat fractal dimension D_{f_throat}，

ln(N_s)＝-D_{f_throat}ln l_b+C (10)

In the above formula, C is an intercept, constant.

7. The method for rapidly analyzing rock permeability based on the throat feature of the scanned image according to claim 1, wherein in step S6, based on the minimum throat diameter, the maximum throat diameter and the fractal dimension of the throat system, the method for analyzing rock sample permeability comprises:

the flow rate of the fluid in a single throat is as follows:

wherein Δ p is the pressure difference across the throat, μ is the viscosity, L (d)_t) The bending length of the throat can be expressed as:

in the formula, L₀Is the length of the direction of flow of the rock sample, D_torIs a fractal dimension of tortuosity of the throat and is based on a fractal dimension D of the rock sample throat_{f_throat}The throat size distribution can be expressed as:

N(d_t) Indicating a throat diameter greater than d_tThe number of pores (2) is obtained by differentiating equation (13) on both sides:

wherein, -dN represents a throat size of d_tTo d_t+ ddt pore size, the flow of fluid in all throats being able to pass through at the smallest throat diameter d_tminAnd maximum throat diameter d_tmaxThe integration between them yields:

according to the darcy formula:

where A is the cross-sectional area in the flow direction of the rock sample, and the core permeability can be expressed as follows according to equations (14) (15):

Images