Reconstruction of seismic data based on SFISTA and curvelet transform

In seismic data processing, the reconstruction and interpolation of missing traces are essential tasks. To overcome the limitations of irregularly sampled seismic data, this paper proposes a seismic data interpolation method combining the smoothing fast iterative soft threshold algorithm (SFISTA) and the curvelet transform; this method uses the curvelet domain as the sparse domain. For comparison, the contourlet transform is used for different sparse domains, and the fast iterative shrinkage-thresholding algorithm (FISTA) is used for different optimization algorithms. Numerical modeling and real data show that the SFISTA method in the curvelet domain can give better reconstruction effects and higher reconstruction accuracy than those in the contourlet domain with the FISTA method; in addition, the SFISTA method in the curvelet domain can be used to reconstruct the missing traces of three-dimensional seismic data.


Introduction
In complex environments, seismic data reconstruction has great significance as a recovery technique. Under external disturbance, irregularly sampled seismic data will affect further geological data processing such as migration imaging and data interpretation. In order to obtain high-quality seismic data, interpolation reconstruction is needed to approximate the original data. In recent years, under the compressive sensing theory, seismic data reconstruction methods based on sparse constraints have become more and more popular. It mainly consists of the sparse transform, measurement matrix, and reconstruction algorithm. The sparse transforms that are often used include the Fourier transform (Zhang et al., 2013;Ciabarri et al., 2014), curvelet transform (Hennenfent et al., 2010;Liu et al., 2015;Zhang et al., 2017;Zhang et al., 2019;Tian and Qin, 2022), contourlet transform (Eslami and Radha, 2006), and seislet transform (Liu W et al., 2016). Because the curvelet transform undergoes multi-scale and multi-direction analysis and can perform the optimal local decomposition of seismic data (Yang et al., 2017), the curvelet transform is employed in this paper as a sparse transform, and the contourlet transform is also used in this paper for comparison analysis.
A classical sparse recovery problem usually requires minimizing the L 0 norm, which is NPhard. The L 1 norm is the approximate of the L 0 norm, which is a convex function, and can be solved by the convex optimization algorithms or tools; so, the L 0 norm is replaced by the L 1 norm for simplicity and effectiveness. The iterative soft threshold algorithm (ISTA) shows great advantages (Daubechies et al., 2004;Mohsin et al., 2015) for a convex optimization algorithm. Gradually, a faster ISTA algorithm (FISTA) has been developed. FISTA is more suitable the synthesis approach to sparse recovery. And for the analysis approach to sparse recovery, Tan et al. (2014) proposed a monotone version of the fast iterative shrinkage-thresholding algorithm, which is the smoothing fast iterative soft threshold algorithm (SFISTA). In this paper, we introduce the SFISTA method-based curvelet transform to the seismic data interpolation reconstruction problem (Tan et al., 2014;Pokala and Seelamantula, 2020;Shen et al., 2020). The theory is given in section 2, and experimental results are given in section 3.

Theory
Assuming that the observed seismic data is y and the downsampling matrix is U, the irregular missing seismic data can be modeled as follows: where n is a randomly generated noise and x is the original seismic data.
The interpolation problem can be expressed as follows: where · 2 is the L 2 norm, · 1 is the L 1 norm, λ is the regularization parameter, and Ψ is the analysis operator. Moreover, Ψ*Ψ I, where Ψ* denotes the adjoint of the operator Ψ and I is the identity matrix. The sparse transform can be expressed as x Ψα, where Ψ is the transform operator and α is the sparse coefficient. Let us denote the following equation: where f(x) is the smooth part.
where g(Ψx) is the non-smooth part that needs to be replaced by the Moreau envelope with an approximately smooth form g μ (Ψx), where μ is the smooth approximate parameter.
where ∇ is the gradient and T λμ is the soft threshold operator; the threshold of T λμ is the multiplication of λ and μ.
In this paper, the process of the smoothing fast iterative soft threshold algorithm (Shen et al., 2020) for optimization is expressed as follows: where k is the number of iterations, andx 0 x 0 and t 0 1.

