Deep learning based restoration of lost sections in Micro-CT core plugs

Lack of petrophysical information is critical for reservoirs development composed of poorly consolidated rocks or for zones bearing wells with core damaged by improper coring operations. The restoration complexity of the digital-core lost sections is associated with the need to consider an enormous amount of data from the existing core image and the necessity to include lithological expert knowledge. That makes deep learning methods a natural choice for solving such problems. We proposed, examined, and compared several deep learning methods convenient for analyzing micro-computed tomography digital core data. It was done under the most simplistic problem statement when the destroyed part (a set of slices) is completely lost. Here, we present the results of comparison interpolation/extrapolation procedures under proposed quality metrics. We discover that the variational autoencoder method can be trained to extract some petrophysical parameters from the digital core plug in an unsupervised manner.


Introduction
Petrophysical parameters are crucial for the characterization of oil and gas reservoirs. The accuracy of the oil field production predicting indicators significantly depends on the amount, reliability, and representativeness of petrophysical data obtained from various sources.
The problem of insufficient petrophysical information arises during the exploration and development of reservoirs composed of poorly consolidated rocks. This situation can also crop out because of incorrect coring operations. Therefore, loose reservoir rocks and intervals with low core recovery are complicated targets for petrophysical examination.
Restoration of the lost or destroyed core sections, depending on the position of the disappeared volume, is the problem of extrapolation or interpolation of the core digital image obtained using micro-computed tomography (Micro-CT). The computational complexity of this task associates with the need to consider an enormous amount of data from the existing parts of the digital core adjusting to the restoring volume. In addition, we need to involve expert knowledge about the potentially admissible core structures, which is problematic. This consideration makes deep learning (DL) methods a natural choice for solving such tasks because of their success in solving various computer vision problems.

Deep learning and data size
The evident problem in the application of DL methods to x-ray Micro-CT data is its amount. For example, the pack of slices generated from the core plug sinogram usually consists of 1000 images IOP Publishing doi: 10.1088/1757-899X/1201/1/012070 2 1000x1000 pixels each. In the single-precision floating-point format, popular in DL applications, the size of a digital slice is about 4 MB, and the core plug is 4 GB.
Efficient processing of data of this size usually requires deeper networks, which in turn will use even more memory to store weights.
For the DL applications, GPGPU is the most popular hardware, but the limitation on the size of its internal memory becomes the main bottleneck. The solution is to use parallel learning strategies on a set of GPGPUs installed in one node. Unfortunately, in this case, arises the problem of synchronization, which can dramatically slow calculations. A tempting alternative is to use the DL-based Super-Resolution (SR) approach. It was successfully applied to several problems [1]. Using DL-based SR, it is possible to diminish the size of the input data, perform the necessary calculations, and flash the result to SR, returning to the original size without a noticeable accuracy lost due to data compression-decompression. To apply this approach, we need to define the metric of the SR, especially for Micro-CT.
Here we have considered only approaches that allow usage of DL at an original image resolution without involving the SR methods or usage parallelly stackable GPGPUs. The proposed methods were trained in a lower image resolution to diminish computational time.

Lost digital core plug sections restoration
Total 3D processing of the digital core plug is the only correct approach for reconstructing the internal porous rock void space. It allows to highlight structural features and generate total vertical and horizontal spatial connectivity between plug's pores. Also, it excludes the orientation error that occurs during independent 2D processing slices belonging to the digital core pack [2]. For the diminishing calculations and simplicity of the results visualization, we will use the classical method, as processing 2D plug slices will significantly reduce the training time and training complexity. Also, all calculations can be executed by one GPGPU with at least 16GB of local memory.
In Figure 1, we present three principal cases of destroyed or lost core sections. Further, we consider the simplest form of the digital segmented core restoration problem when the destroyed part is completely lost. Cases (a) and (b) are slice image prediction problem (interpolation or extrapolation). Obviously, for case (c), this approach is not proper. For this case, we need to create another strategy that will correctly extract and use information from the preserved core pieces, including a method to restore their position and orientation.  We use the binary segmented digital core, as it is much simpler than other image forms. Main essential tasks of petrophysical properties estimation can be conducted on this segmentation [3]. At this stage, it is unclear what is the best metric to evaluate the predicted slice. It is possible to use simple pixel-wise metrics like mean squared error (MSE) or mean absolute error (MAE), learnable perceptual metrics or a learned discriminator from adversarial training [4].
We select MSE and mean intersection over union (IoU), which is applicable due to the binary segmentation of image. IoU proved to be helpful in the evaluation of the quality of objects localization and segmentation.

