A method for calculating effective permeability of digital core based on morphological transformation and graph theory algorithm

Abstract

The present invention proposes a method for calculating the effective permeability of a digital core based on morphological transformation and graph theory algorithms, including: obtaining a three-dimensional digital image of a target digital core sample; extracting effective seepage pores of the three-dimensional digital image based on a morphological refinement method; constructing a topological graph network of the effective seepage pores; obtaining the maximum flow rate and flow distribution based on the minimum seepage resistance according to the topological graph network; and calculating the overall permeability of the target digital core sample according to the maximum flow rate and flow distribution. The present invention can extract the effective pore network in the seepage direction based on the three-dimensional microscopic image of the rock and soil material, apply the graph theory method to characterize the seepage capacity of each pore branch of the material and their mutual influence, and deduce its permeability. It has the characteristics of convenient application and comprehensive characterization.

Description

Technical Field

[0001] The invention belongs to the technical field of digital core analysis and permeability, and in particular relates to a method for calculating effective permeability of a digital core based on morphological transformation and graph theory algorithms. Background Art

[0002] Permeability is an important technical indicator of relevant geotechnical materials in the fields of oil and gas resource exploitation, underground hazardous waste storage shielding projects, and water conservancy facilities construction. Accurately and efficiently obtaining its permeability is a necessary prior technical link for related scientific research and engineering construction. With the maturity and promotion of microscopic imaging equipment and technology, it has become an emerging research method in this field to obtain digital cores based on microscopic digital images, conduct microscopic pore structure characterization analysis, and obtain its permeability because of its advantages such as easy use, multiple analytical functions, and wide research scale.

[0003] Digital cores contain a large number of randomly distributed pores of various shapes, and their geometric shapes and connectivity are the main factors affecting the seepage capacity. On the one hand, they contain a large number of non-conductive, isolated pores or pore networks, which contribute little to the seepage, but increase the calculation time and reduce the accuracy. However, the existing calculation methods based on digital core pore networks rarely accurately extract effective pores that are related to the seepage direction and connected. On the other hand, the branches of the pore network affect each other's seepage capacity. It is necessary to analyze the distribution and tendency of the seepage flow under the interaction as a whole. However, most of the existing work uses real three-dimensional structures or lattice Boltzmann methods to achieve this goal. They are time-consuming to calculate and require high operator experience, otherwise it is difficult to obtain accurate and reasonable calculation results.

[0004] Therefore, it is necessary to propose a digital core permeability calculation method that is reasonable and easy to use, which can extract the effective pores that are related to the seepage direction and are connected, and use mathematical methods to analyze the influence of each branch in the network on the overall seepage capacity of the pore network and the adjacent branches. Summary of the invention

[0005] The purpose of the present invention is to provide a method for calculating the effective permeability of digital cores based on morphological transformation and graph theory algorithms, apply computer morphological refinement methods and graph theory analysis methods to the evaluation of the seepage capacity of digital cores, and realize convenient and accurate acquisition of permeability based on digital cores.

[0006] To achieve the above object, the present invention provides a method for calculating effective permeability of digital core based on morphological transformation and graph theory algorithm, comprising:

[0007] Acquire a three-dimensional digital image of the target digital core sample;

[0008] Extracting effective seepage pores from the three-dimensional digital image based on a morphological thinning method;

[0009] Constructing a topological network of the effective seepage pores;

[0010] According to the topological network, obtaining the maximum flow rate and flow rate distribution based on the minimum seepage resistance;

[0011] The overall permeability of the target digital core sample is calculated according to the maximum flow rate and flow rate distribution.

[0012] Optionally, before extracting the effective seepage pores of the three-dimensional digital image, the method further includes: preprocessing the three-dimensional digital image;

[0013] Preprocessing the three-dimensional digital image includes:

[0014] Performing noise reduction on the three-dimensional digital image;

[0015] Processing the denoised three-dimensional digital image using an image segmentation algorithm to obtain a three-dimensional pore-matrix binary image;

[0016] Based on the fact that mutually connected pore pixels form a same pore connected component, the pore connected component in the three-dimensional pore-matrix binary image is determined.

[0017] Optionally, extracting effective seepage pores from the three-dimensional digital image comprises:

[0018] Determine the required seepage inflow direction and seepage outflow direction, record the surface corresponding to the seepage inflow direction in the three-dimensional pore-matrix binary image as the inflow surface, and record the surface corresponding to the seepage outflow direction as the outflow surface;

[0019] Determine whether the pore-connected components in the three-dimensional pore-matrix binary image are in contact with both the inflow surface and the outflow surface, and select the pore-connected components including pixels on the inflow surface and pixels on the outflow surface to form a preliminary screened three-dimensional pore-matrix binary image;

[0020] A skeleton representing the topological morphology of the seepage path is extracted from the three-dimensional pore-matrix binary image preliminarily screened by using a morphological thinning method; wherein the skeleton is constructed by a pore network composed of connected skeleton pixels, and each skeleton pixel represents a pore pixel in a preset nearby range;

