# Pore Space Reconstruction of Shale Using Improved Variational Autoencoders

Authors: Yi Du; Hongyan Tu; Ting Zhang
Journal: Geofluids
 (2021)

---

## Abstract

Pore space reconstruction is of great significance to some fields such as the study of seepage mechanisms in porous media and reservoir engineering. Shale oil and shale gas, as unconventional petroleum resources with abundant reserves in the whole world, attract extensive attention and have a rapid increase in production. Shale is a type of complex porous medium with evident fluctuations in various mineral compositions, dense structure, and low hardness, leading to a big challenge for the characterization and acquisition of the internal shale structure. Numerical reconstruction technology can achieve the purpose of studying the engineering problems and physical problems through numerical calculation and image display methods, which also can be used to reconstruct a pore structure similar to the real pore spaces through numerical simulation and have the advantages of low cost and good reusability, casting light on the characterization of the internal structure of shale. The recent branch of deep learning, variational auto-encoders (VAEs), has good capabilities of extracting characteristics for reconstructing similar images with the training image (TI). The theory of Fisher information can help to balance the encoder and decoder of VAE in information control. Therefore, this paper proposes an improved VAE to reconstruct shale based on VAE and Fisher information, using a real 3D shale image as a TI, and saves the parameters of neural networks to describe the probability distribution. Compared to some traditional methods, although this proposed method is slower in the first reconstruction, it is much faster in the subsequent reconstructions due to the reuse of the parameters. The proposed method also has advantages in terms of reconstruction quality over the original VAE. The findings of this study can help for better understanding of the seepage mechanisms in shale and the exploration of the shale gas industry.


---

## Body

## 1. Introduction

The demand for oil and gas resources is soaring in the whole world with the fast development of economy. Faced with huge energy demands, the world’s conventional oil and gas production is relatively insufficient, so people begin to pay more attention to unconventional oil and gas resources. Shale oil and shale gas are two kinds of unconventional resources with wide distributions and considerable potentials. However, because of the complex geological conditions of shale reservoirs, it is necessary to figure out the internal structures including the distribution of pore space in shale before any large-scale engineering exploitation [1–3].Shale is a typical porous medium with low porosity, low permeability, and anisotropy, making the exploration and exploitation of shale oil and gas difficult. Some macroscopic properties of shale (such as porosity, permeability, etc.) depend on its microstructures, so it is quite useful to reconstruct a 3D microscopic shale model describing the statistical and topological properties of the pore structure. An important step for building such a model is to reconstruct the internal microstructure of pore space with the characteristics of real shale. The general pore space reconstruction methods are divided into physical experimental methods and numerical reconstruction methods [4–10].Commonly used physical experimental methods are the scanning electron microscope (SEM) [4], the focused ion beam SEM (FIB-SEM) [5], the CT-scanning method [6], and the atomic force microscope SEM (AFM-SEM) [7]. With the help of the above experimental methods and physical equipment, some important parameters and variation characteristics such as the variation in permeability during CO2-CH4 displacement in coal seams [11] can be effectively obtained and studied. Generally speaking, these physical methods normally are fast and their imaging quality mostly is quite good, but the imaging costs are usually quite high and experimental samples are difficult to prepare due to the fragile structure of some porous media like shale, leading to the limitation of physical experimental methods for wide application [12].Some traditional numerical reconstruction methods extract the statistical information (such as porosity and variogram) of the 2D pore structures as constraints to complete reconstruction. For example, the early reconstruction methods such as the simulated annealing method [8], the process-based method [9], and the sequential indicator simulation method [10] use the low-order statistical information for 3D reconstruction. Then, some methods such as the multiple-point statistics (MPS) [13] used for the reconstruction of porous media focus on the reconstruction through the high-order statistical information describing the correlation between multiple points of the pore structures. It is quite obvious that the methods like MPS using high-order information have the advantages in depicting complex pores and the reconstructions are more similar to the realistic pore structures. However, this also makes it a big challenge for the branches of MPS and some improved MPS methods (e.g. direct sampling (DS) [14], filter-based simulation (FILTERSIM) [15], and single normal equation simulation (SNESIM) [13]) to perform reconstruction due to more hardware burdens for CPU and memory as well as a lengthy computational time.Recently, many engineering and scientific researchers have benefited from the fast development of deep learning thanks to its robust extraction capabilities for characteristics, leading to the possibility of deep learning extracting the structural characteristics of porous media for pore space reconstruction [16]. Compared to the traditional reconstruction methods, there are two important advantages introduced by deep learning and its branches. The first one is the training time which can be greatly reduced due to the acceleration of GPU (Graphics Processing Unit) and the reuse of model parameters. Because of the frameworks (e.g., Tensorflow and PyTorch) provided by the public deep learning communities, these deep learning methods can be used by common users even without knowing the complicated assignment and design of parallelization for GPU. The detailed work about GPU is performed using the APIs of the frameworks like the black-box mechanism, which makes it easier for the general public. The second advantage relies on the strong ability of deep learning extracting characteristics from training images (TIs), improving the similarity between reconstructions and the TIs to obtain higher reconstruction quality.As one of the important technologies in deep learning, the autoencoder is an unsupervised learning algorithm whose output can realize the learning and represent of the characteristics of input data. The concept of autoencoder was first put forward in [17]. Hinton et al. [18] derived a fast greedy algorithm based on the autoencoder by using the complementary prior method to generate the deep autoencoder. Vincent et al. [19] added noise to the input data to improve the robustness of the algorithm to form a denoising autoencoder.Based on variational bayes (VB) inference, the variational autoencoder (VAE) model [20] is a deep latent generation model and widely used in image generation, but it often suffers from the complicated expression model leading to the low quality of generated images. In information theory, Fisher information has always been an important tool to describe information behavior in information systems [21]. Fisher introduced Fisher information in the context of statistical estimation, which can describe the behavior of the dynamic system accurately [22]. Recent studies have reported that VAE is difficult to balance the information control between the encoder and the decoder in its structure [23]. For example, when the decoder is too expressive, latent variables generated by the encoder are almost ignored. Therefore, this paper applies the information theory specifically the Fisher information to improve VAE for the reconstruction of shale.Compared to the traditional methods, our method can reuse the parameters obtained in the first training or reconstruction, so the subsequent reconstructions can be largely accelerated without the repeatedly training process existing in the traditional methods, which is a practical contribution for a large-quantity reconstruction mission. In addition, the combination with Fisher information also helps to improve the reconstruction quality of pore space compared to VAE.The sketch of the research technical route of this paper is shown in Figure1. The remainder of this paper is organized below. The main idea of the proposed method is given in Section 2. Section 3 gives a detailed procedure of the proposed method in this paper. Section 4 gives experimental results and comparison with some other methods. Section 5 provides conclusions.Figure 1
The sketch of the research technical route of this paper.


