Quantitative characterization method and device for pore structure of carbonate reservoir

Abstract

The invention relates to a carbonate reservoir pore structure quantitative characterization method and device, and the method comprises the steps: carrying out the pore structure feature extraction of obtained digital images of different scales of a carbonate reservoir, and obtaining corresponding binary images; aiming at the binary image of each scale, taking at least partial region of the binary image as a target region, and searching and screening the binary images of different scales through a similarity matching algorithm based on the pore distribution feature, the pore change feature, the pore fractal feature and the pore shape feature of the target region to obtain a pore key region; and on the basis of a pre-constructed multi-scale comprehensive fractal dimension function and the pore key areas of different scales, the carbonate reservoir multi-scale comprehensive fractal dimension is obtained through calculation. According to the method, the comprehensive quantitative characterization of the multi-scale pore structure is realized, the problem of relevance and fusion characterization among multi-scale images is solved, and data support can be provided for effectively solving the productivity prediction and evaluation of the carbonate reservoir.

Description

Technical Field

[0001] The invention relates to a method and a device for quantitatively characterizing the pore structure of a carbonate reservoir. Background Art

[0002] Carbonate reservoirs are highly heterogeneous, especially the pores, caves, and joints that serve as oil and gas storage spaces, which are diverse, complex, and have pore scale characteristics and heterogeneity at multiple length scales, making the quantitative characterization of carbonate reservoir structure extremely challenging. In order to better understand the characteristics of carbonate reservoir structure, it is necessary to obtain information on the reservoir structure at multiple scales, generally from a few nanometers to millimeters. These characteristics are often controlled by the geometry and connectivity of the reservoir space structure. In order to further improve the production evaluation of carbonate micro-reservoirs, quantitative characterization of reservoir pore structure is crucial to further clarify the structure-activity relationship and performance optimization. Summary of the invention

[0003] In order to enrich the characterization methods of carbonate reservoir pore structure and realize the quantitative characterization of carbonate reservoir pore structure, the present invention proposes a method and device for quantitative characterization of carbonate reservoir pore structure. The technical solution proposed by the present invention is as follows:

[0004] In a first aspect, an embodiment of the present invention provides a method for quantitatively characterizing the pore structure of a carbonate reservoir, comprising:

[0005] Extracting pore structure features from the acquired digital images of carbonate reservoirs at different scales to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by an optical microscope, and microscopic digital images taken by a scanning electron microscope;

[0006] The binary image of each scale is calculated based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics;

[0007] For each scale of binary images, at least a part of the binary image is used as a target area, and based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a similarity matching algorithm is used to search the binary images of different scales to screen out the pore key area;

[0008] Based on the pre-constructed multi-scale comprehensive fractal dimension function and the key pore areas at different scales, the multi-scale comprehensive fractal dimension of the carbonate reservoir is calculated.

[0009] In one or some embodiments, for the binary image of each scale, at least a part of the binary image is used as the target area, and based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a similarity matching algorithm is used to search the binary images of different scales to screen out the pore key area, including:

[0010] For each scale of binary image, perform the following operations:

[0011] Using at least a portion of the binary image as a target area;

[0012] Based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a search area with the highest comprehensive similarity to the target area is selected in binary images of different scales by a similarity matching algorithm to obtain a corresponding target area;

[0013] The Euclidean distances between the target area and the target areas of different scales are determined respectively, and the area with the smallest Euclidean distance is selected to obtain the pore key area.

[0014] In one or some embodiments, based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a search area with the highest comprehensive similarity to the target area is selected in binary images of different scales by a similarity matching algorithm to obtain a corresponding target area, including:

[0015] For each scale of the binary image, gridding the binary image to obtain a grid image;

[0016] respectively calculating the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of each search area in the grid image;

[0017] Based on the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area, and the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the search area, the pore distribution characteristics similarity, pore variation characteristics similarity, pore fractal characteristics similarity and pore shape characteristics similarity of the target area and the search area are calculated respectively, and the comprehensive similarity of the target area and each of the search areas is calculated respectively by the following formula:

[0018]

[0019] Among them, DC(q,t) is the comprehensive similarity between the target area and the search area, q is the target area, t is the search area, S xis the similarity of pore distribution characteristics, γ is the similarity of pore change characteristics, D is the similarity of pore fractal characteristics, G is the similarity of pore shape characteristics, λ a , b , c and λ d are the weights of pore distribution feature similarity, pore change feature similarity, pore fractal feature similarity and pore shape feature similarity respectively;

[0020] In binary images of different scales, search regions having the highest comprehensive similarity to the target region are respectively selected to obtain corresponding target regions.

[0021] In one or some embodiments, the pore distribution characteristics are calculated by the following formula:

[0022]

[0023] Among them, S x (r1, r2) is the pore distribution characteristic, r1 and r2 are any two points in the same scale image, P(r1) is the phase function of point r1, and P(r2) is the phase function of point r2.

[0024] In one or some embodiments, the pore change characteristic is calculated by the following formula:

[0025]

[0026] Among them, γ(r) is the pore change characteristic, r is the equivalent radius of the pore, N(r) is the number of two points with a distance r in the same scale image, and r i is the radius of the ith pore, P(r i ) is r i The phase function value of .

[0027] In one or some embodiments, the pore fractal feature is calculated by the following formula:

[0028]

[0029] Among them, D F is the pore fractal feature, r is the equivalent radius of the pore, N'(r) is the number of pores with a pore diameter greater than r, and a is a constant.

[0030] In one or some embodiments, the pore shape characteristics are calculated by the following formula:

[0031]

[0032] Among them, G T is the pore shape characteristic, S T is the area of ​​the Tth pore region; CT is the perimeter of the Tth pore area.

[0033] In one or some embodiments, the multi-scale integrated fractal dimension of the carbonate reservoir is calculated based on the pre-constructed multi-scale integrated fractal dimension function and the pore key areas of different scales, including:

[0034] Based on the following formula, the pores of different pore size segments in the key pore areas of different scales are calculated to obtain the multi-scale comprehensive fractal dimension of the carbonate reservoir:

[0035]

[0036] Where: d z is the multi-scale comprehensive fractal dimension of carbonate reservoir, B i is the ratio of the pore volume of the critical area of ​​the pores at the i-th scale to the total pore volume, i = 1, 2, ..., m, m is the number of scales, d ij is the fractal dimension of the jth aperture segment in the pore critical region of the i-th scale, j = 1, 2, ..., l, l is the number of aperture segments divided in the pore critical region of the i-th scale, b ij It is the ratio of the pore volume of the jth pore size segment in the pore critical region of the i-th scale to the pore volume of the pore critical region of this scale.

[0037] In one or some embodiments, the method further comprises:

[0038] The maximum circle search method is used to search the radius of the pore area of ​​the binary image, and the radius of the pore area of ​​the single pore structure is determined by the following formula:

[0039]

[0040] Where R is the radius of the pore area, δd is the actual size represented by each pixel, and N r is the number of pixels of the maximum inscribed circle radius of the pore area.

[0041] In one or some embodiments, the method further comprises:

[0042] The area of ​​the pore region in the binary image is determined based on the following formula:

[0043]

[0044] Where S is the area of ​​the pore region, δd is the actual size represented by each pixel, n×n is the number of grids in the pore region, and y ij is the ordinate of the leftmost point in each row, y' ij is the ordinate of the rightmost point in each row, yij -y' ij +1 is the number of pixels contained in each row of the grid, i = 1, 2, ..., n.

[0045] In one or some embodiments, the pore structure features of the acquired digital images of carbonate reservoirs at different scales are extracted to obtain corresponding binary images, including:

[0046] Digital images of carbonate reservoirs at every scale:

[0047] Preprocessing the digital images respectively to obtain preprocessed digital images;

[0048] An iterative image threshold segmentation algorithm is used to identify the pore structure of the preprocessed digital image to obtain a binary image.

[0049] In one or some embodiments, the method of using an iterative image threshold segmentation algorithm to identify the pore structure of the preprocessed digital image to obtain a binary image includes:

[0050] Based on the minimum gray value PF in the preprocessed digital image out | min and maximum gray value PF out | max , the initial threshold T0 is determined by the following formula;

[0051]

[0052] Segmenting the preprocessed digital image into a target portion and a background portion according to the initial threshold, and updating the threshold by the following formula based on the grayscale values ​​of the target portion and the background portion to obtain an updated threshold;

[0053]

[0054]

[0055] Among them, PF out | O is the average gray value of all pixels in the target part, PF out | G is the average gray value of all pixels in the background, PF out (i, j) is the gray value of point (i, j) on the image, N(i, j) is the weight coefficient of point (i, j), T K is the threshold before updating, T K+1 is the updated threshold;

[0056] Segmenting the preprocessed digital image based on the updated threshold value, repeating the process of updating the threshold value until a preset condition is met to obtain an optimal threshold value;

[0057] The preprocessed digital image is segmented based on the optimal threshold to obtain a binary image.

[0058] In one or some embodiments, preprocessing the digital images to obtain preprocessed digital images includes:

