Digital core repository coupled with machine learning as a tool to classify and assess petrophysical rock properties

. To make efficient use of image-based rock physics workflow, it is necessary to optimize different criteria, among which: quantity, representativeness, size and resolution. Advances in artificial intelligence give insights of databases potential. Deep learning methods not only enable to classify rock images, but could also help to estimate their petrophysical properties. In this study we prepare a set of thousands high-resolution 3D images captured in a set of four reservoir rock samples as a base for learning and training. The Voxilon software computes numerical petrophysical analysis. We identify different descriptors directly from 3D images used as inputs. We use convolutional neural network modelling with supervised training using TensorFlow framework. Using approximately fifteen thousand 2D images to drive the classification network, the test on thousand unseen images shows any error of rock type misclassification. The porosity trend provides good fit between digital benchmark datasets and machine learning tests. In a few minutes, database screening classifies carbonates and sandstones images and associates the porosity values and distribution. This work aims at conveying the potential of deep learning method in reservoir characterization to petroleum research, to illustrate how a smart image-based rock physics database at industrial scale can swiftly give access to rock properties.


Introduction
Digital rock analysis recently becomes an important part of the laboratory services in oil and gas industry.Numerical properties accelerate and improve the understanding of the reservoir behavior.
Computed tomography scans from rocks associated with Digital Rock Physics (DRP) analysis receives considerable and extended use in the oil and gas laboratory services.Computed micro-tomography (MCT) scans of rocks record multiple structural information such as the texture and the rock fabric.A 3D image of a rock sample gives access to the actual representation of the mineral phase and the pore space.Rock images with voxel size down to the micrometer resolution allow to extract the topology of the pore space.Once segmented in two (usually pore and solid) or more phases of interest, the MCT images can be used to simulate porous media physical properties such as fluids transport, electrical and geo-mechanical properties as well as processes like enhanced oil recovery simulation [1,2].Ongoing improvements in MCT systems and image analysis software quickly provide increasing amounts of data.Properties computed from rock images accelerate and improve the understanding of the reservoir behavior.
In DRP analysis, key elements are the voxel size and the representative elementary volume (REV) determined during the acquisition phase.A too low image resolution can restrict the pore space reconstruction, causing bad pore space topology determination and rock properties estimation as demonstrated in [3].Al-Raoush [4] and Papadopoulos reveals that the minimal REV for particle size distribution and coordination number is larger than the minimal REV for porosity.This technical restriction results in a permanent search of increasing resolution and volume size of the acquired images.At the same time the constant progress of acquisition technique leads to more and more image acquisition generating dozen terabytes of data needing to be processed afterward.As a consequence, standard image processing techniques become more and more limited and need to be automated to work with an efficient workflow.
Recent advances in high performance computing and machine learning (ML) will probably lead to new and more efficient computations.At the present time, research carried out on deep learning produces various frameworks easily accessible that particularly democratize its uses.Deep neural network provides excellent results for X-ray computed tomography (Wang [5], Würfl et al. [6]).Convolution neural network (CNN) is an important deep learning architecture.It can extract the image features automatically and has high classification accuracy.CNNs have achieved a wide range of applications such as plant classification, face recognition, handwritten Chinese character recognition and so on Mikia et al. [7] and Lopes et al. [8].
Naranjo Leon et al. [9] provided a permeabilityporosity relationship for each rock type, allowing to complement the reservoir characterization in the uncored wells.Chen and Zeng [10] demonstrated the capability of machine learning in rock facies classification from wireline log scalar attributes improving by feature augmentation.Compared to conventional image segmentation methods, machinelearning segmentation could come closer to the ground truth for determining the porosity from noisy MCT images (Berg et al. [11]).Karimpouli and Tahmasebi [12] revealed that a CNN algorithm improves the accuracy of a segmentation comparing with a multiphase thresholding segmentation.Their network also produced valid results for unseen images with a categorical accuracy of about 96%.Araya-Polo et al. [13] used a deep learning architecture to instantaneously predict permeability of clastic rocks from high resolution Scanning Electron Microscopy images.Sudakov et al.
[14] validated a 3D CNN method for predicting permeability of digital Berea sandstone volume subsets.
In this study we pursue tests on the CNN potential for DRP predictions.We provide some proofs of concept and discuss how they could be designed to be integrated in a more global usual analysis workflow.The typical size of a MCT image, a few dozen gigabytes, is huge in comparison to hardware capacity (some gigabytes for GPU) and time consuming.That is why the CNN is relevant to emphasize relevant features with a reasonable time and standard computing resources.
The goal of this study is thus to explore and evaluate the contribution of machine learning as a tool to evaluate geological and petrophysical properties directly from grayscale MCT scan images.The purpose is to deploy an automated workflow directly after the reconstruction of 3D rock images.
The paper is organized as follow: we first introduce the basics of CNN algorithms and the selected image database used for this study.Second, we present the classification of reservoir sedimentary rock types from grayscale MCT images using an adapted version of the Inception-V3 network, named RockClass model.Then an optimized regression-CNN network, named RegPhi model, is used to estimate the total porosity from MCT scan images without segmentation.The porosity is evaluated thanks to an AutoEncoder model realizing an automated segmentation from the grayscale images and associating the total porosity.Finally, the results are discussed.

