Regular Article
Identification of the fractures of carbonate reservoirs and determination of their dips from FMI image logs using Hough transform algorithm
^{1}
Department of Mining, Petroleum and Geophysics, Shahrood University of Technology, Shahrood 3619995161, Iran
^{2}
Department of Petroleum and Chemical Engineering, Science and Research Branch, Islamic Azad University, Tehran 1477893855, Iran
^{*} Corresponding author: kamkarr@yahoo.com
Received:
15
October
2020
Accepted:
29
March
2021
Carbonate reservoirs are of great importance due to having many fractures and the effectiveness of these fractures in oil production. The most effective tools for studying fractures are image logs that capture high resolution images from the well. An example of these images is the FMI tool, which provide important information on the orientation, depth, and type of fracture. Today, the detection of fractures on these logs is done manually, which in the absence of sufficient experience, will encounter errors. The purpose of this study is to identify the reservoir fractures and the dips of the fractures using Canny edge detection algorithm and Hough transform algorithm and image processing operators, so that in the first stage, fractures are identified in Geolog Software and in the second stage, using MATLAB Software, fractures and their dip are interpreted.
© M. Shafiabadi et al., published by IFP Energies nouvelles, 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
In reservoirs with natural fractures, these fractures control the behavior of the reservoir. When fractures are open, they create pathways for hydrocarbon to move into the well, and may even turn low permeability reservoirs into highyield samples. Conversely, when fractures are filled or cemented, they act as a barrier to the movement of hydrocarbons to the well (Haller and Porturas, 1998; Nelson, 1985; Serra, 1989). The core is the main source of information about smallscale fractures. However, the core has several major limitations for studying the fractures: high cost, low recovery at break intervals, and so on. Changing the direction of the brain during brain thinking does not include these limitations against image logs (Khoshbakht et al., 2012; Safinya et al., 1991). Borehole imaging logs is actually a virtual image of a well that has a high resolution which can display the details of the well (Burgess and Peter, 1985). The FMI log provides up to 80% coverage for wells 8.5 inches in diameter. These logs provide important information on the orientation, dip, and type of fracture (Khoshbakht et al., 2012). Torres et al. (1990) Proposed a method based on the Hough transform function to identify fractures. Hall et al. (1996) used the Hough transform method in 3D space to identify fractures. Kherroubi (2008), in an article entitled automatic separation of open fracture traces from the borehole imaging logs, separated fracture traces by morphological operations. Assous et al. (2013), in a report entitled automatic detection of plate characteristics in wells, presented a new algorithm for detecting plate properties in microresistivity wells. Taiebi et al. (2017) automatically identified sinus fractures using Walsh–Hadamard algorithms and the Kmeans clustering method and Hough transform. In their major work, they also relied on edge information to isolate sinusoidal plate features. Shafiabadi et al. (2020) identified the fractures of the reservoir from the FMI image log using Sobel and Canny edge detection algorithms along with comparison of the performance of these two algorithms. The purpose of this study is to identify the reservoir fractures and to find the fractures dips using Hough transform algorithm and image processing operators without using Artificial Intelligence (AI) methods.
1.1 FMI imaging tools
The FMI (Formation MicroImager) was developed in 1991 by Schlumberger (Schlumberger, 1991). It has 4 arms and 4 pads with 24 electrodes on each pad and a total of 192 electrodes. The Formation MicroImager tool is capable of radial microresistivity measurements (vertical resolution: 0.200 (0.5 cm), vertical sampling: 0.100, and depth of investigation: 3000). The image logs of the well provide a cylindrical image. Now, if this element is cut along its axis and opened, flat structures appear as sinusoidal curves (Fig. 1) (Schlumberger, 2002; Serra, 1989).
Fig. 1 Picture of FMI with sine curve of fracture in imaging logs (Schlumberger, 1994; Schwartz, 2006). 
2 Geological setting
The Zagros Sedimentary Basin is part of the southern margin of the Tethys Ocean and one of the most important oil basins in the world (Alavi, 1994). One of the most important zones of the Zagros Belt and Thrust Belt is Dezful embayment that in SW Iran. Dezful embayment is a cavity that developed in the Early Miocene at the front of the Izeh Zone (Bordenave and Hegre, 2005; Kazemi, 2009). Gachsaran Oilfield is one of the biggest oilfields in the Dezful Embayment (Fig. 2) (Bordenave and Hegre, 2010). In this research, two wells from the Gachsaran oil field located in southwest of Iran have been selected for processing and interpretation using Geolog Software and MATLAB.
Fig. 2 Location map of the studied oilfield in the Dezful Embayment of southwest Iran (Sepehr and Cosgrove, 2004). 
3 Image log interpretation
3.1 Open fracture
The craters of these fractures are filled by drilling mud, and if the mud is conductive, the resistance recorded by the image log is much lower than that of the rock background. Therefore, open fracture is seen as a complete or continuous sine wave in the image log (Fig. 3a) (Schlumberger, 2002).
Fig. 3 Example of fractures in wellbore image logs. (a) Open fracture, (b) closed fractures, (c) fault in FMI log in well A and B. 
3.2 Closed fracture
These fractures are more resistant than their surroundings. When the FMI passes through a fractured fracture, due to the diffusion of the flow injected into the well, there is a difference in resistance (dark or light) between the two sides of the fracture plane in the area of fractures. This difference in resistance is called a halo effect, which is a sign of identifying fractures (Fig. 3b) (Schlumberger, 2002).
3.3 Fault
The faults are complications such as fractures, except that the plate displacement is greater in the faults. Although it is important to detect faults in the image logs, this phenomenon is rarely detectable because of the small diameter of the drilled well (Serra, 1989). An example of the fault observed at well B is shown in Figure 3c. One of the reasons for the existence of faults at this depth can be considered as fracture density in this area.
4 Image processing
Image processing is the application of a set of techniques and algorithms to a digital image to analyze, enhance, or optimize image characteristics such as sharpness and contrast. Imageprocessing terminology is mainly derived from the mathematical functions used (e.g., Laplacian and Gaussian filters are used for image enhancement) (Gonzalez and Woods, 1992).
4.1 Histogram equalization
Increasing contrast is one of the important issues in image processing. Histogram Equalization (HE) is one of the common ways to improve contrast in digital images. Histogram equalization increases the contrast of the image to its original state, which means improved image quality and increased processing accuracy (Hum et al., 2014).
4.2 Contrast
Image arithmetic is the implementation of standard arithmetic operations, such as addition, subtraction, multiplication, and division, on images. The addition of two images is the addition of the brightness of the pixels corresponding to the two images. One of the most important uses of adding two images is to add a backdrop to the image. Subtraction two images to reduce the brightness of the pixels corresponding to the two images. Imdivide is used to split two images. Multiplication is used to multiply two images (Gonzalez and Woods, 1992). Multiplication brightens grayscale images and darkens image division.
4.3 Edge detection
Edge detection is actually a set of mathematical operations that can help identify areas of the image where the brightness changes drastically. Edges are usually defined as curved lines. Physical changes appear as color changes and brightness changes as edges in the image. (Umbaugh Scott, 2010).
4.3.1 Canny edge detection
Gradient detection methods reveal edges by searching for the maximum or minimum of the first derivative in the image. The gradient is expressed as (Gonzalez et al., 2004):
The gradient range is obtained from the following equation:
The direction (angle) is calculated as follows:
Among the most important edge detection filters based on the first derivative, we can mention the Roberts, Prewitt, Canny and Sobel filters. The Canny edge detection was created by John F. Canny in 1983 (Canny, 1983). The Canny algorithm is very sophisticated and, of course, contains high accuracy results but will take a long time to implement. The steps in the Canny algorithm are as follows (Gonzalez and Woods, 1992; Liu et al., 2016; Muthukrishnan and Radha, 2011):
Removal of noise, unwanted textural artefacts and unwanted edges, through Gaussian smoothening of the image data.
Finding gradients; after smoothing the image next step is to find the strength of edge by taking the gradient of the image. Sobel operator is used to perform 2D spatial gradient measure on an image. The edges should be marked where the gradients of the image have large magnitudes, finding the gradient of the image by feeding the smoothed image through a convolution operation with the derivative of the Gaussian in both the vertical and horizontal direction.
Nonmaximum suppression; only local maxima should be marked as edges. Finds the local maxima in the direction of the gradient, and suppresses all others, minimizing false edges.
Double threshold. Potential edges are determined by thresholding, instead of using a single static threshold value for the entire image, the Canny algorithm introduced hysteresis thresholding, which has some adaptively to the local content of the image. There are two threshold levels, th, high and tl, low where th > tl. Pixel values above the th value are immediately classified as edges.
Edge tracking by hysteresis. Final edges are determined by suppressing all edges that are not connected to a very strong edge.
The lower the threshold, the more lines become detectable, on the other hand a high threshold may lose weak lines or parts of the lines. σ plays the role of a scale parameter for the edges. Large values of σ produce coarser scale edges and small values of σ produce finer scale edges. Larger values of σ also result in greater noise suppression. This σ should not be confused with standard deviation (Huang and Wang, 2008).
4.4 Hough transform
The Hough transform invented by Richard Duda and Peter Hart in 1972 (Duda and Hart, 1972). The Hough transform is a way of extracting features in image analysis, computer vision, and digital image processing (Duda and Hart, 1972; Gonzalez et al., 2004). Since the boundaries are determined at the edge detection stage, the exact shape of the edge segments that correspond to fractures are found using pattern recognition in order to characterize these curves fully. The Hough transform is an established method for detecting complex patterns of points in binary image data, and has been known to perform well in the presence of noise, extraneous data and occlusions (Illingworth and Kittler, 1988). This edge description is commonly obtained from a feature detecting operator such as the Roberts Cross, Sobel or Canny edge detector and may be noisy, i.e. it may contain multiple edge fragments corresponding to a single whole feature. The idea of the Hough transform is, that every edge point in the edge map is transformed to all possible lines that could pass through that point (Fig. 4). After applying the edge detection algorithm and getting the final results, the next step is to use the Hough algorithm. The input image in the Hough algorithm is the colored image (RGB) obtained from the result of the Canny edge detection algorithm.
Fig. 4 Represents a line in image space as a point in Hough space (Wang and Wang, 2018). 
4.5 Curve fitting
Curve fitting is the process of constructing a curve, or mathematical function, that has the best fit to a series of data points (Fig. 5). Most commonly, one fits a function of the form y = f(x) (Kolb, 1984).
Fig. 5 Apply the Hough algorithm on the detected edge of the closed fracture in well A. 
4.6 Find dip of fractures
The dip is an angle that can have any value between 0 (horizontal plane) to 90° (vertical plane) (Aguilera, 2010). Figure 6 shows how to calculate the dip from the sinusoidal curves in the FMI Log.
Fig. 6 Sinusoidal curve dip calculation method in FMI images log. Dip angle = arc tan (h/d) where h = height of sinusoid and d = borehole diameter (Serra, 1989). 
5 Discussion and results
The field of study consists of FMI images of two wells A and B, for which there are two sections selected from each well. More information on these sections is given in Table 1 and in Figure 7. Figure 7 shows the open and closed fractures and faults in wells A and B along with their σ (σ is the standard deviation of the image data here). With the help of the histogram equalization method, the histogram equals the existing images and is the basis that any image with any histogram can reach the gray pixel frequency equilibrium. We then create contrast images by applying the contrast operators (addition, subtraction, multiplication and division) in the number range shown in Table 2 below. Contrasting gray scale images can cause the edges to intensify or fade. Hence, obtaining the correct sequence of operators and the number intervals will give a better edge representation in the edge recognition algorithm.
Fig. 7 Displays open and closed fractures and fault along with their standard deviation. (a) Closed fracture, (b) open fracture in well A, (c) closed fracture and (d) fault in well B. 
Table of dip and number of pixels of open and closed fractures and fault.
Contrast change interval in imadd, imsubtract, immultiply, imdivide operators.
Figure 8 illustrates the application of contrast operators such as addition, subtraction, and division on fracture and the effect of each of them on the shape. Given the better performance of the addition operator seen in Figure 8, we can obtain different results from the contrast by adding value on closed fracture as shown in Figure 9. To perform the Canny edge detection method on gray scale images, number 4 is placed in the range of 0–1 in the threshold interval, in which there are 6 different states that can be seen in Table 3. These six states create an interval of threshold changes in the edge recognition operator, so that their results in finding the edge differ from those of the other states. On the other hand, the choice of different σ creates new statistical conditions for the edges in each of the 6 states. The study interval for σ is selected from to . Because in the binary image edge is white and the background is black, it is not appropriate to use the multiplication operator. In total, for each case study, 12 000 images were created that the Canny algorithm was able to identify the fracture edges in a certain range of thresholds, σ and operator numbers and are presented in Table 4. Table 4 presents a summary of the threshold and σ and its operators, which due to the high volume of data obtained from this algorithm, we cannot express the whole results in this paper, so some examples are shown in Table.
Fig. 8 Displays changes made after adding and subtracting and division operators on closed fractures. 
Fig. 9 Examples of adding value and changing the edges. 
Six states made for threshold usable in Canny edge detection.
Summary of the appropriate range of threshold and σ and operator parameters.
Figure 10 shows the dip of the fractures obtained with the azimuth of each section. By selecting a range of fracture detection variables using the Hough transform algorithm, we convert the pixels in images adjacent to 3 pixels lines from the binary image space to the coordinate space and create the edges as discontinuous points in the coordinate space we have. Figure 11 shows the application of the Hough transform algorithm to identify fractures on open and closed fractures and fault. Using the existing continuous lines, we can calculate the fracture angle in each case (Tab. 5). The values comparable to the results of the Geolog Software and the geological interpretation are given in Table 5. The final results of fracture identification and dip determination in this study show that the geological results are close to the Hough transform algorithm.
Fig. 10 Displays the dip of fractures in Geolog Software, (a) closed fracture and (b) fracture in well A, (c) closed fracture and (d) fault in well B. 
Fig. 11 The result of applying a Hough transform algorithm. 
Comparison of fracture dip obtained from the proposed algorithm and geologist and Geolog Software.
6 Conclusion
This paper proposed a novel methodology for fracture extraction and fracture dip determination from FMI image logs. This method was performed in three steps. In the first step, the raw data were processed by Geolog Software, and as a result their fractures and dips were determined. In the second step, using the Canny edge detection algorithm, Hough transform algorithm and image contrast operators, their fractures and dips were determined. In the last step, the results of the above two methods were compared with the main results of interpretation. Consequently, the use of Hough transform algorithm in this study to determine the dips of fractures contained the same values of the interpretation results. Finally, it can be concluded that the proposed method for identification and determination of fractures has shown excellent performance.
References
 Aguilera R. (2010) Effect of fracture dip and fracture tortuosity on petrophysical evaluation of naturally fractured reservoirs, J. Can. Petrol. Technol. 49, 9, 69–76. [Google Scholar]
 Alavi M. (1994) Tectonics of the Zagros orogenic belt of Iran: new data and interpretations, Tectonophysics 229, 3, 211–238. [Google Scholar]
 Assous S., Elkington P., Clark S., Whetton J. (2013) Automated detection of planar geological features in borehole images, Geophysics 79, 1, 11–19. [Google Scholar]
 Bordenave M.L., Hegre J.A. (2005) The influence of tectonics on the entrapment of oil in the Dezful Embayment, Zagros Foldbelt, Iran, J. Pet. Geol. 28, 339–368. [Google Scholar]
 Bordenave M.L., Hegre J.A. (2010) Current distribution of oil and gas fields in the Zagros Fold Belt of Iran and contiguous offshore as the result of the petroleum systems, Geol. Soc. Lond. Spec. Publ. 330, 291–353. [Google Scholar]
 Burgess C.J., Peter C.K. (1985) Formation, distribution, and prediction of stylolites as permeability barriers in the Thamama Group, Abu Dhabi. in: Middle East Oil Technical Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Canny J.F. (1983) A variational approach to edge detection, pp. 54–58. [Google Scholar]
 Duda R.O., Hart P.E. (1972) Use of the Hough transformation to detect lines and curves in pictures, Comm. ACM 15, 1, 11–15. [Google Scholar]
 Gonzalez R.C., Woods R.E. (1992) Digital image processing Reading, AddisonWesley, MA. [Google Scholar]
 Gonzalez R.C., Eddins S.L., Woods R.E. (2004) Digital image publishing using MATLAB, Prentice Hall. [Google Scholar]
 Hall J., Ponzi M., Gonfalini M., Maletti G. (1996) Automatic extraction and characterization of geological features and textures from borehole images and core photographs, in: SPWLA 37th Annual Logging Symposium, New Orleans, Louisiana, pp. 1–13. [Google Scholar]
 Haller D., Porturas F. (1998) How to characterize fractures in reservoirs using borehole and core images: Case studies, Geol. Soc. Spec. Publ. Lond. 136, 1, 249–259. [Google Scholar]
 Huang Y., Wang S. (2008) Multilevel thresholding methods for image segmentation with Otsu based on QPSO. in: 2008 Congress on Image and Signal Processing, IEEE, Vol. 3, pp. 701–705. [Google Scholar]
 Hum Y.C., Lai K.W., Mohamad Salim M.I. (2014) Multiobjectives bihistogram equalization for image contrast enhancement, Complexity 20, 2, 22–36. [Google Scholar]
 Illingworth J., Kittler J. (1988) A survey of the Hough transform, Comp. Vision Graph. Image Process 44, 1, 87–116. [Google Scholar]
 Kazemi K. (2009) Seismic imaging of thrust fault structures in Zagros Iranian oil fields, from subsurface and well data, in: 71st EAGE Conference & Exhibition, Amsterdam, The Netherlands. [Google Scholar]
 Kherroubi J. (2008) Automatic extraction of natural fracture traces from borehole images, in: 19th International Conference on Pattern Recognition (ICPR), IEEE, Tampa, pp. 1–4. [Google Scholar]
 Khoshbakht F., Azizzadeh M., Memarian H., Nourozi G.H., Moallemi S.A. (2012) Comparison of electrical image log with core in a fractured carbonate reservoir, J. Pet. Sci. Eng. 86, 289–296. [Google Scholar]
 Kolb W.M. (1984) Curve fitting for programmable calculators, Imtec. [Google Scholar]
 Liu X., Dawn Z., Xu W. (2016) Image edge deduction algorithm based on improved wavelet transform, Int. J. Signal Process. Image Process. Pattern Recogn. 9, 4, 435–442. [Google Scholar]
 Muthukrishnan R., Radha M. (2011) Edge detection techniques for image segmentation, Int. J. Comput. Sci. Inf. Technol. (IJCSIT) 3, 6, 259–267. [Google Scholar]
 Nelson R.A. (1985) Geologic analysis of naturally fractured reservoir, Gulf Publishing Company, Houston, Texas, USA. [Google Scholar]
 Safinya K.A., Lan P.L., Villegas M., Cheung P.S. (1991) Improved formation imaging with extended microelectrical arrays, Soc. Pet. Eng. 22726, 653–664. [Google Scholar]
 Schlumberger (1991) Log Interpretation Principles/Applications, Schlumberger Educational Services. [Google Scholar]
 Schlumberger (1994) FMI fullbore formation microimager, Schlumberger Educational Services, Houston, Texas. [Google Scholar]
 Schlumberger (2002) FMI: borehole geology, geomechanics and 3D reservoir modeling. http://www.slb.com/media/services/evaluation/geology/fmi.pdf. [Google Scholar]
 Schwartz B.C. (2006) Fracture pattern characterization of the Tensleep Formation, Teapot Dome, Wyoming, Master Thesis, Department of Geology and Geography Morgantown, West Virginia, pp. 115–130. [Google Scholar]
 Sepehr M., Cosgrove J.W. (2004) Structural framework of the Zagros FoldThrust Belt, Iran, Mar. Pet. Geol. 21, 7, 829–843. [Google Scholar]
 Serra O. (1989) Formation microscanner image interpretation, Schlumberger Education Services, Houston, TX. [Google Scholar]
 Shafiabadi M., KamkarRouhani A., Ghavami Riabi S.R., Roshandel Kahoo A., Tokhmechi B. (2020) Identification of reservoir fractures on FMI image logs using Canny and Sobel edge detection algorithms, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 76, 11. https://doi.org/10.2516/ogst/2020086. [Google Scholar]
 Taiebi F., Akbarizadeh G.H., Farshidi E. (2017) Detection of reservoir fractures in imaging logs using directional filtering, in: 2017 5th Iranian Joint Congress on Fuzzy and Intelligent Systems (CFIS), IEEE, pp. 150–154. [Google Scholar]
 Torres D., Strickland R., Gianzero M. (1990) A new approach to determining dip and strike using borehole images, in: SPWLA 31th Annual Logging Symposium, Lafayette, Louisiana, pp. 1–20. [Google Scholar]
 Umbaugh S.E. (2010) Digital image processing and analysis, human and computer vision applications with CVIP tools, 2nd edn., CRC Press, Boca Raton, FL. ISBN 9781439802052. [Google Scholar]
 Wang G., Wang W. (2018) The research on edge detection algorithm of lane. EURASIP J. Image Video Process. 2018, 1, 1–9. [Google Scholar]