[0021] Dividing the skeleton into a branch based on a preset division condition, and removing the branch containing the endpoint skeleton pixels that are in contact with the non-inflow surface or the non-outflow surface;

[0022] Determine the range of pore pixels in the three-dimensional pore-matrix binary image initially screened represented by each skeleton pixel in the remaining skeleton, remove the corresponding pore pixels, and obtain the effective seepage pores; wherein the corresponding pore pixels refer to the pore pixels represented by the removed skeleton pixels.

[0023] Optionally, the preset division condition is: the connection condition of common skeleton pixels between each endpoint skeleton pixel and the intersection skeleton pixel;

[0024] Before dividing the skeleton into a branch based on a preset division condition, the following step is further included:

[0025] Determine the mapping relationship between each skeleton pixel and each pore pixel;

[0026] Skeleton pixels are classified, and the number of remaining skeleton pixels in a preset neighborhood of each skeleton pixel is counted. Skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is a first preset value are regarded as endpoint skeleton pixels, skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is a second preset value are regarded as ordinary skeleton pixels, and skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is greater than or equal to a third preset value are regarded as intersection skeleton pixels.

[0027] Optionally, constructing the topological graph network of the effective seepage pores includes:

[0028] Recording the three-dimensional pore-matrix binary image obtained after obtaining the effective seepage pores as a three-dimensional binary image of the effective seepage pores associated with the specified seepage direction, wherein the skeleton represents the path of the effective seepage pores;

[0029] The channel formed by the effective seepage pores in the three-dimensional binary image of the effective seepage pores is called an effective channel, the skeleton pixels corresponding to the formed channel are called effective channel skeleton pixels, and the network nodes of the topological graph network are confirmed based on the effective channel skeleton pixels;

[0030] Constructing a graph network of effective channels according to the topological relationship between skeleton pixels of the effective channels;

[0031] Construct weight indicators of the structural edges of the graph network to complete the construction of the topological graph network; wherein the weight indicators include a first weight indicator and a second weight indicator, the first weight indicator is used to characterize the maximum capacity of the seepage channel represented by each edge, and the second weight indicator is used to characterize the seepage resistance of the seepage channel represented by each edge.

[0032] Optionally, obtaining the maximum flow rate and flow distribution based on the minimum seepage resistance includes:

[0033] S1. The nodes corresponding to the effective channel skeleton pixels representing the pore network seepage input end in the topological map network are regarded as the input nodes of the topological map network, and the nodes corresponding to the effective channel skeleton pixels representing the pore network seepage output end are regarded as the output nodes of the topological map network;

[0034] S2. Calculate the maximum flow rate from the input node to the output node of the topological network using the maximum flow algorithm, record the flow rate and direction of each edge in the network when the flow rate from the input node to the output node is the largest, and record the flow distribution network at this time as G network;

[0035] S3. Construct a flow-increasing network of the G network;

[0036] S4. In the flow-increasing network, find the negative loop about the seepage cost;

[0037] S5. Adjust the flow of each edge in the G network based on the negative loop;

[0038] S6. Repeat S4 and S5 until no negative loop can be found, and determine the maximum flow rate and its flow distribution based on the minimum seepage resistance.

[0039] Optionally, calculating the overall permeability of the target digital core sample includes:

[0040] Using the Ford-Fulkerson algorithm, the minimum cut under the maximum flow based on the minimum seepage resistance is obtained; wherein the minimum cut refers to the set of edges in the network that divides a node into two non-intersecting parts;

[0041] The permeability of the edges in the minimum cut is calculated, and the flow rate of each edge in the minimum cut at the maximum flow under the minimum seepage resistance is used as the contribution of the permeability of each edge in the overall seepage. Based on the permeability of the edges in the minimum cut and the contribution, the overall permeability of the target digital core sample is calculated.

[0042] Optionally, the overall permeability of the target digital core sample is:

[0043]

[0044] Where K is the overall permeability, k t is the permeability of the tth edge in the minimum cut, q t is the flow of the tth edge in the minimum cut, vol t is the maximum capacity of the tth edge in the minimum cut, and T is the total number of edges in the minimum cut.

[0045] The present invention has the following beneficial effects:

[0046] The present invention first obtains a three-dimensional digital image of the target digital core sample, extracts effective seepage pores from the three-dimensional digital image based on the morphological refinement method, and extracts effective pore channels that are connected in the seepage direction from the initial image by using computer image processing and geometric transformation methods, removes non-conductive or "dead-end" pore channels, and improves the accuracy of the calculation work. It also effectively reduces the time consumption of calculating large-scale microscopic images.

[0047] Then, a topological graph network of effective seepage pores was constructed. According to the topological graph network, the maximum flow rate and flow distribution based on the minimum seepage resistance were obtained. According to the maximum flow rate and flow distribution, the overall permeability of the target digital core sample was calculated. The graph theory method was used to characterize and analyze the effective seepage pore network. The established graph network constructed the weights of each edge based on important factors such as the maximum capacity of each pore channel and the seepage resistance calculated by the tortuosity and the roughness of the inner wall, which improved the calculation accuracy. At the same time, the maximum flow rate and its distribution based on the minimum seepage resistance were obtained, which accurately showed the preferred path of seepage under the initial pressure and its flow distribution, and better described the seepage capacity. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] The drawings constituting a part of the present application are used to provide a further understanding of the present application. The illustrative embodiments and descriptions of the present application are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0049] Figure 1 Schematic diagram of a method for calculating effective permeability of a digital core based on morphological transformation and graph theory algorithm according to an embodiment of the present invention;