Examples
In this section, we conduct numerical experiments with different seismic data to demonstrate the reconstruction performance of the SFISTA method in the curvelet domain. The algorithm performance is evaluated by interpolation results, the average amplitude spectrum, single-channel interpolation effect, reconstruction error, and signal-to-noise ratio. Numerical experiments are used to test the method. At last, the paper continues to test the 3D interpolation effect of the interpolation method of the proposed method. The experiments are conducted on a Millet computer running on the Windows 10 operating system and Inter Core m3-6Y30.

Seismic data interpolation in the contourlet domain
The regularization parameters of the FISTA and SFISTA methods are set to 10 −3 . The step size of the FISTA is 1, and the step size of the SFISTA is set to 0.5 ; the number of iterations is set as 500. The two-dimensional data tests have the same parameters.
Part data on the Marmousi2 model (Martin et al., 2006) are selected as test data; the data are shown in Figure 1A. Figure 1B shows the irregular data with 50% of the traces removed randomly, and this part's records are padded with zeros. The sparse domain is the contourlet domain. Figures    Frontiers in Earth Science frontiersin.org reconstruction residuals of the SFISTA method have smaller amplitude ranges than those of the FISTA method. It is obvious that the reconstruction method of SFISTA based on the contourlet transform is more effective. The performance of the proposed method in seismic data interpolation could be evaluated using different qualitative and quantitative analyzing tools. The reconstruction error is the general quantitative evaluation tool used in seismic data interpolation; the reconstruction error is defined to be |x −x|, where x denotes the original data andx denotes the reconstructed seismic data. If the reconstruction error is smaller, the reconstructed seismic data will be closer to the original data. SNR is the signal-to-noise ratio, which is defined as SNR 101og 10 x 2 2 x−x 2 2 . A higher SNR value means that the data have better reconstruction quality. The reconstruction error and signal-to-noise ratio are illustrated in Table 1. Table 1 shows that the SFISTA method has a smaller reconstruction error value and higher SNR. Compared to the two methods, the SFISTA method based on the contourlet transform has a lower reconstruction error and higher SNR, and therefore, the SFISTA method shows better performance.

Seismic data interpolation in the curvelet domain
The part of the Marmousi2 model is shown in Figure 2A. Figure 2B shows the irregular data with 50% of the traces missing randomly. Figures 2C, D show interpolation effects using FISTA and SFISTA methods based on the curvelet transform. Figures 2E, F show the reconstruction errors that correspond to the FISTA and SFISTA methods, in which the reconstruction residuals of the SFISTA method have a smaller magnitude than those of the FISTA method. It is obvious that the reconstruction method of SFISTA in the curvelet domain is more effective.
Comparing Figures 1, 2, the SFISTA method in the curvelet domain has lower reconstruction errors when using the qualitative analyzing method. Comparing Tables 1, 2, the SFISTA method in the curvelet domain has lower reconstruction errors and higher SNR when using the quantitative analyzing method. The SFISTA method shows better performance when using the curvelet domain as the sparse domain.

Coral reef model tests in the contourlet domain
The selected coral reef synthetic data are illustrated in Figure 3A. Figure 3B shows corrupted data with 50% of traces   Figure 3F shows a smaller reconstructed error. Table 3 shows the reconstruction error and SNR, which demonstrates the validity of the SFISTA method in the contourlet domain.

Coral reef model tests in the curvelet domain
The selected Coral reef synthetic data are tested in the curvelet domain; the original data are shown in Figure 4A. Figure 4B shows corrupted data with 50% of traces missing randomly. Figures 4C, D show interpolation results using FISTA and SFISTA methods in the curvelet domain. Figures 4E, F show the reconstructed errors between the reconstruction and original data; Figure 4F shows a smaller reconstructed error. Table 4 shows the reconstruction error and SNR, which demonstrates the validity of the proposed method.
Comparing the detailed values in Tables 3, 4, the reconstruction performance of the SFISTA method in the curvelet domain is better than that in the contourlet domain.

Seismic data interpolation in the contourlet domain
The field data are illustrated in Figure 5A. Figure 5B shows corrupted data with 50% of traces missing randomly. Figures 5C, D show interpolation effects using FISTA and SFISTA methods in the contourlet domain. By comparing the interpolation results of the ellipse regions, we learn that the reconstruction results of the SFISTA method in the contourlet domain are similar to the original data. Figures  5E, F show the reconstructed errors between the reconstruction and original data; Figure 5F shows a smaller reconstructed error. Table 5 shows the reconstruction error and SNR, which demonstrates the validity of the SFISTA method in the contourlet domain.