Machine learning concepts
Convolutional neural networks (CNNs) are a powerful tool widely used in various computer vision problems, like image classification [15], object detection [16], segmentation [17] and image enhancement.As convolutional networks assume that pixels that are close to one another are semantically related, they seem a good candidate to extract physical properties from MCT images.
During the last decade, many network architectures (types of layers, size of layers, interconnexion of layers) have been experimented.In this paper, we use 3 types of networks.
The first one is a customization of the Inception-V3 network as presented in [18] developed by Google to solve the ImageNet Large Scale Visual Recognition Challenge [19].This network has been trained on a very large dataset (1.2 million 2D 299 x 299 images) and is able to classify 1,000 kinds of images.The output of the network is a vector of probabilities for the image to belong to a specific class.
Then we use a network proposed by Sudakov et al. [14] to estimate the porosity from a MCT image.This network has fewer layers than Inception-V3 but is able to accept 3D 100 x 100 x 100 voxel images.Instead of providing a classification, this network is trained to fit a function (porosity in our case).The output of the network is a scalar.
Finally, an AutoEncoder (or encoder-decoder) network is used to produce a pixel-wise output.The key idea is to extract a relatively small set of features (bottleneck) from the input image data and then decode those features into the desired output.Although this kind of CNNs is able to model complicated phenomena due to a large number of parameters, it is still hard to predict the behavior of neural network on "unfamiliar" test examples.There are various techniques to overcome this difficulty: for instance, semi-supervised learning where a mathematical model is incorporated into CNN architecture or loss function as demonstrated in [20 -22].

Dataset
The petroleum reservoir rocks are widely composed of two main lithological classes, sandstones and carbonates [23].The original database is made of a series of reservoir analogue samples including 2 carbonates from Estaillades [24,25] and Savonnieres [26][27][28] formations, and 2 sandstones from Fontainebleau [29,25] and Berea [30][31][32][33] formations.For each sample, the rock type and the rock formation are notified in the generated database which will provide the predicted class of rock type.
See below Table 1 for a description and illustrations of the image collection.
The core plugs are imaged in 3D with a voxel size of 3 to 4 µm.Having the same voxel size ensures a consistent learning of the patterns.This voxel size proposed by most of micro-CT scanners is commonly used for imaging rock porous media.3-µm voxel size enables the internal fine structures to be imaged accurately, though this is the limit between macropore and micropore for the Estaillades sample.Although it remains a challenge to accurately predict porosity from subresolved imaged porous media, we choose to disregard the porosity resolution influence (Saxena et al. [34]).We use only 2-phase segmentation which totally ignores the subresolved porosity.The 2-phase segmentation handed by an expert user reflects the most common arbitrary segmentation and its known sensibility.Sub-resolution porosity involves low contrast and blurred limits between the different phases, which can partly cause a dispersion of the results.These segmented images are taken as references for the training network predicting the porosity.The goal is to try to predict an estimation based on particular subjective expert appreciation, which may differ from the ground truth porosity.For monomineral sample, it is possible to take into account grey levels reported to microporosity thanks to a mean grey segmentation method [25].
Table 1: List of the selected rock images, a collection based on sandstones and carbonates, associated to the CT image information and to the lab and computed basic measurements.Voxaya's software Voxilon [36] is used to extract a dataset of 18 images of size up to 1000 x 1000 x 1000 voxels from those 4 digital plugs.An expert user generates the 3D binary segmented images by selecting the threshold values.The total porosity and absolute permeability are computed from each segmented image.From these large blocks, small 100 x 100 x 100 voxels non-overlapping images are extracted using Python scripting tools [36], leading to a database with more than 36,000 images (18,000 CT and 18,000 segmented).For each segmented image, related pore volume fraction, related permeability and tortuosity are computed, the last two properties not being presented in this paper.For each study, part of the database is used for training and some data is always kept aside for evaluation (training assessment) and testing (model assessment).