Dataset and technical aspects
Proposed restoration methods were tested on a synthetic example -reconstruction of deliberately removed slices from a pack of digital core plug.
There are only a few open datasets of digital core samples obtained from Micro-CT. One of the most comprehensive resources that incorporate data about the internal structure of various porous media is Digitalrocks [5]. For training and testing DL methods, we chose dataset [6] consisting of 11 digital core samples. Examples of segmented sections of each core are presented in Figure 2. The dataset was split into training and validation samples 80% and 20% of its size. We left Buff Berea slices for testing the overall performance of the proposed methods and the capability of their generalization.
To implement and test proposed DL methods we used Python 3 with Tensorflow 2 library. Neural networks (NN) were built with use of tensorflow.keras layers implementation and the provided architectures presented in Keras notation. The calculations were conducted on Nvidia GPUs P100 and V100 provided by Google Colab.

DL slice image prediction
In this work, we consider three different approaches in their basic implementations as we want to know what direction of research is more natural and appropriate for the stated problem.
Since the fact that NN produces continuous outputs, we apply binarization (2.1) on min-max scaled slices as a post-processing technique.
is widely used type of recurrent neural networks (RNNs). Further, we evaluate the convolutional form of LSTM [8]. Inspired by [9] we considered the problem of slice image forecast as next frame prediction where the main goal is to predict frames of dynamic visual scenes using information from the past. We tested different LSTM architectures but failed in training any adequate model. The results were either unsuitable and not applicable for further examination, probably because of underfitting, or the model reproduced only the last slice in the sequence as a prediction (apparently problem of overfitting).
LSTM nets are hard to train due to low parallelization capabilities and large memory consumption. The next frame prediction is usually used to predict deterministic motion on the scene, meaning that such an approach is not appropriate in this case due to the stochastic nature of the core plug slice prediction.