[0050] Figure 2 is a schematic diagram of a morphological thinning operation according to an embodiment of the present invention;

[0051] Figure 3 A schematic diagram of determining branches and nodes of a seepage path according to an embodiment of the present invention;

[0052] Figure 4 This is a schematic diagram of calculating the shape cross-sectional coefficient and seepage resistance of each effective seepage channel branch according to an embodiment of the present invention. DETAILED DESCRIPTION

[0053] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0054] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0055] like Figure 1 As shown, this embodiment proposes a method for calculating effective permeability of digital core based on morphological transformation and graph theory algorithm, including:

[0056] S1. Acquire a three-dimensional digital image of a target digital core sample;

[0057] S2. Extract effective seepage pores from 3D digital images based on morphological refinement method;

[0058] S3. Construct a topological network of effective seepage pores;

[0059] S4. According to the topological network, obtain the maximum flow rate and flow distribution based on the minimum seepage resistance;

[0060] S5. Calculate the overall permeability of the target digital core sample based on the maximum flow rate and flow rate distribution.

[0061] Furthermore, before extracting the effective seepage pores of the three-dimensional digital image, the method further includes: preprocessing the three-dimensional digital image;

[0062] Preprocessing of three-dimensional digital images includes:

[0063] Denoise three-dimensional digital images;

[0064] The denoised three-dimensional digital image is processed using an image segmentation algorithm to obtain a three-dimensional pore-matrix binary image;

[0065] Based on the fact that interconnected pore pixels constitute the same pore connectivity component, the pore connectivity component in the three-dimensional pore-matrix binary image is determined.

[0066] Specifically, in this embodiment, S1 three-dimensional microscopic image acquisition and preprocessing.

[0067] S1.1 Use appropriate methods to obtain a three-dimensional digital image of the target sample.

[0068] S1.2 After applying the image denoising algorithm to perform image denoising, an image segmentation algorithm is selected to obtain its three-dimensional pore-matrix binary image, which is recorded as bw1.

[0069] S1.3 Determine the range of each pore connection component according to the connectivity between pixels. That is, the remaining pore pixels in the 3×3×3 neighborhood of each pore pixel are considered to be interconnected, and the interconnected pore pixels constitute the same pore connection component. A total of C pore connection components are counted in the three-dimensional binary image, and the cth pore connection component is recorded as pore connection component c (c=1, 2, ..., C).

[0070] Furthermore, extracting effective seepage pores from the three-dimensional digital image includes:

[0071] Determine the required seepage inflow direction and seepage outflow direction, record the surface corresponding to the seepage inflow direction in the three-dimensional pore-matrix binary image as the inflow surface, and record the surface corresponding to the seepage outflow direction as the outflow surface;

[0072] Determine whether the pore-connected components in the three-dimensional pore-matrix binary image are in contact with both the inflow surface and the outflow surface, and select the pore-connected components that contain pixels on the inflow surface and pixels on the outflow surface to form a preliminary screened three-dimensional pore-matrix binary image;

[0073] Using the morphological thinning method, a skeleton representing the topological morphology of the seepage path is extracted from the initially screened three-dimensional pore-matrix binary image; wherein the connected skeleton pixels constitute the skeleton of the pore network, and each skeleton pixel represents the pore pixel in the preset nearby range;

[0074] Divide the skeleton into a branch based on a preset division condition, and remove the branch containing the endpoint skeleton pixels that are in contact with the non-inflow surface or the non-outflow surface;

[0075] The range of pore pixels in the three-dimensional pore-matrix binary image initially screened by each skeleton pixel in the remaining skeleton is determined, and the pore pixels represented by the removed skeleton pixels are removed to obtain effective seepage pores.

[0076] Specifically, in this embodiment, S2 extracts effective seepage pores based on morphological transformation.

[0077] S2.1 Determination of seepage direction

[0078] First, the required inflow direction and outflow direction of the seepage are determined. Then, the corresponding inflow and outflow surfaces in the three-dimensional binary image are obtained in the computer image processing software. Specifically, the surface corresponding to the inflow direction of the seepage among the six surfaces of the three-dimensional binary image is recorded as the inflow surface, and the surface corresponding to the outflow direction of the seepage is recorded as the outflow surface.

[0079] S2.2 Obtain the pore connectivity component that is connected to both the inflow surface and the outflow surface

[0080] From the C pore connectivity components, determine whether each pore connectivity component is in contact with both the inflow surface and the outflow surface. That is, select the pore connectivity components that contain pixels on the inflow surface and pixels on the outflow surface. These selected pore connectivity components form a preliminary screened three-dimensional pore-matrix binary image, recorded as bw2.