[0059] The characteristic components of the digital image are extracted based on the following formula:

[0060]

[0061] Among them, v is the feature component, L is the image grayscale, g is the pixel grayscale value, i g is the number of pixels of image gray level g, M and N are the total number of rows and columns of the image matrix corresponding to the digital image respectively;

[0062] Based on the characteristic components, the digital image is subjected to image enhancement processing by the following formula to obtain a preprocessed digital image:

[0063]

[0064] Among them, Z(g) is the enhancement function and E is the enhancement parameter.

[0065] In a second aspect, the present invention provides a device for quantitatively characterizing the pore structure of a carbonate reservoir, comprising:

[0066] A feature extraction module is used to extract pore structure features from the acquired digital images of carbonate reservoirs at different scales, respectively, to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by optical microscope, and microscopic digital images taken by scanning electron microscope;

[0067] A first calculation module is used to calculate the binary image of each scale based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics;

[0068] A search module is used for taking at least a part of the binary image of each scale as a target area, searching the binary images of different scales by a similarity matching algorithm based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, and screening out the key pore area;

[0069] The second calculation module is used to calculate the multi-scale comprehensive fractal dimension of the carbonate reservoir based on the pre-constructed multi-scale comprehensive fractal dimension function and the pore key areas of different scales.

[0070] In a third aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in the first aspect.

[0071] In a fourth aspect, the present invention provides an electronic device, comprising a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other via the communication bus;

[0072] Memory, used to store computer programs;

[0073] The processor is used to implement the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in the first aspect when executing the program stored in the memory.

[0074] Based on the above technical solution, the beneficial effects of the present invention compared with the prior art are as follows:

[0075] The method for quantitative characterization of carbonate reservoir pore structure provided by the embodiment of the present invention lays a data foundation for the fusion and characterization of multi-scale pore structure by extracting characteristic parameters of pore structure of multi-scale pore structure images; calculates the binary image of each scale based on the pre-constructed characteristic function to obtain the corresponding pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics; then based on the determined pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics, it can realize cross-scale image feature area positioning, combine the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics, and realize the search of pore key areas between different scales through the similarity matching algorithm, which can effectively realize the automatic search of pore key areas in the "millimeter-micrometer-nanometer" scale image of carbonate reservoir pores. The present invention constructs a multi-scale comprehensive fractal dimension of carbonate reservoirs, which can reflect the influence of the fractal characteristics of pores in different pore diameter sections on the overall fractal characteristics of carbonate reservoirs, realizes the comprehensive quantitative characterization of multi-scale pore structure, and ultimately solves the correlation and fusion characterization problems between multi-scale images, which can provide data support for effectively solving the production capacity prediction and evaluation problems of carbonate reservoirs.

[0076] Other features and advantages of the present invention will be described in the following description, and partly become apparent from the description, or understood by practicing the present invention. The purpose and other advantages of the present invention are realized and obtained by the structures particularly pointed out in the description, claims and drawings.

[0077] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, preferred embodiments are given below and described in detail with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0079] Figure 1 It is a schematic diagram of the process of the method for quantitative characterization of the pore structure of carbonate reservoirs provided by an embodiment of the present invention;

[0080] Figure 2 Schematic diagram of the relationship between the longest axis of a pore boundary and the shortest axis of the boundary provided by an embodiment of the present invention;

[0081] Figure 3 1 is a result diagram of grid division of the preprocessed digital image F2 provided by an embodiment of the present invention; wherein (a) is a 10*10 grid macro image, (b) is a 10*10 grid meso image, (c) is a 10*10 grid micro image, (d) is a 100*100 grid macro image, (e) is a 100*100 grid meso image, and (f) is a 100*100 grid micro image;

[0082] Figure 4 It is a search graph for key areas of a 10*10 grid image provided by an embodiment of the present invention; wherein (a) is a macro-mesoscopic search graph, (b) is a mesoscopic-microscopic search graph, and (c) is a macro-microscopic search graph;

[0083] Figure 5 It is a search graph for key areas of a 100*100 grid image provided by an embodiment of the present invention; wherein (a) is a macro-mesoscopic search graph, (b) is a mesoscopic-microscopic search graph, and (c) is a macro-microscopic search graph;

[0084] Figure 6 1 is a comparison diagram of the search positioning area and the actual sampling area provided by the embodiment of the present invention; wherein (a) is a macroscopic-mesoscopic search positioning comparison diagram, (b) is a mesoscopic-microscopic search positioning comparison diagram, and (c) is a macroscopic-microscopic search positioning comparison diagram;

[0085] Figure 7 Schematic diagram of the structure of a carbonate reservoir pore structure quantitative characterization device provided by an embodiment of the present invention;

[0086] Figure 8It is a schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0087] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0088] Exemplary embodiments will be described in detail herein, examples of which are shown in the accompanying drawings. When the following description refers to the drawings, the same numbers in different drawings represent the same or similar elements unless otherwise indicated. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. Instead, they are merely examples of devices and methods consistent with some aspects of the present invention as detailed in the appended claims.

[0089] In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer", etc., indicating the orientation or positional relationship, are based on the orientation or positional relationship shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention. In addition, the terms "first", "second", and "third" are used for descriptive purposes only, and cannot be understood as indicating or implying relative importance.

[0090] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "installed", "connected", and "connected" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0091] Carbonate reservoirs are highly heterogeneous, especially the pores, caves, and joints that serve as oil and gas storage spaces, which are diverse, complex, and have pore scale characteristics and heterogeneity at multiple length scales, making the quantitative characterization of carbonate reservoir structures extremely challenging. In order to better understand the characteristics of carbonate reservoir structures, it is necessary to obtain information on reservoir structures at multiple scales, generally from a few nanometers to millimeters. These characteristics are often controlled by the geometry and connectivity of the reservoir space structure. In order to further improve the production evaluation of carbonate micro-reservoirs, quantitative characterization of reservoir structure is crucial to further clarify the structure-activity relationship and performance optimization. Traditional carbonate reservoir pore structure image observation methods are mainly through optical microscopes, cast thin sections, scanning electron microscopes, etc. Traditional image observation methods can directly describe the microscopic pore structure of the reservoir qualitatively and conduct semi-quantitative observations, and are characterized by intuitive observation and mature technology. The inventor found in his work that the existing image-based carbonate reservoir pore structure characterization method still has three main problems: first, the image source is single, and the distribution characteristics of the pore structure are not grasped from the multi-scale image features; second, the correlation characteristics and position matching relationship between the multi-scale images are unclear, and the corresponding relationship between images of different scales is not established; third, the pore structure characteristics of multi-scale and multi-data sources do not reflect the overall structural characteristics of the carbonate reservoir, and it is difficult to achieve a comprehensive quantitative characterization of the multi-scale pore structure, which cannot meet the inventor's expectations. To this end, in order to grasp the different difference characteristics of pore structures of different scales in the image, the inventor has further developed and constructed a multi-scale structural information association and fusion function of the pore structure to achieve a comprehensive quantitative characterization of the multi-scale pore structure. The present invention proposes a carbonate reservoir pore structure quantitative characterization method based on multi-scale image features.

[0092] Invention patent CN201610908816.3 "A method for characterizing the pore structure of carbonate rocks based on rock classification and multi-scale digital cores" A method for characterizing the pore structure of carbonate rocks based on multi-scale digital cores of rock classification, firstly obtain core images of different resolutions by CT scanning technology, then perform binary segmentation of the obtained digital cores, combine the optimized digital cores of different resolutions together, thereby obtaining multi-scale digital cores, evaluate the pore structure parameters of the digital cores and compare them with the experimentally measured data, and finally divide the obtained multi-scale digital cores into regions based on rock bedding, pore structure, etc., and then classify these regions by properties, summarize the classification results and fill them into regions with the same properties. The present invention combines digital cores of different scales, which can not only characterize the physical properties of the medium but also meet its requirements for resolution; at the same time, the present invention classifies the lithology of the obtained multi-scale digital cores, and this method can greatly reduce the difficulty of analysis and calculation, laying the foundation for processing more complex and more heterogeneous cores. This method only realizes the image difference characteristics of pore structure from different CT resolutions, and does not realize multi-scale image analysis research across orders of magnitude.

[0093] Invention patent CN202210159608.3 "Pore and fracture identification method and system based on two-dimensional core scanning image", the present invention belongs to the field of oil and gas reservoir geological exploration technology, specifically relates to a pore and fracture identification method and system based on two-dimensional core scanning image, which aims to solve the problem that the existing technology cannot accurately and efficiently segment pores and fractures, and cannot quantitatively analyze the characteristics of pores and fractures; the method includes scanning and obtaining the initial core plane two-dimensional image, filtering it, and obtaining the first core plane two-dimensional image; then performing image segmentation to obtain the second core plane two-dimensional image; extracting the central coordinates of all pixel points in each pore space, obtaining the centroid of the corresponding pore space, and establishing a pore and fracture identification function; if the value of the pore and fracture identification function is greater than the preset characterization value, the pore space is determined to be a fracture; otherwise, the pore space is determined to be a pore; the present invention can obtain accurate pore and fracture distribution information, and provide effective guidance for reservoir oil and gas exploitation. This method mainly grasps the characteristics of the pore structure from the macroscopic image, and does not consider the pore structure identification and characterization problems of small-scale images such as mesoscopic and microscopic images.