## 2. The Main Idea of the Proposed Method

### 2.1. The Introduction of VAE

In VAE, two models of probability density distribution are established using two neural networks: one is called the encoder, used for the variational inference of original input data to generate the variational probability distribution of the latent variablez; the other is the decoder, reconstructing the data by generating the approximate probability distribution of the original data according to the variational probability distribution of the generated latent variable z.The typical structure of VAE is shown in Figure2, in which “+” and “∗,” respectively, represent the addition and multiplication of elements; X is the input dataset and X′ is the reconstructed result. X usually is quite high-dimensional and used as the training data, whose internal characteristics are learned by the encoder. VAE uses stochastic gradient descent (SGD) to optimize the loss function [20]. SGD is a simple and effective optimization method that updates parameters until losses are within acceptable limits. Although SGD can handle random input, it cannot handle random operation. Hence, VAE needs to “reparameterized” for optimization by introducing an auxiliary parameter ε, which is obtained by sampling from the standard normal distribution N0,1. The latent variable z can be expressed by Equation (1) [24]:
(1)z=µ+σ∗ε,where σ and μ are, respectively, the standard deviation and mean of Gaussian distribution calculated from the encoder. Suppose x is the input vector and x′ is the output vector. The encoder model qφz∣x parameterized by φ is introduced into the encoder network to replace the undetermined real posterior distribution Pθz∣x since the distribution of z is unknowable. Kullback-Leibler (KL) divergence [25] is used to measure the similarity between qφz∣x and Pθz∣x. The constraint parameters φ and θ are repeatedly trained and optimized to minimize the KL divergence, meaning that the evidence lower bound (ELBO) function Lθ,φ;X is maximized, which is defined as follows:
(2)Lθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−DKLqφz∣x∣Pθz,where DKL represents the KL divergence, lb is a base-2 logarithm, Pθx′∣z is the probability distribution of a θ-parameterized decoder, and Pθz represents the probability distribution of the latent variable z. Finally, z is input into the decoder model to obtain the final reconstruction X′. It seems Equation (2) has a penalty item (KL divergence) and a reconstruction item.Figure 2
A typical VAE structural image.


### 2.2. Fisher Information and Shannon Entropy

As can be seen from theLθ,φ;X in Equation (2), only the KL divergence is considered a penalty item for regularization, so it is difficult to balance between the representation of the latent variable z and the likelihood maximization of the reconstructed results [23]. Since the representation of the VAE encoder network is controlled by the divergence error but there is no control exerted on the VAE decoder network, some schemas based on information theory were used to solve this problem; e.g., a constraint item based on Shannon entropy is introduced in ELBO to ensure that latent variables can learn sufficient data characteristics, thus improving the quality of representation learning [26]. However, Shannon information is usually difficult to be processed in calculation, and only approximate substitutes can be obtained in previous studies.For a continuous probability distribution function (PDF), Shannon entropy is a measure of “global characteristics,” which is less sensitive to changes in the distribution of local details of a small area. Fisher information is quite different from Shannon entropy [21]. The former is a measure of the gradient content of the distribution, so it is also quite sensitive to small local details. Let the Shannon entropy and Fisher information of input data X be HX and JX, respectively. Shannon entropy is usually converted into the form of entropy weight [22]:
(3)NX=exp2HX2πexp1,where exp represents the exponential function with base e. When X represents a random vector and Fisher information JX represents a Fisher information matrix, the relation between NX and the trace of JX writes [22, 27]
(4)NX⋅trJX=K,K≥1,where K is a constant and tr∙ represents the trace of a matrix. As shown in Equation (4), Fisher information and Shannon entropy are intrinsically related and have the nature of uncertainty, in which the higher Fisher information is, the lower Shannon entropy is [27], but the calculation of Fisher information is usually easier than Shannon entropy. Hence, Fisher information is considered to be complementary for Shannon entropy and used in this paper.


### 2.3. The Proposed Improved VAE

