A method for automatically calculating a carbonate rock surface porosity based on an imaging logging graph

Abstract

This invention provides a method for automatically calculating the face rate of carbonate rocks based on imaging logging data, comprising: (1) acquiring electrical imaging images, including static and dynamic images, based on unfilled imaging logging data; (2) filling blank bands in the electrical imaging images; (3) performing image preprocessing on the electrical imaging images to remove background noise; (4) calculating the shale grayscale threshold of the static image and obtaining the corresponding shale grayscale threshold of the dynamic image after equalization processing; (5) performing image segmentation on the electrical imaging images and performing face recognition on the electrical imaging images in conjunction with the shale grayscale threshold; (6) marking face pixels in the electrical imaging images, calculating the contribution of each pixel to the face rate based on the grayscale value, calculating the face rate, and obtaining the corrected face rate curve after smoothing processing. The method provided by this invention greatly reduces the workload of manual identification and improves the accuracy of face recognition and the reliability of parameter calculation.

Description

Technical Field

[0001] This invention belongs to the field of well logging technology and relates to a method for automatically calculating porosity, and more particularly to a method for automatically calculating porosity in carbonate rocks based on imaging well logging maps. Background Technology

[0002] Carbonate reservoirs are primarily characterized by secondary pores formed during the post-depositional diagenetic alteration process. Due to the diverse nature of these secondary alterations, the secondary pore structure of carbonate reservoirs is far more complex than that of sandstone and mudstone reservoirs, resulting in highly complex well logging response characteristics.

[0003] Imaging logging data is a crucial tool for identifying pores, providing rich information about the wellbore and its surroundings. It is characterized by its intuitiveness and high accuracy. Since pores appear as "faces" in imaging data, facial parameters reflecting reservoir porosity development can be extracted from the data. Therefore, utilizing imaging logging image information for face recognition is of great significance.

[0004] Currently, various methods for face recognition using imaging logging data have emerged in this field, such as the watershed algorithm for recognizing face contours, wavelet transform modulus maxima for face recognition, and maximum inter-class variance for calculating face thresholds. However, these methods only analyze objective image features and do not incorporate actual formation characteristics, thus presenting several problems in face recognition:

[0005] (1) After acquisition, raw imaging logging data often requires blanking band filling. Due to the large amount of data, most current blanking band filling methods require a significant amount of time to fill the image. Existing repair methods, such as deep learning-based blanking band filling methods, can achieve relatively good results by improving parameters and iterations, but they lack analysis of the characteristics of electrical imaging and understanding of geological information. In addition, there is the traditional inverse distance weighted interpolation method. Although this method is fast, even improved algorithms based on this principle cannot effectively improve the phenomenon of step-like patterns that easily appear in wide blanking bands. The current mainstream filling method is the Criminisi algorithm. Although this algorithm performs well, it may find some low-relevance matching regions during the search for matching regions, thus increasing the time cost.

[0006] (2) Dark areas in dynamic electro-imaging images do not fully represent the characteristics of pores. They are affected by low-resistivity materials, fluids, and lithology. In dynamic images, mud and mud intrusion can make the images appear dark like pores. This leads to an overestimation of the face ratio obtained by simply using classic image segmentation algorithms. Therefore, in terms of face recognition performance, relying solely on algorithms for face recognition without other auxiliary measures to eliminate non-face information and the influence of low-resistivity materials will result in an inaccurate face ratio.

[0007] (3) Since the resistivity of each point obtained from electrical imaging logging is affected by the surrounding rock skeleton or other formation components, the face recognition results include some areas that originally belong to the rock skeleton and other formation components. It is assumed that due to the resistivity measurement, there will be a transition region between the face and these formation components. Thus, the recognized face will include some non-face points, which will lead to an overestimation of the face ratio and a deterioration in the recognition effect.

[0008] (4) Currently, face recognition is generally and solely based on dynamic electro-imaging images. This is because static images represent the absolute scale of the entire well section, while dynamic images represent the local dynamic scale. In comparison, the image differences in dynamic images are more obvious, allowing for a clearer observation of porosity features. However, this leads to face recognition results ignoring the porosity contrast features of the entire well section. Therefore, the current method of using dynamic electro-imaging images alone has the drawback of not being able to accurately reflect the porosity contrast differences across the entire well section.