[0081] S2.3 Determine the branches and nodes of the seepage path, such as Figure 3 As shown;

[0082] First, the morphological thinning method is used to extract the skeleton representing the topological morphology of the seepage path from bw2. The connected skeleton pixels form the skeleton of the pore network, and each skeleton pixel represents the pore pixels that are close (preset nearby range).

[0083] The object of this embodiment is the microscopic images of some geotechnical materials. These materials are considered to contain many irregular pores inside. In other words, there are two types of pixels in the collected microscopic images, one represents pores (so-called empty) and the other represents matrix (so-called solid). Figure 2 As shown in FIG. 1 , the black frame is a microscopic image of the image. After processing, the pixels are divided into those representing pores (gray) and those representing the matrix (white).

[0084] In order to describe the shape of irregular pores, morphological thinning operation is used, which is to extract the skeleton line representing the geometric shape of the object. As shown in the figure below, some red pixels are extracted, and the skeleton pixels composed of them can represent the gray geometric shape. This means that each skeleton pixel (black line segment) represents some pore pixels (gray) nearby. So the proximity here is the range of pore pixels that each skeleton pixel can represent.

[0085] Then, the mapping relationship between each skeleton pixel and each pore pixel is determined, that is, which pore pixels a certain skeleton pixel represents.

[0086] Finally, the skeleton pixels are classified. The number of other skeleton pixels in the 3×3×3 neighborhood of each skeleton pixel is counted, and the skeleton pixel with the number of other skeleton pixels of 1 in the 3×3×3 neighborhood is regarded as the endpoint skeleton pixel, which means that the pore pixel it represents is located at the end of the percolation path; the skeleton pixel with the number of other skeleton pixels of 2 in the 3×3×3 neighborhood is regarded as an ordinary skeleton pixel, which means that the pore pixel it represents is located in the percolation path; the skeleton pixel with the number of other skeleton pixels greater than or equal to 3 in the 3×3×3 neighborhood is regarded as the intersection skeleton pixel, which means that the pore pixel it represents is located at the intersection of the percolation path.

[0087] S2.4 Remove the non-conductive seepage path and its branches

[0088] According to the connection between the common skeleton pixels of each endpoint skeleton pixel and the skeleton pixels of the intersection point, the skeleton is divided into a branch. From all the skeleton pixels, the branches containing the endpoint skeleton pixels that are in contact with the non-inflow surface or the non-outflow surface are eliminated to filter out the "non-conductive" percolation branches.

[0089] Determine the range of pore pixels in bw2 represented by each skeleton pixel, remove the pore pixels represented by the removed skeleton pixels, and obtain the effective seepage pores. The three-dimensional pore-matrix binary image at this time is recorded as bw3. So far, bw3 represents the three-dimensional binary image of the effective seepage pores related to the specified seepage direction, and its skeleton represents the path of the effective seepage pores.

[0090] Furthermore, constructing a topological network of effective seepage pores includes:

[0091] The three-dimensional pore-matrix binary image obtained after the initial screening of the effective seepage pores is recorded as the three-dimensional binary image of the effective seepage pores related to the specified seepage direction, and its skeleton represents the path of the effective seepage pores;

[0092] The channel formed by the effective seepage pores in the three-dimensional binary image of the effective seepage pores is called an effective channel, the skeleton pixels corresponding to the formed channel are called effective channel skeleton pixels, and the network nodes of the topological graph network are confirmed based on the effective channel skeleton pixels;

[0093] Construct a graph network of effective channels based on the topological relationship between skeleton pixels of effective channels;

[0094] Construct weight indicators of the structural edges of the graph network to complete the construction of the topological graph network; wherein the weight indicators include a first weight indicator and a second weight indicator, the first weight indicator is used to characterize the maximum capacity of the seepage channel represented by each edge, and the second weight indicator is used to characterize the seepage resistance of the seepage channel represented by each edge.

[0095] Specifically, in this embodiment, S3 constructs a graph network of effective seepage pores.

[0096] S3.1 Determine the network nodes

[0097] The channel formed by effective seepage pores is called an effective channel, and the skeleton pixels corresponding to the channel formed by effective seepage pores are called effective channel skeleton pixels. The effective channel skeleton pixels in bw3 are divided into three categories: those belonging to effective channel branches, those belonging to effective channel intersections, and those belonging to effective channel endpoints. The distinction method is as follows: count the number of other effective channel skeleton pixels in the 3×3×3 neighborhood of each effective channel skeleton pixel. If the number value is 1, the effective channel skeleton pixel belongs to the effective channel endpoint, that is, the effective channel corresponding to the pixel and its neighborhood represent that fluid seepage can stop. If the number value is 2, the effective channel skeleton pixel belongs to the effective channel branch, that is, the effective channel corresponding to the pixel represents that fluid seepage can pass through. If the number value is greater than 2, the effective channel skeleton pixel belongs to the effective channel intersection, that is, the effective channel corresponding to the pixel represents that multiple fluid seepage can cross and pass through. That is, the effective channel is differentially divided into various branches by the intersection and endpoint of the effective channel.