[0094] Invention patent CN202210096371.9 "A method for extracting digital core structure of low-permeability conglomerate reservoir" obtains the core casting thin section image of the target reservoir, counts the particle size of the mineral components, and determines the number of reservoir core structure types; obtains the digital core of the target reservoir, slices it, draws the grayscale histogram of the digital core voxels of the corresponding slices, and counts the grayscale value range of the digital core voxels of different core structure types; determines the number of thresholds required for digital core structure analysis; determines the size relationship of the voxel grayscale values ​​of different core structure types; uses the maximum class spacing method to perform multi-threshold optimization calculations on digital core slices to obtain the optimal threshold of the digital core of the conglomerate reservoir; based on the optimal threshold result, performs digital core threshold segmentation, obtains digital core voxel lattices of different core structure types, and completes digital core structure extraction. The present invention can realize the accurate extraction of rock chip fillings, and can effectively characterize the three-dimensional structure of digital cores such as skeleton particles, matrix cements, pores and cracks. This method mainly grasps the pore structure recognition of the internal structure of the core from the scale of mesoscopic images, ignoring the pore structure recognition and characterization problems of macroscopic and microscopic images.

[0095] Embodiment 1

[0096] An embodiment of the present invention provides a method for quantitative characterization of the pore structure of a carbonate reservoir, comprising: multi-scale image pore structure feature extraction, cross-scale image feature region positioning, and multi-scale pore structure comprehensive quantitative characterization. By performing feature analysis and data fusion on multi-scale images, quantitative characterization of the pore structure of a carbonate reservoir based on multi-scale image features is achieved, thereby providing more accurate data for production evaluation of carbonate reservoirs.

[0097] Reference Figure 1 As shown, the method specifically includes:

[0098] S101, extracting pore structure features from the acquired digital images of carbonate reservoirs at different scales, respectively, to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by an optical microscope, and microscopic digital images taken by a scanning electron microscope;

[0099] The multi-scale image pore structure feature extraction as described above mainly involves data of macro-scale digital images taken by high-definition optical interlacing, meso-scale digital images taken by optical microscope, and micro-scale digital images taken by scanning electron microscope. Macro-scale digital images are usually at mm level, meso-scale digital images are usually at um level, and micro-scale digital images are usually at nm level.

[0100] S102, calculating the binary image at each scale based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore variation characteristics, pore fractal characteristics, and pore shape characteristics;

[0101] The above characteristic functions include the autocorrelation function S x (r1, r2), variation function γ(r), fractal dimension function and shape characteristic function. In order to characterize the pore distribution characteristics, the autocorrelation function S is constructed. x (r1, r2), this function can reflect the distribution characteristics of multi-scale pores; in order to characterize the pore change characteristics, the present invention constructs the variation function γ(r); in order to characterize the fractal characteristics of the pores, the present invention constructs the fractal dimension function; in order to reflect the pore shape characteristics, the shape characteristic function is established. Based on the constructed autocorrelation function S x (r1, r2), variation function γ(r), fractal dimension function and shape characteristic function are used to calculate the binary image to obtain the corresponding pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics.

[0102] S103, for each scale of the binary image, taking at least a part of the binary image as the target area, searching the binary images of different scales by a similarity matching algorithm based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, and screening out the pore key area;

[0103] S104. Based on the pre-constructed multi-scale comprehensive fractal dimension function and the key pore regions of different scales, the multi-scale comprehensive fractal dimension of the carbonate reservoir is calculated.

[0104] The present invention relates to an experimental data analysis and characterization method in the field of geology, and in particular to a carbonate reservoir pore structure quantitative characterization method based on multi-scale image features. Specifically, the characteristic parameters of the multi-scale image pore structure are first extracted to grasp the image information of the multi-scale image pore structure, and effectively distinguish pores from skeletons. Subsequently, a multi-scale structural information association and fusion function of the pore structure is constructed, and a key area positioning search strategy for pore structure development characteristics is formulated to achieve cross-scale image feature area positioning, and pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics within a certain range are combined, and a similarity matching algorithm is used to achieve different scales. The search for key pore areas between different areas at the same scale can effectively realize the automatic search for key areas in the "millimeter-micrometer-nanometer" scale images of carbonate reservoir pores. Finally, the multi-scale comprehensive fractal dimension of carbonate reservoirs was constructed, which can reflect the influence of the fractal characteristics of pores in different pore diameter segments on the overall fractal characteristics of carbonate reservoirs, realize the comprehensive quantitative characterization of multi-scale pore structure, and finally solve the correlation and fusion characterization problems between multi-scale images, which can provide data support for effectively solving the production capacity prediction and evaluation problems of carbonate reservoirs, and then provide new technical means for the development of research theories and technical models for the development characteristics of pore structure of carbonate reservoirs.

[0105] The multi-scale image pore structure feature extraction involved in the present invention mainly utilizes these three digital images to realize the multi-scale image pore structure feature extraction. The multi-scale image pore structure feature extraction mainly includes three steps: digital image preprocessing, pore structure recognition, and feature parameter calculation.

[0106] In an optional embodiment, the digital images of carbonate reservoirs at different scales acquired in the above step S101 are respectively subjected to pore structure feature extraction to obtain corresponding binary images, including:

[0107] Digital images of carbonate reservoirs at every scale:

[0108] S1011, preprocessing the digital images respectively to obtain preprocessed digital images;

[0109] S1012, using an iterative image threshold segmentation algorithm to identify the pore structure of the preprocessed digital image to obtain a binary image.

[0110] As described above, the traditional image enhancement method for digital image preprocessing is not effective, has a low peak signal-to-noise ratio, is not very applicable, and does not make good use of global information to enhance image details. The present invention, based on the combination of global and local, proposes an image enhancement method that uses a variable gamma value for power law transformation, namely, an improved gamma correction image enhancement method. When the gray value of a pixel point in a digital image is greater than the average brightness, the point is considered to be a bright spot; when the gray value of a pixel point is less than the average brightness, the point is considered to be a dark spot; when the gray value of a pixel point is equal to the average brightness, the point is considered to be a smooth point, and the average brightness of the image is extracted as a feature component. In an optional embodiment, the digital image is preprocessed separately to obtain a preprocessed digital image, including:

[0111] The characteristic components of the digital image are extracted based on the following formula:

[0112]

[0113] Among them, v is the feature component, L is the image grayscale, g is the pixel grayscale value, i g is the number of pixels of image gray level g, M and N are the total number of rows and columns of the image matrix corresponding to the digital image respectively;

[0114] Based on the characteristic components, the digital image is subjected to image enhancement processing by the following formula to obtain a preprocessed digital image:

[0115]

[0116] Among them, Z(g) is the enhancement function and E is the enhancement parameter;

[0117] In order to reduce the enhancement amplitude of pixels that are brighter or darker, and to increase the enhancement amplitude of relatively smooth pixels, the nonlinear transformation function as shown above is used to quantify the enhancement amplitude of pixels. Among them, the value range of Z(g) is (0,1), and E is mainly used to control the slope of the enhancement function Z(g), which is usually 1 and can be adjusted accordingly according to actual conditions. When the gray value is closer to the average value v, Z(g) is closer to 0.5, indicating that the pixel is relatively smooth and the contrast is relatively small; when the gray value is farther away from the average value v, Z(g) is closer to 0 or 1, indicating that the pixel is darker or brighter and the contrast is greater. The grayscale PF expression of the image after improved gamma correction is:

[0118] PF out =cPF in (1.5-Z(g))

[0119] , formula 3;

[0120] Among them, PF inand PF out are the grayscale values ​​of the input and output images respectively. c is a parameter that controls the shape of the conversion curve. It is usually set to 1. If the correction is not ideal, its value can be appropriately increased. The denoising effect of the median filter is relatively ideal. It can filter out noise points with obvious differences in pixel values, and can effectively prevent the blurring of image details such as contours and edges. It is better than the mean filter in suppressing image salt and pepper noise and pulse interference, and can ensure the clarity of the edge, thereby realizing the whole process of digital image preprocessing and converting the original digital image F1 into the preprocessed digital image F2.

[0121] The key to the pore structure recognition as described above is to find the best segmentation threshold of the image and binarize it, and finally perform morphological operation on the binary image to eliminate defects, and use an iterative image threshold segmentation algorithm to perform image segmentation on the pre-processed digital image F2 to obtain a binary image F3. In an optional embodiment, the above step S1012 uses an iterative image threshold segmentation algorithm to perform pore structure recognition on the pre-processed digital image to obtain a binary image, including:

[0122] S10121, based on the minimum grayscale value PF in the preprocessed digital image out | min and maximum gray value PF out | max , the initial threshold T0 is determined by the following formula;