[0009] Therefore, how to provide a method for automatically calculating the face rate of carbonate rocks based on imaging logging maps, reducing the workload of manual identification, and improving the accuracy of face recognition and the reliability of parameter calculation has become an urgent problem to be solved by those skilled in the art. Summary of the Invention

[0010] The purpose of this invention is to provide a method for automatically calculating the face rate of carbonate rocks based on imaging logging maps. This method greatly reduces the workload of manual identification and improves the accuracy of face recognition and the reliability of parameter calculation.

[0011] To achieve this objective, the present invention employs the following technical solution:

[0012] This invention provides a method for automatically calculating the porosity of carbonate rocks based on imaging logging maps, the method comprising the following steps:

[0013] (1) Obtain electrical imaging maps, including static and dynamic maps, based on unfilled imaging logging data;

[0014] (2) Fill the blank bands in the static and dynamic images obtained in step (1);

[0015] (3) Perform image preprocessing on the static and dynamic images obtained in step (2) to remove background noise;

[0016] (4) Calculate the mud gray threshold of the static image obtained in step (3) to remove the influence of mud and mud stripes. After equalization processing, the mud gray threshold of the corresponding dynamic image is obtained.

[0017] (5) Perform image segmentation on the static and dynamic images obtained in step (3), and perform face recognition on the static and dynamic images in combination with the mud grayscale threshold obtained in step (4);

[0018] (6) Mark the face pixels in the static and dynamic images obtained in step (5), calculate the contribution of each pixel to the face ratio based on the gray value, calculate the face ratio and obtain the corrected face ratio curve after smoothing.

[0019] This invention automatically identifies and calculates faces in imaging logging data after filling blank zones. By combining dynamic and static images, it significantly improves the accuracy of face recognition, greatly reduces the workload of manual identification, and minimizes other unknown interferences caused by human factors. In particular, for carbonate reservoirs, the method has high face recognition accuracy and good reliability of parameter calculation, and can accurately identify pore-developing sections in well logging reservoir evaluation.

[0020] Preferably, the blank band filling in step (2) is performed using the Criminisi algorithm optimized based on heuristic information, specifically including the following steps:

[0021] (2.1) Select the effective area of ​​the electro-imaging image and extract the blank band to obtain the label data image. Set the blank band to 0 and the existing area to 1.

[0022] (2.2) Extract the boundary of the blank band, and calculate the boundary two-dimensional gradient, the label data image two-dimensional gradient and the structural information D in sequence, where D takes the absolute value to obtain the global image D;

[0023] (2.3) Calculate the confidence C and priority P of the global image in sequence, search for the repair block where the pixel p with the highest priority is located, and find the area that restricts the existence of the best matching block in the region to be matched by heuristic information. Calculate the color RGB difference SSD between the repair block and the matching block, and take the matching block with the smallest SSD as the best matching block.

[0024] (2.4) Fill the block to be repaired according to the best matching block, update the blank band part of the block to be repaired, and repeat steps (2.3)-(2.4) until the entire area of ​​the blank band is filled.

[0025] In this invention, the correlation between an image and its neighboring images in the Criminisi algorithm is inversely proportional to the distance. If a global search is used, sometimes a matching block that is far away from the block to be repaired and has low correlation may be obtained. That is, an excessively large search range may lead to a large difference between the final repaired image and its neighboring images.

[0026] To address this, the present invention employs a heuristic-based method for finding matching blocks, selecting the source region surrounding the point to be repaired as the matching region. This significantly shortens the search time and ensures that the repaired image is searched within its neighboring image. Specifically, this application uses the Criminisi algorithm optimized based on heuristic information for blank band filling. By adding guidance to the search for matching regions of blank bands, it avoids the time redundancy caused by the algorithm traversing all regions of the image when searching for matching regions, saving a significant amount of time for subsequent face information extraction and significantly accelerating and improving the speed and effectiveness of the filling algorithm.

[0027] Preferably, the formula for calculating the structural information quantity D in step (2.2) is:

[0028]

[0029] In the formula, D(p) represents the amount of structural information, which is used to measure the complexity of the linear structure of the surface; Indicates the direction of the isoluminance lines at pixel p; n p This indicates the direction of the boundary normal of pixel p; α is the pixel value 225.

[0030] Preferably, the formula for calculating the confidence level C in step (2.3) is:

[0031]

[0032] In the formula, C(p) represents the confidence level, which is used to measure the confidence level of the block Ψ to be repaired. p The amount of reliable information in the middle; p is Ψ p The center point; q is Ψ p Points already exist in the middle.

[0033] Preferably, the formula for calculating the priority P in step (2.3) is:

[0034] P(p)=C(p)D(P)

[0035] In the formula, P(p) represents priority.

[0036] Preferably, the search for the best matching block in step (2.3) includes searching the existing area for the matching block whose texture is most similar to the block to be repaired with the highest priority, and using it as the best matching block.

[0037] Preferably, the filling of the block to be repaired in step (2.4) includes copying the corresponding pixels in the best matching block to the unknown pixels in the block to be repaired, thereby converting the unknown pixels into known pixels.

[0038] Preferably, the matching formula between different matching blocks and the block to be repaired with the highest priority is as follows:

[0039]

[0040] In the formula, Indicates a matching block; Ψ p Indicates the block to be repaired; Represents the pixels in the matching block; This indicates the difference between the matching block and the interval of the block to be repaired.

[0041] Preferably, based on the matching criterion of minimizing the sum of the squared differences of pixel gray levels, the formula for calculating the color difference SSD between pixels in the block to be repaired and the matching block is as follows:

[0042]

[0043] In the formula, functions R(i,j), G(i,j), and B(i,j) represent the red, green, and blue primary colors of point (i,j) in the M×M region block, respectively, and M is set to 9.

[0044] Preferably, the formula for obtaining the region where the best matching block exists based on heuristic information is:

[0045]

[0046]

[0047] In the formula, m is the number of rows in the original image; n is the number of columns in the original image; l is the reciprocal of the ratio of the side length of the matching block to the shortest side length of the image; w is the width of the blank band; M is the side length of the block to be repaired centered on the pixel to be repaired p, which is set to 9; L is the side length of the square centered on the pixel to be repaired p.

[0048] Preferably, the image preprocessing in step (3) is performed using an algorithm based on image opening and closing filtering.

[0049] Preferably, the method for calculating the mud grayscale threshold of the static image in step (4) includes: selecting 0.05m as the window length for processing imaging logging data, calculating the global grayscale histogram of the static image, obtaining the mud grayscale range, and using the first peak and valley in the histogram as the mud grayscale threshold.

[0050] Before performing various calculations, this invention requires selecting a window length. The actual resolution of electrical imaging data is 0.002-0.005m, while the resolution of conventional well logging data is 0.1m. Therefore, in order to maintain high data resolution and small errors, this invention specifically selects a window length of 0.05m.

[0051] Furthermore, in static maps, threshold segmentation during the calculation of porosity in carbonate reservoirs can incorrectly classify argillaceous material or argillaceous bands as faces. Moreover, the grayscale values ​​of argillaceous bands are even lower in face locations, and their grayscale value range differs from that of faces. This invention addresses this by establishing a global grayscale histogram of the static map. After analyzing the histogram, the grayscale threshold of argillaceous material in carbonate reservoirs can be clearly identified, thereby reducing the influence of argillaceous bands on porosity calculations.

[0052] Preferably, the method for calculating the mud gray threshold of the dynamic image in step (4) includes: calculating the mud gray threshold of the static image within the length range of each sliding window, and obtaining the mud gray threshold of the corresponding dynamic image after equalization processing.

[0053] Preferably, the image segmentation in step (5) is performed using the maximum inter-class variance algorithm, and the specific calculation formula is as follows:

[0054] g = w0 × w1 × (u0 - u1) 2

[0055] In the formula, g represents the variance; w0 represents the proportion of face pixels to the total number of pixels; w1 represents the proportion of background pixels to the total number of pixels; u0 represents the average grayscale value of face pixels; and u1 represents the average grayscale value of background pixels.

