Automatic Detection of Fractures Based on Optimal Path Search in Well Logging Images

Reservoir fractures are essential locations to gather oil and gas. Recently, imaging logging technology has become a mainstream method for obtaining stratigraphic information. This paper proposed a combined optimal path search strategy to effectively identify and extract the fracture information in well logging images. Specifically, the threshold segmentation of logging images is used to obtain the binary image. In addition, the identification of connected fractures in the logging image is transformed into the optimal path search, and the identification and extraction of reservoir fractures are realized by constructing the optimal path between the two ends of fractures. Finally, an improved ant colony algorithm is applied to filter irrelevant information and extract fractures automatically by recording all the ants’ exploration trajectories in the ant colony. Compared with previous approaches, the proposed method can eliminate irrelevant background features and merely reserve pixels corresponding to fractures. Simultaneously, relative to the conventional strategy, the time consumption is reduced by more than 98%. The findings of this study can help for better extracting fractures automatically and reducing manual workload.


Introduction
Formation fractures are discontinuous profiles widely distributed in different lithologies and gradually formed through diagenesis or tectonic deformation. Naturally fractured reservoirs store over 50% of the known petroleum and gas reserves worldwide [1]. Accordingly, the detection of fractures has attracted plenty of researchers because they are vital indicators of the admissibility of petroleum and gas reservoirs in tight formations [2] and affect reservoir properties, enrichment, and gas reservoir development [3][4][5][6]. For instance, microfractures in shale are considered important fluid transport networks as well as oil and gas migration pathways and are key factors in forming oil and gas reservoirs [7]. So far, there are many technologies for discovering potential fractures in petroleum reservoirs, such as subsurface electromagnetic technology [8], tiltmeter [9], downhole microseismic fracture monitoring [10], and radiotracer diagnosis [11]. Compared with these methods, the imaging logging analysis technology [12] originated in 1986 has become the mainstream method of oil and gas reservoir exploration, which wins the favor of considerable engineers for its intuitive and accurate display of wellbore and stratigraphic structure. Ultrasonic imaging logging [13,14], one of the representative imaging logging technologies, can characterize the geometric features of fractures and distinguish various geological features more clearly. Figure 1 is a schematic diagram of ultrasonic imaging logging. Each pixel in the image corresponds to one arrival time of the ultrasonic signal reflected from the borehole wall. A bright area indicates that the arrival time of echo measured in this area is short. Conversely, a dark area indicates long arrival time or no reflection exists. Therefore, this paper's task is to design a filter which extracts pixels corresponding to fractures as shown in Figure 2.
Due to the complexity of reservoir, it is a timeconsuming process to accurately identify fractures manually. Accordingly, automatic identification of fractures in logging images is of great significance [15], which prompts many researchers to explore proper approaches. Initially, the scholars established mathematical models of the fractures to characterize the characteristics. Changchun et al. regarded the shape of fracture as a sine curve and utilized twodimensional Hough Transform (2D-HT) to detect fractures with fixed patterns in the image [16]. Taiebi et al. used a multiscale technique based on directional filtering and Hough Transform (HT) to extract fractures in the logging images further [17]. Directional filtering can enhance the contrast between the target fracture and the imaging logging background, which is positive to conduct HT for fracture segmentation. Liu et al. used ant colony algorithm to ascertain the edges of fractures. Then, the sinusoidal fracture was found by HT [18]. Jianping et al. improved the HT according to the characteristics of fractures in imaging well logging. For such a sine curve with fixed periods, the initial detection step is divided into two steps: the voting mechanism was used to determine the baseline position of sinusoids; and then the 2D-HT determined the amplitude and initial phase parameters of the sinusoidal curve [19]. Based on the baseline position determined by voting mechanism, Yingming et al. leveraged genetic algorithms to perform nonlinear fitting of sineshaped scattered points in the relevant area to realize pixel extraction of sinusoidal regions such as fracture, bedding, and layer boundaries of imaging logging images [20]. Since squeeze friction in the actual drilling process will cause the borehole wall damaged, the fracture shape is not necessarily in standard sinusoidal shape. Hence, the method based on mathematical models is not effective in fitting and identifying nonstandard sinusoidal fractures.
Recently, with the rapid development of computer vision, researchers have tried to use image segmentation to separate pixels corresponding to fractures from the back-ground of ultrasonic logging images and then realize the extraction of fractures. Wang utilized the valley edge algorithm to determine local dark area to accomplish rock fracture detection without determining a threshold [21]. Xu et al. combined the classification into a cascade system using the K -nearest neighborhood (KNN) classifier, which can classify rock structures and extract the expected features in the ultrasonic well logging images [22]. This method divides the image into superpixel blocks with uniform size and moderate compactness, which provides the basis for subsequent fracture recognition. Sorncharean et al. used grid unit analysis chain and fracture unit verification to eliminate false detection of shadow boundaries on the image, so as to realize the detection of fractures on road images with uneven illumination and strong texture [23]. Zhang used the clustering-minimum spanning tree method to extract the crack area by edge extraction and threshold segmentation and then used the minimum spanning tree method to detect the crack in asphalt pavement [24]. However, in actual images, simply using low-level content information such as color, brightness, and texture of pixels is not enough to generate a good segmentation effect, and it is easy to produce wrong segmentation results. To sum up, all mentioned approaches have the following issues: they fail to identify irregular fractures because of low extraction accuracy, and the identification and extraction efficiency is not efficient. Recognition results often remain some noncrack information.
To address above issues, we have analysed the fractures' semantic information in logging images and utilized the pattern of fracture regions to determine fracture region. Different from the above methods, connectivity of fracture regions in logging images is considered a key standard where pixels corresponding to reservoir fractures connected two ends of an image [25]. On the other hand, heuristic algorithms have played an increasingly vital role in similar optimization problems. For example, Deng et al. [26] proposed multiple strategies for global optimization problems. And Deng et al. [27] further used an improved quantum-inspired differential evolution algorithm to avoid premature convergence and improve the global search ability. Therefore, based on the above inference, this paper proposes a method that combined optimal path search strategy (COPS) to identify and extract fractures in ultrasonic logging images, and the main contributions of this research are as follows: (1) Construct a mathematical model describing this task. And convert fracture recognition to path search (2) Preprocess logging images, and build a path search space by OTSU (OTSU) threshold segmentation algorithm (3) Use the ant colony algorithm (ACA) to search potential paths connecting two ends of fractures, which extracts the entire fracture region (4) Improve ant colony algorithm by designing a parallel searching strategy to accelerate path search dramatically The paper is organized by following steps: first, the proposed methods are presented by detailed computational