[0123]

[0124] S10122, dividing the preprocessed digital image into a target part and a background part according to the initial threshold, and updating the threshold by the following formula based on the grayscale values ​​of the target part and the background part to obtain an updated threshold;

[0125]

[0126]

[0127] Among them, PF out | O is the average gray value of all pixels in the target part, PF out | G is the average gray value of all pixels in the background, PF out (i, j) is the gray value of point (i, j) on the image, N(i, j) is the weight coefficient of point (i, j), T K is the threshold before updating, T K+1 is the updated threshold;

[0128] S10123, segmenting the preprocessed digital image based on the updated threshold value, repeating the process of updating the threshold value until a preset condition is met to obtain an optimal threshold value;

[0129] If T is not reached K =T K+1 , then the above S10122-S10123 is cycled. When T K =T K+1 The cycle ends and the optimal threshold is obtained for segmentation.

[0130] S10124. Segment the preprocessed digital image based on the optimal threshold to obtain a binary image.

[0131] Pore ​​structure segmentation image by phase function For a two-phase system considering only pores and skeletons, when When located in the pore, the value is 1. When in the skeleton phase, the value is 0.

[0132] As described above, the characteristic parameter calculation is performed by using the maximum circle search method to search the radius of the pore area (white area) in the binary image F3. If the number of pixels with the maximum inscribed circle radius of the pore area is N r , the radius of the pore area is represented by R. In an optional embodiment, the method further includes:

[0133] S10131, using the maximum circle search method to perform radius search on the pore area of ​​the binary image, and determining the pore area radius R of the single pore structure by the following formula:

[0134]

[0135] Where R is the radius of the pore area, δd is the actual size represented by each pixel, and N r is the number of pixels of the maximum inscribed circle radius of the pore region. The longest axis and the shortest axis of the boundary of the single pore structure can also be determined. The longest axis of the boundary is defined as the Euclidean distance between the two farthest points on the outer boundary of the pore region. The line segment connecting these points is called the longest axis of the boundary of the cell cavity. The short axis of the boundary is defined as the line segment perpendicular to the longest axis of the boundary, denoted as d. The method for identifying the outer boundary of the pore region can refer to the description in the prior art, which will not be repeated here.