[0098] At this point, the effective channel skeleton pixels representing each effective channel branch form a corresponding effective channel skeleton pixel cluster. Then a total of I effective channel branches are obtained, and a total of I effective channel skeleton pixel clusters are corresponding. That is, the i-th effective channel skeleton pixel cluster ske i Represents the i-th valid channel branch chan i , (i=1,2,…,I), S3.2 Preliminary construction of the graph network structure of effective seepage pore channels

[0099] Construct a graph network of valid channels based on the topological relationship between valid skeleton pixels. Set corresponding nodes in the graph network for valid skeleton pixels belonging to valid channel intersections and valid channel endpoints. Then construct corresponding undirected edges based on the connectivity between valid channel endpoints and intersections. That is, if there are valid channel branch skeleton pixels connecting the valid channel endpoints or intersections corresponding to two nodes, then construct an undirected edge between the two nodes; otherwise, no undirected edge is constructed between the two nodes. The i-th valid channel branch chan i The corresponding undirected edge is denoted as edge1 i .

[0100] At this point, the graph network nodes corresponding to the intersections of the valid channels are recorded as intersection nodes, the graph network nodes corresponding to the endpoints of the valid channels are recorded as intersection endpoints, and each undirected edge corresponds to a branch of each valid channel in the graph.

[0101] S3.2 The weight index 1 of the network structure edge is constructed to describe the maximum capacity of the seepage channel represented by each edge. The maximum capacity of each branch is used as the weight index 1 of the corresponding edge. In order to reduce the time consumption of processing complex pore networks, the cross-section of the pore seepage channel corresponding to each branch is assumed to be circular.

[0102] S3.2.1 Calculate the i-th valid channel branch chan i Total length, written as length i .

[0103] S3.2.2 Determine ske i The intersection of the effective channels in contact with each other or any one of the intersections, denoted as ske i，start . Then, from ske i Of all the pixels, calculate their i Along the way i，start The distance is selected as length i ×1 / 4, length i ×2 / 4 and length i ×3 / 4 three pixels, respectively denoted as ske i，1 、ske i，2and ske i，3 .

[0104] S3.2.3 Calculate ske i，1 、ske i，2 and ske i，3 The distances to the nearest matrix pixels are calculated as the average of these three distance values ​​d i .

[0105] S3.2.4 Calculation of chan i The maximum capacity of vol i

[0106] vol i =length i ×d i 2 ×Π / 4 (1)

[0107] So far, the maximum capacity of the i-th valid channel branch is denoted as vol i , is regarded as the weight 1 of the corresponding edge in the graph network.

[0108] S3.3 Construction of weight index 2 of the edge of the graph network structure, which is used to describe the seepage resistance of the seepage channel represented by each edge.

[0109] The seepage resistance of each effective seepage channel branch is evaluated based on its cross-sectional area and its change, cross-sectional shape, path length, etc., thereby constructing the weight index of the graph network structure edge.

[0110] S3.3.1 Calculate the surface area of ​​each effective channel. Determine all matrix voxels at the matrix-pore interface, that is, the matrix pixels that are in contact with the pore voxels. Determine the closest valid skeleton pixels corresponding to these matrix pixels, and assign the closest valid skeleton pixels to the skeleton. i The matrix pixels of are regarded as the surface of the i-th effective channel. The surface area of ​​the three-dimensional volume composed of them is regarded as the surface area of ​​the i-th effective channel, denoted as As i .

[0111] S3.3.2 Calculate the shape cross-sectional coefficient and seepage resistance of each effective seepage channel branch, such as Figure 4 shown.

[0112] The cross-sectional shape coefficient C of the i-th effective channel i The calculation formula is:

[0113] C i =vol i ×length i / As i 2 (2)

[0114] Count the i-th valid channel skeleton pixel cluster ske i There are K skeleton pixels in total. Count the distance between the skeleton pixels and the nearest matrix pixel. i The distance between the kth skeleton pixel and its nearest matrix pixel is denoted as d i，k (k=1,2,…,K). Then the seepage resistance w of the i-th effective channel is i for:

[0115]

[0116] At this point, w i and As i As the edgeG i The weight of .

[0117] Furthermore, obtaining the maximum flow rate and flow distribution based on the minimum seepage resistance includes:

[0118] Step 1. The nodes corresponding to the effective channel skeleton pixels representing the percolation input end of the pore network in the topological map network are regarded as the input nodes of the topological map network, and the nodes corresponding to the effective channel skeleton pixels representing the percolation output end of the pore network are regarded as the output nodes of the topological map network;

[0119] Step 2. Use the maximum flow algorithm to calculate the maximum flow from the input node to the output node of the topology network, record the flow size and direction of each edge in the network when the flow from the input node to the output node is the largest, and record the flow distribution network at this time as G network;

[0120] Step 3. Construct a flow-increasing network of the G network;

[0121] Step 4. In the flow-increasing network, find the negative loops about the seepage cost;

[0122] Step 5. Adjust the flow of each edge in the G network based on the negative loop;