Mathematical Model.
In the stratum, different structural planes show certain differences on account of different morphological and physical characteristics. In fractured reservoirs, the expansion zone generated by structural fractures is where a large amount of oil and gas transfer and accumulate. Consequently, the method designed in this paper is mainly aimed at the identification and extraction of such fractures in ultrasonic logging images. This type of fracture is a fracture zone surrounding the wellbore and formed along the borehole wall. The differences mainly manifest in the color depth, shape difference, width change of the band-shaped sinusoid, and the irregularity of the rock spot and texture characteristics around the curve [28]. The fracture presents a three-dimensional ellipse shape on the wellbore wall, and if the cylindrical well wall is cut along its axis, the three-dimensional elliptical fracture will be a sinusoidal curve after unfolding on the two-dimensional image [29][30][31][32], as shown in Figure 3. The function of the fracture can be described by Equation (1) [33]: in the formula, y 0 is the baseline position of the sinusoidal crack curve, A is the amplitude of the curve, β is the initial phase, and ω is the angular velocity. For reservoir fractures distributed along the circumference of the well wall, the shape of fractures can be approximately considered to be sinusoidal. However, in practice, due to the mechanical vibration in drilling process and the influence of geological activities, the fractures spreading in the logging images may not show a standard sinusoidal shape, and sometimes, it is even far from a sinusoidal form. Hence, in the process of image processing, if a sinusoidal shape is regarded as a prior condition of the fracture morphology, the algorithm will have a high probability to generate false recognition of fractures or fail to fit the shape of fractures well.
According to the fractures distributed around the borehole wall, a feature that is obviously different from other interference information is that the fractures connect the left and right sides in the logging image. Therefore, we suppose the ultrasonic logging image is G, containing an ordered pair ðV, EÞ, where V is a vertex set composed of all pixels and E is a set of disordered point pairs composed of elements in V, denoted as edge set, whose elements are edges, and a same point pair can appear multiple times in E, denoted as G = ð V, EÞ. As a result, for the identification and extraction of fractures in the logging image, the problem is transformed into finding a subset V subset of vertices in graph G that can connect the regions on both sides of the graph; this V subset can correspond to the area where the fracture is located.

