Identification of reservoir fractures on FMI image logs using Canny and Sobel edge detection algorithms

Because of the significant impact of fractures on production in hydrocarbon reservoirs, identification of these phenomena is a very important issue. Image logs are one of the best tools for revealing and studying fractures in reservoir and researcher can get lots of information about geological features in wells, by studying and analyzing these logs. In this research, two approaches have been used to determine the fractures in two wells A and B located in one of the oil fields in southwest of Iran. In the first approach, using Geolog software (version-7), after processing and correction of raw image log data, the number, position, dip, extension, layering, density and expansion of fractures have been identified. In the second approach, considering that the fractures in FMI images have edges, the Canny and Sobel filters as edge detection operators in image processing have been used to detect fractures in these images.


Introduction
Fractures are one of the important complications in reservoir heterogeneity related to fluid flow and are more characteristic of large carbonate oil fields (Rezaie, 2008). The term fracture refers to any physical breakdown or discontinuity in the rock that is above the threshold of rock resistance. Fractures include joints and faults (Zahrabzadeh, 2009). Natural fractures control the behavior of the reservoir. When the fractures are open, they create pathways for the movement of hydrocarbons to the well, and may even convert low permeability tanks into high-yielding samples. Conversely, when the fractures are closed or cemented, they act as a barrier against the movement of hydrocarbons towards the well. This dual behavior of fractures makes it necessary to identify them (Javid et al., 2012). Direct and indirect methods for identifying fractures include the use of seismic sections, petrophysical logs, well testing, flower blasting, and descriptions of cores and image logs (Tingay et al., 2008). Seismic sections have low resolution for detection of fractures but normally identify large-scale faults (Tokhmechi, 2009). Cores are the main source of information for small scale fractures. In addition, the cores for study of the fractures have three major limitations: high cost, low recycle at broken intervals and core shift during core sampling, while such limitations do not exist in image logs. The image logs present high-resolution virtual images of well walls that are capable of displaying delicate wall phenomena (Kulander et al., 1987). Image analysis began with the start of its use in 1986. Since then, many studies have begun to identify fractures on these logs. Torres et al. (1990) proposed a method based on Hough transform to detect the amplitude and phase of the fracture sine curve and the gradient corresponding to the fractures. This method is not the perfect response for interpretation and analysis of the fractures (Torres et al. 1990). Ye and Baviler (1998) proposed a method based on Fourier transform to separate the padding and heterogeneous plates of each other. However, this method works weak for all the carbonate reservoirs (Ye and Baviler, 1998). Wang (2005) used the Valley-edge algorithm to find fractures in the rocks without using threshold values. By this way, he detected the edges of the Charge Coupled Device (CCD) images taken from rocks. A multi-scale technique was used in the Wang's algorithm for detecting the edges (Wang, 2005). Wang et al. (2007) provided a method based on edge detection and backup machines (Wang et al., 2007). Assous et al. (2014) represented a new algorithm for detection of plane features on micro-resistance images of a well in a report titled "Automatic detection of plane features in image logs". To separate the sinusoidal plane features, they also relied on information of the edges in the main part of their work (Assous et al., 2014). Considering the importance of fractures and problems in image interpretation, the purpose of this study is to introduce and apply image processing techniques and edge detection algorithms in image processing called Canny and Sobel algorithms to identify fractures in image logs. The purpose of this study is to use the Canny and Sobel algorithms to find the edges of the fracture and finally compare the two methods to find the best algorithm to do this.