VAE
Variational Autoencoder (VAE) is a generative model based on variational approximations [10]. The idea of application generative DL models for studying core plug porous media is not new [11]. But we did not find any information concerning the application of VAE for this purpose. VAE describes the latent space in a probabilistic distribution manner rather than producing points in the latent space. VAE learns a probability distribution, which makes this approach attractive considering the quasi-stochastic nature of slices images sequence. The key idea is to make the autoencoder latent space similar to some prior distribution p(z) and then generate new data using it. The classical VAE represents an autoencoder model, where the encoder maps data sample x in its latent space representation z by producing parameters of learned approximate posterior distribution q(z|x) and sampling from it. The decoder then maps z to parameters of data likelihood distribution p(x|z) In VAE, we regularize the latent space with Kullback-Leibler (KL) divergence, making it similar to a known prior p(z) which is usually multivariate Gaussian N (0, I). Negative loglikelihood represents the reconstruction loss. Then VAE loss (2.3) is the sum of reconstruction loss and the prior regularization.
where D KL is the Kullback-Leibler divergence.
Generating new data samples with VAE comprises picking a random point in latent space, which corresponds to sampling ξ from some prior distribution and then passing it to the decoder (generator) to produce a sample in data space.
A significant property of VAE is the ability for smooth interpolation of data samples in latent space. It happens because the KL term forces each dimension of the posterior distribution to be concentrated around the origin with a nonzero variance, resulting in overlapping between the distribution of distinct data points. At the same time, the reconstruction term encourages the decoder to produce the intended output. That means that close points in the latent space should generate neighboring data samples in the sense of the chosen reconstruction loss function, which IOP Publishing doi:10.1088/1757-899X/1201/1/012070 5 usually can not be achieved when using a plain autoencoder model. This property makes VAE a manifold learning algorithm.
VAE is straightforward for implementation in different versions of model architectures. We used 2D convolution layers and dense layers to produce the parameters of the posterior distribution. We build an extensive network for faster convergence to test this approach, but the architecture should be carefully designed as large networks are prone to overfitting. The total architecture is shown in Table 1. Conv2DTranspose, 128 3x3 kernels, ReLU activation, strides 2 conv t 1 conv t 3 Conv2DTranspose, 64 3x3 kernels, ReLU activation, strides 2 conv t 2 conv t 4 Conv2DTranspose, 32 3x3 kernels, ReLU activation, strides 2 conv t 3 conv Conv2D 1 3x3 kernels, sigmoid activation, strides 1 conv t 4 resize Bilinear image resize to input width and height (output) conv We trained VAE on 250x250px size slice images with ADAM optimizer, learning rate 10 −4 and batch size 4 for 100 epochs. The example of VAE reconstruction for a sample from the training split (Berea core plug) is shown in Figure 3. generate images with sharp edges and fine details [12]. To mitigate or overcome this problem, one may use a learnable reconstruction similarity function as proposed in VAEGAN [13] or a special multistage architecture with a refinement network to improve the coarse output from VAE [12].
Latent space analysis may provide some insights into how VAE encodes the multidimensional input information into a vector. Also, it answers the question -whether the independent components of the vector have some interpretable meaning, if it is disentangled or not. Usually, digital core plugs are needed for the evaluation of petrophysical parameters. We selected the porosity P per slice (2.4) as the most straightforward in both guises as a physical and graphical parameter for calculation ( Figure 4). Then we mapped each slice to its mean encoding (z = µ, where µ is the mean of learned Gaussian posterior distribution).
To find latent variables similar to porosity, we used correlation distance c(u, v) (2.5) between min-max scaled values of porosity and each latent variable independently.
where u and v is the mean of elements of vectors u and v, x · y is the dot product of x and y. We found that there is more than one latent variable similar to the porosity (Figure 4). VAE highlighted porosity in slices latent space. However, the representation is entangled between different latent variables. Checking VAE on the test core plug showed it didn't generalize well and possessed some overfitting. Nevertheless, we didn't consider only online slice generation, and we fine-tuned VAE for ten epochs on test core plug slices, except 20, to check if one can use VAE to interpolate IOP Publishing doi:10.1088/1757-899X/1201/1/012070 7 these missing slice images. After the fine-tuning, we encoded all slices of the test core plug with their mean and plotted the 128-dimensional latent space points with PCA and t-SNE ( Figure  5). We can see that core plug latent space points form a direction prediction curve ordered according to slices sequence. It is significant that the unseen slices in the process of training also belong to this curve. Now it is possible to convert the interpolation/extrapolation problem from the image (slice) domain to one-dimensional vector function domain g(s) = l, where s is a slice index, l is a latent vector and g(s) is a function we need to interpolate/extrapolate using produced slices encodings. We used Akima spline [14] to interpolate the encodings for missing slices. The prediction results for the first six slices are shown in Figure 6, the resulting mean IoU is 0.71857774, and MSE is 0.32730734.
(a) (b) Figure 6. Example of VAE interpolation application. The first row is ground-truth, the second is reconstruction, the third is prediction. (a) slices, (b) slice patches.
Providing an essentially larger dataset of core plugs makes it possible to use LSTM for encodings interpolation and extrapolation.
It seems that the proposed technique with well generalized VAE can be used for interpolation/extrapolation of a large number of slices. However, it is still unclear whether the predicted slices will have numerically estimated petrophysical parameters such as permeability close to the real ones. Possible solutions are controlling the physics directly by conditional VAE (CVAE) or improving the disentangling ability of the encoder network. This can be achieved, for example, with a weighted KL term in VAE loss function L V AE = L rec + βL KL [15] where β > 1, thought this modification degrades the curve formation of core plug latent representation. There also exist more sophisticated methods for solving next frame prediction by VAE application [16].

