Caving Depth Classification by Feature Extraction in Cuttings Images

The estimation of caving depth is of particular interest in the oil industry. During the drilling process, the rock classification problem is studied to analyze the concentration of cuttings at the vibrating shale shakers through the classification of caving images. To date, depth estimation based on caving rock images has not been treated in the literature. This paper presents a new depth caving estimation system based on the classification of caving images through feature extraction. To extract the texture descriptors, the cutting images are first mapped on a common space where they can be easily compared. Then, textural features are obtained by applying a multi-scale and multi-orientation approach through the use of Gabor transformations. Two different depth classifiers are developed; the first separates the textural features by using a soft decision based on the Euclidean distance, and the second performs a hard decision classification by applying a thresholding procedure. A detailed mathematical formulation of the developed classifiers is presented. The developed estimation system is verified using real data from rock cutting images in petroleum wells. Several simulations illustrate the performance of the proposed model using real images from a wellbore in a Colombian basin. The correct classification rate of a database containing 17 depth estimates is 91.2%. Resumen La estimacion de la profundidad de la que provienen los derrumbes que usualmente se presentan en las caras del pozo o tambien llamados cavings es de gran interes en la industria petrolera. Durante el proceso de perforacion de un pozo, el problema de clasificacion de rocas ha sido estudiado con el fin de analizar la concentracion de recortes o ripios de perforacion en las zarandas vibratorias a traves de la clasificacion de imagenes de cavings. Sin embargo, la estimacion de la profundidad de los derrumbes basada en la utilizacion de imagenes de los mismos no ha sido tratada en la literatura. Este articulo presenta un nuevo modelo para la estimacion de la profundidad de derrumbes a traves de extraccion de caracteristicas. Para la extraccion de estas caracteristicas o descriptores de textura, imagenes de recortes son transformadas en un espacio comun, el cual permite su comparacion. Luego, las caracteristicas se obtienen aplicando la transformacion de Gabor, un enfoque que se caracteriza por proporcionar un analisis multi-escala y multi-orientacion. Se desarrollaron dos clasificadores, el primero separa las caracteristicas de textura usando un enfoque basado en la norma Euclideana y el segundo basado en decisiones por umbral. La formulacion matematica detallada de los clasificadores desarrollados se presenta en este articulo. El sistema de estimacion desarrollado se evalua usando datos reales de imagenes de derrumbes pertenecientes a un pozo petrolero. Simulaciones muestran el rendimiento del modelo propuesto usando imagenes reales de un pozo perteneciente a una cuenca Colombiana. La correcta clasificacion para una base de datos de imagenes que contiene 17 clases o profundidades es de 91.2%.