In this section, Fisher information is used to balance the encoder and decoder of the VAE model. According to Equation (4), it can be found that there is a relationship between Fisher information JX and Shannon entropy weight NX. Shannon entropy is usually difficult to be processed in calculation [28], so Fisher information instead of Shannon entropy is used to calculate information and two new penalty items of Fisher information are added for Lθ,φ;X. The improved method of VAE is called information variational autoencoder (IVAE) hereafter, which not only includes the KL divergence as a penalty item but also includes Fisher information as two penalty items in the encoding stage and the decoding stage, respectively, i.e. IVAE maximizes Lθ,φ;X under the constraint of Fisher information. The new Lθ,φ;X of IVAE writes
(5)Lθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−DKLqφz∣x∣Pθz‐λz∣trJz‐Fz∣⏟PI‐1‐λx′∣trJx′‐Fx′∣⏟PI‐2.Compared withLθ,φ;X shown in Equation (2), Equation (5) adds two more penalty items (PI-1 and PI-2) that adjust the encoder and decoder, respectively: PI-1 controls Fisher information of the encoder network, and PI-2 controls Fisher information of the decoder network. λz and λx′ are adjustment coefficients. Both Fx′ and Fz are positive constants, representing the expected Fisher information values in the decoder and the encoder, respectively, meaning that Fisher information in the decoder and encoder can be controlled by Fx′ and Fz. The larger Fx′ and Fz are, the more likely the distribution is estimated by the θ- and φ-parameterized model; otherwise, the smaller Fx′ and Fz show that the influence of Shannon information is enhanced. Fisher information estimation can be directly calculated according to its definition [28].The optimization procedures of the model are discussed, respectively, for the encoder network and the decoder network in the following section. Take the encoder network as an example. If the Fisher information is only considered in the encoder network, setλx′=0 in Equation (5). Meanwhile, the KL divergence is considered, but the reconstruction item in Equation (5) is not considered. The ELBO of the encoder network, denoted as Leθ,φ;X, contains the error item of KL divergence and information error item:
(6)Leθ,φ;X=−DKLqφz∣x∣Pθz−λztrJz−Fz,where the posterior distribution qφz∣x obeys the normal distribution and Pθz obeys the standard normal distribution. Therefore, the KL divergence in Equation (6) is calculated as [20] follows:
(7)DKLqφz∣x∣Pθz=−121+lbσ2−μ2−σ2.According to the definition,trJz in Equation (6) can be obtained [25]:
(8)trJz=∫z∂∂zqφz∣x2=1σ2.Equations (7) and (8) are substituted into Equation (6), and then Leθ,φ;X writes
(9)Leθ,φ;X=121+lbσ2−μ2−σ2−λz1σ2−Fz.Similarly, Fisher information is used in ELBO at the decoder network, denoted asLdθ,φ;X. Now, the reconstruction item in Equation (5) is considered, but the KL divergence is not considered. Then, set λz=0 in Equation (5), and Ldθ,φ;X writes
(10)Ldθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−λx′trJx′−Fx′.As shown in Equation (2), the original VAE only considers KL divergence and a reconstruction item, while IVAE includes a reconstruction item and three penalty items (KL divergence, PI-1, and PI-2) related to the Fisher information and KL divergence, as shown in Equation (5). IVAE balances both the likelihood estimation and the dependence between input data and latent variables through the new penalty items, making the new model include KL divergence and Fisher information together so as to improve the reconstruction quality.


## 2.1. The Introduction of VAE

In VAE, two models of probability density distribution are established using two neural networks: one is called the encoder, used for the variational inference of original input data to generate the variational probability distribution of the latent variablez; the other is the decoder, reconstructing the data by generating the approximate probability distribution of the original data according to the variational probability distribution of the generated latent variable z.The typical structure of VAE is shown in Figure2, in which “+” and “∗,” respectively, represent the addition and multiplication of elements; X is the input dataset and X′ is the reconstructed result. X usually is quite high-dimensional and used as the training data, whose internal characteristics are learned by the encoder. VAE uses stochastic gradient descent (SGD) to optimize the loss function [20]. SGD is a simple and effective optimization method that updates parameters until losses are within acceptable limits. Although SGD can handle random input, it cannot handle random operation. Hence, VAE needs to “reparameterized” for optimization by introducing an auxiliary parameter ε, which is obtained by sampling from the standard normal distribution N0,1. The latent variable z can be expressed by Equation (1) [24]:
(1)z=µ+σ∗ε,where σ and μ are, respectively, the standard deviation and mean of Gaussian distribution calculated from the encoder. Suppose x is the input vector and x′ is the output vector. The encoder model qφz∣x parameterized by φ is introduced into the encoder network to replace the undetermined real posterior distribution Pθz∣x since the distribution of z is unknowable. Kullback-Leibler (KL) divergence [25] is used to measure the similarity between qφz∣x and Pθz∣x. The constraint parameters φ and θ are repeatedly trained and optimized to minimize the KL divergence, meaning that the evidence lower bound (ELBO) function Lθ,φ;X is maximized, which is defined as follows:
(2)Lθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−DKLqφz∣x∣Pθz,where DKL represents the KL divergence, lb is a base-2 logarithm, Pθx′∣z is the probability distribution of a θ-parameterized decoder, and Pθz represents the probability distribution of the latent variable z. Finally, z is input into the decoder model to obtain the final reconstruction X′. It seems Equation (2) has a penalty item (KL divergence) and a reconstruction item.Figure 2
A typical VAE structural image.


## 2.2. Fisher Information and Shannon Entropy