[0123] Step 6. Repeat steps 4 and 5 continuously, update the G network, until no negative loop can be found, and determine the maximum flow rate and its flow distribution based on the minimum seepage resistance.

[0124] Specifically, in this embodiment, S4 calculates the maximum flow rate and flow distribution based on the minimum seepage resistance.

[0125] After obtaining the pore topology graph network, the maximum seepage flow rate between the input and output ends is obtained, and the flow distribution in the graph network is continuously adjusted to reduce the seepage resistance consumption, thereby obtaining the maximum flow and its distribution with the minimum seepage resistance.

[0126] S4.1 Determine the nodes corresponding to the percolation input and output in the graph network. The nodes corresponding to the valid channel skeleton pixels representing the percolation input of the pore network are regarded as the input nodes of the graph network. The nodes corresponding to the valid channel skeleton pixels representing the percolation output of the pore network are regarded as the output nodes of the graph network.

[0127] S4.2 uses the maximum flow algorithm to calculate the maximum flow from the input node to the output node of the graph network, as well as the flow distribution of each edge at this time. Specifically, the Boykov-Kolmogorov algorithm is applied under the condition that the flow of each edge does not exceed its maximum capacity. By advancing forward and backward, the flow in the graph is continuously adjusted to find the maximum flow from the input node to the output node. The requirements during the calculation are as follows:

[0128] (1) The flow rate of each edge does not exceed its maximum capacity;

[0129] (2) Except for the edges connected to the input or output nodes, no requirements are made on the flow direction of the other edges. When the flow rate is not zero, the edges connected to the input node can only flow away from the input node. When the flow rate is not zero, the edges connected to the output node can only flow toward the output node.

[0130] (3) Except for the input and output nodes, the sum of the input flows of all other nodes is always equal to the sum of the output flows;

[0131] (4) The sum of the input flows of each input node is always equal to the sum of the output flows of each output node;

[0132] (5) The maximum flow from the input node to the output node is obtained based on the input flow of each input node and the maximum time.

[0133] Record the flow size and direction of each edge in the network when the flow from the input node to the output node is the largest, and record the flow distribution graph network at this time as G.

[0134] S4.3 Constructing the flow-increasing network D of G f That is, when the flow between the input / output nodes of the current graph network is the largest, the flow and maximum capacity of each effective channel branch are compared to determine the variable range in future flow adjustment. Specifically, the edge G i The flow rate is recorded as flow i At this time, 0≤flow i ≤vol i .

[0135] Then construct a flow-increasing graph network with nodes distributed in the same way as G, denoted as D f . D f The edges in are all set as bidirectional edges. fThe flow in all edges is 0. Then compare G and D f , determine D f The marking method of each edge in D f Middle and edgeG i The edges with the same connected nodes and the same connection direction are denoted as edgeD1 i , D f Middle and edge1 i The connected nodes are the same, but the edges with opposite connection directions are denoted as edgeD2 i .

[0136] Now, according to the relationship between the flow rate of each edge in G and the maximum capacity, determine the flow-increasing network D f The flow of each edge in: If flow i =0, then edgeD1 i The flow rate becomes vol i , then edgeD2 i The flow rate is still 0; if 0 <flow i <vol i , then edgeD1 i The flow rate becomes vol i -flow i , then edgeD2 i Flow i ;if flow i =vol i , then edgeD1 i The flow rate becomes vol i -flow i , then edgeD2 i Flow i .

[0137] After determining the flow of each edge, the flow increase network D is constructed. f Except for the input or output nodes, all nodes are connected by bidirectional edges.

[0138] S4.4 In the flow increasing network D f , find negative loops with respect to the seepage cost.

[0139] First, calculate D f The percolation cost of each edge in D is calculated. In this step, two edges connecting two nodes (non-input nodes or non-output nodes) are considered to be independent edges. f Taking any edge in the equation as an example, its theoretical seepage cost = the flow rate of the edge × the seepage resistance of the edge.

[0140] Then, determine D fNegative loops in . A loop refers to several nodes connected in series by edges, where each node has only one input edge and one output edge, that is, these edges form a loop. If the percolation method of a certain edge is consistent with the overall flow direction of the loop, its actual percolation cost on the loop is equal to its theoretical percolation cost; if the percolation method of a certain edge is opposite to the overall flow direction of the loop, its actual percolation cost on the loop is equal to the opposite of its theoretical percolation cost. Calculate the sum of the actual percolation costs of all edges on the loop. If the actual percolation costs of all edges on the loop are negative, the loop is a negative loop. In actual calculations, in order to reduce the iterative calculation time, only loops with no less than 20 edges are considered.

[0141] S4.5 combined with flow enhancement ring D f The negative loop in G adjusts the flow of each edge in G.

[0142] Determine D f The negative loop in G sorts out the edges at the corresponding positions in G and updates the flow of the edges at the corresponding positions in G. f The flow direction of the corresponding edge on the negative loop in G is the same, then the flow of this edge in G is equal to the existing flow plus D f The flow of the corresponding edge on the negative circuit in G; if an edge in G is connected to D f The flow direction of the corresponding edge on the negative loop in G is opposite, so the flow rate of this edge in G is equal to the existing flow rate minus D f The flow of the corresponding edge on the negative circuit in . This updates the flow distribution and total flow in the G network.