Introduction
A crucial problem in the oil industry is the collapse of the wellbore during the drilling process.The collapse of the wellbore occurs when pieces of rock fall from the walls of the borehole, which is usually known as caving (Aldred, et al., 1999).The costs associated with the drilling process are in the millions, and between 15% and 30% of the total resources dedicated to an oil project are designated to cover losses, including material and drilling equipment and continuity of the drilling process; these losses are known as non-productive time (Aldred, et al., 1999).Hence, classifying the rocks associated with the collapse of the wellbore is an important area of research in this industry.
The problem of rock classification during the drilling process is investigated by analyzing the cutting concentration at the vibrating shale shakers by classifying caving images (Marana, et al., 2009).The detection of a high caving concentration indicates problems such as the collapse of the well borehole walls (Castillo & Moos, 2000).The cutting concentration has been addressed using several supervised classifiers, such as the Optimum-Path Forest (OPF), Artificial Neural Networks with Multi Layer Perceptron (ANN-MLP), Support Vector Machines (SVM) and Bayesian Classifiers (BC) (Marana, et al., 2009).Artificial neural networks and SVM are the most used and result in up to 93.5% correct classification when the rocks texture is analyzed (Singh, Singh, Tiwary, & Sarkar, 2010), (Marmo, Amodio, Tagliaferri, Ferreri, & Longo, 2005).The OPF technique offers a lower computational complexity than other methods (Marana, et al., 2009).The Bayesian Classifier is the most traditional classification method, and its principal advantages include its simplicity, better performance when more data are available, and self-correction, which means that the result changes for different data (Islam, Wu, Ahmadi, & Sid-Ahmed).However, those classification methods have only a few well-established classes.Despite the diverse research into cutting concentration, research on depth estimation using caving images is lacking in the literature (Galvis, Arguello, & Tarazona, 2011).
This paper presents a new system for estimating caving depth based on features extracted from rock images.Given a set of textural rock features, the estimation system divides the caving depth into a set of depth ranges previously established as classes in the classification system.The depth classifier aims at separating the textural features by using a soft decision based on the Euclidean distance and a hard decision by using a threshold constant.The developed estimation system is verified using real data from rock cutting images in petroleum wells.Additionally, the system for estimating the depth of caving can be extended to other areas, such as mining, where problems associated with caving in a specific region of the wellbore can be avoided or addressed.This is the first work applying a classification system to estimate the depth of the caving.The principal contributions of this paper include the establishment of an extensive mathematical model for the rock classification problem and the verification of the developed model using images from real caving rocks.
An important set of caving rocks, called cuttings, are rocks whose depth origin is previously known.Each cutting rock image is considered a mosaic of different textural regions.To extract these textural regions, the cutting images are first mapped onto a common space where they can be easily compared.This mapping reduces the dimensionality of the images by retaining only the relevant information.A two-dimensional (2D) Gabor transformation is used to map the rock images as this transform emphasizes the textural features.Texture plays a particularly important role in the composition of natural images, and its analysis and classification are areas of continued research (Kachanuban & Udomhunsakul, 2007).
Depth classification is a multi-step process.In the first stage, a large database of rock images is acquired from a real wellbore in a Colombian basin.The sensed gray-scale images are divided into training and testing sets.This function has been recognized as a very useful tool in computer vision and image processing, especially for texture analysis (Shen, Jia, & Chen, 2011) (Clausi & Jernigan, 2000).The Gabor transformation provides a multi-scale and multi-orientation representation of the underlying images (Clausi & Jernigan, 2000).The features are extracted from the Gabor transformations and through the co-occurrence matrix (Marana, et al., 2009) (Tou, Tay, & Lau, 2009).The set of features is represented as a matrix of features F C for the training set and F X for the test set.Specifically, the κ-extracted features correspond to the mean µ, the standard deviation σ, the contrast α, the homogeneity ß, the energy ɣ, and the correlation ρ. each image in a set is divided in subimages, and then, the Gabor transformations are obtained.
In the classification stage, a multiclass approach is applied.In this approach, a number of binary classifiers is developed with the set of image features in the training set F C .Then, given a caving image from the test set, its set of features F X is used to test the classifier.To make the decision about the class or depth to which the caving belongs, the classifier uses two methods: one method based on rules that make the decision using a hard thresholding approach and another method based on the ℓ 2 -norm.The number of binary classifiers depends on the number of classes or, in this case, the number of known depths D. The mathematical model for classification and the approaches presented here correctly classify depths up to 91% of the time.
The paper is organized as follows.Section 2 describes the data and process methodology, including the division of the data, the feature extraction step, the Gabor transformation used, and the classification.Section 3 presents the experimental setup.Section 4 presents the results of the simulations using data from a real wellbore in a Colombian basin.Let c dj and xd ℓj be the vector representations of C dj and Xd ℓj , such that c dj = Vec(C dj ) and xd ℓj = Vec(Xd ℓj ).More specifically, the vector representation of the j th subimage is given by (c dj ) ℓ = (C dj ) (ℓ-rP)r´ for ℓ = 0, …, P 2 -1 and r = [ ℓ -P ].