As can be seen from theLθ,φ;X in Equation (2), only the KL divergence is considered a penalty item for regularization, so it is difficult to balance between the representation of the latent variable z and the likelihood maximization of the reconstructed results [23]. Since the representation of the VAE encoder network is controlled by the divergence error but there is no control exerted on the VAE decoder network, some schemas based on information theory were used to solve this problem; e.g., a constraint item based on Shannon entropy is introduced in ELBO to ensure that latent variables can learn sufficient data characteristics, thus improving the quality of representation learning [26]. However, Shannon information is usually difficult to be processed in calculation, and only approximate substitutes can be obtained in previous studies.For a continuous probability distribution function (PDF), Shannon entropy is a measure of “global characteristics,” which is less sensitive to changes in the distribution of local details of a small area. Fisher information is quite different from Shannon entropy [21]. The former is a measure of the gradient content of the distribution, so it is also quite sensitive to small local details. Let the Shannon entropy and Fisher information of input data X be HX and JX, respectively. Shannon entropy is usually converted into the form of entropy weight [22]:
(3)NX=exp2HX2πexp1,where exp represents the exponential function with base e. When X represents a random vector and Fisher information JX represents a Fisher information matrix, the relation between NX and the trace of JX writes [22, 27]
(4)NX⋅trJX=K,K≥1,where K is a constant and tr∙ represents the trace of a matrix. As shown in Equation (4), Fisher information and Shannon entropy are intrinsically related and have the nature of uncertainty, in which the higher Fisher information is, the lower Shannon entropy is [27], but the calculation of Fisher information is usually easier than Shannon entropy. Hence, Fisher information is considered to be complementary for Shannon entropy and used in this paper.


## 2.3. The Proposed Improved VAE

In this section, Fisher information is used to balance the encoder and decoder of the VAE model. According to Equation (4), it can be found that there is a relationship between Fisher information JX and Shannon entropy weight NX. Shannon entropy is usually difficult to be processed in calculation [28], so Fisher information instead of Shannon entropy is used to calculate information and two new penalty items of Fisher information are added for Lθ,φ;X. The improved method of VAE is called information variational autoencoder (IVAE) hereafter, which not only includes the KL divergence as a penalty item but also includes Fisher information as two penalty items in the encoding stage and the decoding stage, respectively, i.e. IVAE maximizes Lθ,φ;X under the constraint of Fisher information. The new Lθ,φ;X of IVAE writes
(5)Lθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−DKLqφz∣x∣Pθz‐λz∣trJz‐Fz∣⏟PI‐1‐λx′∣trJx′‐Fx′∣⏟PI‐2.Compared withLθ,φ;X shown in Equation (2), Equation (5) adds two more penalty items (PI-1 and PI-2) that adjust the encoder and decoder, respectively: PI-1 controls Fisher information of the encoder network, and PI-2 controls Fisher information of the decoder network. λz and λx′ are adjustment coefficients. Both Fx′ and Fz are positive constants, representing the expected Fisher information values in the decoder and the encoder, respectively, meaning that Fisher information in the decoder and encoder can be controlled by Fx′ and Fz. The larger Fx′ and Fz are, the more likely the distribution is estimated by the θ- and φ-parameterized model; otherwise, the smaller Fx′ and Fz show that the influence of Shannon information is enhanced. Fisher information estimation can be directly calculated according to its definition [28].The optimization procedures of the model are discussed, respectively, for the encoder network and the decoder network in the following section. Take the encoder network as an example. If the Fisher information is only considered in the encoder network, setλx′=0 in Equation (5). Meanwhile, the KL divergence is considered, but the reconstruction item in Equation (5) is not considered. The ELBO of the encoder network, denoted as Leθ,φ;X, contains the error item of KL divergence and information error item:
(6)Leθ,φ;X=−DKLqφz∣x∣Pθz−λztrJz−Fz,where the posterior distribution qφz∣x obeys the normal distribution and Pθz obeys the standard normal distribution. Therefore, the KL divergence in Equation (6) is calculated as [20] follows:
(7)DKLqφz∣x∣Pθz=−121+lbσ2−μ2−σ2.According to the definition,trJz in Equation (6) can be obtained [25]:
(8)trJz=∫z∂∂zqφz∣x2=1σ2.Equations (7) and (8) are substituted into Equation (6), and then Leθ,φ;X writes
(9)Leθ,φ;X=121+lbσ2−μ2−σ2−λz1σ2−Fz.Similarly, Fisher information is used in ELBO at the decoder network, denoted asLdθ,φ;X. Now, the reconstruction item in Equation (5) is considered, but the KL divergence is not considered. Then, set λz=0 in Equation (5), and Ldθ,φ;X writes
(10)Ldθ,φ;X=Eqφz∣xlbPθx′∣z⏟Reconstructionitem−λx′trJx′−Fx′.As shown in Equation (2), the original VAE only considers KL divergence and a reconstruction item, while IVAE includes a reconstruction item and three penalty items (KL divergence, PI-1, and PI-2) related to the Fisher information and KL divergence, as shown in Equation (5). IVAE balances both the likelihood estimation and the dependence between input data and latent variables through the new penalty items, making the new model include KL divergence and Fisher information together so as to improve the reconstruction quality.


## 3. The Procedure of the Proposed Method

The steps of the proposed method are as follows:Step 1.
Input the TI of real shale and initialize all parameters.Step 2.
Use SGD to update the parameters.Step 3.
Use the designed deep neural network to iteratively fit the training data, and then the probability distribution and information content of the training data are obtained. Optimize the encoder network according to Equation (6).Step 4.
The random sampling results fromN0,1 as well as σ and μ obtained from the encoder network are substituted into Equation (1) to calculate the latent result z. Once the training error reaches the thresholds, the network parameters of the encoder are saved.Step 5.
Take the latent resultz as the input of the decoder network, which is iteratively optimized according to Equation (10) to obtain z. The final decoding result from decoder is the final reconstruction result, and the network parameters of the decoder are saved.The flowchart of the above procedures is shown in Figure3.Figure 3
The flowchart of the proposed method.