Imaging logs
In 1985, Schlumberger introduced the Formation Micro Scanner (FMS) service. The FMS sonde consists of four orthogonal imaging pads each containing 16 microelectrodes that are in direct contact with the borehole wall during the recording. The tool works by emitting a focused current from the four pads into the formation. Electrical imaging device FMI (Formation Micro Imager) is the second generation of electrical imaging devices that was produced in 1991 by Schlumberger Petroleum Company. This device was designed to increase the wall cover of the well and to obtain a more coherent image than the FMS (Javid et al., 2012). The structure of the microprocessor is such that the pads and flaps are each mounted on 24 copper electrodes in two rows. The device consists of 192 electrodes with an effective diameter of 2 inches (Fig. 1). This machine, with 80% borehole coverage, produces high quality images. It is capable of analyzing very thin layers and detecting the type of fine structures such as cross-layers and stylolites, very fine fractures and distribution analysis of complex reservoirs (Schlumberger, 2004).
In general, an image log is taken cylindrically from the wellbore; for interpretation and observation, the image is opened and flattened and its geographic directions are determined (Fig. 2). Then, the plate structures that cut the wall of the well horizontally, in the form of sine curves with the following general equation (Rider, 1996): where A = R tan u and B ¼ p 2 À b. In equation (1), R is the well radius, u is the dip angle, b is the azimuth and C is the sinusoidal position in the image.
The processes performed on the images are divided into two categories: A: The fundamental corrections necessary for electrical image logs include depth correction, speed correction, correction of broken and disabled pads, and so on. B: Improving image quality processing, which can be equalization and image normalization.
By interpreting and processing these images, parameters such as layering, fracture, stylolite, and fault can be detected (Serra, 1984(Serra, , 1989.

Geological setting
The study area is one of the oil fields in southwest of Iran. It is situated in Zagros belt, which extends in northwestsoutheast direction. This fold-thrust belt extends about 2000 km from southeastern Turkey through northern Syria and Iraq to western and southern Iran. The study field is located in Dezful embayment of Zagros basin in southwest of Iran (Alavi, 2004). Geological history of this basin includes long time subsidence and deposition, which has interrupted by short time uplift. Folding process of this basin occurred in Miocene and Pliocene and continued until now that led to forming long anticlines (Mottie, 1995). In this research, image logs or data as well as other relevant data and information from two wells within Gachsaran oil field, located in southwest of Iran, have been selected (Fig. 3). Processing and interpretation of these data have been carried out in Geolog software. The Gachsaran oil field is an elongated doublyplunging anticlinal structure with dimension of 70 km length and 6-15 km width, and it is considered as one of the most important productive oil fields in the Oligocene-Lower Miocene carbonate horizons and the Middle Cretaceous carbonate horizons ( Fig. 4) (Mottie, 1995).

Methodology
As a result of the processing the FMI images from wells A and B, the sinusoidal curves of the closed fractures are displayed as clear lines and open fractures as dark lines and induced fractures due to drilling operations as  (Schlumberger, 2005).

Fig. 2. The collision plates with well walls (fractures or layers)
are represented in the image logs as sinusoidal curves (Serra, 1989). perpendicular to the layering (Fig. 5). In the image logs, faults are complications such as fractures, with the difference that the displacement of the plates in the faults is more. Although fault detection is very important in image logs, but due to the low diameter of the well, this phenomenon is rarely identifiable (Rider, 1996). As shown in Figure 6, the fracture planes are seen in the image (sine lines dark and light); They are selected by the arrows in the image as an open fracture (dark lines) or a closed fracture (the light lines) and the observed fault which are identified as blue, red, and pink sinuous curves with their dip and azimuth, respectively. Due to the dissolution of the well wall caused by the drilling fluid, the wall of the well may breakout. Therefore, the loosening of the well wall prevents full attachment of the pad to the wall, and at all these intervals, there is a black image. An example of a wellbore breakout is seen in well B at a depth of 2365-2368 m in the following figure that is confirmed by the caliper logs referred to by arrows (Fig. 7). After studying these two wells, a total   Table 1, the results obtained from the main interpretation of the wells made with the help of Geolog software are shown for comparison. As can be seen from these two interpretation results shown in Table 2, no difference is seen between the open fractures in well B and the closed fractures in well A. However, the number of open fractures in well A identified by these two interpretation results vary by "1" difference. Similarly, a difference of "2" exists between the number of the closed fractures in well B are seen between these two interpretation results. After identification of fractures and considering that the number of fractures in each region is high, we draw these fractures  in the stereo diagram to obtain their statistical maxima. In this study, the total trend of fractures was investigated by plotting the rose diagram. With regard to the general trend of fractures, it can be concluded that the dominant fractures (open and closed) in well A have northeast-southwest trends, and a general east-west trend can be seen in well B (Fig. 8). Figure 9 demonstrates the general trend of the faults from the FMI images taken from wells A and B. The pink sinusoidal curves in this figure displayed by the arrow indicate the plane of fault detected at wells A and B. The fault dip in well A is 19°with azimuth 16°from north to east. However, in well B, the fault dip is 69°with azimuth 357°from east to west. According to Figure 9, the general trend of the fault in well A is approximately northwest-southeast and the trend of the fault in well B is approximately northeast-southwest.

Edge detection techniques
Edge detection means the process of detecting and locating sharp brightness discontinuities in an image. Discontinuities are sudden changes in the brightness of pixels that determine the boundary of the objects in the image (Clerk Maxwell, 1868;Gonzalez and Woods, 2002). Revealing the edges in images with noise is a difficult task because both the noise and the edges contain high frequencies.
Trying to reduce the noise of the images by melting the image will cause the edges to fail. Edge detector operators that reduce noise also reduce the accuracy of the edges that are visible due to changes in the image (Maini and Aggarwal, 2011). Generally, different edge detection methods are classified into two categories: gradient-based methods and Laplacian-based methods.

Gradient-based techniques
Edge detection methods can be classified into directional and non-directional or gradient-based operators. Directional operators use two masks and two convolutions. While non-directional use single mask and convolution but they are sensitive to noise due to gradient nature of the operators. Gradient-based detection methods reveal edges by searching for the maximum or minimum of the first derivative in the image (Gonzalez and Woods, 1992).

Laplacian-based techniques
In the Laplacian-based detection method, the location of the zero crossing of the second derivative of the image is searched for and defined as follows: Among the most important edge detection filters that are the basis of their first derivative work, we can refer to the Roberts, Prewitt, Canny and Sobel filters, and the second derivative filters can also be used to describe the Marr-Hildreth and LOG filter (Gonzalez and Woods, 2002).

Sobel edge detection operator
The Sobel edge detection method is introduced by Sobel (1970). The Sobel operator also uses the first derivative to find the edges of the image. This edge is one of the weakest edges, and can only detect edges with high light intensity variations. The Sobel technique performs a 2-D spatial gradient quantity on an image, and thus, highlights regions of high spatial frequency that correspond to edges (Sobel, 1970). In general, it is used to find the estimated absolute gradient magnitude at each point in n input gray-scale image. In conjecture, at least the operator consists of a pair of 3 Â 3 complication kernels as given in the tables shown in Figure 10. The gradient magnitude is given by Gonzalez and Woods (1992): and typically, an approximate magnitude is computed using: The angle of orientation of the edge (relative to the pixel grid) giving rise to the spatial gradient is given by: uðx:yÞ ¼ arctanðGx=GyÞ:

Canny edge detection operator
In industry, the Canny edge detection technique is one of the standard edge detection techniques. It was first created by John Canny in his Master's thesis at MIT in 1983, and still outperforms many of the newer algorithms that have been developed (Owens, 1997). Canny algorithm is very sophisticated and, of course, with high precision results that

Results and discussion
In this section, the Canny and Sobel edge-detection algorithms are implemented in MATLAB environment to find the edges of fractures, and therefore, to identify the fractures on the FMI image logs. The first step in this research is to introduce and define the input image in the MATLAB software. The input image is a dynamic and normalized image originally generated from the FMI raw material data. The existence of vertical white lines on the image that has been made after correction in the Geolog software, causes the discontinuities of the images to appear, and poses a unified fragmentation pixel accumulation. Therefore, using the Geolog software, these lines are removed from the image and the image is continuously added. Since the color images are composed of three Red, Green and Blue (RGB) color matrices, these color images are converted to a gray-scale images to reduce the volume of operations, so that we can display the image with only one matrix. On the FMI imaging screens, this filter is applied throughout the two wells A and B (Fig. 11).

Removal of noise from images
Noise contained in image is smoothed by convolving the input image I (i, j) with Gaussian Filter G. Mathematically, the smooth resultant image is given by An approach to noise reduction in the image logs is to use the Non-Local Minimum (NLM) filter. The NLM filter is a common filter for processing medical images. The window size of this filter varies depending on the brightness of the Well A Well B image. The NLM filter is implemented by applying equations (7)-(10) (Ergin and Kilinc, 2014): where, I(q) represents a 2D image, p is the desired pixel in the image to perform the filtering operation, q is the pixel in the neighboring g of the research window, and w (p, q) is also called the weight function, which in the relation, the production of its components is shown: where d is the Euclidean distance between the p and q pixels; h is called the continuous kernel, which controls the reduction of the exponential function; and z(p) is the normalization constant calculated by the following equation.

Finding gradients
In this step we detect the edges where the change in grayscale intensity is maximum. Required areas are determined with the help of gradient of images. Sobel operator is used to determine the gradient at each pixel of smoothened image (Gonzalez and Woods, 2002): : These Sobel masks are convolved with smoothed image and giving gradients in i and j directions: Therefore, edge strength or magnitude of gradient of a pixel is given by The direction of gradient is given by

Non-maximum suppression
Non maximum suppression is carried out to preserve all local maxima in the gradient image, and deleting everything else this results in thin edges. For a pixel M(i, j): Firstly round the gradient direction u nearest 45°, then compare the gradient magnitude of the pixels in positive and negative gradient directions i.e. If gradient direction is east then compare with gradient of the pixels in east and west directions say E (i, j) and W (i, j) respectively. If the edge strength of pixel M (i, j) is largest than that of E (i, j) and W (i, j), then preserve the value of gradient and mark M (i, j) as edge pixel, if not then suppress or remove.

Hysteresis thresholding
Thresholding is the first step in most image processing systems and the simplest method of image segmentation. At this point, the image field is separated from the objects by finding a threshold (Huang and Wang, 2008). r is the standard deviation parameter of the gaussian filter and its value is fixed as per the degree of noise present in the input image (Bland and Altman, 1996;Gonzalez and Woods, 1992) are used. For a pixel M(i, j) having gradient magnitude G following conditions exists to detect pixel as edge: If G < t low than discard the edge. If G > t high keep the edge. If t low < G < and t high and any of its neighbors in a 3 Â 3 region around it have gradient magnitudes greater than t high , keep the edge. If none of pixel (x, y)'s neighbors have high gradient magnitudes but at least one falls between t low and, t high search the 5 Â 5 region to see if any of these pixels have a magnitude greater than thigh. If so, keep the edge. Else, discard the edge.
In this research, we use two thresholds (high and low thresholds) as follows: (a) (b) (c) Fig. 11. Conversion of (a) the FMI log to (b) the gray-scale image, and (c) the black and white (binary) image.
The reason for using this Terscheld range in this study was that the lower the thresholds such as Th = [0.1, 0.2] applied to the image, the more edges would appear other than the fractures. The higher the thresholds such as Th = [0.7, 0.9] are applied to the image, the edge of the sinusoidal curve of the fractures is removed. As a result, the most moderate Terscheld in this study (Th = [0.3, 0.4]) was selected so that neither the many edges of the fractures were removed nor the many edges other than the fractures appeared.

The edge detection result and analysis
Since color plays an important role in the analysis of images and color images are more informative than the gray ones, thus, in this study, the Canny edge algorithm operations are applied to the color images to achieve a more favorable result. Figure 12 shows the results of the closed fracture of FMI images in wells A and B at along with the changes made the Canny algorithm operation and its edges. Figure 13 shows the results of the Canny algorithm operation, along with its changes on the borehole breakout observed at a depth of 2432-2434 m in well B. Figure 14 shows the results of the closed fracture of FMI images in wells A and B at along with the changes made the Sobel algorithm operation and its edges. Figure 15 shows the results of the Sobel algorithm operation, along with its changes on the borehole breakout observed at a depth of 2432-2434 m in well B. Figure 16 shows the results of the comparison of applying the Canny algorithm with the Sobel algorithm on the fault in the FMI images obtained from well B. Every edge detection method has its own constraints and merits, and to find an appropriate edge (c) (b) (a) detection algorithm for a particular case, we must implement these edge detection algorithms to get the desired results and to find the appropriate one. The Canny algorithm, due to the more detailed presentation of fractures, in comparison with the Sobel algorithm, has a better performance on the FMI images, and has been capable of revealing the sinusoidal curves of the fault and open and closed fractures to a large extent. The comparison results of Canny and Sobel edge detection algorithms are shown in Table 2.

Conclusion
Image logs are among the powerful tools that are used in the study of hydrocarbon reservoirs to identify fractured zones in the walls of wells. Two approaches were taken to identify fractures on FMI logs. In the first approach, after loading the FMI data from two wells, situated in an oil field in southwest of Iran, in the Geolog software, various corrections and processing steps were made to the images. As a result, two static and dynamic images were obtained that indicate the types of fractures and breakouts in the two wells. In the second approach taken in this study, we used the Canny and Sobel edge detection algorithms so that we could show the performance of these edge detectors on the FMI log to identify fractures. The Sobel algorithm indicates very weak margins and discontinuities. It also contains false margins, however, the Canny algorithm shows clean, almost continuous, and real edges strongly. The result shows Canny operator is optimum even for noisy images. This operator fills the gap between strong and weak edges of the image. Compared to other edge detection techniques, this operator is less fooled by specious noise. In Sobel edge detection method, it looks for both horizontal and vertical edges independently then represents them together. As a result, the proposed method can be a suitable solution for identifying fractures in electrical image logs.