Gabor Representation
The vectors c dj and xd ℓj are first transformed to extract their texture features.A two-dimensional Gabor function is used because it provides a multi-scale and multi-orientation representation of the underlying signals (Clausi & Jernigan, 2000).In the spatial domain, a 2D Gabor function is a sinusoidal plane wave modulated by a Gaussian function.More specifically, the 2D Gabor function is given by (1) where u0 is the frequency and θ is the anti-clockwise rotation of the Gaussian function and the plane wave.The values σx and σy represent the size of the Gaussian function or the scale in the x and y dimensions, respectively.
The Gabor function in (1) can be expressed in matrix notation as: (2) where Wf u θ v represents the Gabor filter bank at frequency u and rotation v.
Using W, a discrete set of transformations is obtained.This set with different frequencies and orientations is used to extract features from the subimages.Selection of the frequencies emphasizes the intermediate frequency band as the most significant information about a texture often appears in the middle frequency channels (Zhang, Tan, & Ma, 2002), (Chang & Kuo, 1993).The frequencies are given by (3) where the set of frequencies is fu = { FL, FH } and P is the size of the subimages.However, the literature indicates that a finer quantization of orientation (larger number of rotations) is needed.The restriction on the decision about the number V of rotations, θ, is based on the computational efficiency.
Given a training vector c dj , the respective training Gabor transformation can be written in matrix notation as (4) where u and v indicate the frequency and the rotation, respectively.
Analogously, the Gabor transformation of the test vectors xd ℓ j can be written as ( 5) Figure 1 shows the representation of the data sets and the process used to obtain the respective transformations.

Feature Extraction
Two features, the mean and the standard deviation of the classification system, can be directly estimated using the Gabor transformation in ( 4) and ( 5); the index C represents the set of images to which the transformation belongs.These features are calculated for each subimage of the training set using the transformation in ( 4) as (6) (7) where P is the dimension of the subimage and u is a P 2 -long one-valued vector.Equivalently, the mean and standard deviation are calculated for the test images by using the respective in ( 5) as (8) (9) The energy, contrast, homogeneity and correlation are also calculated for the subimages, in addition to the mean and the standard deviation given in ( 6), ( 7), ( 8) and ( 9).The co-occurrence matrix is defined as a matrix or distribution of co-occurring values at a given offset.The use of this matrix is based on the assumption that the texture information on the image vector or is contained in the overall or average spatial relationships of these vectors (Haralick, Shanmugan, & Dinstein, 1973).The co-occurrence matrix is defined as G = Q ( ) where Q is an operator that defines the position of two pixels relative to each other.Let the elements of G be gij.Then, gij represents the number of times that pixel pairs with intensities i and j are present in an image in the position specified by Q.
The size of this co-occurrence matrix is determined by the number of possible intensity levels in the image.For an 8-bit (256 possible levels) image, G will be 256x256.To reduce computation load, a frequently used approach is to quantize the intensities into a few bands to keep the size of the matrix G treatable.In the case of the 256 intensities, a quantization on B levels, where L = 8 results in a co-occurrence matrix of size B x B (Gonzalez & Woods, 2008).The co-occurrence matrix associated with the training image vector is given by .
(10) Accordingly, the co-occurrence matrix for the image vector is Using the co-occurrence matrices, the energy, contrast, homogeneity and correlation are calculated.The angular second moment or energy for the training set is calculated by .
(12) The contrast or the measure of the amount of local variation is given by .
(13) The direct measure of the local homogeneity is calculated by ( 14) Finally, the correlation or the measure of linear dependency ρ of neighboring image pixels is given by ( 15) where µ k , µ r , σ k and σ r are the means and standard deviations of g c and g c , respectively.
Similarly, textural features are calculated for the testing image vector .The 6 characteristics for the training set in ( 6), ( 7) and ( 12)-( 15) are combined into the matrix F C shown in equation ( 16) for the H sub-images and the UV Gabor transformations of each of those sub-images.Similarly, the matrix F X is constructed using the characteristics for the testing set.The columns of this matrix represent the H subimages, and the rows represent the filtered images UV for each of the k features.A flowchart of the procedure used to obtain the matrix of features of any image and to make a final decision about its classification is presented in Figure 2. The matrices of features F C and F X represent the texture at determined frequencies U and orientations V after filtering by the Gabor transformation.Some of the frequency-orientation combinations will be more representative than others, taking into account that the real frequency-rotation of the texture is adjusted to the combination given by the Gabor function.
The analysis of the data allows an image to be characterized using only a small set of the images filtered by Gabor filters.Thus, a portion of the data extracted from and can be used for the classification while the other portion is noise, making the classification more difficult.For this reason, a classification based on distances is used in this model, following heuristic rules to use the most representative features.Two approaches for the estimation of the depth d by the classification of a caving are proposed.