## 4. Experimental Results and Analyses

In this section, the experiments were implemented with an Intel Core i7-9700k 4.1 GHz CPU, 16GB memory, and GeForce RTX2070s GPU with 8GB video RAM to evaluate the effectiveness of IVAE in the pore space reconstruction of shale by the metrics of porosity, variogram curves, permeability, distribution of pores, and CPU/GPU/memory performance. The experimental software framework for IVAE is Tensorflow-GPU.


### 4.1. Training Data

To estimate the effectiveness of the proposed method in the pore space reconstruction of shale, the real shale volume data with the resolution of 64 nanometers obtained by nano-CT scanning were used as the test data for the following tests. Figure4 shows three cross-sections (64×64 pixels) of the shale image with two facies: grain and pore space.Figure 4
Some shale cross-section images from real volume data.
(a)
#1 shale image(b)
#2 shale image(c)
#3 shale imageNote that in 2D conditions the basic image unit is pixel while in 3D conditions the unit of 3D images is voxel. Before applying any reconstruction methods, a 3D cube with64×64×64 voxels was used as the TI, extracted from the original shale volume data. Figures 5(a)–5(c) are the exterior (64×64×64 voxels), cross-sections (X=32, Y=32, and Z=32) and pore space of the TI (porosity=0.2690).Figure 5
TI of shale: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)


### 4.2. Reconstructions Using IVAE, VAE, SNESIM, and DS

In the section, IVAE, VAE, and two typical reconstruction methods—SNESIM and DS—were, respectively, used to reconstruct shale for comparison. The main parameters of IVAE are initialized as follows: learning rate is 0.001, the number of training epochs is 4000, batch size is 64,λz=λx′=0.1, and Fx′=Fz=5. Figure 6 is the reconstruction by IVAE. Figures 7–9 are, respectively, the reconstructed images using VAE, SNESIM, and DS. All the reconstructed images are 64×64×64 voxels. It can be found that the results of the four reconstruction methods all have similar long-connected pore spaces with the TI.Figure 6
A reconstructed image using IVAE: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 7
A reconstructed image using VAE: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 8
A reconstructed image using SNESIM: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 9
A reconstructed image using DS: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)


### 4.3. Comparison of Porosity

Porosity of shale indicates the shale’s ability to store fluids and is one of the characteristics for evaluating reconstruction quality. The definition of porosity 𝜙 is as follows:(11)ϕ=VpV,where V is the total volume and Vp is the pore volume.The porosities of the TI and the images shown in Figures6–9 are displayed in Table 1. It is seen that the porosity of the IVAE-reconstructed image is closest to that of the TI. To obtain an average performance in porosity, IVAE, VAE, SNESIM, and DS were performed for another 20 times to, respectively, achieve 20 reconstructed images (64×64×64 voxels). As shown in Table 2, the porosity of 3D images reconstructed by IVAE is closest to that of the TI, and the variance is smallest, indicating a high reconstruction quality and low fluctuation.Table 1
Porosities of the TI and the reconstructed images (shown in Figures6–9) using IVAE, SNESIM, VAE, and DS.
TISNESIMVAEDSIVAEPorosity0.26900.21010.23340.15430.2665Table 2
Average and variance of the porosities of 20 reconstructed shale images using IVAE, SNESIM, VAE, and DS.
SNESIMVAEDSIVAEAverage of porosity0.21780.23820.16820.2673Variance of porosity0.00240.00230.00270.0013


### 4.4. Comparison of Variogram

Variogram is widely used to evaluate the variability of spatial structures in different direction [29], which is defined as follows:
(12)γh=12EZx+h−Zx2,where E is the mathematical expectation, h is the lag distance between two positions: x+h and x, and Zx means the property value at the position x. If the variogram curves of two geological bodies are similar in a certain direction, it means that the spatial structures of the two geological bodies are similar in this direction; otherwise, it means that their structures are quite different. In this section, the variogram curves of TI and average variograms of 20 reconstructed images by IVAE, SNESIM, VAE, and DS were, respectively, calculated in the three directions of X, Y, and Z, as shown in Figure 10. In general, the variogram curves of IVAE are closest to those of the TI.Figure 10
The variogram curves of TI and average variogram curves of 20 reconstructed images by IVAE, SNESIM, VAE, and DS in three directions.
(a)
X direction(b)
Y direction(c)
Z directionIn order to quantitatively measure the difference of variogram curves between the TI and the reconstruction of each method, a difference degree (DD) is defined as follows:
(13)DDTI,re=∑h=1nXh−xh2,where TI and re, respectively, represent the TI and the reconstruction method and Xh and xh represent the variogram value of the TI and the reconstructed image when the distance is h. A smaller DD means a smaller difference between the TI and the reconstruction. DDs between the TI and the reconstructions of IVAE, SNESIM, VAE, and DS measured by variogram differences in the X, Y, and Z directions are shown in Table 3. It seems that reconstructed images by IVAE are closest to the TI.Table 3
Difference degrees (DDs) between the TI and the reconstructions of IVAE, SNESIM, VAE and DS measured by variogram differences in the X, Y, and Z directions.
DirectionsDifference degrees (DDs)DD (TI, SNESIM)DD (TI, DS)DD (TI, VAE)DD (TI, IVAE)X0.05690.07760.04450.0112Y0.07450.05790.04530.0159Z0.04710.02010.06450.0056


### 4.5. Permeability Estimation Using the Lattice Boltzmann Method