Image Preprocessing.
According to the principle of ultrasonic logging images, fractures are darker in the image and appear dark gray sinusoidal shapes. Factors like the unnatural etching of the strata, extrusion deformation caused by pressure, and the damage to the borehole wall by the drill bit during the drilling process will leave a geometric shape similar to the fracture or a similar gray value on the ultrasonic logging image. Therefore, before performing fracture identification and extraction, it is necessary to perform image segmentation to filter unnecessary interference information and enhance the fracture characteristics at the same time.
As mentioned in introduction, based on the OTSU algorithm, this paper performs preliminary threshold segmentation on the original ultrasonic logging image to extract effective fracture area information. The OTSU algorithm, also known as the maximum between-class variance method, is an algorithm to determine threshold. This threshold is used to perform fixed threshold binarization of the image to maximize the variance between classes. According to the gray characteristics of the image, it is divided into two parts: the background and the foreground. The segmentation with the largest variance between the classes means the smallest probability of misclassification. Suppose the size of a single ultrasonic logging image is X × Y, grayscale range R = ½0,255, R ∈ N * , and the probability of the gray level i is Define the fracture area C t , the background area C b , and a gray-scale threshold T; then assign all the pixels in the logging image to C t and C b according to this threshold, the probabilities of C t and C b are θ t and θ b , respectively; so the mean gray value of the fracture area μ t and the mean value of the background area μ b are According to the literature [34], this paper adopts an easily calculated interclass variance evaluation function σ 2 c as the standard of the threshold segmentation, which is defined as follows: among them, μ T = ∑ 255 i=0 ip i and μðtÞ = ∑ T i=0 ip i . Finally, find the best threshold to transformed into As shown in Figure 4, the ultrasonic logging image is transformed into a binary image after threshold segmentation. Compared with the maximum entropy threshold segmentation algorithm [35], the OTSU algorithm has a better extraction effect in response to ultrasonic imaging logging fracture identification and minimizes the interference of background noise.

Optimal Path Search Based on OTSU Threshold
Segmentation Algorithm and Ant Colony Algorithm. As mentioned above, the search of fracture area in logging images can be transformed into seeking the connected area on the left and right sides of the connection image G = ðV, EÞ. Therefore, the search space is established according to the binary image segmented by a threshold value. In a binary image, the gray value of the potential fracture area is 0, which is a passable region; the gray value of the background area is 1, which is an obstacle region. In this paper, the serial number method [36] is used to establish the grid set corresponding to the binary ultrasonic logging image, and the effect is shown in Figure 5.
After rasterizing the image and establishing the search space, it is necessary to search the area connecting the left and right sides to determine the fracture area. For path search using the grid method, the search efficiency is positively related to the size of the space, so an efficient and stable search algorithm is needed. Compared with other algorithms, the ACA has strong robustness and adaptability and has achieved good results in solving path planning [37]. Accordingly, the ACA is used to search for the fracture area.
The ACA is a swarm intelligence algorithm that simulates the foraging behavior of ants in nature. By releasing pheromones on the foraging path, the ant colony will tend to walk on the path with higher pheromone concentration and releasing pheromones simultaneously. Pheromone can attract more ants and form a feedback loop. Eventually, the entire ant colony will find the most suitable path. The search steps of fractures area based on ACA are as follows: (1) Suppose the number of ants in the ant colony is m, the number of pixels in the potential fracture area is n, and the distance between pixels V i and V j in the fracture area is d ij ði, j = 1, 2,⋯,nÞ, the pheromone concentration between any target pixel points V i and V j in the same fracture region at time t is τ ij ðtÞ, and the initial pheromone concentration is τ ij ð0Þ = τ 0 . The pixel points V io and V jo are randomly selected on the left and right sides of the fracture area as the starting and ending points of the path search (2) Let P k ij ðtÞ be the probability of the ant k in the ant colony moving from the pixel point i to the pixel point j; the calculation formula is among them, η ij ðtÞ is the heuristic function, allow k k = ð1, 2,⋯,kÞ is the set of pixels to be searched by ant k, a is the pheromone importance factor, and b is the importance factor of the heuristic function. Through the roulette wheel selection method, according to the transition probability of the remaining pixels, the ant k will go to the next pixel (3) When an ant has completed a traversal and there is no way to go, the pheromone concentration τ ij ðtÞ on the path between the traversed pixels will be updated, as shown in Δτ k ij = Q L k , the kth ant goes from pixel i to pixel j, 0, others, where Δτ k ij represents the pheromone released by the ant k on the connection path between the pixel V i and V j , ρ denotes pheromone volatile factor, Δτ ij represents the sum of the pheromone concentration accumulated by all ants on the connecting path between the pixel V i and V j , Q is a constant, and L k is the length of the path that the ant k passes 4 Journal of Sensors (4) When the maximum number of iterations is reached, the path search ends and the search path is printed out Through the optimal path search (OPS) algorithm based on the OTSU algorithm and ACA, the connected areas at both ends of the connecting image G can be traversed by different ants in the ant colony. When all ants' footprints are recorded, the fracture areas in the ultrasonic logging image can be identified and extracted. Simultaneously, since there are no ants traversing the nonfracture areas, these areas can be filtered out, and only information about the fracture areas can be retained.