Rock classification
Convolution neural network (CNN) can automatically extract image features and presents high classify accuracy.
The goal is to prove that the rock types and formations classification can be realized using existing technologies for image classification.

Network architecture
Here we choose to use the pre-trained Inception version 3 (Inception-V3) model specialized in 2D image classification by feature extraction, provided by the TensorFlow framework.This network is quite optimized and contains 48 layers (see Figure 1 below).A specialized algorithm for training is designed to overcome the limitations of usual algorithms.

Transfer Learning Concept
All the power of the inception approach lives in the Transfer Learning concept.The inception network model is a complex network, so training the model directly from the beginning would cost at least a few days.However, using the method of Transfer Learning, the parameters of the first layers are kept unchanged and only the last layer from the network is adapted to the use case and trained.The last layer is a Softmax classifier, a mathematical function which outputs a probability distribution [37].The porosity for a patch 100 x 100 x 100 is individually estimated.To evaluate the porosity for the whole rock sample one could compute porosity for every patch and average the result.It is replaced by a layer with as many neurons as there are classes to choose from.The network has to be retrained to update the weights of this layer using a back propagation algorithm and the cross entropy loss function.In our case, we use 4 classes: Berea, Fontainebleau, Estaillades and Savonnieres.

Training
As input data, we use 4,000 grayscale 2D slices (1,000 per class) resized to be compatible with the network (initially 1000 x 1000 resized to 299 x 299) and rescaled between 0 and 1.0.Training duration is about 6 minutes on a 32 cores Intel Xeon(R) CPU E5-2667 v3 @ 3.20GHz system.

Results
For testing, 1,200 slices (300 slices per class) from images were analyzed and the computation took 64 seconds.Eventually the data for testing reveals a final accuracy of 100% for the selected rock typing classification.Recognition between the 2 lithological rock types (sandstone and carbonate) is 100% and recognition of the 4 rock formations is also 100%.

Porosity estimation Porosity from CT scan without segmentation
In this paper we try to predict porosity from grayscale 3D images.Araya et al. [13] already used grayscale SEM image as input for DL-based permeability predictions avoiding the segmentation, the highest sensitive step linked to the resolution, and speeding-up processing workflows.
Following that work, we focus on the Berea sandstone sample for training and prediction imaged in 3D.For the network architecture we take inspiration from Sudakov et al. [14].Using this network, the authors try to predict permeability by Linear Regression and CNN from segmented images.Our model RegPhi is implemented using the Keras framework [38] and is described in Table 2.

Training
In our study we use non-overlapping sub-blocks (100 x 100 x 100 voxels) from 2 Berea CT images.We have 8,000 of these sub-blocks, 3,500 being used for training while keeping 600 blocks for validation.

Results
Figure 2 shows the distribution of the relative error on porosity prediction compared to the reference.     ������ −  ��������� / ���������  These predictions are realized on a set of 1,000 subblocks of 100 x 100 x 100 voxels from grayscale Berea sandstone images without any segmentation or improvement process.We use a REV-independent approach in order to provide prediction on large amount of images.Here we can see encouraging results about the porosity prediction, the median of the series is under 15% despite a mean relative error of 18% (see Figure 3).

Network architecture
In this work, we propose an AutoEncoder network that segments pores in 3D rock images (see Figure 4 for a detailed description of the network architecture).
The idea is to compute the high-level features from the input 3D volume that already contain information about porosity in the training example and then decode those features into the segmentation mask, where every value corresponds to the probability of the voxel being a pore.To make sure that the bottleneck layer captures information about porosity, we predict it as an additional pathway in the decoder.We believe that this type of guidance helps the network to cope well with the segmentation task since it already takes the porosity into account while decoding.
Theoretically speaking, we should have trained network with the loss function for porosity and without.Nevertheless we could in general argue that some guidance for the bottleneck is shown to be a promising technique in deep learning.For instance, variational autoencoders are good example how this type of guidance improves the robustness of classification problem (Pu et al. [39]).
Figure 4: Network architecture with 12 convolutional layers, bottleneck, and 12 corresponding upsampling layers.The layers are organized in 4 blocks with 3 layers each.On top of the encoder part, we show the number of features after each layer.On the bottom of the encoder, we show the tensor size.With dashed lines, we illustrate skip connections [40], where encoder features are copied to the corresponding decoder layers.This helps to bring small details that were lost after downscaling.The decoder part has two pathways, the first one recovers the segmentation mask and the second one computes the porosity.