Classification based on a hard thresholding
In this approach, the goal is to find the subset of feature vectors that determines the texture of an image.Each feature matrix of the training classes is compared with the feature matrix of a test image.A sum of the differences of the matrices for small distances is calculated.To calculate those distances and to make a decision related to the similarity among vectors, a threshold is included in the model.This threshold τ is estimated using the training set of images.The distances in this approach are given by: (17) for h = 0, …, H -1; ℓ = 0, …, UVH -1 and .Then, the maximum value for similar features present in is used for the classification, and this selection can be represented as (18) To make the decision, the training set class with the greater number of similar feature sets is established as the class to which the test image belongs.This decision is given by the comparison among the r C of the possible classes and is defined as: class = max c r c (19)

Classification based on the Euclidean Norm (ℓ 2 )
Instead of the sum presented with the previous approach, this method uses the distances between vectors.To calculate those distances, is calculated as: (20) for h = 0, …, H -1; ℓ = 0, …, UVH -1 and .
The matrix contains the minimum distances between the features in the training and test images, and those minimum distances characterize each subimage and image, respectively.To use the most representative distances, a distance per class is obtained by r c =min h min ℓ (S c ) ℓ,h (21) For comparison among classes, the distances between the classes are evaluated and the class with the minimum distance is selected as the class to which the test image belongs.This decision is given by the evaluation of the distances r C as: class = min c r c (22) Finally, weight can be applied in the model to give more importance to some of the features.To apply this weight, a vector ω with the weights is applied to the feature vector before the comparisons in equations ( 17) and ( 20).Then, the vector ω can be expressed as a function of the weights, ω Given the images of the test set, the accuracy of each of the approaches is calculated as the sum of the successful classifications over the number of evaluations, which consists of the evaluation of the L images in the test set in this case.The equations for the calculation of the performance are ( 23)

Experimental setup
To evaluate the efficiency of the presented model, rock samples extracted from a wellbore in a Colombian basin were used.Cuttings from depths between 1371.6 and 3825.54meters (4500-12550 ft) were acquired via realtime monitoring by the Colombian Petroleum Institute (ICP).To facilitate the process of acquisition and reduce uncertainty, specimens were created in an ICP laboratory containing samples of cuttings and caving materials.Images of some samples are presented in Figure 3.For the preparation of the specimens, the rock samples were placed in a metallic mold lined with vacuum grease or acrylic.A mixture of resin and hardening in a ratio of 1 to 2 grams was added and allowed to dry.When the specimen was completely hardened, it was extracted from the mold.Finally, the surface was polished until the visible area of the sample surface was maximized.Figure 4 shows the process used to generate specimens.A set of images was acquired using a reflected-light microscope.Figure 5 presents some of these images, and it should be noted that different textures are present.Training set C and testing set X were thus randomly created.To simulate the classification process, a set of D = 17 classes was chosen, each representing a depth in a well as illustrated in Table 1.Then, the mathematical model presented in section 2 was applied.A subsection of each of the images in C and X was extracted to eliminate the blurred borders that were present.Then, each image C d and Xd ℓ was divided into 12 non-overlapping subimages.
The preprocessing was performed using a bank of Gabor filters W, as presented in equation ( 2).For each subimage of C d and Xd ℓ , a set of U = 14 frequencies and V=6 rotations were applied, for a total of UV = 84 Gabor filters.For each of the transformations and , k = 6 texture descriptors were calculated.An analysis of the descriptors among the different D classes allowed determining the difficulty of separation.This is the main reason to use this approach instead of traditional classification methods, given that the classes are not linearly separable.To show the non-separable characteristic, the following figures present the histograms and some examples of the features.The histograms in Figure 6 show the concentration of the descriptors in the same value ranges for two randomly selected classes.It is important to clarify that the domain does not significantly change between classes and that the overlapping behavior is similar for any pair of selected classes.
Figure 7 shows samples of two randomly selected features from the data set, and again, the overlap is similar for any pair of selected classes.Each data point is labeled according to the classes it belongs to.Here, the classes are depicted in green and blue, and the goal is to use these data as a training set to classify a new observation according to the feature classes.Note that class overlap is present in the projected space.
In evaluating the computational complexity of the algorithms, we consider two stages: feature extraction and classification.The two stages can be performed efficiently, requiring Ơ(UVH + D 2 ) operations per image, where H is the number of subimages (12 in our experiments), U and V are the frequency and rotation of the Gabor filter bank, respectively, and D is the number of classes.showing that the classes are not linearly separable.

