MDA GAN: Adversarial-Learning-based 3-D Seismic Data Interpolation and Reconstruction for Complex Missing

Our approach is based on GAN, and in V Appendix we add a basic explanation for GAN and the current underlying framework for GAN-based image inpainting in computer vision.

The major disadvantage of current 2D seismic interpolation is that reconstruction 2D slices one by one leads to discontinuities in brutal stacked section[12] when there are complex cases such as large scale or continuous missing in shot gathers. To solve this problem, we extend the seismic interpolation and reconstruction to three dimensions.

The generator is a 3D CNN structure, and the input and output are 3D data of the same size. This article describes the structure of the generator and its components detailed in II-B.

The 3D seismic data have obvious anisotropy, where the seismic reflection axis is approximately perpendicular to the inline and crossline slices, and a clear seismic layered texture can be observed in these two directions, while the timeline direction shows irregular texture. Therefore, we use two 2D discriminators $\mathcal{D}^{m}_{\theta_{D}^{\text{xi}}}$ and $\mathcal{D}^{m}_{\theta_{D}^{\text{t}}}$, inline and crossline share discriminator $\mathcal{D}^{m}_{\theta_{D}^{\text{xi}}}$, and timeline direction uses $\mathcal{D}^{m}_{\theta_{D}^{\text{t}}}$ discriminator. 2D discriminators are employed for each directional slice to maintain the anisotropy of the restored seismic data, and we also use a 3D discriminator $\mathcal{D}^{t}_{\theta_{D}}$ to ensure the continuity of the results in 3D space.

The framework consists of a 3D generator, two 2D discriminators and a 3D discriminator. Accordingly, we use three adversarial losses $\mathcal{L}_{\mathcal{G}^{t}}$, $\mathcal{L}_{\mathcal{G}^{m}_{\text{xi}}}$ and $\mathcal{L}_{\mathcal{G}^{m}_{\text{t}}}$. In combination with the reconstruction (V-B), the loss of the supervised generator can be expressed as equation (1).

where $\mathcal{L}_{\text{rec}}$ is the reconstruction loss, which we will describe in detail in II-D, $\lambda$ is a scaling factor to adjust the reconstruction loss to the same order of magnitude as the adversarial loss. Each adversarial loss is expressed as equation (2).

where $\textbf{I}_{m}$ is the missing input data. This loss tries to make the generator fool the discriminator so that the generated hidden variables are approximated to the real data.

Each discriminator determines the authenticity of the generated tensor or matrix, using the cross-entropy loss for supervision, equation (3).

For reconstructing high-quality seismic data, our approach employs adversarial learning in multiple dimensions, hence the term Multi-Dimensional Adversarial (MDA) GAN. The overall framework of MDA GAN is shown in Fig.1-(a).

The discriminator network uses the structure of a standard encoder. The discriminator consists of five convolutional layers, and each layer employs the LeakyReLU[37] activation function, and also incorporates the spectral normalization[44] to stabilize the training of the GAN.

Most of the current reconstruction losses are based on $L_{1}$ and $L_{2}$ [45, 46, 47, 48], where $L_{1}$ provides stable gradients and $L_{2}$ has the more stable solution, but both loss functions have the common disadvantage of unsmoothed gradients in GANs (equations (10) and (13)). Generator networks generally use Tanh as the final activation layer to keep the network’s value domain at $[-1,1]$, while ensuring training stability[37]. Next we derive the gradient expressions for $L_{1}$ and $L_{2}$ under the Tanh activation function.

First, we derive the gradient of $L_{1}$ loss back propagation to backbone.

The most reasonable relationship between $\partial\mathcal{L}_{L_{1}}/\partial\hat{x}_{t}$ and $\hat{y}_{t}$ should be the one that triggers the minimum gradient (minimum absolute value) when $\hat{y}_{t}=y_{t}$. As $\hat{y}_{t}$ and $y_{t}$ become farther apart, the greater the gradient they cause should be, and this increasing relationship should be linear, i.e., smooth, its slope should not change with $\hat{y}_{t}$. We define this optimal assumption as the equation (11).

Fig. 4 shows a sharp fluctuation in the neighborhood of $L_{1}$ loss when $\hat{y}_{t}=y_{t}$. Although it provides the correct gradient direction, which leads to it being trainable, but its values are coarse and the incremental relationship is nonlinear, so it obviously does not satisfy the above optimality assumption.

Equations (12) and (13) illustrate the gradient form of $L_{2}$ back propagation to backbone.

In Fig. 5, the three curves before reaching the extreme value point execute a relatively reasonable trend, although it is not smooth, the magnitude of their absolute values is proportional to the distance from the point $\hat{y}_{t}=y_{t}$. However, when they are in the position after the extreme point, i.e. the area marked by the red line segment, they show an unreasonable trend and the absolute values that should be large are becoming reduced. The $L_{2}$ loss forms a parabola between $\hat{y}_{t}=y_{t}$ and $\hat{y}_{t}=1$. Although this can still be trained and does not affect the propagation of the gradient, we still hope to find a more reasonable gradient induced by the loss.