The permeability of shale means the ability of allowing fluid to pass through shale and often is related to porosity, geometry of pores in the direction of liquid penetration, and other factors. A 3D 19-velocity model of the Lattice Boltzmann Method (LBM), called D3Q19, was used to compute the permeabilities of reconstructed shale. There are altogether 19 velocity vectors defined in this model. The conditions of no-slip velocity are reached by the bounce-back scheme. Since the reconstructed results are cubic in our experiments, two parallel faces are open and the other four are sealed when supposing the fluids are perpendicularly passing through the two open parallel faces. Hence, there are altogether three directions for the computation of permeabilities using the Darcy’s law. The detailed procedures and equations can be referred to [29, 30].The data of the TI and the reconstructed images were, respectively, used as the input files of LBM simulation to calculate the permeabilities of those models with the size of64×64×64 voxels. As shown in Table 4, the permeability of the TI and the average permeabilities of 20 reconstructed images using IVAE, SNESIM, VAE, and DS in three directions were computed by LBM. The permeability values of the reconstructed images using IVAE are quite close to those of the TI, displaying the good reconstruction quality of IVAE.Table 4
Permeability of the TI and average permeabilities of 20 reconstructed images in three directions using IVAE, SNESIM, VAE, and DS.
DirectionPermeability (10-3 μm2)TISNESIMDSVAEIVAEX6.1235.5916.5826.2105.901Y5.2346.3545.8196.0215.721Z7.3188.7286.7899.1296.980


### 4.6. Distribution and Numbers of Pores

The diameter of a pore is defined as follows:(14)Diameter=6Vπ3,where V is the volume of a pore. The 3D shale images reconstructed by IVAE, SNESIM, VAE, and DS are imported into the software Avizo [31] to analyze the number of pores and the pore size. Table 5 shows the pore number in the TI and the average pore numbers in the reconstructions of these methods. Table 6 shows the maximum, minimum and average pore diameters in the TI and the reconstructions of IVAE, SNESIM, VAE, and DS. From Tables 5 and 6, it is seen that IVAE has shown the best reconstruction quality since the reconstruction of IVAE has the most similar pore number and average pore diameter with those of the TI.Table 5
The pore number in the TI and the average pore numbers in the reconstructions of IVAE, SNESIM, VAE, and DS.
TISNESIMDSVAEIVAEPore number206184178195204Table 6
The maximum, minimum ,and average pore diameters in the TI and the reconstructions of IVAE, SNESIM, VAE, and DS.
TISNESIMDSVAEIVAEAverage (voxel)5.875.386.215.695.83Min (voxel)3.384.934.793.613.45Max (voxel)23.3222.8318.9223.5222.88


### 4.7. Comparison of Reconstruction Time and CPU/GPU/Memory Usage

Since the pore space reconstruction of shale is quite CPU-intensive and normally uses large memory, the reconstruction time and CPU/GPU/memory usage of IVAE, SNESIM, VAE, and DS are compared in this section. As for the computational hardware platforms, these four methods are different. SNESIM and DS are typical CPU-based reconstruction methods, while IVAE and VAE can be performed by both CPU and GPU. Table7 shows the time consumption and the average usage of CPU, GPU, and memory by IVAE, SNESIM, VAE, and DS for 20 reconstructions.Table 7
The time consumption and average usage of CPU, GPU, and memory by IVAE, SNESIM, VAE, and DS for 20 reconstructions.
SNESIMDSVAEIVAEAverage CPU usage (%)23.133.737.539.8Average GPU usage (%)NoneNone62.865.0Average memory usage (MB)53273843984627Time for the first reconstruction (s)22391861167286627Time for the other 19 reconstructions (s)42602353292501490In Table7, the time for 20 reconstructions is divided into two parts: the time for the first reconstruction and the time for the other 19 reconstructions. The former means the time used for the first reconstruction of the traditional reconstruction methods (SNESIM and DS) and deep learning methods (VAE and IVAE). The latter is the time consumed for the other 19 reconstructions by IVAE, SNESIM, VAE, and DS.Deep learning methods normally benefit from the ability of saving model parameters during training. Therefore, the reconstruction process after the first model training that is finished in the first reconstruction only needs to reuse the saved parameters and takes much less time to reconstruct images. Hence, as typical deep learning methods, VAE and IVAE require much longer training time to set up the training model before reconstruction, usually consuming more time in the first reconstruction than the time used for the subsequent reconstruction. On the contrary, SNESIM and DS need to rescan the TI from scratch for each reconstruction because they only store training modes in memory. When a reconstruction process is over, training data of SNESIM and DS in memory will be cleared, resulting in the repeated processes of scanning the TI and extracting characteristics from the TI in the subsequent reconstruction.Since VAE and IVAE can save and reuse the parameters after the first reconstruction, each subsequent reconstruction only needs to input the parameters and complete the reconstruction very soon. Hence, as shown in Table7, although VAE and IVAE spend more time on the first reconstruction than SNESIM, the time for the subsequent reconstructions decreases largely. From an overall point of view, VAE and IVAE have great advantages in the speed of multiple reconstructions over the other CPU-based reconstruction methods. As for the overall performance of time consumption and average use of CPU/GPU/memory, VAE and IVAE are evenly matched and each one has its own strengths.


## 4.1. Training Data

To estimate the effectiveness of the proposed method in the pore space reconstruction of shale, the real shale volume data with the resolution of 64 nanometers obtained by nano-CT scanning were used as the test data for the following tests. Figure4 shows three cross-sections (64×64 pixels) of the shale image with two facies: grain and pore space.Figure 4
Some shale cross-section images from real volume data.
(a)
#1 shale image(b)
#2 shale image(c)
#3 shale imageNote that in 2D conditions the basic image unit is pixel while in 3D conditions the unit of 3D images is voxel. Before applying any reconstruction methods, a 3D cube with64×64×64 voxels was used as the TI, extracted from the original shale volume data. Figures 5(a)–5(c) are the exterior (64×64×64 voxels), cross-sections (X=32, Y=32, and Z=32) and pore space of the TI (porosity=0.2690).Figure 5
TI of shale: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)