Self-supervised DL optical flow net
Optical flow estimation is an ill-posed inverse problem solution of evaluating an object's motion patterns between images of a scene for different time steps. Alternatively, it can be presented as a calculation of the motion of image intensities or pixels [17]. DL methods for optical flow estimation gained rapid development recently and considered being the state-of-the-art methods tested on many benchmarks, for example KITTI-2015 [18]. Optical flow is applicable in many video processing problems; automatic navigation tasks like visual odometry, object detection, tracking, and segmentation [19]. It is also one of the principal research areas in video frames interpolation or next frame prediction.
We describe dense optical flow f low I t I t+k (x, y) = [δx, δy] between two slices I t and I t+k as a per-pixel point displacements δx and δy. This pixel interconnection I t+k (x, y) = I t (x+δx, y+δy) makes it is possible to interpolate slices in-between by interpolating intermediate displacements in optical flow and warping corresponding slices where t < r < t + k , I t is the reference slice, I t+k is the target slice, I t+r is an intermediate slice between the reference and target slices, w(I, f low) is some warping procedure and g(r; a, b) is interpolated value at r between points a and b.
The lack of ground truth requires training the network in an unsupervised or self-supervised manner. We propose to use the reconstruction error between the warped frame I t+k and the target frame I t+k as the loss function (2.7). (2.7) To train the network by means of stochastic gradient descent methods w must be differentiable.
We will use realization of differentiable sampling [20] from tensorflow graphics.image.transformer.sample module with a bilinear kernel. For simplicity, we set g(x; a, b) as a linear interpolation. Here we do not consider an occlusion problem, as we suppose that the number of pixels in the reference slice is always enough to warp it into the target slice and vice versa. In Table 2 we present total architecture. We trained optical flow DL net on 250x250px size slice images with ADAM optimizer, learning rate 10 −4 and batch size 4 for 100 epochs. We trained this network on pairs (I t , I t+k ), where k = 21, so we can compare interpolation results with VAE on the same stack of slices. The example of the warped image and the estimated optical flow visualized with [21] for a pair of samples from the training set (Berea core plug) are shown in Figure 7. Results for test core plug slices interpolation with optical flow presented in Figure 8, the mean IoU is 0.7173934 and MSE is 0.32103884. We can conclude that linear interpolation is an incorrect choice. Pixel movements look unnatural, meaning that optical flow training needs additional regularization, possibly encourage small displacements. Another problem with the traditional optical flow approach is that unseen structures or formations can not pop up during direct interpolation.

Results
As the quality measurement baseline for presented methods, we use the naïve constant first frame interpolation -assignment of the last slice to all predicted I t+r predicted := I t+k ground−truth . Figure 9 present results of interpolation for the most distant slice (number 10 from 20 slices long interpolation task). Table 3 contains quality metrics for the same slice image.  Results of VAE and optical flow DL methods are slightly different in MSE. VAE produced two times better results in terms of porosity error in comparison with the naïve approach. We conclude that MSE is an inappropriate metric for problems of this type. IoU quality metric proved to be better than MSE. However, appropriate metrics need to be developed based on the graphic properties controlling results of future physical processes simulations on a digital core plug.
According to Table 3, VAE has the best score in all metrics. Because of the heuristic nature of deep learning methods, further research needs to be conducted with other DL approaches and architectures. Obtained results show that DL and especially VAE are efficient for lost or destroyed core body restoration.

Conclusion
The problem of lost Micro-CT scanned data restoration is poorly covered in the existing informational sources. We consider that restoring lost cross-sections in the data domain is a more promising area for research than predicting petrophysical parameters of the lost core plug volume. Restoration of void space and porous media matrix allows researchers to apply a wide range of numerical methods for assessing different petrophysical parameters.
This article proposes DL-based solutions to this problem. It shows that the developed VAE method is the best among those considered in the pixel-wise error of the predicted slices and their petrophysical properties. We need to emphasize the inappropriateness of pixel-wise quality metrics to this specific task, which requires the development of physics-aware metrics and loss functions for DL.
We hope further research improves the quality of lost data restoration and the actual depth of interpolation/extrapolation volume.