Both the $L_{1}$ and $L_{2}$ losses are deficient in the form of the gradient passed to the backbone. We give the optimal gradient form assumption in equation (11). We want $\varepsilon$ to be directly used as the gradient form of the loss to backbone, and the resulting loss will satisfy the optimal gradient assumption, so we integrate it. Then the loss function we expect can be expressed by the equation (14).

The loss of all voxels should be averaged during training, so that the loss expression is (15).

Its form is similar to cross-entropy, so we named it Tanh Cross-Entropy (TCE). This loss causes the gradient at training as in equation (16).

In Fig. 6, TCE shows the smoothed gradient of its feedback to backbone. When $\hat{y}_{t}=y_{t}$, the gradient is $0$ and the distance between $\hat{y}_{t}$ and $y_{t}$ is linear with the gradient size, which is the optimal loss we want to obtain.

Although $L_{1}$ and $L_{2}$ can give correct feedback to the generator in the gradient direction, their gradient values are coarse and they do not reasonably reflect the distance between GT and the predicted value, while the gradient passed from TCE to the generator backbone is the residual of GT and the predicted value, which perfectly captures the change of the distance between the output and the target.

The over-arching training framework follows the Fig. 1. The optimizer uses Adam[49]. In general, the training of the discriminator should be ahead of the generator, which is more conducive to producing high-quality data. We use the TTUR strategy for training [50], setting the learning rate of discriminators $\eta_{d}=0.0004$ and the generator $\eta_{g}=0.0001$, finally alternating the training of the generator and discriminators.

We designed three control groups and one experimental group using the current most common UNet as the baseline model (Group 1), UNet has excellent performance in interpolating seismic data, which has been demonstrated in several works ([18, 20, 22] etc.). Group 2 is the replacement of the backbone of Group 1 with the generator we proposed in Fig. 1. Group 2 does not apply adversarial learning and therefore does not require tanh activation function and TCE loss, this group can verify the validity of the generator and the FSM compared to Group 1. Group 3 incorporates multidimensional adversarial learning, which differs from the experimental group in that it uses $L_{1}$ as the reconstruction loss. This group can be compared with Group 2 to verify the effectiveness of multidimensional adversarial and with the experimental group to verify the effectiveness of TCE loss. Group 4 is the experimental group, which uses exactly the method proposed in this paper. The settings of the four groups of experiments are shown in Table I.

The current experimental approach of most data-driven seismic interpolation studies is to train a new set of weights for each survey [5, 17, 18, 19, 20, 21, 22, 25, 22, 27, 33, 34], but retraining the neural network is extremely time-consuming and resource-wasting, resulting in a much higher time and resource cost than the theory-driven approach. Furthermore, the traces of the gathers are often highly similar to each other, so the neural network weights obtained from training for a specific survey may be overfitted, and this overfitting tends to show extremely high performance for the test set divided by the same survey, so this way of dividing the data set may not objectively reflect the performance of the model.

In this study, we want to train a general model that can be applied to the reconstruction of most seismic data without having to retrain new weights for different surveys, which is the current trend in AI development, CV and NLP scientists are currently working on building general models for various tasks[51, 52]. We have used multiple surveys (both pre-stack and post-stack) to jointly train one model in the expectation that it will generalize to more data, we have made the final model publicly available on Github¹¹1Available: https://github.com/douyimin/MDA_GAN.

The usage of the data is shown in Table II.

All training data are cropped to $128\times 128\times 128$, because the generator is fully convolutional, so the validation and test sets can be of arbitrary size. During training, 30%-99% of slices are randomly removed in the training cuboid as input data, the original data are used as supervised information, and the whole process does not require manual annotation. The training batch size is 4, and it is trained on two RTX3090Ti with 500000 steps. The model with the highest SSIM (details in III-C) metric on the validation set was selected as the final model and tested.

Next, we present and analyze the qualitative and quantitative results for the test set in each survey separately, and then present the combined quantitative metrics on the test set for all surveys. The trace of test data deletion in this paper’s experiments were all performed using python’s random library, and all random seeds were ’2022’.

All experiments in this paper were evaluated quantitatively using Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity (SSIM) metrics[53], and all data were normalized to $[0,1]$ prior to evaluation. To facilitate comparison and reproduction, the quantitative metrics are calculated directly using the public library scikit-image’s implementation²²2Available: https://scikit-image.org of PSNR and SSIM, and all hyperparameters in SSIM are default.

In SEG C3, we stacks gather into 3D, as shown in Table II 3/5 of the data were involved in training, 1/5 were used for validation, and we used the last 1/5 for testing. Certainly, there are other surveys involved in the training.

We present it using one of the $192\times 240\times 64$ grid points (Fig. 7), which is obtained by interpolating (resize) the original data $168\times 205\times 31$.