Combined Optimal Path Search Strategy.
Through strategies mentioned above, the fracture area in the ultrasonic logging image can be correctly identified and extracted. However, due to the large resolution of a single logging image, the algorithm could converge slowly and easily fall into a locally optimal solution if the optimal path search is performed directly. Therefore, this paper proposes the COPS given the morphological characteristics of fracture regions and ultrasonic logging images.
The current common ACA application is designed for the shortest path search. The pheromone concentration Δ τ k ij is calculated according to the total length of the ant's path, which means that the shorter path ant colonies select to travel, the higher the pheromone concentration will be left. Then, the overall ant colonies tend to the shortest path. The object of this paper is the overall search of the fracture area; that is, we hope the ant colony can traverse the entire fracture area of the logging image as much as possible. To   (11) to The kth ant goes from pixel i to pixel j, 0, others, where d ij represents the distance between the pixel points V i and V j passed by the ant. This formula uses local information of the path that ants pass through to calculate the released pheromone concentration. Therefore, it can increase the probability of the ant going to the different pixels in fracture areas and as many pixels as possible can be traversed. It can be seen from Equation (1) that the approximate position of the sinusoidal fracture curve in the logging image can be determined by the baseline y 0 . Thus, a specific path search range can be determined by seeking the baseline position, and the scale of the path search can be reduced. For any point s 1 on the sine curve, another point s 2 on the curve can be determined, which meets Xðs 1 Þ − Xðs 2 Þ = T/2, where T is the period of the sine curve, and the midpoint S 0 = ðx 0, y 0 Þ between s 1 and s 2 must fall on the baseline y 0 . By searching for such point pairs and using the voting accumulation mechanism, the ordinate information of all midpoints is counted to determine the baseline position. Baseline positioning is performed on the binarized ultrasonic logging image in Figure 5, and the area where the baseline is located is cropped to get the target search area G s , and the result shown in Figure 6 is obtained.
In order to further improve the efficiency of ACA and accelerate the convergence speed, the ACA in this paper is improved by adopting a strategy searching single area path parallelly. For the target search region G s , the fracture area is connected through left and right. Therefore, G s is divided into l different subregions G si ði = 1, 2,⋯,lÞ along the longitudinal direction, and the fracture area within each subregion G si is still connected from left to right. The improved ACA is utilized to seek the path in each subarea G si simultaneously, and then all the subfracture areas searched in the subarea are spliced to complete the identification and extraction of the fractures in the logging image. The calculation process is shown in Figure 7. Using the strategy of searching a single region path parallelly, the region that the algorithm needs to traverse simultaneously is only 2Δy/Yl of the original region, which can accelerate the search speed and make the algorithm converge rapidly. Since the search range is reduced and each area is searched by an independent ant colony, it can also effectively improve the situation that the ant colony falls into optimal local results.

Experiment Platform.
For verifying the performance of the optimal path search strategy proposed in this paper, we Split image

Splitce path
Search each region parallel  Journal of Sensors tested the algorithm in real logging images. Figure 8 depicts the overall structure of the ultrasonic image well logging instrument. The well logging tool consists of rotary ultrasonic transducer driven by motor, main control board receiving ultrasonic signal and conducting signal processing, sensor drive board driving ultrasonic transducer, and cable for communication with a host computer. 250 points are collected evenly by a transducer rotating one circle, and a logging image is drawn according to the arrival time and maximum amplitude of the echo. The test object is the ultrasonic logging image G Z from 4423 meters to 4559 meters downhole of the Zhanjiang production oil well. Images are provided by the Oilfield Technology Research Institute of China Oilfield Services Co., Ltd. (COSL). Logging image G Z is cut into 72 subimages with a resolution of 352 × 352. The image processing algorithm runs on a hardware computing platform based on an i9-9900k processor with 64GB memory, and the programming language is MATLAB R2020a. Figure 9 shows the operation of the circumferential ultrasonic imaging logging tool in the Zhanjiang production well.