Simulations and results
The classification was performed following the mathematical model in section 2. To calculate the performance of this method, the correct classification of each of the binary classifiers must be considered.The general performance calculated using equation ( 23) for each of the approaches-hard thresholding and ℓ2 -norm-is presented in Figure 8.The performance results obtained using real data and the two approaches are 91.3% and 87.5%, respectively, for the hard thresholding and ℓ2-norm approach.
To highlight the influence of the different features, the performance was calculated using different numbers and combinations of features.Figure 9 shows the performance using 4,5 and 6 features, and each of the bars represents a different combination of features.The comparison shows that the performance increases with the number of features and becomes more precise for k = 4,5,6.The best performance is obtained using the complete set of features.Given the variability in the performance, especially when 4 features were used, an analysis was performed to determine whether some features are dominant over the others following the weighting method in the mathematical model.Simulations involving the total number of features and assigning weights were performed by applying the vector ω.Improvements in the performance obtained using the hard thresholding classification resulted in 92%, 93% and 94% success when the weight values ω 2 , ω 4 and ω 5 were respectively applied to the standard deviation, contrast and homogeneity with greater consideration compared with the other features.The weight values considered in the analysis were ω 2 = 1.2, ω 4 = 1.6 and ω 5 = 1.6.Different settings were tested in the hard thresholding classification, which includes a threshold value τ. Figure 10 shows the variation in classification accuracy for 13 different values of τ for the proposed κ features and distance measures.Figure 10 shows the independence of the performance from the ℓ2norm distance, represented by the red line, and the best performance achieved was 91.3% correct classification using the hard thresholding approach and a τ value of 0.008.In Figure 11b, some of the classifier errors are present in binary classifiers involving the same class.For instance, the errors in classes 2 -16, 11 -16 and, 12 -16 indicate that some classes as the 16 in this case present very similar features compared to other classes, which makes the classification more difficult.