[0056] Preferably, the face recognition in step (5) includes: calculating the average gray value of the static image within the length of each sliding window, using gray level 128 as the minimum limit and gray level 224 as the maximum limit, defining the gray value range for face recognition using the static image, and using the dynamic image for face recognition in the remaining gray value range.

[0057] In this invention, the Otsu's algorithm is better for image segmentation when the gray values ​​of an image are not concentrated on one side, or when there is a significant difference between the foreground and background. Equalization, on the other hand, amplifies subtle differences in the image, which is why dynamic images are commonly used for face recognition. However, when the resistivity of formation porosity and formation components is not significantly different—that is, when the high-resistivity formation components in the static image do not obscure the pore resistivity characteristics—using dynamic images for face recognition can cause faces in the static image to transform into formation components. Therefore, this invention chooses to use static images to calculate the face rate when the average gray value is below 224. Furthermore, when the average gray value of the static image is above 128 within the sliding window, the low-resistivity features of faces do not obscure the resistivity characteristics of formation components, while equalization would disrupt this balance. Therefore, using static images to recognize faces in this range is more accurate.

[0058] This invention uses the maximum inter-class variance algorithm for image segmentation. First, the segmentation threshold of faces in the image is determined. Then, the segmentation threshold is compared with the gray value of the pixel. The point with the smaller value is the face. Finally, based on the characteristic that different gray values ​​of each pixel correspond to different resistivity, the gray values ​​of face pixels within the sliding window are normalized.

[0059] Preferably, the formula for calculating the contribution of each pixel to the face ratio in step (6) is:

[0060]

[0061] In the formula, K represents the contribution of the corresponding gray level to the face rate; T represents the maximum inter-class variance threshold; min g The minimum grayscale value of the region is set as the grayscale threshold for muddy texture; G represents the grayscale value of the current pixel, ranging from T to min. g between.

[0062] Preferably, the formula for calculating the face rate in step (6) is:

[0063]

[0064] In the formula, SPOR represents the face ratio after correction; sum(K) represents the total contribution of faces in the region, and the contribution of a single point is at most 1; sum(i) represents the total number of pixels in the region.

[0065] In this invention, the gray value of each pixel corresponds to the resistivity at the current location. According to the principle of measuring resistivity, the resistivity is affected by the resistivity of the surrounding area. Therefore, a specific correction formula is used to correct the identified face pixels, thereby reducing the above-mentioned influence.

[0066] Preferably, the smoothing process in step (6) includes a five-point triple smoothing process.

[0067] Compared with the prior art, the present invention has the following beneficial effects:

[0068] This invention automatically identifies and calculates faces in imaging logging data after filling blank zones. By combining dynamic and static images, it significantly improves the accuracy of face recognition, greatly reduces the workload of manual identification, and minimizes other unknown interferences caused by human factors. In particular, for carbonate reservoirs, the method has high face recognition accuracy and good reliability of parameter calculation, and can accurately identify pore-developing sections in well logging reservoir evaluation. Attached Figure Description

[0069] Figure 1 This is a flowchart of the method for automatically calculating the porosity of carbonate rocks based on imaging logging provided by the present invention;

[0070] Figure 2 This is a schematic diagram illustrating the selection of the best matching region when obtaining the best matching block based on heuristic information in Example 1;

[0071] Figure 3 This is the global gray-scale density histogram of well A in Example 1;

[0072] Figure 4 This is the Criminisi electro-imaging fill effect diagram optimized based on the heuristic information method in Example 1;

[0073] Figure 5 The results are the logging curves and porosity curves of well A from 4028 to 4108m in Example 1. Detailed Implementation

[0074] The technical solution of the present invention will be further illustrated below through specific embodiments. Those skilled in the art should understand that the embodiments described are merely illustrative of the present invention and should not be construed as limiting the invention in any way.

[0075] Example 1

[0076] This embodiment provides a method for automatically calculating the porosity of carbonate reservoirs based on imaging logging maps. Taking well A in a carbonate reservoir of an oilfield as an example, ... Figure 1 As shown, the method includes the following steps:

[0077] (1) Obtain electrical imaging maps, including static and dynamic maps, based on unfilled imaging logging data;

[0078] (2) The Criminisi algorithm, optimized based on heuristic information, is used to fill blank bands in the static and dynamic graphs obtained in step (1), specifically including the following steps:

[0079] (2.1) Select the effective area of ​​the electro-imaging image and extract the blank band to obtain the label data image. Set the blank band to 0 and the existing area to 1.

[0080] (2.2) Extract the boundary of the blank band, and calculate the boundary two-dimensional gradient, the label data image two-dimensional gradient and the structural information D in sequence, where D takes the absolute value to obtain the global image D;

[0081] The formula for calculating the structural information content D is as follows:

[0082]

[0083] In the formula, D(p) represents the amount of structural information, which is used to measure the complexity of the linear structure of the surface; Indicates the direction of the isoluminance lines at pixel p; n p This indicates the direction of the boundary normal of pixel p; α is the pixel value 225;

[0084] (2.3) Calculate the confidence C and priority P of the global image in sequence, search for the repair block where the pixel p with the highest priority is located, and find the area that restricts the existence of the best matching block in the region to be matched by heuristic information. Calculate the color RGB difference SSD between the repair block and the matching block, and take the matching block with the smallest SSD as the best matching block.

[0085] The confidence level C is calculated using the following formula:

[0086]

[0087] In the formula, C(p) represents the confidence level, which is used to measure the confidence level of the block Ψ to be repaired. p The amount of reliable information in the middle; p is Ψ p The center point; q is Ψ p Points already exist in the middle;

[0088] The formula for calculating the priority P is:

[0089] P(p)=C(p)D(P)

[0090] In the formula, P(p) represents priority;

[0091] The search for the optimal matching block involves searching the existing region for the matching block whose texture is most similar to the highest priority repair block, and selecting this as the optimal matching block. The matching formula between different matching blocks and the highest priority repair block is as follows:

[0092]

[0093] In the formula, Indicates a matching block; Ψ p Indicates the block to be repaired; Represents the pixels in the matching block; This indicates the difference between the matching block and the block to be repaired.

[0094] Based on the matching criterion of minimizing the sum of the squared differences of pixel gray levels, the formula for calculating the color difference SSD between pixels in the block to be repaired and the matching block is as follows:

[0095]

[0096] In the formula, the functions R(i,j), G(i,j), and B(i,j) represent the red, green, and blue primary colors of the point (i,j) in the M×M region block, respectively, and M is set to 9;

[0097] The formula for obtaining the region where the best matching block exists based on heuristic information is:

[0098]

[0099]

[0100] In the formula, m is the number of rows in the original image; n is the number of columns in the original image; l is the reciprocal of the ratio of the matching block's side length to the shortest side length of the image; w is the width of the blank band; M is the side length of the block to be repaired centered on the pixel to be repaired p, set to 9; L is the side length of the square centered on the pixel to be repaired p (see...). Figure 2 );

[0101] (2.4) Copy the corresponding pixels in the best-matching block to the unknown pixels in the block to be repaired, thereby converting the unknown pixels into known pixels. Update the blank band portion of the block to be repaired. Repeat steps (2.3)-(2.4) until the entire area of ​​the blank band is filled (see...). Figure 4 );

[0102] (3) The static and dynamic images obtained in step (2) are preprocessed using an image opening and closing filtering algorithm to remove background noise.

