Conditional Denoising Diffusion Probabilistic Model for Ground-roll Attenuation

For conventional c-DDPM, the target data distribution should be either clean data or ground roll data. Once one of them is obtained by the c-DDPM, the other one can be obtained by subtracting it from the noisy data. In the conventional c-DDPM, we typically set the clean data $\mathbf{x}_{0}$ as the target distribution and the noisy data $\mathbf{y}$ as the condition. Correspondingly, the ground roll $\mathbf{z}_{0}$ can be obtained by subtracting the generated clean data from the noisy data. The forward process for the c-DDPM is defined on a Markov chain within $T$ timesteps:

where $\mathbf{x}_{0}$ is sampled from the training datasets $q(\mathbf{x}_{0}|\mathbf{y})$, $\mathbf{x}_{1},..,\mathbf{x}_{T}$ are the latent variables fixed on the Markov chain. $q(\mathbf{x}_{t}|\mathbf{x}_{t-1})$ is the pre-defined Gaussian translation and obtains $\mathbf{x}_{t}$ by adding random noise to $\mathbf{x}_{t-1}$:

where $\beta_{t}\in(0,1)$ is a predefined variance schedule that increases with timestep $t$. The property of the Markov chain enables:

A more explicit form of Eq. 3 is:

where ${\alpha_{t}}=1-\beta_{t}$ and ${\bar{\alpha}}_{t}=\prod_{s=1}^{t}{\alpha}_{s}$. At timestep $T$, $\mathbf{x}_{t}$ converges to a prior distribution, i.e, a standard normal distribution. Thus, the forward process is a diffusion process that gradually converts a clean image to a pure noise image.

The reverse process is also defined on a Markov chain, which converts pure noise $\mathbf{x}_{T}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ to the data distribution $\mathbf{x_{0}}$, under the guidance of $\mathbf{y}$:

where the Gaussian translation $p_{\theta}(\mathbf{x}_{t-1}|\mathbf{x}_{t},\mathbf{y})$ has learned mean $\bm{\mu}_{\theta}$ and fixed variance $\mathbf{\Sigma}_{\theta}={\sigma_{t}}^{2}$:

For the detailed inference of the formulas, see [36]. The optimization function of the network is:

After training the network, we can gradually generate the clean data $\hat{\mathbf{x}_{0}}$ with the sampling step in c-DDPM:

where $\hat{\mathbf{x}}_{t}$ is the estimated latent variable in the sampling process.

One limit of the conventional c-DDPM is that the ground roll $\hat{\mathbf{z}}_{0}$ obtained by subtracting the predicted clean data $\hat{\mathbf{x}}_{0}$ from the noisy data $\mathbf{y}$ usually suffer from some leakage of reflections because of the inaccurate prediction of $\hat{\mathbf{x}}_{0}$. Thus, we design two target distribution, clean data and ground roll, in the improved c-DDPM given the noisy data as condition. The forward and reverse processes for the improved c-DDPM are shown in Fig. 1. In detail, the forward process in the improved c-DDPM is

where $\mathbf{x}_{0}$ and $\mathbf{z}_{0}$ refer to the clean data and ground roll, respectively. The reverse process is formulated as:

Correspondingly, the optimization function becomes

where $\boldsymbol{\epsilon}_{\theta,\mathbf{x}}(\mathbf{x}_{t},\mathbf{z}_{t},\mathbf{y},t)$ and $\boldsymbol{\epsilon}_{\theta,\mathbf{z}}(\mathbf{x}_{t},\mathbf{z}_{t},\mathbf{y},t)$ represents the estimated noise for sampling $\mathbf{x}_{t-1}$ and $\mathbf{z}_{t-1}$ respectively. Note that, both the noise are jointly estimated from $\mathbf{x}_{t}$ and $\mathbf{z}_{t}$ with the network. Once the training process is finished, we can gradually generate the clean data $\hat{\mathbf{x}}_{t}$ and the ground roll $\hat{\mathbf{z}}_{t}$ with the noisy data as condition:

The proposed training process and sampling process of the improved c-DDPM are shown in Algorithm 1 and Algorithm 2, respectively.

