Automatic fracture–vug identification and extraction from electric imaging logging data based on path morphology

Petroleum Science, Dec 2018

We present a path morphology method to separate total rock pore space into matrix, fractures and vugs and derive their pore structure spectrum. Thus, we can achieve fine pore evaluation in fracture–vug reservoirs based on electric imaging logging data. We automatically identify and extract fracture–vug information from the electric imaging images by adopting a path morphological operator that remains flexible enough to fit rectilinear and slightly curved structures because they are independent of the structuring element shape. The Otsu method was used to extract fracture–vug information from the background noise caused by the matrix. To accommodate the differences in scale and form of the different target regions, including fracture and vug path, operators with different lengths were selected for their recognition and extraction at the corresponding scale. Polynomial and elliptic functions are used to fit the extracted fractures and vugs, respectively, and the fracture–vug parameters are deduced from the fitted edge. Finally, test examples of numerical simulation data and several measured well data have been provided for the verification of the effectiveness and adaptability of the path morphology method in the application of electric imaging logging data processing. This also provides algorithm support for the fine evaluation of fracture–vug reservoirs.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs12182-018-0282-6.pdf

Automatic fracture–vug identification and extraction from electric imaging logging data based on path morphology

Petroleum Science February 2019, Volume 16, Issue 1, pp 58–76 | Cite as Automatic fracture–vug identification and extraction from electric imaging logging data based on path morphology AuthorsAuthors and affiliations Xi-Ning LiJin-Song ShenWu-Yang YangZhen-Ling Li Open Access Original Paper First Online: 13 December 2018 107 Downloads Abstract We present a path morphology method to separate total rock pore space into matrix, fractures and vugs and derive their pore structure spectrum. Thus, we can achieve fine pore evaluation in fracture–vug reservoirs based on electric imaging logging data. We automatically identify and extract fracture–vug information from the electric imaging images by adopting a path morphological operator that remains flexible enough to fit rectilinear and slightly curved structures because they are independent of the structuring element shape. The Otsu method was used to extract fracture–vug information from the background noise caused by the matrix. To accommodate the differences in scale and form of the different target regions, including fracture and vug path, operators with different lengths were selected for their recognition and extraction at the corresponding scale. Polynomial and elliptic functions are used to fit the extracted fractures and vugs, respectively, and the fracture–vug parameters are deduced from the fitted edge. Finally, test examples of numerical simulation data and several measured well data have been provided for the verification of the effectiveness and adaptability of the path morphology method in the application of electric imaging logging data processing. This also provides algorithm support for the fine evaluation of fracture–vug reservoirs. KeywordsPath morphology Image automatic identification Electric imaging logging Fracture–vug reservoir  Edited by Jie Hao 1 Introduction Fractures and vugs are the main reservoir space and seepage channels in fracture–vug reservoirs (Yang et al. 2011). Due to their complexity and heterogeneity, their distribution mechanism has a significant impact on detailed reservoir evaluation (Li et al. 2015; Wang et al. 2013). Currently, identification and quantitative evaluation of fractures and vugs have become the focus for exploration and development of fracture–vug reservoirs (Wu et al. 2013; Li et al. 2017a, b, c; Thachaparambil et al. 2013). Geologists have worked deeply on pore structure characterization (Bian et al. 2014; Gunter et al. 2014; Chen and Heidari 2014; Holden et al. 2014; Gao et al. 2015). Many articles have been published about the acoustic and electrical parameters of fracture–vug rocks (Wang et al. 2011; Tang et al. 2013; Xu et al. 2014; Chen et al. 2014; Zhang et al. 2015; Pan et al. 2017). Both conventional logging and electric imaging logging have been applied as the main well logging methods for fracture–vug identification and extraction. Though conventional logging has been well developed (Lyu et al. 2016; Zazoun 2013; Shalaby and Islam 2017; Amartey et al. 2017), it cannot achieve accurate results because of the low vertical resolution. Electric imaging logging, which is characterized by its high vertical resolution and borehole coverage, can display fractures and vugs intuitively. It has been a key step for interpreters to identify and describe fracture–vug information from conductivity images (Lai 2011). To calculate the fracture–vug parameters accurately, logging image segmentation plays an important role in separating fractures and vugs from the formation matrix. However, various segmentation algorithms have different effects on extraction accuracy (Tang 2013; Cui et al. 2016; Liu et al. 2017; Xie et al. 2017). For automatic fracture–vug identification and extraction, it is necessary and important to study the image segmentation method suitable for the characteristics of the electric imaging images. Until now, the main electrical imaging logging software, including GeoFrame from the Schlumberger corporation, eXpress from the Atlas corporation, LEAD from CNPC logging and Logview from GEOTECH, requires interpreters to extract the fracture or vug based on a human–computer interaction platform. Thus, it depends on human subjectivity and experience, which results in poor efficiency and involves a lot of labor. Domestic and international scholars have performed much research on the automatic recognition and quantitative evaluation of electrical imaging logging. For example, marker-controlled image segmentation has been used to mark the outline of necessary geological features and compute geological parameters (Delhomme 1992). The Hough transform was used to extract sinusoids from the electric imaging images and automatically calculate the formation occurrence of stratification and fracture (Hall et al. 1996). To trace fractures in the reservoir, a Gauss–Laplace operator was used to remove horizontal geologic features such as the stratigraphic boundary of the enhanced electric imaging image (Chitale et al. 2004). Wavelet transforms were used to estimate the fracture density (Tokhmechi et al. 2009). MATLAB was used to identify fractures automatically, detect sinusoids in the images and calculate the fracture parameters (Cornet 2013). In addition, some scholars draw mathematical morphology into the image processing of electrical imaging logging and achieve good results (Xiao et al. 2015; Xavier et al. 2015; Li et al. 2017a, b, c). With the limitation of the complex geologic environment and logging conditions, the methods above cannot automatically identify fracture–vug information completely and achieve the required quantitative interpretation accuracy. In particular, they cannot accurately describe irregular fractures and vugs. In recent years, in the field of image processing, path morphology, by constructing flexible structuring elements, has been proposed to identify elongated and curved subimages (Heijmans et al. 2005; Talbot and Appleton 2007). The method has been used for medicine, road extraction, etc. (Wang 2014), but not for geophysical logging. With the complexity of the image distribution, this research aims to realize separation of fractures and vugs from a complex geological background and obtain their quantitative parameters automatically. According to the basic theory of path morphology and path operators suitable for a curved structure, we apply path morphology operators to electrical imaging images and develop a novel automatic identification and extraction method for fractures and vugs. By comparing and analyzing the conventional data, Logview software analysis data and core scanning data, we verify the effectiveness and high precision of the path morphology method in electric imaging logging data processing. 2 Response characteristics of the borehole geological structures in the electric imaging images Since the conductivity difference between fracture, vug and matrix with mud invasion, high conductivity fractures and vugs appear dark, while matrix appears bright in the electric imaging images. It provides favorable conditions for efficiently extracting fractures and vugs. However, other geological structures, such as a stratification plane, argillaceous stripe, induced fracture, will always have similar behavior to fractures and vugs. Therefore, before introducing the path morphological algorithm, it is necessary to summarize the characteristics of the borehole geological structures in the electric imaging images used for extracting fractures and vugs. (1) Stratigraphic interfaces and stratification planes: These often appear as a group of continuous complete sinusoidals that are parallel to each other. Moreover, they have uniform thickness, gentle dip direction and angle (Fig. 1a). Open image in new window Fig. 1 Response characteristics of typical borehole geological structures in electric imaging images. a Stratigraphic interfaces and stratification planes. b Argillaceous stripe. c Fracture. d Vug. e Drilling tool fracture. f Fracturing fracture. g Stress release fracture   (2) Argillaceous stripes: The electric imaging images display bigger thickness continuously and completely, and they are always parallel to the stratigraphic interfaces and stratification planes (Fig. 1b).   (3) Fractures: These usually show as low resistivity strips, shaped like approximate sinusoids with extended length and a lack of completeness and regularity in the image (Fig. 1c).   (4) Vug: The solution pores on the small scale appear as approximate elliptical or speckly, while a cavity on a large scale appears with short extension as an irregularity clot, slice or dark stripe (Fig. 1d).   Note that in the electrical images belonging from a fracture–vug reservoir, irregular fractures always overlap the solution pores or cavities. (5) Drilling tool fracture: These are caused by drilling tools with fine radial extensions. They often occur as in pinnate or echelon arrangements (Fig. 1e).   (6) Fracturing fracture: These always possess small radial extensions, large longitudinal extensions, and appear symmetrical with an approximately 180° delay (Fig. 1f).   (7) Stress release fracture: They often have uniform thickness and regular arrangement (Fig. 1g).   3 Path morphology principle and algorithm in electric imaging logging data Based on mathematical morphology (Serra 1988; Soille and Talbot 2001), the path morphology with a more flexible structuring element group provides an effective detection method for the identification of elongated and curved structures rather than just straight lines (Heijmans et al. 2005). We apply the path morphology method to accomplish the automatic fracture–vug extraction from the electric imaging images by constructing the adjacency relation and determining the length of the path operator. 3.1 Path morphology fundamental principles3.1.1 Adjacency relation Suppose E is a point set representing pixel locations in a two-dimensional (2-D) image domain. A direction, i.e., an adjacency relation, among the pixels in E is denoted by the symbol “\(\mapsto\).” The set E and its adjacency relation compose a directed graph called the adjacency graph (Fig. 3). If \(x \mapsto\, y\) indicates that there exists an asymmetric directed path from x to y, y is called the successor of x and x the predecessor of y. If the adjacency relation “\(\mapsto\)” is known, an arbitrarily subset X in the image domain set E is defined as $$\delta (X) = \{ y \in E|x \mapsto\, y\;{\text{for}}\;{\text{some}}\;x \in X\} ,$$ (1) where \(\delta (X)\) is a successor set satisfying the adjacency relation and comprises all pixel points that have a predecessor in X. Similarly, \(\tilde{\delta }(X)\) is a predecessor set satisfying the adjacency relation and contains all the pixel points that have a successor in X. As shown in Fig. 2, \(b_{1}\), \(b_{2}\), \(b_{3}\) are successors of a, and \(a_{1}\), \(a_{2}\), \(a_{3}\) are predecessors of b, so we may express successor or predecessor set as \(\delta (\{a\}) = \{ b_{1} ,b_{ 2} ,b_{ 3} \}\), \(\tilde{\delta }(\{b\}) = \{ a_{1} ,a_{2} ,a_{3} \}\). Open image in new window Fig. 2 Adjacency relation. Hollow points represent the image domain set E, while solid points represent subsets of E 3.1.2 Path operator For \(\varvec{a} = \{ a_{1} ,a_{ 2} , \ldots ,a_{L} \}\), if \(a_{n} \mapsto a_{n + 1}\), or equivalently \(a_{n + 1} \in \delta (\{ a_{n} \} )\) (n  = 1, 2,…, L − 1), \(\varvec{a} = \{ a_{1} ,a_{ 2} , \ldots ,a_{L} \}\) is defined as a \(\delta\)-path of length L. \(\sigma\)(a) represents elements of the path a existing in E; in other words, \(\sigma\)(a) contains all pixel points which are located on the path a, $$\sigma (a_{1} ,a_{ 2} , \ldots ,a_{L} ) = \{ a_{1} ,a_{ 2} , \ldots ,a_{L} \}$$ (2) The operator \(\alpha_{L} (X)\) is defined as the union of all \(\delta\)-paths of length L in X belonging to image domain set, $$\alpha_{L} (X){ = } \cup \{\sigma (\varvec{a})|\varvec{a}\in \varPi_{L} (X ) \}$$ (3) where \(\varPi_{L}\) represents the \(\delta\)-path set of length L in a subset X. Accessibly, the operator \(\alpha_{L}\) has the algebraic properties of the mathematical morphology opening, including increasing monotonicity, antiextensivity and idempotence (Heijmans et al. 2005), so \(\alpha_{L}\) is defined as the path opening operator, and L is defined as the length of \(\alpha_{L}\). In this paper, in order to extract fractures and vugs from electric imaging images regarded as subset X, we need to look for the union of all pixel points which satisfy \(\delta\)-path of length L (\(\alpha_{L} (X)\)). 3.1.3 Adjacency graph available for electric imaging image characteristics According to the different vertical and horizontal extensions for the different target geological bodies, we need to design the vertical and horizontal adjacency graph available for the electric imaging image characteristics (Fig. 1). The adjacency graph measures their angles, while the operator length L measures their extension and connectivity. According to the different characteristics of borehole geological structures, we designed the vertical and horizontal adjacency graph (Fig. 3). Open image in new window Fig. 3 Adjacency graph available for electric imaging image characteristics. a Vertical adjacency graph. b Horizontal adjacency graph. Black dots represent pixel points in 2-D image domain, and the arrow represents the adjacency relation (1) Vertical adjacency graph   For fractures with large vertical extension, we design a vertical adjacency graph to extract them from the images. We define each pixel of the set E in the 2-D image domain satisfying this adjacency relationship: starting from the chosen pixel to its adjacency pixels in the 0°, 45°, 90° and 135° directions, so the set E and its adjacency relation compose a directed graph, called the vertical adjacency graph (Fig. 3a). We place each pixel of the electric imaging image I in the vertical adjacency graph and perform the path opening. It can keep the high-angle fracture extended from 0° to 135°, the drilling tool fracture and the fracturing fracture mainly in the vertical direction (Fig. 3a). Then, we determine the corresponding path lengths L as threshold values to separate the high-angle fracture from the drilling tool and fracturing fractures. (2) Horizontal adjacency graph   Similarly, for target geologic structures with a large horizontal adjacency extension, we design a horizontal adjacency graph to distinguish them. We define each pixel of the set E satisfying this adjacency relationship: starting from the chosen pixel to its adjacency pixels in the − 45°, 0° and 45° directions, so the set E and its adjacency relation compose a directed graph called the horizontal adjacency graph (Fig. 3b). We place each pixel of the electric imaging image I into the horizontal adjacency graph and perform the path opening. It can keep the low-angle fracture extended from − 45° to 45°, stratigraphic interfaces, stratification planes, argillaceous stripes and stress release fractures in the horizontal direction. Then, we determine the corresponding path lengths L as threshold values to separate low-angle fractures from other fractures. Specifically, the path lengths of the vug are much smaller than of the tectonic or induced fractures, so we determine smaller path lengths L as the threshold value to extract them. Note that the number of the \(\delta\)-paths grows exponentially with the operator length L for the common periodic adjacency graph (Fig. 3). For example, for the horizontal adjacency graph shown in Fig. 3b, except for at the boundary, there exist irreversible \(3^{L - 1}\)\(\delta\)-paths of length L starting from an arbitrary pixel, and the path opening can be understood as the union performing the morphological opening with the \(3^{L - 1}\) structuring elements. However, these structuring elements are not translational to each other. This is the difference between path opening and morphological opening. Obviously, it is inefficient to compute all these operations. We apply an efficient decomposition algorithm (Talbot and Appleton 2007) to realize the morphological path opening algorithm applicable to electric imaging logging as follows. 3.2 Path opening algorithm for fracture–vug extraction Primarily, we need to convert the electric imaging image X into a binary image B. Take the vertical adjacency graph as an example in Fig. 3a, and let B(x) denote the value of the pixel point x on the binary image B and (x1, x2) denote the coordinates of x. We calculate two values at each pixel x: the longest path length \(l^{ - } [x]\) traveling upward from pixel x (not including x itself) and the longest path length \(l^{ + } [x]\) traveling downward from pixel x. If B(x) = 1, the maximum path length through x is \(l[x]= l^{ - } [x]{ + }l^{ + } [x] - 1\). If B(x) = 0, set \(l[x]{ = 0}\). We use the recursive formulas (4) and (5) to compute \(l^{ - } [x]\) and \(l^{ + } [x]\): $$l^{ - } [x] = 1 + \hbox{max} (l^{ - } [(x^{1} + 1,x^{2} )],l^{ - } [(x^{1} - 1,x^{2} + 1)],l^{ - } [(x^{1} ,x^{2} + 1)],l^{ - } [(x^{1} + 1,x^{2} + 1)])$$ (4) $$l^{ + } [x] = 1 + \hbox{max} (l^{ + } [(x^{1} - 1,x^{2} )],l^{ + } [(x^{1} - 1,x^{2} - 1)],l^{ + } [(x^{1} ,x^{2} - 1)],l^{ + } [(x^{1} + 1,x^{2} - 1)])$$ (5) Analogously, for the horizontal adjacency graph in Fig. 3b, according to recursive Eqs. (6) and (7), we compute \(l^{ - } [x]\) and \(l^{ + } [x]\) for each pixel x: $$l^{ - } [x] = 1 + \hbox{max} (l^{ - } [(x^{1} + 1,x^{2} - 1)],l^{ - } [(x^{1} + 1,x^{2} )],l^{ - } [(x^{1} + 1,x^{2} + 1)])$$ (6) $$l^{+} [x] = 1 + \hbox{max} (l^{+} [(x^{1} - 1,x^{2} - 1)],l^{+} [(x^{1} - 1,x^{2} )],l^{+} [(x^{1} - 1,x^{2} + 1)])$$ (7) For details of the recursive computation of \(l^{ - } [x]\) and \(l^{ + } [x] ,\) see Heijmans et al. (2005) and Talbot and Appleton (2007). Searching the path length L of the operator \(\alpha_{L} (X)\) for each borehole geological structure as the appropriate threshold value, the pixels x that satisfy \(l[x] < L\) are removed, and the set comprising all the pixels for \(l[x] \ge L\) is the result of the path opening operation. The image of the maximum path length is generated by substituting the maximum path length l[x] for B(x). The minimum value between two peaks in the frequency distribution of the maximum path length l[x] for all pixels corresponds to the threshold value of the path length L for each geological structure along the borehole. To identify these different structures, we perform the path opening operation to the electric imaging images by their appropriate threshold values L. According to the adjacency graph (Fig. 3a), we implement the opening transform to separate fractures in the electric imaging logging as shown in Fig. 4. In Fig. 4a, the black dots indicate the target pixels set X contained in image domain set E, and the light gray boundary including dots and arrows indicates the adjacency graph. In Fig. 4b, the green, red and blue dots represent isolated points, planar structures and curved structures, respectively. We, respectively, ascertain L = 2 and L = 6 as the threshold values of the path length for separating the three kinds of structures. As shown in Fig. 4b, we compute the path opening \(\alpha_{2} (X)\) of a set X and remove the green isolated points (background noise), while \(\alpha_{6} (X)\) is computed for removing the red planar structures (vugs) and ultimately reserving the blue curved structures (fractures). Consequently, for fracture–vug recognition from electric imaging images, we determine the suitable adjacency graph and the corresponding path length of the operator \(\alpha_{L}\) to suppress the noise and extract various geologic bodies. Open image in new window Fig. 4 Path opening operation for extracting fractures and vugs 4 Noise suppression and fracture–vug extraction from simulated data based on path morphology In this section, we apply the path operator \(\alpha_{L}\) to process the simulated data of the electric imaging images. Placing each pixel point from the electric imaging image into the vertical and horizontal adjacency graph, we determine the appropriate path length L as the threshold value to split the different geological images into subimages with respective maximum path length. The path opening method can achieve noise suppression and automatic fracture–vug extraction and efficiently characterize the fracture–vug reservoirs. 4.1 Quantitative estimation of the fracture–vug parameters in the conductivity images In electric imaging logging, the original conductivity data are acquired from the electrode array of instruments scanning the borehole wall, which is then transferred to a ground system that processes and shows these data as a 2-D conductivity image covering the full borehole. We used the FMI instrument belonging to the MAXIS-500 imaging logging series developed by Schlumberger and designed the simulated formation models comprising of fracture–vug information and calculated their parameters. There are four main plates and four folding plates, including 192 microelectrodes acquiring 192 conductivity curves in the FMI instrument. They provide higher vertical resolution and borehole coverage, i.e., 0.508 cm and 80%, respectively, in a 20.32 cm borehole. After filling the blank strips of the electric imaging logging images, the number of horizontal sampling points along the borehole will reach almost 250. After pretreatment, the conductivity data from 192 plates are calibrated to a static image with the corresponding color distributed in the borehole direction from 0° to 360° (Gaillot et al. 2007). To display fractures and vugs more completely, we need to convert the static image into a dynamic image. We select 2ft as the length of sliding window and 1/2 window length as step length and repeat static color calibration in the window for the whole well. For fracture identification, the image shows a curved interface similar to a sinusoid when the fracture plane intersects the borehole. Combining the geometric relationship between the fracture plane and borehole (Fig. 5), we determine the fracture spatial position along the pixel edge, fit the fracture centerline with a polynomial function and calculate the fracture occurrence including the dip angle and direction. The dip angle α of the fracture obeys the geometric relationship $$\alpha = \text{arctg}(H/d)$$ (8) where H is the vertical distance between the highest and lowest points of the polynomial and d is the borehole diameter. Open image in new window Fig. 5 a Geometric relationship between the fracture planes intersecting the borehole. b Corresponding planar distribution Moreover, θ represents the dip direction, i.e., the azimuth located at the lowest point of the polynomial. For the vug identification, we need to mark vugs and extract their edges from the binary image. A vug is considered an ellipse fitted with the least squares method. Therefore, we can acquire the aspect ratio, area and other parameters by computing the major LENmax and minor LENmin axes of the fitted ellipse (Fig. 6). Open image in new window Fig. 6 Schematic of the solution pore or cavity fitting with the ellipse $${\text{Aspect}}\,{\text{ratio}}\,A_{\text{spect}} : A_{\text{spect}} = {\text{LEN}}_{\hbox{min} } /{\text{LEN}}_{\hbox{max} }$$ (9) $${\text{Area}}\,S:S = \pi {\text{LEN}}_{\hbox{min} } {\text{LEN}}_{\hbox{max} } /4$$ (10) In terms of the aspect ratio, the rock matrix is closer to 1.0, while a fracture with a small aperture and long extension is far less than 1.0. Moreover, the aspect ratio of the vug, whose shape is nearly circular, is between the rock matrix and fracture. For the fine reservoir evaluation, we display the plane porosity, which is the respective ratio of the area of the fracture, vug or rock matrix to the total area of the electric imaging image per unit depth, by a 2-D waveform diagram in ascending order of the aspect ratio. In other words, we call the diagram a fracture–vug pore structure spectrum (Fig. 7). Open image in new window Fig. 7 Schematic of the fracture–vug pore structure spectrum distribution based on pore aspect ratio (the blue dotted line represents the fractures; the black one represents the vugs; the red one represents the matrix pores) To calculate the fracture–vug parameters automatically and performing the opening operation of the binary image based on the path operator from Fig. 3, we can acquire the target geological structures such as fractures and vugs, whose interiors are filled with the morphological region-filling algorithm. After thinning the binary images, we obtain the fracture centerline fitted with a polynomial function for extracted fractures, while we detect their morphological boundaries fitted with elliptic functions for the extracted vugs. 4.2 Noise suppression and fracture–vug extraction based on the model data According to the measured features of the electric imaging logging instrument, we build two formation response models, including the fracture model and fracture–vug model. Figures 8 and 9 are shown as binary images switched from the conductivity response images. Here, the vertical depth is 1 m and the number of horizontal sampling points distributed from 0° to 360° is 250. It is well known that the real conductivity images always comprise areas with different degrees of noise, mainly caused by the electronic parts, circuits, collision between instrument and borehole, and others, so we add 10% salt and pepper noise to the conductivity models for the binary images. Open image in new window Fig. 8 a Fracture formation model with 10% salt and pepper noise, including one high-angle fracture, three low-angle fractures and background noise. The ordinate represents the vertical relative depth (m), while the abscissa represents the azimuth. b Separated noise from the fracture formation model. c Frequency distribution of the maximum path length l(x) for each pixel x (the red line represents the threshold value between noise and fracture). d Low-angle fractures extraction. e Centerlines of the three low-angle fractures. f Polynomial fit of the centerlines of low-angle fractures. g Low-angle fracture occurrence (circumferential and radical direction represents dip direction and dip angle, respectively). h High-angle fractures extraction. i Centerlines of the high-angle fractures. j Polynomial fit of the centerline of the high-angle fracture. k High-angle fracture occurrence Open image in new window Fig. 9 a Fracture–vug formation model with 10% salt and pepper noise including four fractures and seven pores. The ordinate represents the vertical relative depth (m), whereas the abscissa represents the azimuth. b Separated noise from the fracture–vug formation model. c Frequency distribution of the maximum path length l(x) for each pixel x (the red, purple, blue and orange dotted line represent the threshold values dividing into noise, small pores, medium pores, large pores and fracture, respectively). d Pore extraction. e Pore edge detection. f Ellipse fit of pore edges. g Fractures extraction. h Fracture centerlines. i Polynomial fits of the fracture centerlines 4.2.1 Noise suppression and fracture extraction for the fracture formation model We design a conductivity response model of the fracture formation composed of one high-angle fracture, three low-angle fractures and background noise; particularly, the high-angle fracture and two low-angle fractures are parallel to each other (Fig. 8a). The extension length of the fracture connected horizontally is much greater than of the isolated noise. As shown in Fig. 8c, implementing Fig. 3b as an adjacency graph for opening operation, the results indicate obviously that the path length of the fracture (L > 8) is greater than the background noise (0 ≤ L ≤ 8). Using the maximum path length of the background noise as the threshold value, we can remove the noise (Fig. 8b) and reserve fractures at L > 8. To perform the opening operation to the process denoising binary image, the low-angle fractures are extracted using the horizontal adjacency graph, whereas the high-angle fractures are extracted using the vertical adjacency graph (Fig. 3). To acquire fracture parameters, we pick up the fracture centerlines fitted with polynomials. The result shows that the path morphology method not only effectively separates the high and low-angle fractures, but also commendably keeps the connectivity of each fracture (Fig. 8). 4.2.2 Fracture and pore separation for the fracture–vug formation model The fracture–vug formation model comprises four fractures connected horizontally at various angles, especially the two low-angle fractures which cut each other. We assume the shape of the simulated pores on different scales is circular, and their diameters are, respectively, the length of 12, 20 and 35 sampling intervals (Fig. 9a). It is well known that the extension length of the fracture is much greater than of the pores. To separate the fractures and vugs at different scales, we perform the horizontal adjacency graph for the opening operation and determine a suitable path length of the operator \(\alpha_{L}\) as the corresponding threshold value for each geological body. Obviously, the path length of these pores is between fractures and background noise. According to Fig. 9c, there are four peak ranges on frequency distribution of l(x) and the minimum value between two peaks as the threshold value of the path length L for vugs on different scales. Taking the path opening length L = 8 as the threshold value, those pixels belonging to noise (where l(x) < 8) are removed after the path opening operation \(\alpha_{8} (X)\). Similarly we perform \(\alpha_{13} (X)\) to the denoising image by threshold value (L = 13), and these pixels belonging to small-scale pores (where 8 ≤ l(x) < 13) are filtered out. Then, using the path lengths L = 26 and L = 40 as the threshold values, executing \(\alpha_{ 2 6} (X)\) and \(\alpha_{ 4 0} (X)\) to the previous processed image in turn, we can successively remove medium-scale pores (where 13 ≤ l(x) < 26) and large-scale pores (where 26 ≤ l(x) < 40) and ultimately isolate fractures at l(x) ≥ 40. The result shows clearly that the path morphology method effectively extracts fractures and pores from the background noise (Fig. 9). Particularly, for extracted pores, we detect their morphological boundaries fitted with an elliptic function, whereas we acquire the fracture centerlines fitted with polynomial functions for the extracted fracture. The specific parameters of fracture and pore are shown in Tables 1 and 2, respectively. Table 1 Pore parameters of fracture–vug formation model Pore (scale) Aspect ratio Size, cm2 Small 0.90 13.2 Medium 0.95 38.4 Large 0.98 125.4 Table 2 Fracture occurrence of fracture–vug formation model Fracture Dip angle, ° Dip direction, ° Blue line 48 88 Red line 35 85 Black line 23 268 Green line 48 88 5 Filling blank strips and matrix separation in electric imaging images Before applying the path morphology to identify fractures and vugs, we need to deal with the original electric imaging images, including the singular spectral interpolation (Recht 2011; Li et al. 2017a, b, c; Cai et al. 2008; Li et al. 2017a, b, c) for the blank strips and the Otsu method (Otsu 1979; Mala and Sridevi 2016) for matrix separation. 5.1 Singular spectrum analysis and interpolation reconstruction on conductivity response During the electric imaging logging measurement, images appear as blank strips owing to the loss of conductivity data in the electric imaging logging data because the electrode plates do not cover the full borehole (Fig. 10). In this section, we reconstruct the full borehole 2-D electric imaging images with the singular spectral analysis method. Open image in new window Fig. 10 Schematic of the electric imaging logging data matrix (“*” represents practical electric imaging logging data and “○” missing data) (1) Convert the electric imaging images to 2-D data in the frequency domain by a Fourier transform.   (2) Implement the singular spectrum interpolation for every frequency slice of the 2-D data.   For each frequency slice \(S = \left\{ {s_{1} ,s_{2} , \ldots ,s_{N} } \right\}\), construct the Hankel matrix: $$M = \left[ {ll_{1} ,ll_{2} , \ldots ,ll_{K} } \right] = \left[ {\begin{array}{*{20}c} {s_{1} } & {\begin{array}{*{20}c} {s_{2} } & {s_{3} } \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{\lambda } } \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{K} } \\ \end{array} } \\ {s_{2} } & {\begin{array}{*{20}c} {s_{3} } & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {s_{\lambda } } & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {s_{K} } & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } \\ {\begin{array}{*{20}c} {s_{3} } \\ {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{\lambda } } \\ {s_{\lambda } } & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{K} } \\ {s_{K} } & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } \\ {s_{\lambda } } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{K} } \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} \\ \end{array} } & {\begin{array}{*{20}c} {\mathinner{\mkern2mu\raise1pt\hbox{.}\mkern2mu \raise4pt\hbox{.}\mkern2mu\raise7pt\hbox{.}\mkern1mu}} & {s_{N} } \\ \end{array} } \\ \end{array} } \right]$$ (11) and $$ll_{i} = \left[ {s_{i} ,s_{i + 1} , \ldots ,s_{i + \lambda - 1} } \right]^{\text{T}} ,{1 \le \lambda \le N,1 \le i \le K,}$$ (12) where \(ll_{i}\) denotes the delay vectors for the 1-D conductivity signal in the frequency domain, λ the length of the delay vector and K the number of delay vectors. Performing singular value decomposition for the Hankel matrix, we construct the reduced-rank Hankel matrix after cutting off small singular values. $$M^{ * } = U\left[ {\Sigma_{j} } \right]V^{H},\, j \le r$$ (13) where \(\Sigma_{j}\) are the first j larger singular values, r denotes the rank of the original \(\lambda \times K\) Hankel matrix and U and V denote the λ unitary matrix and K unitary matrix satisfying the singular value decomposition of the Hankel matrix, respectively. By performing the inverse Hankel transform for the low-rank Hankel matrix \(M^{ * }\), the elements in the first column and the last row of the new matrix compose the reconstructed signal. (3) To reconstruct the 2-D electric imaging logging data in the spatial domain, we perform an inverse Fourier transformation of the 2-D conductivity signal in the frequency domain.   Figure 11b shows the singular spectrum interpolation that is used to fill the blank strips in Fig. 11a. The result indicates that the filled image makes up for the information that the instrument cannot detect, especially formations with strong heterogeneity. The interpolation method can not only suppress noise and improve image quality, but also reconstruct complete electric imaging logging images. Open image in new window Fig. 11 a Original electric imaging logging image. b Electric imaging image interpolated by singular spectrum interpolation. c Fracture–vug image separated by the Otsu method 5.2 Automatic matrix segmentation from electric imaging images based on the Otsu method In this section, we apply the Otsu method to automatically separate fractures and vugs from the rock matrix. The fundamental principle of the Otsu method is that we determine the optimal threshold to divide the gray level histogram of the images into two classes whose variance is the largest. As shown in Fig. 12, we define the image F whose size is Q × R pixels, or equivalently, \(F = \{ f(x,y)|x \in \{ 1,2, \ldots ,Q\} ,y \in \{ 1,2, \ldots ,R\} \}\). Let \(f(x,y)\) be the gray value located in \((x,y)\). We test k at the gray level range and determine k as the optimal segmentation threshold when between-class variance \(\sigma^{2} (k)\) is the maximum. Open image in new window Fig. 12 Schematic diagram of automatic segmentation based on the Otsu method Figure 11c shows how the Otsu method is utilized to automatically segment Fig. 11b. Comparing Fig. 11c, b, we see that the Otsu method is an efficient way to separate fractures and vugs from the filled electric imaging image. 6 Analysis and processing of practical logging data We use the novel path morphology method to extract fractures and vugs from the fracture–vug reservoir in the wells A, B and C (as shown in Figs. 13, 14, 15, respectively). Wells A and C are carbonate formations, whereas well B is a clastic formation. To implement the automatic fracture–vug extraction and identification based on the opening transformation, after filling the blank strips of original electric imaging logging images and applying conventional logging curves for shale correction, we need to apply the Otsu method for the threshold segmentation that targets geological structures such as fractures and vugs, which are separated from the matrix background in the real conductivity images. Open image in new window Fig. 13 Fracture extraction and occurrence characterizations of the fracture reservoir in well A. a Original electric imaging image. b Binary electric imaging image after removing the matrix. c Stratification and background noise extracted from b. d Fracture extracted from b. e Polynomial fit of the fracture centerlines. f Fracture occurrence Open image in new window Fig. 14 Pore extraction of the pore reservoir and comparison with the thin section in well A. a Original electric imaging image. b Vugs extracted from a by path opening operation. c Edge detection of b. d Ellipse fit of c. e the thin section located in 2704 m. f Solution pores extracted from e Open image in new window Fig. 15 Fracture–vug extraction and comparison with the multi-threshold segmentation in well B. a Original electric imaging image. b Fracture extraction based on multi-threshold segmentation. c Vug extraction based on multi-threshold segmentation. d Threshold selection between fractures and vugs. e Fracture extraction based on path opening operation. f Vug extraction based on path opening operation The original electric imaging images show how the high-angle fractures cutting the stratification are extremely developed and complexly distributed, so we determine the formation as the fracture reservoir (in Fig. 13). According to the vertical and horizontal adjacency graph, as shown in Fig. 3, the three high-angle fractures are completely extracted from the noise and stratification, so the path opening operation has the advantage of being efficient in detecting fractures in all directions, especially for groups of fractures with a large curvature. To acquire the fracture parameters, we pick the fracture centerlines fitted with a polynomial. Figure 14 shows the result of vug extraction from pores by a path opening operation in contrast to the result of thin section analysis. Figure 14e shows that intercrystalline pores are mainly developed in this reservoir. Based on the path opening operation, we adequately separate the vugs from the noise background and accurately detect their edges (Fig. 14c, d). Comparing the plane porosity of vugs extracted from Fig. 14b by path opening operation with those from thin section (Fig. 14f), the plane porosity in Fig. 14b, e is 7.5% and 6.8%, respectively, with a relative error of less than 10%. The result shows that it is reasonable to use the path opening operation for the segmentation of the electric imaging images. We perform the path opening operation to process the electric imaging image from the fracture–vug reservoir and compare the extraction result with that of multi-threshold segmentation (Fig. 15). We can distinguish the fracture and vug subimages extracted by multi-thresholding (Ge et al. 2015), according to the fracture roundness > 3.5 and the length–width ratio > 1.5 (Fig. 15b, c). We determine the minimum value between the first peak and the second peak in Fig. 15d as the threshold value of the path length L, and the extracted fractures and vugs are shown in Fig. 15e, f based on the path opening operation. As shown in Fig. 15, both the path opening operation and the multi-threshold segmentation have achieved good segmentation results and verify the validity of the path opening algorithm for separating fractures and vugs. In particular, with respect to the multi-threshold segmentation method, the path opening operation can better keep the integrity of fractures while extracting vugs. We use the path opening algorithm to deal with the electric imaging image from the fracture–vug reservoir of the well B in the clastic formation and compare the result with the Logview software analysis (Fig. 16). The whole-bore electric imaging image shows that the formation develops horizontal stratification and vugs, but high conductivity fractures. In contrast to the results of the Logview software analysis, the vugs surrounded by the pocks automatically picked up by the software are consistent with the vugs separated from the horizontal stratification by the path opening operation, and the plane porosity calculated by the two methods is coherent, which prove the reliability of the path opening algorithm in vug extraction. Open image in new window Fig. 16 Comprehensive logging interpretation of the vug reservoir in a clastic formation, well B. The first column: depth; the second column: conventional logging curves; the third column: the whole-bore electric imaging image; the fourth column: matrix extraction; the fifth column: the horizontal stratification extraction; the sixth column: vug extraction; the seventh column: the vug subimage with pocks picked up automatically by the Logview software; the eight column: comparison of the plane porosity between the two methods (the red line represents the extraction results based on the path opening operation, and the black represents that based on the Logview software) For oil-bearing strata in a carbonate formation (well C), we compare the interpretations of conventional logging and core data by analyzing and evaluating the fracture–vug reservoir after performing the path opening operation (Fig. 17). Because the reservoir develops secondary pores and is relatively high in uranium, it causes the natural gamma radiation to be slightly higher than in the adjacent layer. The effective porosity calculated by conventional logging interpretation is 2.3%–4.4%. The whole-bore electric imaging logging image and core scanning image indicate that fractures develop very well with vugs. The path opening operation is applied to detect fractures in all directions and to separate the fractures from the vugs. On the fracture–vug pore structure spectrum, the fractures with a small aspect ratio show a high amplitude peak on the left; the high amplitude peaks in the middle indicate that ellipsoid vugs are developed in the formation; the peak area on the right, whose aspect ratio is almost 1, represents the matrix proportion. Given the threshold value of the corresponding aspect ratio for fracture, vug and matrix, their porosity percentage can be calculated by integrating the area within the threshold. Calibrating porosity with the logging data in the pure matrix layer, we compare the cumulative porosity curve with the core test porosity (blue points), where the average relative error was 5.0%, which also indicated the effectiveness of the path opening operator in extracting fractures and vugs. As shown in Figs. 16 and 17, the path opening algorithm can achieve good results in automatically extracting fractures and vugs whether in carbonate rocks or clastic rocks. Open image in new window Fig. 17 Logging comprehensive interpretation of the fracture–vug reservoir in carbonate formation, well C. The first column: depth; the second column: conventional logging curves; the third column is the fullbore electric imaging image; the fourth column: matrix extraction; the fifth column: the fracture extraction; the sixth column: vug extraction; the seventh column: the fracture–vug pore structure spectrum; the eighth column: porosity accumulation curve (blue points represent core test porosity); the ninth column is the core scanning image 7 Conclusion On the basis of the fracture–vug response characteristics of the electric imaging logging data, we have developed a novel fracture–vug identification and extraction method based on the path morphology. To extract multi-scale fractures and vugs from the electric imaging logging images automatically, we have designed an adjacency graph available for linear and curved structures such as different angle fractures and determined a suitable path length of the operator \(\alpha_{L}\) as the corresponding threshold value for each geological feature. We have successfully applied the new method to deal with simulated data and practical electric imaging logging data and draw the following conclusions: (1) The path morphology method provides an efficient way to identify fractures and vugs, and we calculate their parameters automatically by path operators tracing the linear and curved structures available for fracture distribution with different linear structures.   (2) The new method has the advantage of being precise enough to distinguish fractures at different angles by applying the vertical adjacency graph for high-angle fractures and horizontal adjacency graph for low-angle fractures.   (3) The path morphology method can be implemented to effectively suppress random noise on the electric imaging logging images and automatically identify fractures and vugs. Furthermore, the fracture–vug pore structure spectrum denoted by the pore aspect ratio not only describes the distribution of multi-scale pore types, but also performs better in dividing reservoir types and quantitative evaluation of the pore structure.   In this paper, the new fracture–vug extraction method based on the path morphology is still at the exploratory stage. At present, the path length of operator is chosen at the minimum point between two peaks belonging to fractures and vugs in the frequency distribution of the maximum path length. However, it is difficult to choose the point in practical applications, so we need to explore searching automatically for threshold values. Furthermore, we need to study deeply the effect of different adjacency relations and path opening operators on the extraction of argillaceous stripes, stratification and high resistivity fractures and thus summarize the extraction templates and applicable conditions for different borehole geological structures. Notes Acknowledgements This study was granted access to projects supported by the National Major Fundamental Research Program of China “On basic research problems in applied geophysics for deep oil and gas fields” (Grant Number 2013CB228605), CNPC Science and Technology Project (Grant Number 2016A-3303) and CNPC Logging Project (Grant Number 2017E-15). References Amartey EO, Akiti TT, Armah T, et al. Integrating gamma log and conventional electrical logs to improve identification of fracture zones in hard rocks for hydrofracturing: a case study from Ghana. Appl Water Sci. 2017;7(3):1091–8.  https://doi.org/10.1007/s13201-016-0450-z.CrossRefGoogle Scholar Bian HL, Ju G, Mao ZQ, et al. Pore structure effect on reservoir electrical properties and well logging evaluation. Appl Geophys. 2014;11(4):374–83.  https://doi.org/10.1007/s11770-014-0462-0.CrossRefGoogle Scholar Cai JF, Candes EJ, Shen ZW. A singular value thresholding algorithm for matrix completion export. Soc Ind Appl Math. 2008;20(4):1956–82.  https://doi.org/10.1137/080738970.CrossRefGoogle Scholar Chen H, Heidari Z. Pore-scale evaluation of dielectric measurements in formations with complex pore and grain structures. In: 55th SPWLA annual logging symposium. 2014; Abu Dhabi, Paper NNNN.Google Scholar Chen XL, Tang XM, Qian YP. Propagation characteristics of multipole acoustic logging in cracked porous tight formations. Chin J Geophys. 2014;57(9):2961–70.  https://doi.org/10.6038/cjg20140921 (in Chinese).CrossRefGoogle Scholar Chitale DV, Quirein J, Perkins T, et al. Application of new borehole imager and technique to character secondary porosity and net-to-gross in vugular and fractured carbonate reservoirs in Permian basin. In: 45th SPWLA annual logging symposium, Netherlands. Paper RR (2004)Google Scholar Cornet J. Fracture detection and analysis from image log raw data. Norway: Norwegian University of Science and Technology; 2013.Google Scholar Cui WP, Yang YQ, Zhang X. Calculation of visual width of fracture in electric imaging logging. Well Logging Technol. 2016;40(5):633–6.  https://doi.org/10.16489/j.issn.1004-1338.2016.05.020 (in Chinese).CrossRefGoogle Scholar Delhomme JP. A quantitative characterization of formation heterogeneities based on borehole image analysis. In: 33th SPWLA annual logging symposium. 1992; Oklahoma, Paper T.Google Scholar Gaillot P, Brewer T, Pezard P, et al. Borehole imaging tools—principles and applications. Sci Drill. 2007;5:1–4.  https://doi.org/10.2204/iodp.sd.5.07S1.2007.CrossRefGoogle Scholar Gao JS, Wang DS, Liu JY, et al. Theory study on fracture–vug reservoirs resistivity based on dual and triple porosity models. Well Logging Technol. 2015;39(3):277–82.  https://doi.org/10.16489/j.issn.1004-1338.2015.03.003 (in Chinese).CrossRefGoogle Scholar Ge XM, Fan YR, Li JT, et al. Pore structure characterization and classification using multifractal theory—an application in Santanghu basin of western China. J Pet Sci Eng. 2015;127:297–304.  https://doi.org/10.1016/j.petrol.2015.01.004.CrossRefGoogle Scholar Gunter GW, Spain DR, Viro EJ, et al. Winland pore throat prediction method-a proper retrospect: New examples from carbonates and complex systems. 55th SPWLA Annual Logging Symposium. 2014; Abu Dhabi, Paper, KKK.Google Scholar Hall J, Ponzi M, Gonfalini M, et al. Automatic extraction and characterization of geological features and textures from borehole images and core photographs. In: 37th SPWLA annual logging symposium. 1996; Paper, CCC.Google Scholar Heijmans H, Buckley M, Talbot H. Path openings and closings. J Math Imaging Vis. 2005;22(2):107–19.  https://doi.org/10.1007/s10851-005-4885-3.CrossRefGoogle Scholar Holden A, Lehmann C, Ryder K, et al. Integration of production logs helps to understand heterogeneity of Mishrif reservoir in Rumaila. In: 55th SPWLA annual logging symposium. 2014; Abu Dhabi, Paper, GGG.Google Scholar Lai FQ. Micro-resistivity image logging processing and interpretation methods research. Qingdao: China University of Petroleum (East China); 2011 (in Chinese).Google Scholar Li C, Huang JP, Li ZC, et al. Regularized least-squares migration of simultaneous-source seismic data with adaptive singular spectrum analysis. Pet Sci. 2017a;14(1):61–74.  https://doi.org/10.1007/s12182-016-0134-1.CrossRefGoogle Scholar Li P, Liu ZC, Yang N, et al. Application of fine target reservoir imaging in fractured–vuggy carbonate reservoir, Yubei area. Geophys Prospect Pet. 2015;54(4):443–51.  https://doi.org/10.3969/j.issn.1000-1441.2015.04.011 (in Chinese).CrossRefGoogle Scholar Li XN, Shen JS, Li ZL, et al. Characterization of reservoir fracture and vug parameters by conductivity image of FMI based on multi-scale mathematical morphology method. J China Univ Pet (Edition of Natural Sciences). 2017b;41(1):69–77.  https://doi.org/10.3969/j.issn.1673-5005.2017.01.008 (in Chinese).CrossRefGoogle Scholar Li ZL, Shen JS, Li S, et al. Singular spectral interpolation of blank strips in formation micro-scanner conductivity image and separation of porosity between fracture and karst cave. Well Logging Technol. 2017c;41(1):33–40.  https://doi.org/10.16489/j.issn.1004-1338.2017.01.006 (in Chinese).CrossRefGoogle Scholar Liu RL, Xie F, Xiao CW, et al. Extracting fracture–vug plane porosity from electrical imaging logging data using dissection of wavelet-transform-based image. Chin J Geophys. 2017;60(12):4945–55.  https://doi.org/10.6038/cjg20171233 (in Chinese).CrossRefGoogle Scholar Lyu WY, Zeng LB, Liu ZQ, et al. Fracture responses of conventional logs in tight-oil sandstones: a case study of the Upper Triassic Yanchang Formation in southwest Ordos Basin, China. AAPG Bull. 2016;100(9):1399–417.  https://doi.org/10.1306/04041615129.CrossRefGoogle Scholar Mala C, Sridevi M. Multilevel threshold selection for image segmentation using soft computing techniques. Soft Comput. 2016;20(5):1793–810.  https://doi.org/10.1007/s00500-015-1677-6.CrossRefGoogle Scholar Otsu N. A threshold selection method from gray-level histogram. Syst Man Cybern IEEE Trans. 1979;9(1):62–6.CrossRefGoogle Scholar Pan BZ, Yuan MX, Fang CH, et al. Experiments on acoustic measurement of fractured rocks and application of acoustic logging data to evaluation of fractures. Pet Sci. 2017;3:520–8.CrossRefGoogle Scholar Recht B. A simpler approach to matrix completion. J Mach Learn Res. 2011;12:3413–30.Google Scholar Serra J. Image analysis and mathematical morphology. New York: Academic Press; 1988.Google Scholar Shalaby MR, Islam MA. Fracture detection using conventional well logging in carbonate Matulla Formation, Geisum oil field, southern Gulf of Suez, Egypt. J Pet Explor Prod Technol. 2017;3:1–13.  https://doi.org/10.1007/s13202-017-0343-1.CrossRefGoogle Scholar Soille P, Talbot H. Directional morphological filtering. IEEE Trans Pattern Anal Mach Intell. 2001;23(11):1313–29.  https://doi.org/10.1109/34.969120.CrossRefGoogle Scholar Talbot H, Appleton B. Efficient complete and incomplete path openings and closings. Image Vis Comput. 2007;25(4):416–25.  https://doi.org/10.1016/j.imavis.2006.07.021.CrossRefGoogle Scholar Tang JW. Image logging data processing and fracture identification methods research. Xian: Xi’an University of Science and Technology; 2013 (in Chinese).Google Scholar Tang XM, Qian YP, Chen XL. Laboratory study of elastic wave theory for a cracked porous medium using ultrasonic velocity data of rock samples. Chin J Geophys. 2013;56(12):4226–33.  https://doi.org/10.6038/cjg20131225 (in Chinese).CrossRefGoogle Scholar Thachaparambil MV, Wu RC, Wang D, et al. Automated extraction of fracture network from 3D poststack seismic data for fracture and fault characterization and modeling. Lead Edge. 2013;32(4):380–4.  https://doi.org/10.1190/tle32040380.1.CrossRefGoogle Scholar Tokhmechi B, Memarian H, Rasouli V, et al. Fracture detection from water saturation log data using a Fourier–wavelet approach. J Pet Sci Eng. 2009;69(1–2):129–38.  https://doi.org/10.1016/j.petrol.2009.08.005.CrossRefGoogle Scholar Wu G, Yang H, Li H, et al. Controlling effects of the Ordovician carbonate pore structure on hydrocarbon reservoirs in the Tarim Basin, China. Pet Sci. 2013;10(3):282–91.  https://doi.org/10.1007/s12182-013-0277-2.CrossRefGoogle Scholar Wang H, Sun SZ, Yang H, et al. The influence of pore structure on P- & S-wave velocities in complex carbonate reservoirs with secondary storage space. Pet Sci. 2011;8(4):394–405.  https://doi.org/10.1007/s12182-011-0157-6.CrossRefGoogle Scholar Wang R, Wang Z, Shan X, et al. Factors influencing pore-pressure prediction in complex carbonates based on effective medium theory. Pet Sci. 2013;10(4):494–9.  https://doi.org/10.1007/s12182-013-0300-7.CrossRefGoogle Scholar Wang S. Study on the extraction methods of high resolution remote sensing images. Nanjing: Nanjing University of Science and Technology; 2014 (in Chinese).Google Scholar Xavier A, Guerra CE, Andrade A. Fracture analysis in borehole acoustic images using mathematical morphology. J Geophys Eng. 2015;12(3):492–501.  https://doi.org/10.1088/1742-2132/12/3/492.CrossRefGoogle Scholar Xiao XL, Jin XJ, Zhang X, et al. Fracture identification based on information fusion of conventional logging and electrical imaging logging. Oil Geophys Prospect. 2015;50(3):542–7.  https://doi.org/10.13810/j.cnki.issn.1000-7210.2015.03.023 (in Chinese).CrossRefGoogle Scholar Xie F, Xiao CW, Liu RL, et al. Multi-threshold de-noising of electrical imaging logging data based on the wavelet packet transform. J Geophys Eng. 2017;14(4):900–8.  https://doi.org/10.1088/1742-2140/aa6ad3.CrossRefGoogle Scholar Xu S, Su YD, Chen XL, et al. Numerical study on the characteristics of multipole acoustic logging while drilling in cracked porous medium. Chin J Geophys. 2014;57(6):1999–2012.  https://doi.org/10.6038/cjg20140630 (in Chinese).CrossRefGoogle Scholar Yang H, Sun SZ, Cai L, et al. A new method of formation evaluation for fractured and caved carbonate reservoirs: a case study from the Lundong area, Tarim Basin, China. Pet Sci. 2011;8(4):446–54.  https://doi.org/10.1007/s12182-011-0162-9.CrossRefGoogle Scholar Zazoun RS. Fracture density estimation from core and conventional well logs data using artificial neural networks: the Cambro-Ordovician reservoir of Mesdar oil field, Algeria. J Afr Earth Sci. 2013;83(1):55–73.  https://doi.org/10.1016/j.jafrearsci.2013.03.003.CrossRefGoogle Scholar Zhang G, Li N, Guo HW, et al. Fracture identification based on remote detection acoustic reflection logging. Appl Geophys. 2015;12(4):473–81.  https://doi.org/10.1007/s11770-015-0522-0.CrossRefGoogle Scholar Copyright information © The Author(s) 2018 Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Authors and Affiliations Xi-Ning Li12Jin-Song Shen12Email authorWu-Yang Yang3Zhen-Ling Li41.State Key Laboratory of Petroleum Resources and ProspectingChina University of PetroleumBeijingChina2.CNPC Key Laboratory of Geophysical ProspectingChina University of PetroleumBeijingChina3.North-west Division of Exploration and Production Research InstituteCNPCLanzhouChina4.Huabei Division of China Petroleum Logging CorporationXi’anChina


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs12182-018-0282-6.pdf

Xi-Ning Li, Jin-Song Shen, Wu-Yang Yang, Zhen-Ling Li. Automatic fracture–vug identification and extraction from electric imaging logging data based on path morphology, Petroleum Science, 2018, 58-76, DOI: 10.1007/s12182-018-0282-6