We fold the prestack gather into 3D, and intercept $512\times 256\times 128$ grid points in the test data as a sample for qualitative analysis, as shown in Fig. 10.

Fig. 11 shows the qualitative interpolation results of MDA GAN and UNet, and Table V shows the corresponding quantitative results.

At 50% discrete deficiency, both MDA GAN and UNet are able to perform the interpolation task well. When the ratio comes to 75%, UNet starts to be overwhelmed, while our method still shows high reliability. With more complex 90% and 95% missing, UNet invalidates, while our method still yields promising results. In the case of consecutive missing cases, all interpolation results of UNet are failures. Fig. 11 shows that the qualitative results of our method in complex cases are significantly higher than those of UNet, which is currently the most used, and the same conclusion can be obtained in the quantitative results in Table V.

This work focuses on the interpolation and reconstruction of complex cases. The above experiments validate the reconstruction performance of our method under complex missing for the set of prestack gather. The purpose of this subsection is to verify the performance of our method with complex surveys rich in anomalies such as faults and salt bodies.

The missing in original engineering often occurs in shot gathers, and anomalous bodies also appear in brutal stack sections [12] formed by shot gathers, but the post-stack data containing anomalous bodies is more abundant in the current public data, so we use the post-stack as the test data in this subsection.

Finally we calculate the SSIM and PSNR metrics jointly for all the test sets shown in Table II. Because the size of the test parts differs for each survey, we calculate only the proportional discrete missing cases, the results are shown in Table VI. For a more visual analysis of the performance of each group of experiments, we organized Table VI into Fig. 14.

In Fig. 14 (a), the MDA-based metrics are remarkably higher than the other two groups, indicating the effectiveness of multidimensional adversarial. While using $L_{1}$ and TCE losses do not differ significantly in terms of SSIM, with $L_{1}$ slightly higher than TCE, which is one of the reasons why L1 is currently widely used as a reconstruction loss, but the experiments from above show that the qualitative results of TCE losses are more advanced (Fig. 8 and 9), and it can provide smoother and continuous reconstruction results. In (b), the PSNR of MDA GAN (TCE) was higher than the other groups. Similar to (a), there was no significant difference between TCE and $L_{1}$, and TCE was slightly higher than $L_{1}$ in terms of PSNR metrics. This may be because PSNR is more sensitive to pixel differences, while SSIM prefers to show structural differences in regions.

Generative Adversarial Network (GAN) are proposed by Goodfellow et al.[54]. It consists of a generator network $\mathcal{G}_{\theta_{G}}$ and a discriminator network $\mathcal{D}_{\theta_{D}}$, where the task of the generator is to yield images $x\in\mathcal{R}^{N_{w}\times N_{h}}$ with a latent noise prior vector, $z\in\mathcal{R}^{d}$ as input $z$ is sampled from a known distribution, i.e. latent vector $z$, $z\sim\mathcal{U}[-1,1]^{d}$ [54]. The task of the discriminator is to distinguish real images from generated ones. The generator and the discriminator play a zero-sum game in which the two networks learn from each other and the data obtained by the generator becomes increasingly close to the real data, so that the desired data can be generated. The objective function of GANs can be expressed by the equation (19).

where $p_{data}(x)$ represents the distribution of the training set and $p_{z}(z)$ represents the normal distribution of the noise, $\mathcal{D}_{\theta_{D}}(\mathcal{G}_{\theta_{G}}(z))$ represents the generated data. $\mathcal{G}_{\theta_{G}}$ expects $\mathcal{D}_{\theta_{D}}(\mathcal{G}_{\theta_{G}}(z))$ to be as large as possible, when $V(\mathcal{G}_{\theta_{G}},\mathcal{D}_{\theta_{D}})$ becomes smaller, so equation is minimized for $\mathcal{G}_{\theta_{G}}$. $\mathcal{D}_{\theta_{D}}(x)$ is enhanced with the discriminator performance, when $\mathcal{D}_{\theta_{D}}(\mathcal{G}_{\theta_{G}}(z))$ becomes smaller and $V(\mathcal{G}_{\theta_{G}},\mathcal{D}_{\theta_{D}})$ becomes larger, so equation is maximized for $\mathcal{D}_{\theta_{D}}$.

The key to applying GAN to image inpainting is to join discriminator and adversarial loss on the base of encoder regression [28, 29, 30, 31]. The existing image inpainting work is for 2D images, its loss can be expressed as equation (20), the basic framework is shown in Fig 15.

where $\mathcal{L}_{\text{adv}}$ is adversarial loss, commonly used CE-GAN[54], LS-GAN [55], W-GAN[56], etc., $\lambda\mathcal{L}_{\text{rec}}$ is the reconstruction loss (regression loss), $L_{1}$ and $L_{2}$ are generally used [28, 29, 30, 31].

Most of the current GAN-based image inpainting follows this framework.