We first test the proposed methods DDPM-1c and DDPM-2c on the synthetic data example. Besides, we also apply the LTF and U-Net to the synthetic data as comparison. Fig. 5 shows the 2D slice of the SEAM Arid model. To provide the benchmark for comparison, we can directly obtain the clean data and the ground roll by applying the elastic finite-difference modeling method to the modified velocity model. We use the $v_{p}$ and $\rho$ model shown in Fig. 5a and 5c but set $v_{s}$ to zero to simulate the clean data, while we retain the near-surface part above 600 m of Fig. 5 but make the model beneath 600 m homogeneous to simulate the ground roll. The 25-Hz Ricker wavelet is used to simulate these seismic data. The simulated clean data and ground roll are shown in Fig. 6b and Fig. 6c, respectively. The simulated data has a total of 640 traces and 224 time samples per trace, with a time sampling interval of 0.004 s. We first apply automatic gain control (AGC) to these simulated data. Then, we add the weighted ground-roll data to the clean data to obtain the noisy data shown in Fig. 6a. The weighting factor for the ground roll can be used to adjust the SNR of noisy data. We can see from Fig. 6a that the coherent ground-roll noise corrupts the shallow reflections at near offsets and deeper reflections at far offsets. We use our proposed method (DDPM-1c, DDPM-2c), LTF and U-Net to remove the ground roll from the noisy data. Fig. 6d and 6e show the clean data and ground roll obtained by the LTF method, respectively. We can obviously see some ground-roll noise in the denoising result. The clean data and ground roll corresponding to the U-Net method are shown in Fig. 6f and 6g. We can still see that the ground roll is not removed completely, especially for the part pointed by the red box. Besides, the energy leakage of reflections is observed in Fig. 6g. The clean data predicted by using DDPM-1c is shown in Fig. 6h. As pointed by the red box, the ground-roll energy leakage is relatively weak. The separated ground-roll data by using DDPM-1c is shown in Fig. 6i. We can see that the energy leakage of reflections is less visible compared with LTF and U-Net. Fig. 6j and 6k show the predicted clean data and ground roll with DDPM-2c. We can see that the predicted clean data is comparable to the simulated clean with very slight contamination from ground roll. As pointed by the red box, the ground roll is completely removed, and the masked reflections are recovered well. In addition, there is obviously less remaining reflection energy in the predicted ground roll compared with the result of DDPM-1c. The comparison between these figures tell us that the proposed DDPM-2c shows better performance in ground-roll attenuation than LTF, U-Net and DDPM-1c.

We apply the F-K plane analysis to further evaluate the performance of the ground-roll attenuation methods. The F-K spectra corresponding to the data in Fig. 6a and Fig. 6b are shown in Fig. 7a and 7b, respectively. From the spectra, we observe that the clean data and ground roll have overlap in the frequency band of 20 - 35Hz. Fig. 7c displays the F-K spectra of the predicted clean data in Fig. 6d, containing some residual ground-roll energy. Fig. 7d, 7e, and 7f show the F-K spectra corresponding to the predicted clean data in Fig. 6f, 6h, and 6j, respectively. The F-K spectra corresponding to DDPM-2c has less residual ground roll and looks closer to that of the simulated clean data.

Moreover, we utilize amplitude spectra to assess the performance of these methods for attenuating ground roll in noisy data. The amplitude spectra of the true clean data and the predicted clean data in Fig. 6 are illustrated in Fig. 8a. We compute the spectra of these data over the same time and space window. Obviously, the U-Net result has an increased energy from 10 to 40 Hz, arising from the inaccurate amplitude of the predicted clean data. On the other hand, the clean data predicted by LTF, DDPM-1c and DDPM-2c have consistent frequencies with the true clean data. Furthermore, the amplitude spectra of the true noise data and the separated ground-roll data in Fig. 6 are shown in Fig. 8b. We can see that amplitude spectra corresponding to DDPM-2c shows the best consistency with the target spectra. Although the prediction of DDPM-1c slightly degrades, it is still better than the LTF and U-Net methods.

For a quantitation assessment, we then use four evaluation metrics, including Mean absolute error (MAE), mean square error (MSE), peak signal-to-noise ratio (PSNR) and structural similarity index measure (SSIM) to measure the difference between the ground truth and predicted results. You are referred to Appendix A for more details about these evaluation metrics. The lower the values of MAE and MSE, and the higher the values of PSNR and SSIM, the closer the predicted data is to the target. Considering that the DDPMs has the inherent generative diversity, we perform the sampling process for the DDPM-based methods five times and compute the mean value of the metrics. The evaluation results for the clean data and ground roll predicted by different methods are summarized in Table. I and II, respectively. We can see from Table. I that our proposed methods (DDPM-1c and DDPM-2c) perform better than the LTF and U-Net methods in predicting clean data. Meanwhile, Table. II illustrates that DDPM-2c has superior performance in predicting ground roll.

We then apply the c-DDPM method to the field data with no needs to re-train the model. Fig. 9a shows a shot gather from the Line 001 dataset after AGC and down-sampling. The trace number is 640, and the recording length is 224 with a time sampling of 0.0047 s. The ground-roll noise is predominantly located at near offset and seriously contaminates the reflection events. The denoising result of the LTF method is shown in Fig. 9b. We can see that the reflections are still masked by strong ground-roll residuals. The ground roll corresponding to the LTF method is shown in Fig. 9c. There is obvious leakage of reflections at the masked area of the ground roll and reflections. The predicted clean data using U-Net is shown in Fig. 9d. As shown by the red box, the reflections are partially recovered, but the ground roll still obviously contaminates the seismic data. Fig. 9e shows the corresponding ground roll. We can see the ground roll is more complete, while there are some leakage of reflections as pointed by the red arrows. The clean data predicted by DDPM-1c and DDPM-2c are shown in Fig. 9f and Fig. 9h, respectively. As pointed by the red box, the reflections are recovered better than the LTF and U-Net methods. The predicted ground roll by using DDPM-1c and DDPM-2c are shown in Fig. 9g and Fig. 9i, respectively. As pointed out by the red arrows, we can see the slight reflection energy in the result of DDPM-1c, while the result of DDPM-2c has no visible reflections. These figures show that the DDPM-2c has better performance on the field data than the LTF, U-Net and DDPM-1c methods.