Conclusions
The classification of caving depth estimation using hard thresholding and the ℓ2-norm has been developed in this work.The developed methods and a complete mathematical model are then used to classify caving depths by features extracted from cutting images.The multi-scale and multi-orientation representation obtained using Gabor filters permit a much more complex characterization of the texture due to the independence of the rotation used in the acquisition process.The performed analysis precisely extracts the texture features for different rotations, which allows an optimal comparison among the images in a training set and a testing set.
The size of the images and the number of Gabor filters applied are important because these parameters allow a complete characterization of the texture.However, as a consequence, additional information is included.This information is noise and should be removed for the classification process.The approaches in this work treat this problem by searching a reduced and a representative set.Additionally, simulations using different numbers of feature combinations reveal variation in the performance of these methods.However, none of the combinations with fewer features than the established κ achieved a greater performance.
The performance achieved with the hard thresholding classification is 91.3%, whereas the ℓ2-norm based classification achieved 87.5% success.These values demonstrate the applicability of the methods and techniques developed in this work.In addition, these methods can be applied directly in rock classification problems.
The training set is denoted as C=[C 1 , …, C d , …, C D ], and the testing set is denoted as X=[Xd 1 , …, Xd L ] where D and L indicate the number of images in each group and, traditionally, D » L. The images C d and Xd ℓ are N M matrices representing the d th training and testing images, respectively.Figure 1 shows the representation of the data sets.The rock images in the training set belong to rocks that have a known depth d, and therefore, they are called cuttings.Conversely, images in the testing set belong to regions of caving, whose depth d must be estimated.Each image C d and Xd ℓ is again divided into H non-overlapping P P subimages to work with small image sections and non-repeated information.The next step involves the extraction of the textural features from the training and test subsets.The selected method of texture analysis for feature extraction is critical to the success of the texture classification.A multichannel Gabor function is used for the feature extraction.

Figure 1 .
Figure 1.Data representation in the model.There are two sets of images; each image in a set is divided in subimages, and then, the Gabor transformations are obtained.
Let C d and Xd ℓ be N M matrices representing cutting and caving images, where d and dℓ represent the image depth, and let C and X be the ensemble of images C d and Xd ℓ , respectively.Thus, C=[C 1 , …, C d , …, C D ] and X=[Xd 1 , …, Xd ℓ , …, Xd L ].Note that C represents the training image set, X represents the test image set, and D is the number of depths.These values are associated with the number of classes in the system.Note also that C 1 is the image from the first depth and that C d is the image from the depth d.The images C d and Xd ℓ are divided into H non-overlapping P×P subimages, such that H={h 1 }{h 2 }, where h 1 =[N/P] and h 2 =[M/P].The j th subimage of C d is denoted by C dj , and the j th subimage of Xd ℓ is denoted by Xd ℓj .The image C d can be expressed as a function of its subimages: C d =[C d1 , C d2 , …, C dh ].Similarly, the image Xd ℓ can be expressed as a function of its subimages: Xd ℓ =[Xd ℓ1 , …, Xd ℓH ]. Figure 1 shows this representation of the data sets.

Figure 2 .
Figure 2. Flowchart of the process used to estimate the matrix of features.

Figure 3 .
Figure 3. Real cuttings samples selected for creation of specimens.

Figure 4 .
Figure 4. Specimen preparation.The rock samples are placed into a mold lined with acrylic, and a resin-hardener is added and allowed to dry.The specimen is extracted from the mold when it is consistent and polished to obtain the maximum visible area of the samples.

Figure 5 .
Figure 5. Rock images for classification.Note the different textures present in the set.

Figure 6 .
Figure 6.Comparison of histograms for different classes.The descriptors for different classes are concentrated in the same value ranges.

Figure 7 .
Figure 7. Scatter plot of depth features for two randomly chosen features represented by green and blue.There is a class overlap in the projected spaceshowing that the classes are not linearly separable.

Figure 8 .
Figure 8. Overall performance evaluated using hard thresholding and ℓ2norm distances.

Figure 9 .
Figure 9.Comparison of the performance using different values of k.

Figure 10 .
Figure 10.Performance for different τ values represented by the blue bars.The performance using the ℓ2-norm distances is shown in red and is independent of the τ values.The performance of the binary classifiers is presented in Figure 11; values of 0 (blue) above the diagonal represent classification errors and values of 1 (red) represent a successful classification.The point highlighted in Figure 11a represents a classification error in the binary classifier among classes 2 and 4 for the hard thresholding classification.

Figure 11 .
Binary classifiers using the (a) heuristic and the (b) ℓ2-norm decisions.In both figures, the 0 (blue) values represent errors and the 1 (red) values represent successful classifications.