All Tables
Comparison of fracture dip obtained from the proposed algorithm and geologist and Geolog Software.
All Figures
Fig. 1 Picture of FMI with sine curve of fracture in imaging logs (Schlumberger, 1994; Schwartz, 2006). 

In the text 
Fig. 2 Location map of the studied oilfield in the Dezful Embayment of southwest Iran (Sepehr and Cosgrove, 2004). 

In the text 
Fig. 3 Example of fractures in wellbore image logs. (a) Open fracture, (b) closed fractures, (c) fault in FMI log in well A and B. 

In the text 
Fig. 4 Represents a line in image space as a point in Hough space (Wang and Wang, 2018). 

In the text 
Fig. 5 Apply the Hough algorithm on the detected edge of the closed fracture in well A. 

In the text 
Fig. 6 Sinusoidal curve dip calculation method in FMI images log. Dip angle = arc tan (h/d) where h = height of sinusoid and d = borehole diameter (Serra, 1989). 

In the text 
Fig. 7 Displays open and closed fractures and fault along with their standard deviation. (a) Closed fracture, (b) open fracture in well A, (c) closed fracture and (d) fault in well B. 

In the text 
Fig. 8 Displays changes made after adding and subtracting and division operators on closed fractures. 

In the text 
Fig. 9 Examples of adding value and changing the edges. 

In the text 
Fig. 10 Displays the dip of fractures in Geolog Software, (a) closed fracture and (b) fracture in well A, (c) closed fracture and (d) fault in well B. 

In the text 
Fig. 11 The result of applying a Hough transform algorithm. 

In the text 