[0136] As mentioned above, the pore area represents the area of ​​the pixel set in the white closed area of ​​the image, which is represented by S. The pore area is divided into n×n grids, and the leftmost point of each row of the grid is P(x i ,y ij ), the rightmost point is Q(xi ,y' ij ), then the pixel points contained in this row of grid are y ij -y' ij +1, in an optional embodiment, the method further includes:

[0137] S10132, determining the area of ​​the pore region in the binary image based on the following formula:

[0138]

[0139] Where S is the area of ​​the pore region, δd is the actual size represented by each pixel, n×n is the number of grids in the pore region, and y ij is the ordinate of the leftmost point in each row, y' ij is the ordinate of the rightmost point in each row, y ij -y' ij +1 is the number of pixels contained in each row of the grid, i = 1, 2, ..., n.

[0140] The pore area represents the side length of the area formed by the outermost pixel set in the white closed area in the image, which is represented by C. The perimeter of the pore area is the length of the calculated chain code. The characteristic parameter calculation as described above also includes:

[0141] By traversing the pore area contour point by point, the contour length calculated statistically is the perimeter of the pore area. Specifically, initialize the perimeter C = 0, traverse all chain code elements of a single contour, and if the chain code value is an odd number, then If the chain code value is an even number, then C = C + 1. For each contour chain code, the above operation is performed, and the sum of all contour perimeters is the total pore area perimeter C. For macroscopic pores, the corresponding pore area perimeter is 2C / π.

[0142] The cross-scale image feature region localization as described above mainly includes two steps: feature function construction and pore key region search.

[0143] The characteristic function construction mentioned above, in order to characterize the pore distribution characteristics, construct the autocorrelation function S x (r1, r2), which can reflect the statistical characteristics of multi-scale pores, and the autocorrelation function S x (r1, r2) is defined as the probability that any two points in the same scale image are distributed in the same pore.

[0144] In an optional embodiment, the autocorrelation function S x (r1, r2) represents the correlation between any two points in the same scale image. For a two-phase system that only considers pores and skeletons, the pore distribution characteristics are calculated by the following formula:

[0145]

[0146] Among them, S x (r1, r2) is the pore distribution feature, r1 and r2 are any two points in the same scale image, specifically any two points with a distance of r in the same scale image, P(r1) is the phase function of point r1, and P(r2) is the phase function of point r2. When the distance r = 0, the autocorrelation function S x The value of (r1, r2) is the porosity of the image; when the distance r increases to a certain value, the autocorrelation function tends to be stable.

[0147] In order to reflect the degree of change of the regionalized variable within a certain distance range in a certain direction, a variogram γ(r) is constructed. The variogram of the same scale image is closely related to its structure. The variogram γ(r) is an indispensable function for evaluating the structural properties of the image. In an optional embodiment, the pore change characteristics are calculated by the following formula:

[0148]

[0149] Among them, γ(r) is the pore change characteristic, r is the equivalent radius of the pore, N(r) is the number of two points with a distance r in the same scale image, and r i is the radius of the ith pore, P(r i ) is r i The phase function value of .

[0150] In order to characterize the fractal characteristics of pores, a fractal dimension function D is constructed. F , the perimeter and area of ​​each pore are measured in the binary image F2, and twice the slope of the double logarithmic plot of the area and perimeter is the fractal dimension function D F According to the fractal description, the relationship between the pore size distribution fractal and the fractal dimension can be obtained. In an optional embodiment, the pore fractal feature is calculated by the following formula:

[0151]

[0152] Among them, D F is the fractal characteristic of pores, r is the equivalent radius of the pore, N'(r) is the number of pores with a pore size greater than r, and a is a constant, usually taking the natural logarithm.

[0153] In order to reflect the pore shape characteristics, the shape characteristic function G is established. T Since the shape of the pore area in carbonate reservoirs is very complex, in order to make the simplified pores have the same geometric characteristics as the real pores, it is necessary to evaluate the shape characteristics of the real pore area, that is, to obtain the shape characteristic function G of the pore area TIn an optional embodiment, the pore shape characteristic is calculated by the following formula:

[0154]

[0155] Among them, G T is the pore shape characteristic, S T is the area of ​​the Tth pore region; C T is the perimeter of the Tth pore area.

[0156] The pore key area search described above combines the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics within a certain range, and uses a similarity matching algorithm to realize the pore key area search between digital images of different scales and between different areas of digital images of the same scale. Let q be the target area, t be the search area, S x represents the similarity of pore distribution characteristics between the search area and the target area, γ represents the similarity of pore change characteristics between the search area and the target area, D represents the similarity of pore fractal characteristics between the search area and the target area, and G represents the similarity of pore shape characteristics between the search area and the target area. The present invention takes the pore key area search between digital images of different scales as an example for explanation.

[0157] In an optional embodiment, for the binary image of each scale described in step S103, at least a part of the binary image is used as the target area, and based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a similarity matching algorithm is used to search the binary images of different scales to screen out the pore key area, including:

[0158] For each scale of binary image, perform the following operations:

[0159] S1031, taking at least a part of the binary image as a target area, and selecting search areas with the highest comprehensive similarity to the target area in binary images of different scales through a similarity matching algorithm based on pore distribution characteristics, pore change characteristics, pore fractal characteristics, and pore shape characteristics of the target area, to obtain a corresponding target area;

[0160] In the search for key pore areas at different scales, first, the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristic parameters are manually set, and the system automatically matches the target area q on the specified scale image. Or by selecting the specified area image as the target area q. The pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area are calculated respectively by equations 9 to 12. Subsequently, the regional image of the corresponding scale is selected in the images of different scales as the search area t. For each scale image, taking the macroscopic digital image as an example, the macroscopic digital image, the mesoscopic digital image and the microscopic digital image are gridded to obtain the corresponding grid image. The regional image of the corresponding scale is selected as the search area t in the grid image corresponding to the mesoscopic scale and the digital image corresponding to the microscopic scale, respectively. Specifically, each grid in the grid image of each scale can be traversed as the search area. The pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics in the same search area are combined to calculate the comprehensive similarity between the search area and the target area, and the first target area with the highest comprehensive similarity to the target area is screened in the mesoscopic grid image, and the second target area with the highest comprehensive similarity to the target area is screened in the microscopic grid image. The first target area and the second target area are the target areas corresponding to the macroscopic digital image. Similarly, the target area corresponding to the mesoscopic digital image and the target area corresponding to the microscopic digital image can be obtained.

[0161] S1032, respectively determining the Euclidean distances between the target area and the target areas of different scales, selecting the area with the smallest Euclidean distance, and obtaining the pore key area.

[0162] The pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics in the same target area are combined and the Euclidean distance L qt’ The size of is used to sort the distances of target areas at different scales, and the closest distance is the key pore area that needs to be searched and located.

[0163] Specifically, for each scale of digital image, the Euclidean distance between its target area and each target area is calculated, and the target area with the smallest Euclidean distance to the target area is selected to obtain the pore key area corresponding to the digital image of this scale. The Euclidean distance is used to measure the comprehensive similarity DC, and the Euclidean distance L from the target area q to the target area t' is calculated. qt’ , the Euclidean distance L qt’ Represents the actual distance between two regions. This distance combines the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics. Euclidean distance L qt’ The expression is:

[0164]

[0165] In the above formula, q1, q2, ..., q M Respectively represent the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area, t'1, t'2, ..., t' M They represent the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area respectively, and M represents the number of target areas.

[0166] In an optional embodiment, the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area described in step S1031 above are used to select search areas with the highest comprehensive similarity to the target area in binary images of different scales through a similarity matching algorithm to obtain the corresponding target area, including:

[0167] S10311, for each scale of the binary image, gridding the binary image to obtain a grid image;

[0168] S10312, respectively calculating the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of each search area in the grid image;

[0169] In the different scale images, the regional images of the corresponding scales are selected as the search area t, and specifically, each grid in the grid image can be used as the search area. The pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of each search area are calculated respectively by equations 9 to 12.

[0170] S10313, based on the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area, and the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the search area, respectively calculate the pore distribution characteristics similarity, pore variation characteristics similarity, pore fractal characteristics similarity and pore shape characteristics similarity of the target area and the search area, and respectively calculate the comprehensive similarity of the target area and each of the search areas by the following formula:

[0171]

[0172] Among them, DC(q,t) is the comprehensive similarity between the target area and the search area, q is the target area, t is the search area, S x is the similarity of pore distribution characteristics, γ is the similarity of pore change characteristics, D is the pore fractal characteristics, G is the similarity of pore shape characteristics, λ a , b ,c and λ d are the weights of pore distribution feature similarity, pore change feature similarity, pore fractal feature similarity and pore shape feature similarity respectively; wherein, the pore distribution feature similarity, pore change feature similarity, pore fractal feature similarity and pore shape feature similarity are all calculated by using the similarity calculation formula, which will not be described here;

[0173] S10314. Select search areas with the highest comprehensive similarity to the target area in binary images of different scales to obtain corresponding target areas.

[0174] For each scale image, taking the macro-scale digital image as an example, the macro-scale digital image, the meso-scale digital image and the micro-scale digital image are grid-divided to obtain the corresponding grid image. The regional images of the corresponding scale are selected as the search area t in the grid image corresponding to the meso-scale and the digital image corresponding to the micro-scale, respectively. Specifically, each grid in the grid image of each scale can be traversed as the search area. The pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics in the same search area are combined to calculate the comprehensive similarity between the search area and the target area, and the first target area with the highest comprehensive similarity to the target area is selected in the meso-scale grid image, and the second target area with the highest comprehensive similarity to the target area is selected in the micro-scale grid image. The first target area and the second target area are the target areas corresponding to the macro-scale digital image. Similarly, the target area corresponding to the meso-scale digital image and the target area corresponding to the micro-scale digital image can be obtained.

[0175] As mentioned above, for the comprehensive quantitative characterization of multi-scale pore structure, the pore volume calculation of the digital image F2 assumes that the image is a rectangular parallelepiped with a two-dimensional image as the bottom plane and the maximum grayscale of the pixel as the height. Each pixel can be regarded as a column formed by its own bottom plane as the bottom and its own grayscale as the height. The sum of the volumes of all columns can be used as the skeleton volume, and the remaining volume is the pore volume. The pore volume ratio can reflect the degree of influence of the fractal characteristics of pores in different pore sizes on the overall fractal characteristics of carbonate reservoirs. Therefore, the fractal dimensions of pores of different scales are weighted and calculated using the pore volume of pores in each pore size segment, and the result is the multi-scale comprehensive fractal dimension d of carbonate reservoirs. z In an optional embodiment, the multi-scale integrated fractal dimension of the carbonate reservoir is calculated based on the pre-constructed multi-scale integrated fractal dimension function and the pore key areas of different scales in step S104, including:

[0176] S1041. Calculate the pores of different pore diameters in the key pore areas of different scales based on the following formula to obtain the multi-scale comprehensive fractal dimension of the carbonate reservoir:

[0177]

[0178] Where: d z is the multi-scale comprehensive fractal dimension of carbonate reservoir, B i is the ratio of the pore volume of the pore key area of ​​the i-th scale to the total pore volume. The total pore volume refers to the total pore volume in the digital image of this scale, i = 1, 2, ..., m, m is the number of scales, d ij is the fractal dimension of the jth aperture segment in the pore critical region of the i-th scale, j = 1, 2, ..., l, l is the number of aperture segments divided in the pore critical region of the i-th scale, b ij It is the ratio of the pore volume of the jth pore size segment in the pore critical region of the i-th scale to the pore volume of the pore critical region of this scale.

[0179] The method for quantitative characterization of carbonate reservoir pore structure based on multi-scale image features provided by the embodiment of the present invention mainly realizes the following functions and effects: first, by preprocessing the multi-scale pore structure image, the image information of the multi-scale image pore structure can be grasped, thereby effectively distinguishing pores from skeletons, realizing the extraction of characteristic parameters of the pore structure, and laying a data foundation for the fusion and characterization of the multi-scale pore structure; second, by constructing the multi-scale structural information association and fusion function of the pore structure, a key area positioning search strategy for pore structure development characteristics is formulated, which can realize the cross-scale image feature area positioning, and the pore distribution characteristics, pore change characteristics, and pores within a certain range are combined. The fractal characteristics and pore shape characteristics are combined, and the similarity matching algorithm is used to realize the search of key pore areas between different scales and different areas of the same scale. It can effectively realize the automatic search of key areas in the "millimeter-micrometer-nanometer" scale image of carbonate reservoir pores. Third, the multi-scale comprehensive fractal dimension of carbonate reservoirs is constructed, which can reflect the influence of the fractal characteristics of pores in different pore diameter sections on the overall fractal characteristics of carbonate reservoirs, and realize the comprehensive quantitative characterization of multi-scale pore structure. Finally, the correlation and fusion characterization problems between multi-scale images are solved, which can provide data support for effectively solving the production capacity prediction and evaluation problems of carbonate reservoirs. Compared with traditional methods, a quantitative characterization method of carbonate reservoir pore structure based on multi-scale image features is scientific, reliable and practical.

[0180] Embodiment 2

[0181] In order to more clearly illustrate the method for quantitatively characterizing the pore structure of carbonate reservoirs provided in an embodiment of the present invention, the present invention applies the method for quantitatively characterizing the pore structure of carbonate reservoirs to quantitatively characterize the pore structure of a certain carbonate reservoir, as follows.

[0182] A method for quantitative characterization of carbonate reservoir pore structure based on multi-scale image features specifically includes: multi-scale image pore structure feature extraction 1, cross-scale image feature region positioning 2, multi-scale pore structure comprehensive quantitative characterization 3. By performing feature analysis and data fusion on multi-scale images, quantitative characterization of carbonate reservoir pore structure based on multi-scale image features is achieved, providing more accurate data for carbonate reservoir production evaluation.

[0183] As mentioned above, the multi-scale image pore structure feature extraction 1 mainly involves data of macro-scale digital images taken by high-definition optical interlacing, meso-scale digital images taken by optical microscope, and micro-scale digital images taken by scanning electron microscope. The macro-scale digital images are usually at mm level, the meso-scale digital images are usually at um level, and the micro-scale digital images are usually at nm level. Therefore, the multi-scale image pore structure feature extraction 1 involved in the present invention mainly uses these three digital images to realize the pore structure feature extraction of multi-scale images. The multi-scale image pore structure feature extraction 1 mainly includes three steps: digital image preprocessing 1.1, pore structure identification 1.2, and characteristic parameter calculation 1.3.

[0184] As described above in digital image preprocessing 1.1, the traditional image enhancement method has insignificant effect, low peak signal-to-noise ratio, weak applicability, and does not make good use of global information to enhance image details. The present invention, based on the combination of global and local, proposes an image enhancement method that uses a variable gamma value for power law transformation, namely, an improved gamma correction image enhancement method. When the gray value of a pixel point in a digital image is greater than the average brightness, the point is considered to be a bright point; when the gray value of a pixel point is less than the average brightness, the point is considered to be a dark point; when the gray value of a pixel point is equal to the average brightness, the point is considered to be a smooth point, and the average brightness of the image is extracted as a feature component. The calculation formula of the feature component is as shown in Formula 1 above.

[0185] In order to reduce the enhancement amplitude of pixels that are brighter or darker and increase the enhancement amplitude of relatively smooth pixels, a nonlinear transformation function is used to quantify the enhancement amplitude of the pixels. The expression of the enhancement function Z(g) is as shown in Formula 2 above.

[0186] When the gray value is closer to the average value v, Z(g) is closer to 0.5, indicating that the pixel is relatively smooth and the contrast is relatively small; when the gray value is farther away from the average value v, Z(g) is closer to 0 or 1, indicating that the pixel is darker or brighter and the contrast is greater. The grayscale PF expression of the image after improved gamma correction is as shown in formula 3 above. The denoising effect of median filtering is ideal. It can filter out noise points with obvious pixel value differences and effectively prevent image details such as contours and edges from blurring. It is better than mean filtering in suppressing image salt and pepper noise and pulse interference, and can ensure edge clarity, thereby realizing the entire process of digital image preprocessing 1.1, converting the original digital image F1 into the preprocessed digital image F2.

[0187] As mentioned above, the key to pore structure recognition 1.2 is to find the best segmentation threshold of the image and binarize it, and finally perform morphological operation on the binary image to eliminate defects. The iterative image threshold segmentation algorithm is used to segment the digital image F2 to obtain the binary image F3. The iterative image threshold segmentation algorithm steps can refer to the above steps S10121-S10124, and the best threshold can be obtained for segmentation. After completing the threshold binarization segmentation of the image, the target boundary formed is often very uneven. Mathematical morphology is used to smooth the object contour, wherein the open operation is used to remove pepper noise, and the closed operation is used to remove sandstone noise. If the open and closed operations are used continuously, the foreground and background noise in the binary image can be significantly improved. The pore structure segmentation image is obtained by the phase function For a two-phase system considering only pores and skeletons, when When located in the pore, the value is 1. When in the skeleton phase, the value is 0.

[0188] As described above, the feature parameter calculation 2.1 uses the maximum circle search method to search the radius of the pore area (white area) in the binary image F3. If the number of pixels of the maximum inscribed circle radius of the pore area is N r , the radius of the pore area, represented by R, the radius of the pore area of ​​the single pore structure is calculated by the above formula 7.

[0189] The longest axis of the boundary is defined as the Euclidean distance between the two farthest points on the outer boundary of the pore region. If the longest axis of the boundary is a useful descriptor, it is best applied to boundaries with a single pair of farthest points. When there is more than one pair of farthest points, they should be close to each other and are the main parameters that determine the shape of the boundary. The line segment connecting these points is called the longest axis of the cell cavity boundary. The short axis of the boundary is defined as the line segment perpendicular to the longest axis of the boundary, denoted by d. The relationship between the longest axis of the pore boundary and its short axis is shown in Figure 2 As shown,

[0190] The pore area refers to the area of ​​the area formed by the set of pixels in the white closed area in the image, which is represented by S. The pore area is divided into n×n grids, and the leftmost point in each row is P(x i ,y ij ), the rightmost point is Q(x i ,y' ij ), then the pixel points contained in this row of grid are y ij -y' ij +1, the area S of the pore region is calculated by the above formula 8.

[0191] The pore area represents the side length of the area formed by the outermost pixel set in the white closed area in the image, which is represented by C. The perimeter of the pore area is the length of the calculated chain code. By traversing the contour of the pore area point by point, the length of the contour calculated by statistics is the perimeter of the pore area. Initialize the perimeter C = 0, traverse all the chain code elements of a single contour, if the chain code value is an odd number, then If the chain code value is an even number, then C = C + 1. For each contour chain code, the above operation is performed, and the sum of all contour perimeters is the total pore area perimeter C. For macroscopic pores, the corresponding pore area perimeter is 2C / π.

[0192] The cross-scale image feature region localization 2 as described above mainly includes two steps: feature function construction 2.1 and key region search 2.2.

[0193] As mentioned above, the characteristic function is constructed in 2.1. In order to characterize the pore distribution characteristics, the autocorrelation function S is constructed. x (r1, r2), which can reflect the statistical characteristics of multi-scale pores, and the autocorrelation function S x (r1, r2) is defined as the probability that any two points in the same scale image are distributed in the same pore. The autocorrelation function S x (r1, r2) represents the correlation between any two points in the same scale image. For a two-phase system considering only pores and skeletons, its expression is shown in Equation 9 above.

[0194] In order to reflect the degree of change of the regionalized variable within a certain distance range in a certain direction, the variogram γ(r) is constructed. The variogram of the same scale image is closely related to its structure. The variogram γ(r) is an indispensable function for evaluating the structural properties of the image. The expression of the variogram γ(r) is shown in Equation 10 above.

[0195] In order to characterize the fractal characteristics of pores, a fractal dimension function D is constructed. F , the perimeter and area of ​​each pore are measured in the binary image F2, and twice the slope of the double logarithmic plot of the area and perimeter is the fractal dimension function D FThe specific calculation formula is the above formula 11.

[0196] In order to reflect the pore shape characteristics, the shape characteristic function G is established. T Since the shape of the pore region of carbonate reservoirs is very complex, in order to make the simplified pores have the same geometric characteristics as the real pores, it is necessary to evaluate the shape characteristics of the real pore region, that is, to construct the shape characteristic function G of the pore region as shown in the above formula 12. T .

[0197] As mentioned above, the key area search 2.2 combines the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics in a certain search area, and uses the similarity matching algorithm to realize the pore key area search between different scales and between different areas of the same scale. Let q be the target area, t be the search area, S x represents the similarity of pore distribution characteristics between the search area and the target area, γ represents the similarity of pore change characteristics between the search area and the target area, D represents the similarity of pore fractal characteristics between the search area and the target area, and G represents the similarity of pore shape characteristics between the search area and the target area. The comprehensive similarity DC(q,t) between the two regional images can be calculated by using the above formula 14. The target area t' corresponding to the digital image of each scale is screened by the comprehensive similarity DC(q,t).

[0198] The Euclidean distance is used to measure the comprehensive similarity DC, and the Euclidean distance L from the target area q to the destination area t' is calculated. qt’ , the Euclidean distance L qt’ Represents the actual distance between two regions. This distance combines the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics. Euclidean distance L qt’ It can be calculated by the above formula 13.

[0199] In the search for key pore areas at different scales, first, the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristic parameters are set manually, and the system automatically matches the target area q on the specified scale image. Or, by selecting the specified area image as the target area q. The pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area are calculated respectively by equations 9 to 12. Then, the regional image of the corresponding scale is selected in the images of different scales as the search area t. Finally, the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics within a certain range, that is, the above-mentioned collection area, are combined, and the Euclidean distance L is used. qt’ The size of is used to sort the distances of target areas at different scales, and the closest distance is the key pore area that needs to be searched and located.

[0200] As mentioned above, the multi-scale comprehensive quantitative characterization of pore structure 3, the pore volume calculation of the digital image F2 assumes that the image is a rectangular parallelepiped with a two-dimensional image as the bottom plane and the maximum grayscale of the pixel as the height. Each pixel can be regarded as a column formed by its own bottom plane as the bottom and its own grayscale as the height. The sum of the volumes of all columns can be used as the skeleton volume, and the remaining volume is the pore volume. The pore volume ratio can reflect the degree of influence of the fractal characteristics of pores in different pore sizes on the overall fractal characteristics of carbonate reservoirs. Therefore, the fractal dimensions of pores of different scales are weighted and calculated with the pore volume of the pores in each pore size segment as the weight, and the result is the multi-scale comprehensive fractal dimension dz of the carbonate reservoir. The specific calculation formula is shown in the above formula 15.

[0201] As described above, the cross-scale image feature region localization 2 first divides the pre-processed digital image F2 into grids of 10*10 and 100*100 respectively. Figure 3 As shown, (a) is a 10*10 grid macro image, (b) is a 10*10 grid meso image, (c) is a 10*10 grid micro image, (d) is a 100*100 grid macro image, (e) is a 100*100 grid meso image, and (f) is a 100*100 grid micro image.

[0202] In the 10*10 grid image key area search, the search for the corresponding area in the macroscopic image is determined based on the mesoscopic image. Figure 4 As shown in the figure, first, the mesoscopic image is selected as the target area, then the regional image matrix [1,1] in the macroscopic image is used as the search area, and the two regional images are subjected to comprehensive similarity calculation. After the search of the regional image matrix [1,1] in the macroscopic image is completed, the regional image matrix [1,2] in the macroscopic image is used as the search area, and so on, until the search of the regional image matrix [10,10] in the macroscopic image is completed. The search results are shown in FIG. Figure 4 As shown in (a). In the key area search of the 10*10 grid image, the search for the corresponding area in the mesoscopic image is determined based on the microscopic image. First, the microscopic image is selected as the target area, and then the regional image matrix [1,1] in the mesoscopic image is used as the search area, and the two regional images are subjected to a comprehensive similarity calculation. After completing the search for the regional image matrix [1,1] in the mesoscopic image, the regional image matrix [1,2] in the mesoscopic image is used as the search area, and so on, until the search for the regional image matrix [10,10] in the mesoscopic image is completed. The search results are shown in FIG. Figure 4As shown in (b). In the key area search of the 10*10 grid image, the search for the corresponding area in the macro image is determined based on the micro image. First, the micro image is selected as the target area, and then the regional image matrix [1,1] in the macro image is used as the search area, and the two regional images are subjected to a comprehensive similarity calculation. After completing the search of the regional image matrix [1,1] in the macro image, the regional image matrix [1,2] in the macro image is used as the search area, and so on, until the search of the regional image matrix [10,10] in the macro image is completed. The search results are shown in FIG. Figure 4 (c). Figure 4 As can be seen from (a), the maximum value of the comprehensive similarity is 0.6104, and the corresponding search matrix is ​​the 5th row and 7th column, indicating that the regional image matrix [5,7] in the macroscopic image has the highest comprehensive similarity with the mesoscopic image. Figure 4 As can be seen from (b), the maximum value of the comprehensive similarity is 0.5824, and the corresponding search matrix is ​​the 10th row and 10th column, indicating that the regional image matrix [10,10] in the mesoscopic image has the highest comprehensive similarity with the microscopic image. Figure 4 As can be seen in (c), the maximum value of the comprehensive similarity is 0.6099, and the corresponding search matrix is ​​the 7th row and 9th column, indicating that the regional image matrix [7,9] in the macro image has the highest comprehensive similarity with the micro image.

[0203] In the 100*100 grid image key area search, the search for the corresponding area in the macroscopic image is determined based on the mesoscopic image. Figure 5 As shown in the figure, first, the mesoscopic image is selected as the target area, then the regional image matrix [1,1] in the macroscopic image is used as the search area, and the two regional images are subjected to comprehensive similarity calculation. After the search of the regional image matrix [1,1] in the macroscopic image is completed, the regional image matrix [1,2] in the macroscopic image is used as the search area, and so on, until the search of the regional image matrix [100,100] in the macroscopic image is completed. The search results are shown in FIG. Figure 5 As shown in (a). In the key area search of the 100*100 grid image, the search for the corresponding area in the mesoscopic image is determined based on the microscopic image. First, the microscopic image is selected as the target area, and then the regional image matrix [1,1] in the mesoscopic image is used as the search area, and the two regional images are subjected to a comprehensive similarity calculation. After completing the search for the regional image matrix [1,1] in the mesoscopic image, the regional image matrix [1,2] in the mesoscopic image is used as the search area, and so on, until the search for the regional image matrix [100,100] in the mesoscopic image is completed. The search results are shown in FIG. Figure 5As shown in (b). In the key area search of the 100*100 grid image, the search for the corresponding area in the macro image is determined based on the micro image. First, the micro image is selected as the target area, and then the regional image matrix [1,1] in the macro image is used as the search area, and the two regional images are subjected to a comprehensive similarity calculation. After completing the search of the regional image matrix [1,1] in the macro image, the regional image matrix [1,2] in the macro image is used as the search area, and so on, until the search of the regional image matrix [100,100] in the macro image is completed. The search results are shown in FIG. Figure 5 (c) as shown.

[0204] from Figure 5 As can be seen from (a), the maximum value of the comprehensive similarity is 0.8220, and the corresponding search matrix is ​​row 81 and column 35, indicating that the regional image matrix [81,35] in the macroscopic image has the highest comprehensive similarity with the mesoscopic image. Figure 5 As can be seen from (b), the maximum value of the comprehensive similarity is 0.8216, and the corresponding search matrix is ​​row 17 and column 77, indicating that the regional image matrix [17,77] in the mesoscopic image has the highest comprehensive similarity with the microscopic image. Figure 5 As can be seen in (c), the maximum comprehensive similarity value is 0.7644, and the corresponding search matrix is ​​the 85th row and 49th column, indicating that the regional image matrix [85,49] in the macro image has the highest comprehensive similarity with the micro image.

[0205] In order to verify the accuracy of the method in locating the cross-scale image feature area, the search results of the key area of ​​the 10*10 grid image, the search results of the key area of ​​the 100*100 grid image and the area where the real extracted image is located are compared. The comparison diagram of the search location area and the real sampling area is shown in the figure below. Figure 6 shown. Figure 6 (a) The sampling area corresponding to the macroscopic scale is searched through the mesoscopic image. Figure 6 (b) The sampling area corresponding to the mesoscopic scale is searched through microscopic images. Figure 6 (c) The microscopic image is used to search for the sampling area corresponding to the macroscopic scale. The darkest area in the figure represents the actual sampling area, the darker area represents the 10*10 grid search positioning area, and the light area represents the 100*100 grid search positioning area. Figure 6 As can be seen in (a), when searching for the sampling area corresponding to the macroscopic scale through the mesoscopic image, a 100*100 grid is used for searching, and the estimated positioning area is closer to the actual sampling area. Figure 6As can be seen in (b), when searching for the sampling area corresponding to the mesoscopic scale through microscopic images, a 100*100 grid is used for searching, and the estimated positioning area is closer to the actual sampling area. Figure 6 As can be seen in (c), when searching the sampling area corresponding to the macro scale through the microscopic image, the 100*100 grid is used for search, and the estimated positioning area is closer to the actual sampling area. This is mainly because the actual sampling size area is smaller, and the 100*100 grid division is more accurate than the 10*10 grid division, and the search area after division is closer to the actual area. In the search target image using the macroscale image as the target area, the reliability of using the mesoscale image as the target area is higher than that of the microscale image as the target area, which is mainly because the mesoscale is closer to the pore structure characteristics of the macroscale than the microscale.

[0206] Embodiment 3

[0207] The embodiment of the present invention provides a carbonate reservoir pore structure quantitative characterization device, referring to Figure 7 As shown, including:

[0208] The feature extraction module 201 is used to extract pore structure features from the acquired digital images of carbonate reservoirs at different scales, respectively, to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by optical microscope, and microscopic digital images taken by scanning electron microscope;

[0209] A first calculation module 202 is used to calculate the binary image of each scale based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics;

[0210] A search module 203 is used for searching binary images of each scale, taking at least a part of the binary image as a target area, searching binary images of different scales by a similarity matching algorithm based on pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area, and screening out key pore areas;

[0211] The second calculation module 204 is used to calculate the multi-scale comprehensive fractal dimension of the carbonate reservoir based on the pre-constructed multi-scale comprehensive fractal dimension function and the pore key areas of different scales.

[0212] The implementation principle and technical effect of the carbonate reservoir pore structure quantitative characterization device provided in the embodiment of the present invention are similar to those of any of the aforementioned method embodiments, and will not be repeated here.

[0213] Embodiment 4

[0214] An embodiment of the present invention provides a computer-readable storage medium having a computer program stored thereon. When the program is executed by a processor, the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in any one of the aforementioned method embodiments is implemented.

[0215] The computer-readable storage medium may be included in the device / apparatus described in the above embodiment; or it may exist independently without being assembled into the device / apparatus. The above computer-readable storage medium carries one or more programs, and when the above one or more programs are executed, the method according to the embodiment of the present invention is implemented.

[0216] According to an embodiment of the present invention, the computer-readable storage medium may be a non-volatile computer-readable storage medium, for example, may include but is not limited to: a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination thereof. In the present invention, the computer-readable storage medium may be any tangible medium containing or storing a program, which may be used by or in combination with an instruction execution system, apparatus, or device.

[0217] Embodiment 5

[0218] An embodiment of the present invention provides an electronic device, referring to Figure 8 As shown, it includes a processor 111, a communication interface 112, a memory 113 and a communication bus 114, wherein the processor 111, the communication interface 112 and the memory 113 communicate with each other through the communication bus 114.

[0219] Memory 113, used for storing computer programs;

[0220] The processor 111 is used to implement the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in any one of the above method embodiments when executing the program stored in the memory 113.

[0221] The implementation principle and technical effect of the electronic device provided by the embodiment of the present invention are similar to those of any of the aforementioned method embodiments, and will not be repeated here.

[0222] The above-mentioned memory 113 can be an electronic memory such as a flash memory, an EEPROM (electrically erasable programmable read-only memory), an EPROM, a hard disk or a ROM. The memory 113 has a storage space for program codes for executing any method steps in the above-mentioned method. For example, the storage space for program codes may include various program codes for implementing various steps in the above method respectively. These program codes can be read from or written into one or more computer program products. These computer program products include program code carriers such as hard disks, compact disks (CDs), memory cards or floppy disks. Such computer program products are usually portable or fixed storage units. The storage unit may have a storage segment or storage space arranged similarly to the memory 113 in the above-mentioned electronic device. The program code can be compressed, for example, in an appropriate form. Generally, the storage unit includes a program for executing the method steps according to an embodiment of the present invention, that is, a code that can be read by, for example, a processor 111, and when these codes are run by an electronic device, the electronic device performs the various steps in the method described above.

[0223] In this document, relational terms such as first and second, etc. are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprises" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. The orientation or position relationship indicated by the terms "upper", "lower", etc. is based on the orientation or position relationship shown in the accompanying drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operate in a specific orientation, and therefore cannot be understood as a limitation on the present invention.

[0224] Each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. It should be noted that the embodiments of the present invention and the features in the embodiments can be combined with each other without conflict. The present invention is not limited to any single aspect, any single embodiment, or any combination and / or permutation of these aspects and / or embodiments. Each aspect and / or embodiment of the present invention can be used alone, or in combination with one or more other aspects and / or other embodiments.

[0225] Finally, it should be noted that the above-described embodiments are only specific implementations of the present invention, which are used to illustrate the technical solutions of the present invention, rather than to limit them. The protection scope of the present invention is not limited thereto. Although the present invention is described in detail with reference to the above-described embodiments, ordinary technicians in the field should understand that any technician familiar with the technical field can still modify the technical solutions recorded in the above-described embodiments within the technical scope disclosed by the present invention, or can easily think of changes, or make equivalent replacements for some of the technical features therein; and these modifications, changes or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention, and should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims

1. A method for quantitatively characterizing the pore structure of a carbonate reservoir, characterized in that: include: Extracting pore structure features from the acquired digital images of carbonate reservoirs at different scales to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by an optical microscope, and microscopic digital images taken by a scanning electron microscope; The binary image of each scale is calculated based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics; For each scale of binary images, at least a part of the binary image is used as a target area, and based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a similarity matching algorithm is used to search the binary images of different scales to screen out the pore key area; Based on the pre-constructed multi-scale comprehensive fractal dimension function and the key pore areas at different scales, the multi-scale comprehensive fractal dimension of the carbonate reservoir is calculated.

2. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: For each scale of the binary image, at least a part of the binary image is used as the target area, and based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, the binary images of different scales are searched by a similarity matching algorithm to screen out the pore key area, including: For each scale of binary image, perform the following operations: Using at least a portion of the binary image as a target area; Based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a search area with the highest comprehensive similarity to the target area is selected in binary images of different scales by a similarity matching algorithm to obtain a corresponding target area; The Euclidean distances between the target area and the target areas of different scales are determined respectively, and the area with the smallest Euclidean distance is selected to obtain the pore key area.

3. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 2, characterized in that: Based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, a search area with the highest comprehensive similarity to the target area is selected in binary images of different scales by a similarity matching algorithm to obtain a corresponding target area, including: For each scale of the binary image, gridding the binary image to obtain a grid image; respectively calculating the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of each search area in the grid image; Based on the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the target area, and the pore distribution characteristics, pore variation characteristics, pore fractal characteristics and pore shape characteristics of the search area, the pore distribution characteristics similarity, pore variation characteristics similarity, pore fractal characteristics similarity and pore shape characteristics similarity of the target area and the search area are calculated respectively, and the comprehensive similarity of the target area and each of the search areas is calculated respectively by the following formula: Among them, DC(q,t) is the comprehensive similarity between the target area and the search area, q is the target area, t is the search area, S x is the similarity of pore distribution characteristics, γ is the similarity of pore change characteristics, D is the similarity of pore fractal characteristics, G is the similarity of pore shape characteristics, λ a , b , c and λ d are the weights of pore distribution feature similarity, pore change feature similarity, pore fractal feature similarity and pore shape feature similarity respectively; In binary images of different scales, search regions having the highest comprehensive similarity to the target region are respectively selected to obtain corresponding target regions.

4. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The pore distribution characteristics are calculated by the following formula: Among them, S x (r1, r2) is the pore distribution characteristic, r1 and r2 are any two points in the same scale image, P(r1) is the phase function of point r1, and P(r2) is the phase function of point r2.

5. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The pore change characteristics are calculated by the following formula: Among them, γ(r) is the pore change characteristic, r is the equivalent radius of the pore, N(r) is the number of two points with a distance r in the same scale image, and r i is the radius of the ith pore, P(r i ) is r i The phase function value of .

6. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The pore fractal characteristics are calculated by the following formula: Among them, D F is the pore fractal feature, r is the equivalent radius of the pore, N'(r) is the number of pores with a pore diameter greater than r, and a is a constant.

7. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The pore shape characteristics are calculated by the following formula: Among them, G T is the pore shape characteristic, S T is the area of ​​the Tth pore region; C T is the perimeter of the Tth pore area.

8. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The multi-scale comprehensive fractal dimension of the carbonate reservoir is calculated based on the pre-constructed multi-scale comprehensive fractal dimension function and the pore key areas of different scales, including: Based on the following formula, the pores of different pore size segments in the key pore areas of different scales are calculated to obtain the multi-scale comprehensive fractal dimension of the carbonate reservoir: Where: d z is the multi-scale comprehensive fractal dimension of carbonate reservoir, B i is the ratio of the pore volume of the critical area of ​​the pores at the i-th scale to the total pore volume, i = 1, 2, ..., m, m is the number of scales, d ij is the fractal dimension of the jth aperture segment in the pore critical region of the i-th scale, j = 1, 2, ..., l, l is the number of aperture segments divided in the pore critical region of the i-th scale, b ij It is the ratio of the pore volume of the jth pore size segment in the pore critical region of the i-th scale to the pore volume of the pore critical region of this scale.

9. The method for quantitatively characterizing the pore structure of carbonate reservoirs according to claim 1, characterized in that: The method further comprises: The maximum circle search method is used to search the radius of the pore area of ​​the binary image, and the radius of the pore area of ​​the single pore structure is determined by the following formula: Where R is the radius of the pore area, δd is the actual size represented by each pixel, and N r is the number of pixels of the maximum inscribed circle radius of the pore area.

10. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The method further comprises: The area of ​​the pore region in the binary image is determined based on the following formula: Where S is the area of ​​the pore region, δd is the actual size represented by each pixel, n×n is the number of grids in the pore region, and y ij is the ordinate of the leftmost point in each row, y' ij is the ordinate of the rightmost point in each row, y ij -y' ij +1 is the number of pixels contained in each row of the grid, i = 1, 2, ..., n.

11. The method for quantitative characterization of pore structure of carbonate reservoirs according to claim 1, characterized in that: The pore structure features of the obtained digital images of carbonate reservoirs at different scales are extracted to obtain corresponding binary images, including: Digital images of carbonate reservoirs at every scale: Preprocessing the digital images respectively to obtain preprocessed digital images; An iterative image threshold segmentation algorithm is used to identify the pore structure of the preprocessed digital image to obtain a binary image.

12. The method for quantitative characterization of carbonate reservoir pore structure according to claim 11, characterized in that: The method of using an iterative image threshold segmentation algorithm to identify the pore structure of the pre-processed digital image to obtain a binary image includes: Based on the minimum gray value PF in the preprocessed digital image out | min and maximum gray value PF out | max , the initial threshold T0 is determined by the following formula; Segmenting the preprocessed digital image into a target portion and a background portion according to the initial threshold, and updating the threshold by the following formula based on the grayscale values ​​of the target portion and the background portion to obtain an updated threshold; Among them, PF out | O is the average gray value of all pixels in the target part, PF out | G is the average gray value of all pixels in the background, PF out (i, j) is the gray value of point (i, j) on the image, N(i, j) is the weight coefficient of point (i, j), T K is the threshold before updating, T K+1 is the updated threshold; Segmenting the preprocessed digital image based on the updated threshold value, repeating the process of updating the threshold value until a preset condition is met to obtain an optimal threshold value; The preprocessed digital image is segmented based on the optimal threshold to obtain a binary image.

13. The method for quantitatively characterizing the pore structure of a carbonate reservoir according to claim 11, characterized in that: The preprocessing of the digital images to obtain preprocessed digital images comprises: The characteristic components of the digital image are extracted based on the following formula: Among them, v is the feature component, L is the image grayscale, g is the pixel grayscale value, i g is the number of pixels of image gray level g, M and N are the total number of rows and columns of the image matrix corresponding to the digital image respectively; Based on the characteristic components, the digital image is subjected to image enhancement processing by the following formula to obtain a preprocessed digital image: Among them, Z(g) is the enhancement function and E is the enhancement parameter.

14. A device for quantitatively characterizing the pore structure of a carbonate reservoir, characterized in that: include: A feature extraction module is used to extract pore structure features from the acquired digital images of carbonate reservoirs at different scales, respectively, to obtain corresponding binary images; wherein the digital images at different scales include macroscopic digital images taken by high-definition optical interlacing, mesoscopic digital images taken by optical microscope, and microscopic digital images taken by scanning electron microscope; A first calculation module is used to calculate the binary image of each scale based on a pre-constructed characteristic function to obtain corresponding pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics; A search module is used for taking at least a part of the binary image of each scale as a target area, searching the binary images of different scales by a similarity matching algorithm based on the pore distribution characteristics, pore change characteristics, pore fractal characteristics and pore shape characteristics of the target area, and screening out the key pore area; The second calculation module is used to calculate the multi-scale comprehensive fractal dimension of the carbonate reservoir based on the pre-constructed multi-scale comprehensive fractal dimension function and the pore key areas of different scales.

15. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the program is executed by a processor, the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in any one of claims 1 to 13 is implemented.

16. An electronic device, characterized in that: It includes a processor, a communication interface, a memory and a communication bus, wherein the processor, the communication interface and the memory communicate with each other via the communication bus; Memory, used to store computer programs; The processor is used to implement the method for quantitatively characterizing the pore structure of a carbonate reservoir as described in any one of claims 1 to 13 when executing the program stored in the memory.

Images