## 4.2. Reconstructions Using IVAE, VAE, SNESIM, and DS

In the section, IVAE, VAE, and two typical reconstruction methods—SNESIM and DS—were, respectively, used to reconstruct shale for comparison. The main parameters of IVAE are initialized as follows: learning rate is 0.001, the number of training epochs is 4000, batch size is 64,λz=λx′=0.1, and Fx′=Fz=5. Figure 6 is the reconstruction by IVAE. Figures 7–9 are, respectively, the reconstructed images using VAE, SNESIM, and DS. All the reconstructed images are 64×64×64 voxels. It can be found that the results of the four reconstruction methods all have similar long-connected pore spaces with the TI.Figure 6
A reconstructed image using IVAE: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 7
A reconstructed image using VAE: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 8
A reconstructed image using SNESIM: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)Figure 9
A reconstructed image using DS: (a) exterior; (b) cross-sections (X=32, Y=32, and Z=32); (c) pore space.
(a)(b)(c)


## 4.3. Comparison of Porosity

Porosity of shale indicates the shale’s ability to store fluids and is one of the characteristics for evaluating reconstruction quality. The definition of porosity 𝜙 is as follows:(11)ϕ=VpV,where V is the total volume and Vp is the pore volume.The porosities of the TI and the images shown in Figures6–9 are displayed in Table 1. It is seen that the porosity of the IVAE-reconstructed image is closest to that of the TI. To obtain an average performance in porosity, IVAE, VAE, SNESIM, and DS were performed for another 20 times to, respectively, achieve 20 reconstructed images (64×64×64 voxels). As shown in Table 2, the porosity of 3D images reconstructed by IVAE is closest to that of the TI, and the variance is smallest, indicating a high reconstruction quality and low fluctuation.Table 1
Porosities of the TI and the reconstructed images (shown in Figures6–9) using IVAE, SNESIM, VAE, and DS.
TISNESIMVAEDSIVAEPorosity0.26900.21010.23340.15430.2665Table 2
Average and variance of the porosities of 20 reconstructed shale images using IVAE, SNESIM, VAE, and DS.
SNESIMVAEDSIVAEAverage of porosity0.21780.23820.16820.2673Variance of porosity0.00240.00230.00270.0013


## 4.4. Comparison of Variogram

Variogram is widely used to evaluate the variability of spatial structures in different direction [29], which is defined as follows:
(12)γh=12EZx+h−Zx2,where E is the mathematical expectation, h is the lag distance between two positions: x+h and x, and Zx means the property value at the position x. If the variogram curves of two geological bodies are similar in a certain direction, it means that the spatial structures of the two geological bodies are similar in this direction; otherwise, it means that their structures are quite different. In this section, the variogram curves of TI and average variograms of 20 reconstructed images by IVAE, SNESIM, VAE, and DS were, respectively, calculated in the three directions of X, Y, and Z, as shown in Figure 10. In general, the variogram curves of IVAE are closest to those of the TI.Figure 10
The variogram curves of TI and average variogram curves of 20 reconstructed images by IVAE, SNESIM, VAE, and DS in three directions.
(a)
X direction(b)
Y direction(c)
Z directionIn order to quantitatively measure the difference of variogram curves between the TI and the reconstruction of each method, a difference degree (DD) is defined as follows:
(13)DDTI,re=∑h=1nXh−xh2,where TI and re, respectively, represent the TI and the reconstruction method and Xh and xh represent the variogram value of the TI and the reconstructed image when the distance is h. A smaller DD means a smaller difference between the TI and the reconstruction. DDs between the TI and the reconstructions of IVAE, SNESIM, VAE, and DS measured by variogram differences in the X, Y, and Z directions are shown in Table 3. It seems that reconstructed images by IVAE are closest to the TI.Table 3
Difference degrees (DDs) between the TI and the reconstructions of IVAE, SNESIM, VAE and DS measured by variogram differences in the X, Y, and Z directions.
DirectionsDifference degrees (DDs)DD (TI, SNESIM)DD (TI, DS)DD (TI, VAE)DD (TI, IVAE)X0.05690.07760.04450.0112Y0.07450.05790.04530.0159Z0.04710.02010.06450.0056


## 4.5. Permeability Estimation Using the Lattice Boltzmann Method

The permeability of shale means the ability of allowing fluid to pass through shale and often is related to porosity, geometry of pores in the direction of liquid penetration, and other factors. A 3D 19-velocity model of the Lattice Boltzmann Method (LBM), called D3Q19, was used to compute the permeabilities of reconstructed shale. There are altogether 19 velocity vectors defined in this model. The conditions of no-slip velocity are reached by the bounce-back scheme. Since the reconstructed results are cubic in our experiments, two parallel faces are open and the other four are sealed when supposing the fluids are perpendicularly passing through the two open parallel faces. Hence, there are altogether three directions for the computation of permeabilities using the Darcy’s law. The detailed procedures and equations can be referred to [29, 30].The data of the TI and the reconstructed images were, respectively, used as the input files of LBM simulation to calculate the permeabilities of those models with the size of64×64×64 voxels. As shown in Table 4, the permeability of the TI and the average permeabilities of 20 reconstructed images using IVAE, SNESIM, VAE, and DS in three directions were computed by LBM. The permeability values of the reconstructed images using IVAE are quite close to those of the TI, displaying the good reconstruction quality of IVAE.Table 4
Permeability of the TI and average permeabilities of 20 reconstructed images in three directions using IVAE, SNESIM, VAE, and DS.
DirectionPermeability (10-3 μm2)TISNESIMDSVAEIVAEX6.1235.5916.5826.2105.901Y5.2346.3545.8196.0215.721Z7.3188.7286.7899.1296.980