[0143] S4.6 Repeat the above steps

[0144] Repeat steps S4.4 and S4.5, update network G, and find its flow-increasing network D f， Find its negative loop, thereby determining its negative loop. Repeat this operation until no negative loop can be found. Thus, the maximum flow based on the minimum seepage resistance and its flow distribution are determined. In order to reduce the iterative calculation time, only loops with no less than 20 edges are considered.

[0145] Furthermore, calculating the overall permeability of the target digital core sample includes:

[0146] The Ford-Fulkerson algorithm is used to obtain the minimum cut under the maximum flow based on the minimum percolation resistance; the minimum cut refers to the set of edges in the network that divides a node into two non-intersecting parts;

[0147] The permeability of the edges in the minimum cut is calculated, and the flow rate of each edge in the minimum cut at the maximum flow under the minimum seepage resistance is used as the contribution of the permeability of each edge in the overall seepage. Based on the permeability and contribution of the edges in the minimum cut, the overall permeability of the target digital core sample is calculated.

[0148] Specifically, in this embodiment, S5 calculates the overall permeability based on the minimum cut.

[0149] S5.1 Obtain the minimum cut set of the maximum flow with minimum seepage resistance

[0150] Use the Ford-Fulkerson algorithm to obtain the minimum cut under the maximum flow of the graph network based on the minimum percolation resistance. The minimum cut refers to the set of edges in the network that divide the node into two non-intersecting parts. By finding the minimum cut, the cut set of edges is obtained, that is, the flow of these edges is equal to the maximum flow, so they are the bottleneck of the overall flow. Specifically, in this situation, there are T edges in the minimum cut, and the tth edge in the minimum cut is recorded as edge3 t (t=1,2,…,T), edge3 t The flow rate at this time is recorded as q t , the maximum capacity is vol t .

[0151] S5.2 Overall permeability calculation

[0152] The permeability of the edges in the minimum cut is calculated, and the permeability of the entire percolation network is inferred from this.

[0153] First, the Hagen-Poiseuille principle is applied to calculate the permeability of each edge in the minimum cut. Get edge3 t The permeability of the corresponding branch is k t Then, based on the flow rate of each edge in the minimum cut at the maximum flow under the minimum seepage resistance, the contribution of the permeability of each edge to the overall seepage is taken, and the overall permeability K is calculated:

[0154]

[0155] At this point, the effective permeability of the pattern to be measured is obtained.

[0156] This embodiment uses computer image processing and geometric transformation methods to extract the effective pore channels in the seepage direction from the initial image, remove the non-conductive or "dead-end" pore channels, improve the accuracy of the calculation work, and effectively reduce the calculation time of processing large-scale microscopic images.

[0157] The graph theory method is applied to characterize and analyze the effective seepage pore network. The established graph network constructs the weight of each edge based on important factors such as the maximum capacity of each pore channel and the seepage resistance calculated by the tortuosity and the roughness of the inner wall, which improves the calculation accuracy. At the same time, the maximum flow rate and its distribution based on the minimum seepage resistance are obtained, which accurately shows the preferred path of seepage under the initial pressure and its flow distribution, and better describes the seepage capacity.

[0158] Easy to use. Compared with existing methods such as the lattice Boltzmann method or numerical simulation based on real three-dimensional models, there are few parameters that need to be manually adjusted or involved when processing different digital cores. Ordinary technicians who do not have rich experience in this field can also complete the work flow and obtain reasonable calculation results.

[0159] The above are only preferred specific implementations of the present application, but the protection scope of the present application is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed in the present application should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A digital core effective permeability calculation method based on morphological transformation and graph theory algorithm, characterized in that: include: Acquire a three-dimensional digital image of the target digital core sample; Extracting effective seepage pores from the three-dimensional digital image based on a morphological thinning method; Constructing a topological network of the effective seepage pores; According to the topological network, obtaining the maximum flow rate and flow rate distribution based on the minimum seepage resistance; Obtaining the maximum flow rate and flow distribution based on the minimum seepage resistance includes: S1. The nodes corresponding to the effective channel skeleton pixels representing the pore network seepage input end in the topological map network are regarded as the input nodes of the topological map network, and the nodes corresponding to the effective channel skeleton pixels representing the pore network seepage output end are regarded as the output nodes of the topological map network; S2. Calculate the maximum flow rate from the input node to the output node of the topological network using the maximum flow algorithm, record the flow rate and direction of each edge in the network when the flow rate from the input node to the output node is the largest, and record the flow distribution network at this time as G network; S3. Construct a flow-increasing network of the G network; S4. In the flow-increasing network, find the negative loop about the seepage cost; S5. Adjust the flow of each edge in the G network based on the negative loop; S6. Repeat S4 and S5 until no negative loop can be found, and determine the maximum flow rate and flow distribution based on the minimum seepage resistance; The overall permeability of the target digital core sample is calculated according to the maximum flow rate and flow rate distribution.