Test Results.
In order to further verify the recognition accuracy of the proposed algorithm for different forms of fractures and different background noise interference, five typical fracture regions in G Z were selected for testing and compared with the results of the common threshold segmentation algorithm and the artificially marked fracture regions as ground truth. The results are shown in Figure 10. As shown in Figure 10, whether it is a sinusoidal fracture, a horizontal fracture, or a fracture with obvious background noise, the combined optimal search strategy algorithm can eliminate the residual background noise and more accurately preserve fracture information while compared with the traditional direct threshold algorithm.
Confusion matrix, which is shown in Figure 11, is an index for evaluating classification performance. According to pixels' classification in ground truth and segmentation results of two algorithms, we calculate average "true posi-tive" (TP), "false positive" (FP), "true negative" (TN), and "false negative" (FN) in the above images to draw the confusion matrices. Consequently, precision and recall which determine performances of classification can be calculated as  Table 1 shows the corresponding F1 score, mean Intersection over Union (IoU), precision, and recall generated by two mentioned methods. Obviously, toward fracture segmentation in well logging images, the IoU can be increased from 20.72 to 43.21, and the F1 score can be increased by 0.26. Precision and recall are commonly contradictory when assessing the same segmentation result, and higher recall from Maximum Entropy Threshold represents that it classifies more pixels as fractures, but its precision is reduced. In summary, all results demonstrate that the fracture segmentation generated by the proposed method is closer to the fracture information marked manually.  Figure 12 shows the time-consuming comparison of fracture identification and extraction of the above five typical logging images using conventional ACA and the algorithm proposed in this paper to perform path search on the computing platform of this paper. Both colony sizes are 20, and the number of iterations is 100 while pheromone importance factor a is 1, importance factor of the heuristic function b is 7, and pheromone volatile factor ρ is 0.3. It can be seen from Figure 12 that the algorithm proposed in this paper has a significant efficiency improvement.
For exploring the factor affecting segmentation further, we conducted comparative tests by adjusting iterations as the number of iterations determines whether the ant colony could traverse all possible paths. We use parameters men-tioned above with different numbers of iteration to carry out experiments. The average time consumed and IoUs are shown in Table 2. By observing the outcomes, while 100 is determined as the iteration number, the proposed method can acquire a balance between time consumed and visual results. Accordingly, we performed 100 iterations when using our method.

Discussion
In this paper, we propose an effective approach coping with fracture recognition and extraction. Relative to conventional methods, i.e., threshold segmentation and Hough Transform, not only can the method identify fracture regions accurately, but it can also extract pixels corresponding to fracture regions to show intact shape of fracture clearly.
Compared with the conventional ant colony algorithm, one advantage is apparent acceleration in the rate of path search by parallel search strategy as search space is reduced by 4 times. By analysing the principle of our method, another reason of acquiring remarkable results based on our method is that ant colonies tend to search paths with high concentration of pheromone. This indicates that when an ant reaches the end of the path, the rest of ant colonies are likely to follow it to avoid invalid searches, which have a positive impact on filtering irrelevant information. However, due to the same reason that ants tend to search paths with high concentration of pheromone, all fractures may not be extracted completely when more than one fracture exists in logging images.
On the other hand, based on prior knowledge of fracture connectivity, the proposed method cannot extract the complete fracture regions if actual fractures in logging images are discontinuous and disconnected, which is intractable in our completed study. Therefore, we will focus on the task of identifying and extracting multiple fractures in the next stage.

Summary and Conclusions
In this paper, the COPS algorithm is proposed which filters background interference information and automatically extracts fractures in ultrasonic logging images. The method's key is to convert the target recognition in the fracture area into optimal path search and search the path according to the fracture's feature that connected on both sides of the image. The main remarks of this study are as follows: (1) Establish a mathematical model of the fracture, determine the reference position of the fracture through the voting accumulation mechanism, and appropriately cut the area where the crack is located according to the reference position (2) The images obtained by cutting are further cut along the vertical direction to obtain each subregion with partial fracture information, and the ACA was used to search and extract the fractures in each subspace (3) The fractures extracted from each subspace are spliced to obtain a complete fracture Experiments and analysis have been conducted to illustrate our method's advantages and potential disadvantages. In the next stage, we will concentrate on improving multiple fractures extraction in a well logging image.

Data Availability
The well logging data used to support the findings of this study were supplied by Ao Qiu under license and so cannot be made freely available. Requests for access to these data should be made to Wei Zhang (weizhang@uestc.edu.cn).

Conflicts of Interest
The authors declare that they have no conflict of interest.