Seismic data interpolation in the curvelet domain
The field data are illustrated in Figure 6A. Figure 6B shows incomplete seismic data with 50% of traces missing randomly. Figures 6C, D show interpolation effects using FISTA and SFISTA methods in the curvelet domain. By comparing the interpolation results of the ellipse regions, we learn that the reconstruction data on the SFISTA method in the curvelet domain are similar to the original data. Figures 6E, F show the reconstructed errors between the reconstruction and original data; Figure 6F shows a smaller reconstructed error.   Table 6 shows the reconstruction error and SNR results of the two algorithms in the curvelet domain. Comparing the two methods, SFISTA performs better with a lower error and higher SNR. Figure 7A plots the interpolation single trace of the missing trace with zero values. In Figures 7, the black, pink, and blue lines correspond to the original data and those obtained by FISTA and SFISTA methods, respectively. Figure 7B shows the zoomed traces from the transparent gray window of Figure 7A. A detailed comparison reveals that the blue line is similar to the black line. We can, therefore, affirm that the proposed SFISTA method can better restore the significant features of the useful signal than the FISTA method.
Comparing the detailed values in Tables 5, 6, the SFISTA method in the curvelet domain is better than that in the contourlet domain.

FIGURE 8
Average amplitude spectrum from original data, interpolated with FISTA and SFISTA methods.
Frontiers in Earth Science frontiersin.org The performance of these methods could be investigated more by comparing the average amplitude spectrum (Mafakheri et al., 2022). The average amplitude spectrum presents the original data and those obtained by FISTA and SFISTA methods, as shown in Figure 8 by the black, pink, and blue lines, respectively. The SFISTA method in the curvelet domain gives a closer spectrum to that of the original data, particularly in the range of 5-30 Hz. In this case, our method has better performance.

Three-dimensional seismic data tests
The experimental results of three sets of 2D data show that SFISTA based on the curvelet transform shows good performance. Next, this method will be tested for 3D seismic data reconstruction. The 3D data (size: 64 × 64 × 64) are from the software package of MathGeo 2020 (https://gitee.com/sevenysw/MathGeo2020). The 3D discrete curvelet transform method comes from the reference of Ying et al. (2005).
The original data are shown in Figure 9A. Figure 9B shows corrupted data with 51% of the traces removed randomly; the iteration number of the FISTA and SFISTA methods is 1000. Figures 9C, D show the interpolation results using the FISTA and SFISTA methods in the curvelet domain. Figures 9E, F show reconstructed errors obtained by the FISTA and SFISTA methods. The quantitatively recovered reconstruction error and SNR are shown in Table 7; Table 7 illustrates that the SFISTA method is better than the FISTA method. Figure 10A plots the interpolation single trace from Figures  9C, D. The black, pink, and blue lines represent the original data and those obtained by FISTA and SFISTA methods, respectively. Figure 10B shows the zoomed traces from the transparent gray window of Figure 10A. A detailed comparison reveals that the blue line is similar to the black line. We can, therefore, affirm that the proposed SFISTA method can better restore the significant features of the useful signal than the FISTA approach.
The performance of these methods could be investigated more by comparing the average amplitude spectrum. The average amplitude spectrum of the original data and FISTA and SFISTA methods are shown in Figure 11 by the black, pink, and blue lines, respectively. Interestingly, in this case, the reconstruction data obtained by our method are very similar to the original data.

Conclusion
In this paper, we proposed a new seismic data reconstruction method combining the curvelet transform and the SFISTA method.  Comparing the results obtained when using the contourlet and curvelet domains as the sparsedomain, it can be concluded that the optimization algorithm in the curvelet domain has better performance. In the same sparse domain, comparing the FISTA and SFISTA methods, it can be concluded that SFISTA shows better performance. The seismic data reconstruction effects of SFISTA based on the curvelet transform have been demonstrated by quantitative and qualitative comparisons with several sets of data. The proposed method can be used for 2D and 3D seismic data reconstruction.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.