Computational framework
The input data of the network is a set of grayscale CT image volumes of dimension 100 x 100 x 100 voxels.Larger images are segmented into these patches, so that our network can deal with input data of any shape.
The basic ingredients for the encoders and decoders are residual blocks.To decrease resolution, we employ strided convolutions with stride 2, so that our network is fully convolutional.The choice of this architecture is well supported in the literature [41,42] as it allows to speed up the training process and deals with "vanishing gradients", which is a common problem in deep learning.The residual blocks have a very simple structure and allow direct pass-through of the batch normalized input, see Figure 4.
In the encoder pathway, four groups of three residual blocks are chained together, 12 blocks in total.In the first three groups, every third block reduces the patch resolution via strided convolution while increasing feature depth, with the overall goal of gradually reducing dimensionality.Others keep spatial resolution the same while increasing the number of features.In the last group, both blocks decrease spatial resolution such that the output shape is 3 x 3 x 3 x 192, see Figure 4.The feature output at the representation level is concatenated.This is the final output of the encoder, and the bottleneck of the network.
After passing the bottleneck, the low-dimensional representation is decoded again by a chain of residual layers.The latent variables enter two decoder pathways.One is the path that predicts the porosity and another one outputs the segmentation mask.Decoder pathways that lead to the segmented 3D volume use transposed convolutions to exactly revert the encoder on the corresponding level.However, the only link between them is through the latent representation and skip connections, see Figure 5.To prevent overfitting, we apply dropout layers on the bottleneck with 0.2 probability of each node to be discarded.
Figure 5: Single residual block of the encoder and decoder networks.After batch normalization, a first path leads through a (possibly strided) convolution or upsampling layer and a leaky ReLU.A second path either keeps the input or passes it through a strided (transposed) convolution in case it needs to be resampled.Both paths are added together to produce the final output.In the decoder block we skip connect the corresponding encoder features to the output of the decoder.The idea is that it is much easier for such blocks to learn the identity transformation, or perform only small modifications to the input [42], which helps the encoder-decoder paths to gradually add details.

Loss function
In the classical supervised machine learning used to train the model, we need training data together with a label or target.By observing the random variable  and its label  in the training set, the supervised learning tries to fit the model to the training samples ( � ,  � ), . . ., ( � ,  � ).
Test data has similar structure as training data, but the network never sees it.
The training process can be considered as learning the mapping  to predict label  * from the new instance ( * ,  * ) of the test data via ( * )  =   * .
A typical cost function for image classification and segmentation tasks is cross-entropy with logits.
In order to predict the porosity from the bottleneck, we use the standard mean squared error (MSE): ���  = ||ɸ �   ɸ � || � (2) where ɸ � denotes the target porosity and ɸ � is the predicted porosity computed at the bottleneck level.

Training the CNN model
From the training data, we leave aside 10% for a validation set whereas the rest is used for training.Several images are also completely held back and used only for testing.We implement the network using TensorFlow in Python3, and train it on an Intel Xeon(R) CPU E5-2667 v3 @ 3.20GHz system with one Nvidia Quadro M6000.Weights are initialized using the same strategy as for residual networks [40].Stochastic optimization using the Adam optimizer [43] took roughly two days, after which loss remained stable.We train with batch size 10 and learning rate 1-e4.
Reconstruction of a single pathway during evaluation requires about 0.5 -1.5 seconds for the 3D volume of size 100 x 100 x 100.The complete segmentation of a 1000 x 1000 x 1000 voxels image takes about 15 -20 minutes.

Results
Figures 6 and 7 show the source CT image and the resulting segmentations using Hysteresis segmentation and the AutoEncoder network.Globally, AutoEncoder segmentations tend to be sharper than the binary thresholding but also thinner (one voxel difference on the boundaries).Thus, porosity computed on the AutoEncoder images will always be smaller than the reference obtained by Hysteresis segmentation.

Computing resources
The most time-consuming task is the training part of the networks.GPU runs (using one Nvidia Quadro M6000 graphics card) are 15 faster than CPU runs (on a 32 cores Intel Xeon workstation).