2. The method for calculating effective permeability of digital cores based on morphological transformation and graph theory algorithm according to claim 1 is characterized in that: Before extracting the effective seepage pores of the three-dimensional digital image, the method further includes: preprocessing the three-dimensional digital image; Preprocessing the three-dimensional digital image includes: Performing noise reduction on the three-dimensional digital image; Processing the denoised three-dimensional digital image using an image segmentation algorithm to obtain a three-dimensional pore-matrix binary image; Based on the fact that mutually connected pore pixels form a same pore connected component, the pore connected component in the three-dimensional pore-matrix binary image is determined.

3. The method for calculating effective permeability of digital core based on morphological transformation and graph theory algorithm according to claim 2 is characterized in that: Extracting effective seepage pores from the three-dimensional digital image includes: Determine the required seepage inflow direction and seepage outflow direction, record the surface corresponding to the seepage inflow direction in the three-dimensional pore-matrix binary image as the inflow surface, and record the surface corresponding to the seepage outflow direction as the outflow surface; Determine whether the pore-connected components in the three-dimensional pore-matrix binary image are in contact with both the inflow surface and the outflow surface, and select the pore-connected components including pixels on the inflow surface and pixels on the outflow surface to form a preliminary screened three-dimensional pore-matrix binary image; A skeleton representing the topological morphology of the seepage path is extracted from the three-dimensional pore-matrix binary image preliminarily screened by using a morphological thinning method; wherein the skeleton is constructed by a pore network composed of connected skeleton pixels, and each skeleton pixel represents a pore pixel in a preset nearby range; Dividing the skeleton into a branch based on a preset division condition, and removing the branch containing the endpoint skeleton pixels that are in contact with the non-inflow surface or the non-outflow surface; Determine the range of pore pixels in the three-dimensional pore-matrix binary image initially screened represented by each skeleton pixel in the remaining skeleton, remove the corresponding pore pixels, and obtain the effective seepage pores; wherein the corresponding pore pixels refer to the pore pixels represented by the removed skeleton pixels.

4. The method for calculating effective permeability of digital cores based on morphological transformation and graph theory algorithm according to claim 3 is characterized in that: The preset division condition is: the connection condition of common skeleton pixels between each endpoint skeleton pixel and the intersection skeleton pixel; Before dividing the skeleton into a branch based on a preset division condition, the following step is further included: Determine the mapping relationship between each skeleton pixel and each pore pixel; Skeleton pixels are classified, and the number of remaining skeleton pixels in a preset neighborhood of each skeleton pixel is counted. Skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is a first preset value are regarded as endpoint skeleton pixels, skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is a second preset value are regarded as ordinary skeleton pixels, and skeleton pixels whose number of remaining skeleton pixels in the preset neighborhood is greater than or equal to a third preset value are regarded as intersection skeleton pixels.

5. The method for calculating effective permeability of digital core based on morphological transformation and graph theory algorithm according to claim 3 is characterized in that: Constructing the topological graph network of the effective seepage pores includes: Recording the three-dimensional pore-matrix binary image obtained after obtaining the effective seepage pores as a three-dimensional binary image of the effective seepage pores associated with the specified seepage direction, wherein the skeleton represents the path of the effective seepage pores; The channel formed by the effective seepage pores in the three-dimensional binary image of the effective seepage pores is called an effective channel, the skeleton pixels corresponding to the formed channel are called effective channel skeleton pixels, and the network nodes of the topological graph network are confirmed based on the effective channel skeleton pixels; Constructing a graph network of effective channels according to the topological relationship between skeleton pixels of the effective channels; Construct weight indicators of the structural edges of the graph network to complete the construction of the topological graph network; wherein the weight indicators include a first weight indicator and a second weight indicator, the first weight indicator is used to characterize the maximum capacity of the seepage channel represented by each edge, and the second weight indicator is used to characterize the seepage resistance of the seepage channel represented by each edge.

6. The method for calculating effective permeability of digital cores based on morphological transformation and graph theory algorithm according to claim 1 is characterized in that: Calculating the overall permeability of the target digital core sample includes: Using the Ford-Fulkerson algorithm, the minimum cut under the maximum flow based on the minimum seepage resistance is obtained; wherein the minimum cut refers to the set of edges in the network that divides a node into two non-intersecting parts; The permeability of the edges in the minimum cut is calculated, and the flow rate of each edge in the minimum cut at the maximum flow under the minimum seepage resistance is used as the contribution of the permeability of each edge in the overall seepage. Based on the permeability of the edges in the minimum cut and the contribution, the overall permeability of the target digital core sample is calculated.

7. The method for calculating effective permeability of digital cores based on morphological transformation and graph theory algorithm according to claim 6 is characterized in that: The overall permeability of the target digital core sample is: Where K is the overall permeability, k t is the permeability of the tth edge in the minimum cut, q t is the flow of the tth edge in the minimum cut, vol t is the maximum capacity of the tth edge in the minimum cut, and T is the total number of edges in the minimum cut.

Images