[0103] (4) Select 0.05m as the window length for processing imaging logging data, calculate the global grayscale histogram of the static image obtained in step (3), and obtain the shale grayscale range (see Figure 3 The first peak in the histogram is used as the mud gray threshold to remove the influence of mud and mud stripes. The mud gray threshold of the static image is calculated within the length of each sliding window. After equalization, the mud gray threshold of the corresponding dynamic image is obtained.

[0104] (5) The maximum inter-class variance algorithm is used to segment the static and dynamic images obtained in step (3) respectively. The average gray value of the static image is calculated within the length of each sliding window. Gray level 128 is used as the minimum limit and gray level 224 is used as the maximum limit. The gray value range for face recognition using the static image is specified, and the remaining gray value range is used for face recognition using the dynamic image.

[0105] The specific calculation formula for the Otsu's inter-class variance algorithm is as follows:

[0106] g = w0 × w1 × (u0 - u1) 2

[0107] In the formula, g represents the variance; w0 represents the proportion of face pixels to the total number of pixels; w1 represents the proportion of background pixels to the total number of pixels; u0 represents the average grayscale value of face pixels; u1 represents the average grayscale value of background pixels.

[0108] (6) Mark the face pixels in the static and dynamic images obtained in step (5), calculate the contribution of each pixel to the face ratio based on the gray value, calculate the face ratio, and obtain the corrected face ratio curve after five-point cubic smoothing (see Figure 5 );

[0109] The formula for calculating the contribution of each pixel to the face ratio is as follows:

[0110]

[0111] In the formula, K represents the contribution of the corresponding gray level to the face rate; T represents the maximum inter-class variance threshold; min g The minimum grayscale value of the region is set as the grayscale threshold for muddy texture; G represents the grayscale value of the current pixel, ranging from T to min. g between;

[0112] The formula for calculating the face rate is:

[0113]

[0114] In the formula, SPOR represents the face ratio after correction; sum(K) represents the total contribution of faces in the region, and the contribution of a single point is at most 1; sum(i) represents the total number of pixels in the region.

[0115] Therefore, this invention automatically identifies and calculates faces in imaging data after filling blank zones based on imaging logging data. By combining dynamic and static images, it significantly improves the accuracy of face recognition, greatly reduces the workload of manual identification, and reduces other unknown interferences caused by human factors. In particular, for carbonate reservoirs, the method has high face recognition accuracy and good reliability of parameter calculation, and can accurately identify pore-developing intervals in well logging reservoir evaluation.

[0116] The applicant declares that the above description is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Those skilled in the art should understand that any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention fall within the protection and disclosure scope of the present invention.

Claims

1. A method for automatic calculation of carbonate rock surface porosity based on imaging logging maps, characterized in that, The method comprises the following steps: (1) obtaining an electrical imaging diagram according to un-filled imaging logging data, including a static diagram and a dynamic diagram; (2) filling the static diagram and the dynamic diagram obtained in step (1) respectively with blank bands; (3) pre-processing the static diagram and the dynamic diagram obtained in step (2) respectively to remove image background noise; (4) calculating a shale gray threshold of the static diagram obtained in step (3) to remove shale and shale bands, and obtaining a shale gray threshold of the dynamic diagram after equalization processing of the shale gray threshold; (5) performing image segmentation on the static diagram and the dynamic diagram obtained in step (3) respectively, and performing face hole recognition on the static diagram and the dynamic diagram respectively in combination with the shale gray threshold obtained in step (4); the face hole recognition comprises: calculating an average gray value of the static diagram in each sliding window length range, taking gray level 128 as the minimum limit and gray level 224 as the maximum limit, defining a gray value range for face hole recognition using the static diagram, and using the dynamic diagram for face hole recognition in the remaining gray value range; (6) marking the static diagram and the dynamic diagram obtained in step (5) respectively with face hole pixels, calculating the contribution of each pixel point to the face hole rate based on the gray value, calculating the face hole rate, and performing smoothing processing to obtain a corrected face hole rate curve.

2. The method of claim 1, wherein, The blank band filling in step (2) is performed by using a Criminisi algorithm optimized based on heuristic information, and specifically comprises the following steps: (2.1) selecting an effective area of the electrical imaging diagram and extracting a blank band to obtain a label data image, setting the blank band to 0 and the existing area to 1; (2.2) extracting the boundary of the blank band, and sequentially calculating the two-dimensional gradient of the boundary, the two-dimensional gradient of the label data image and the structure information D, wherein D takes an absolute value, to obtain a global image D; (2.3) sequentially calculating the confidence C and the priority P of the global image, searching for a repair block in which a pixel point p with the maximum priority exists, searching for a region in which a best matching block exists in a to-be-matched region by a heuristic information method, calculating the color RGB difference SSD between the to-be-repaired block and the matching block, and taking the matching block with the minimum SSD as the best matching block; (2.4) filling the to-be-repaired block according to the best matching block, updating the blank band part of the to-be-repaired block, and repeating steps (2.3)-(2.4) until all areas of the blank band are filled.

3. The method of claim 2, wherein, The calculation formula of the structure information D in step (2.2) is: In the formula, represents the structural information amount, which is used to measure the complexity of the linear structure of the surface; represents the isophote line direction of the pixel point p; represents the boundary normal direction of the pixel point p; is the pixel value 225.

4. The method of claim 2, wherein, The calculation formula of the confidence C in step (2.3) is: In the formula, Indicates confidence level, used to measure the block to be repaired. The amount of reliable information in the middle; p is The center point; q is Points already exist in the middle.

5. The method of claim 2, wherein, The calculation formula of the priority P in step (2.3) is: In the formulae, Indicates priority.

6. The method of claim 2, wherein, The search for the best matching block in step (2.3) comprises searching for a matching block with the most similar texture to the to-be-repaired block with the maximum priority in the existing region as the best matching block.

7. The method of claim 2, wherein, The filling of the to-be-repaired block in step (2.4) comprises copying the pixel points in the best matching block to the unknown pixel points of the to-be-repaired block, so as to convert the unknown pixel points into known pixel points.

8. The method of claim 6, wherein, The matching formula between different matching blocks and the to-be-repaired block with the maximum priority is: In the formula, represents a matching block; represents a block to be repaired; represents a pixel point in the matching block; represents a difference between the matching block and the block to be repaired.

9. The method of claim 8, wherein, Based on the matching criterion of the minimum sum of pixel gray square differences, the color difference SSD between the to-be-repaired block and the pixels in the matching block is calculated as follows: In the formula, the functions , , respectively represent the red, green, and blue primary colors of the point in the M x M region block, and M is set to 9.

10. The method of claim 2, wherein, The formula for obtaining the best matching block existence area based on heuristic information is: In the formula, m is the number of rows of the original image; n is the number of columns of the original image; l is the reciprocal of the ratio of the matching block side length to the shortest side length of the image; w is the width of the blank band; M is the side length of the block to be repaired with the pixel point p as the center, and is set to 9; and L is the square side length with the pixel point p as the center.

11. The method of claim 1, wherein, The image preprocessing in step (3) is performed by using an algorithm based on image opening and closing filtering.

12. The method of claim 1, wherein, The mud gray scale threshold calculation method of the static image in step (4) comprises the following steps: 0.05 m is selected as the window length for processing the imaging logging data; the global gray scale histogram of the static image is calculated; the mud gray scale range is obtained; and the first peak and valley in the histogram are taken as the mud gray scale threshold.

13. The method of claim 1, wherein, The mud gray scale threshold calculation method of the dynamic image in step (4) comprises the following steps: the mud gray scale threshold of the static image is calculated in each sliding window length range; and the mud gray scale threshold of the corresponding dynamic image is obtained after equalization processing.

14. The method of claim 1, wherein, The image segmentation in step (5) is performed by using the maximum inter-class variance algorithm, and the specific calculation formula is: In the formula, g represents a variance; represents a proportion of the number of face pixels to the total number of pixels; represents a proportion of the number of background pixels to the total number of pixels; represents an average value of the gray scale of face pixels; represents an average value of the gray scale of background pixels.

15. The method of claim 1, wherein, The contribution degree calculation formula of each pixel point to the face ratio in step (6) is: In the formula, K represents the contribution degree corresponding to the face rate of the gray level; T represents the maximum inter-class variance threshold value; The minimum value of the area gray level is set as the argillaceous gray level threshold value; G represents the gray value of the current pixel point, and the range is between T .

16. The method of claim 1, wherein, The calculation formula of the face ratio in step (6) is: In the formula, SPOR represents the corrected aperture ratio; represents the total sum of the regional pixel point contribution degrees, and the single-point contribution degree is highest as 1; represents the total number of regional pixel points.

17. The method of claim 1, wherein, The smoothing processing in step (6) comprises five-point cubic smoothing processing.

Images