Rock classification
In this study, 4 classes of lithological rock formations are distinguished with a highly successful recognition rate.This work will be extended to numerous classes of reservoir rock facies.All authors [4,10,12,14] referring to machine learning mentioned the importance of the quality of the dataset and required a massive amount of data.This database generation with various available reservoir rock formations imaged at several voxel sizes will be a crucial step.
Rock classification using Inception-V3 (RockClass) proves to be a very promising tool.It can be used for pre-calibrating geological-oriented workflow in digital rock analysis proposing a suitable image processing and analysis workflow.The classification rate could be used as a quality control indicator for validating the image quality before any numerical computations.Then a rock class recognition is interesting to offer a fast estimation of the numerical properties already processed from similar rock class.
To further improve digital rock classification, additional physical attributes as mineralogy and pore fabric could probably enhance RockClass capability in detailing geological rock types.For example, pore fabric is directly related to the hydraulic radius size distribution, which is the key parameter for the permeability computation.

Porosity estimation
The porosity is computed on 1,000 sub-blocks of 100 x 100 x 100 voxels by 3 porosity modules: Hysteresis segmentation, RegPhi an AutoEncoder networks.
The porosities estimated from Hysteresis segmentation are used as reference values.They are computed from segmented images obtained by binary thresholding.Porosity estimation obtained by Hysteresis segmentation is benchmarked and validated in several studies with diverse partners (IFPEN, Geosciences Montpellier, [2]).
Figure 8 illustrates the porosity distribution obtained by these 3 methods on 1,000 juxtaposed blocks of size 100 x 100 x 100 voxels belonging to the same 1000 x 1000 x 1000 voxels image.The reference porosity values indicate that the porosity ranges from 9 to 36% with a mean value of 21.5%.This large distribution is mostly due to variation of the pore-solid ratio along the series of small REV of 100 x 100 x 100 voxels with a high voxel size.The RegPhi and AutoEncoder results present a wider porosity distribution, respectively, from 8 to 41% with a mean value of 18.6% and from 8 to 37% with a mean value of 18.8%.The AutoEncoder model is considered being a more reliable network for predicting porosity values close to the values obtained by conventional image analysis.
RegPhi model allows to estimate the porosity without segmentation whereas AutoEncoder generates a segmented image to assess the porosity.Both networks seem to predict similar porosity distributions and from 2.7% (AutoEncoder) to 2.9% (RegPhi) under evaluate the porosity compared to the reference one.This promising result validates improbable capabilities of both models to evaluate porosity directly from grayscale images.An intensified training is required to enhance their performance.
RegPhi and AutoEncoder networks are both trained on these similar subs-volumes of Berea grayscale image.This similarity of results demonstrates how efficient and accurate segmented images are generated by the AutoEncoder model in this study.From these automated segmented images, series of numerical petrophysical properties could be swiftly computed.A complete digital petrophysical workflow can thus be proposed, here for Berea sandstone.

Conclusion
In this study we aim at demonstrating the potential of neural network algorithms to make digital rock analysis simple.We work on implementing an automated workflow from MCT images to digital petrophysical measurements without operator intervention.
The numerical petrophysics workflow is based on a 3-step approach.First, a rock classification from grayscale MCT images using RockClass model is investigated in the two main lithological rock reservoirs types, carbonates and sandstones.Then RegPhi model is used to predict the total porosity from grayscale MCT images.
An AutoEncoder network allows to automatically generate the pore space segmented image and then evaluate the total porosity.
The results of this study are not limited to 4 reservoir rock classes and porosity estimation, but the proposed workflow can be adapted for any CT rock types in order to deliver a complete series of numerical analysis.This work was supported by the ERC Starting Grant ``Light Field Imaging and Analysis'' (LIA 336978, FP7-2014).

Figure 2 :
Figure 2: Distribution of absolute error for porosity prediction in percentage from the same Berea sandstone image.

Figure 3 :
Figure 3: Corresponding box plot of absolute error for porosity prediction in percentage from the same Berea sandstone image.

Figure 6 :Figure 7 :
Figure 6: Comparison between Hysteresis segmentation (Reference) and AutoEncoder model methods: original CT image in grayscale, segmented image with black pixels for pore and greyish for solid.

Figure 8 :
Figure 8: Porosity distribution computed by image analysis (reference), RegPhi and AutoEncoder models for Berea Sandstone sub-blocks.