We also show the f-k spectra of the original seismic data and the denoising results of different methods in Fig. 10. The amplitude spectra of the original seismic data and denoising results are shown in Fig. 11a. We can see that the LTF does not attenuate the ground roll with frequency from 8 Hz to 15 Hz effectively. The amplitude spectra of the original seismic data and the predicted ground roll in Fig. 9 are shown in Fig. 11b. We can see that DDPM-1c and DDPM-2c predicts the ground roll with frequency from 10 to 20 Hz very well.

It is more challenging when the energy of ground roll is much stronger than the reflections in the seismic data. Thus, we perform AGC on the noisy input in the synthetic and field data examples, to balance the amplitudes across the whole data. AGC adjust the amplitudes of seismic data in a window according to the mean energy within the window. The reflections overlapped by strong ground-roll noise will have smaller gain than the reflections unaffected by the ground-roll noise. Thus, the predicted reflections in the strongly-affected area tends to be weaker than the predicted reflections in the unaffected area. This can be seen in Fig. 9f and 9h.

There are two possible solutions to this issue:

1) Improve the network so that it can work well when the ground roll energy is dominant in the noisy data. DDPMs are deep generative models originally made for image synthesis. The network architecture for a range of 0-255 RGB of an image should be carefully modified for the high-resolution task such as seismic processing.

2) Reconstruct the clean data instead of separating it from the noisy data after AGC. Both the DDPM-1c and DDPM-2c are required to separate the reflection energy from the noisy condition after AGC. The network can learn to reconstruct the contaminated area with continuous reflection energy by lessening the constraint of the condition. In other words, the c-DDPM should possess more diversity to construct the contaminated area.

The preparation of training dataset is an important component of the DL method. We applied the convolution modeling to generate the ground roll. We configure the source wavelet frequency for the convolution modeling to fall within a certain range, ensuring the trained model generalizes effectively across both synthetic and field data examples. The training dataset used above is named as ”dataset1”. To test how the frequency of the ground roll influence the generation quality, we enlarge the frequency range of the ground roll and generate a new training dataset, named as ”dataset2”. The new dataset is then used to train the DDPM-2c with the same training parameters as before. After training, we then apply the trained model to the noisy field data in Fig. 9a. We can see that the residual ground-roll noise with higher frequency are also removed, but the direct arrivals are also missing. Additionally, there are more leakage of reflections into the predicted ground roll when using dataset2 for training compared with the scenario for dataset1. This is because the model trained on dataset2 tends to generate the events with higher frequency (which usually refer to the reflections), especially in the area where the reflections and the ground roll heavily overlap.

We then compute the amplitude spectra of the predicted clean data and ground roll corresponding to dataset2 to further demonstrate the influence of the frequency band of the ground roll for training on the data generation. The amplitude spectra of the predicted clean data using dataset2 exhibits an obvious decrease after 20Hz, while the amplitude spectra of the ground roll increases after 20Hz. The experiment result demonstrates that the frequency band of the training dataset affects the generation quality of the ground roll and the clean data. Thus, we should carefully prepare the training dataset with a reasonable frequency band to improve the predicted clean data and ground roll.

We implemented a conventional sampling process of the DDPM in the c-DDPM. This sampling process is time-consuming because many timesteps are required. In order to accelerate the sampling, DDIM defines both the forward and reverse processes on a non-Markov chain instead of a Markov chain [41]. This non-Markov chain allows DDIM to sample within significantly reduced timesteps, making it much more computationally efficient. However, the reduced timesteps inevitably degrades the quality of the generation. Therefore, we employ the conventional sampling process of the DDPM instead of the sampling process of the DDIM. Here, we provide a brief illustration of the sampling process of the DDIM for reference. The sampling step in DDIM is

where $\tau_{i}$ is the sampled timesteps from the original range 0 … $\tau$, $\alpha$ is the pre-defined noise schedule, and $\sigma_{\tau_{i}}^{2}=\eta\cdot\sqrt{(1-\alpha_{\tau_{i-1}})/(1-\alpha_{\tau_{i}})}\sqrt{(1-\alpha_{\tau_{i}}/\alpha_{\tau_{i-1}})}$. When $\eta=0$, there is no random noise added in every sampling step, thus the sampling process is fixed. When $\eta=1$, the sampling process is the same to the DDPM, only with reduced timesteps.