## 4.6. Distribution and Numbers of Pores

The diameter of a pore is defined as follows:(14)Diameter=6Vπ3,where V is the volume of a pore. The 3D shale images reconstructed by IVAE, SNESIM, VAE, and DS are imported into the software Avizo [31] to analyze the number of pores and the pore size. Table 5 shows the pore number in the TI and the average pore numbers in the reconstructions of these methods. Table 6 shows the maximum, minimum and average pore diameters in the TI and the reconstructions of IVAE, SNESIM, VAE, and DS. From Tables 5 and 6, it is seen that IVAE has shown the best reconstruction quality since the reconstruction of IVAE has the most similar pore number and average pore diameter with those of the TI.Table 5
The pore number in the TI and the average pore numbers in the reconstructions of IVAE, SNESIM, VAE, and DS.
TISNESIMDSVAEIVAEPore number206184178195204Table 6
The maximum, minimum ,and average pore diameters in the TI and the reconstructions of IVAE, SNESIM, VAE, and DS.
TISNESIMDSVAEIVAEAverage (voxel)5.875.386.215.695.83Min (voxel)3.384.934.793.613.45Max (voxel)23.3222.8318.9223.5222.88


## 4.7. Comparison of Reconstruction Time and CPU/GPU/Memory Usage

Since the pore space reconstruction of shale is quite CPU-intensive and normally uses large memory, the reconstruction time and CPU/GPU/memory usage of IVAE, SNESIM, VAE, and DS are compared in this section. As for the computational hardware platforms, these four methods are different. SNESIM and DS are typical CPU-based reconstruction methods, while IVAE and VAE can be performed by both CPU and GPU. Table7 shows the time consumption and the average usage of CPU, GPU, and memory by IVAE, SNESIM, VAE, and DS for 20 reconstructions.Table 7
The time consumption and average usage of CPU, GPU, and memory by IVAE, SNESIM, VAE, and DS for 20 reconstructions.
SNESIMDSVAEIVAEAverage CPU usage (%)23.133.737.539.8Average GPU usage (%)NoneNone62.865.0Average memory usage (MB)53273843984627Time for the first reconstruction (s)22391861167286627Time for the other 19 reconstructions (s)42602353292501490In Table7, the time for 20 reconstructions is divided into two parts: the time for the first reconstruction and the time for the other 19 reconstructions. The former means the time used for the first reconstruction of the traditional reconstruction methods (SNESIM and DS) and deep learning methods (VAE and IVAE). The latter is the time consumed for the other 19 reconstructions by IVAE, SNESIM, VAE, and DS.Deep learning methods normally benefit from the ability of saving model parameters during training. Therefore, the reconstruction process after the first model training that is finished in the first reconstruction only needs to reuse the saved parameters and takes much less time to reconstruct images. Hence, as typical deep learning methods, VAE and IVAE require much longer training time to set up the training model before reconstruction, usually consuming more time in the first reconstruction than the time used for the subsequent reconstruction. On the contrary, SNESIM and DS need to rescan the TI from scratch for each reconstruction because they only store training modes in memory. When a reconstruction process is over, training data of SNESIM and DS in memory will be cleared, resulting in the repeated processes of scanning the TI and extracting characteristics from the TI in the subsequent reconstruction.Since VAE and IVAE can save and reuse the parameters after the first reconstruction, each subsequent reconstruction only needs to input the parameters and complete the reconstruction very soon. Hence, as shown in Table7, although VAE and IVAE spend more time on the first reconstruction than SNESIM, the time for the subsequent reconstructions decreases largely. From an overall point of view, VAE and IVAE have great advantages in the speed of multiple reconstructions over the other CPU-based reconstruction methods. As for the overall performance of time consumption and average use of CPU/GPU/memory, VAE and IVAE are evenly matched and each one has its own strengths.


## 5. Summary and Conclusions

The properties of shale such as low porosity, low permeability, and complex inner structures are the main reasons challenging the exploitation of shale reservoirs. The establishment of a 3D pore space model of shale can help to analyze the characteristics of shale, improving the producing efficiency of shale resources. Main conclusions, including the advantages and disadvantages (or limitations) of our method, are as follows:(1)
A real 3D shale cube acquired by nano-CT scanning was used as the TI, providing the necessary real structural information of pore space for pore space reconstruction, so the reconstructed structures have similar structures with the real shale(2)
Traditional numerical reconstruction methods need to repeatedly scan the TI to extract the statistical information in reconstruction, leading to consuming more time in repeated reconstructions. The proposed method shows the advantage in shale reconstruction in terms of both speed and quality by reusing the saved models after the first reconstruction or training(3)
Our method combines VAE and Fisher information together for the pore space reconstruction of shale, using new penalty terms related to Fisher information in the encoder and the decoder, respectively, so our method performs better than the original VAE in the reconstruction quality, also verified by the real experiments(4)
Our method still has some disadvantages or limitations. First, compared to traditional reconstruction methods, our method has higher CPU/memory usage and consumes much more time in the first reconstruction. Second, our method is established on the framework of Tensorflow-GPU, so a GPU is necessary for our method, increasing the hardware cost of reconstruction



---

*Source: 5545411-2021-03-18.xml*