Fast and Robust Low-Rank Approximation for Five-Dimensional Seismic Data Reconstruction

Five-dimensional (5D) seismic data reconstruction becomes more appealing in recent years because it takes advantage of five physical dimensions of the seismic data and can reconstruct data with large gap. The low-rank approximation approach is one of the most effective methods for reconstructing 5D dataset. However, the main disadvantage of the low-rank approximation method is its low computational efficiency because of many singular value decompositions (SVD) of the block Hankel/Toeplitz matrix in the frequency domain. In this paper, we develop an SVD-free low-rank approximation method for efficient and effective reconstruction and denoising of the seismic data that contain four spatial dimensions. Our SVD-free rank constraint model is based on an alternating minimization strategy, which updates one variable each time while fixing the other two. For each update, we only need to solve a linear least-squares problem with much less expensive QR factorization. The SVD-based and SVD-free low-rank approximation methods in the singular spectrum analysis (SSA) framework are compared in detail, regarding the reconstruction performance and computational cost. The comparison shows that the SVD-free low-rank approximation method can obtain similar reconstruction performance as the SVD-based method but with a large computational speedup.


I. INTRODUCTION
Data reconstruction is extremely important during the entire seismic processing chain due to the fact that our acquired data are never complete and regular in spatial dimensions. More specifically, abundant physical and economic restrictions happen throughout the acquisition phase resulting in data incompletion, e.g., feathering, near-offset gap, cross-line under-sampling in towed marine streamer acquisition and obstacles, terrain restrictions in land acquisition [1], [2], [4]- [8], [10]- [12], [16], [18]- [20]. Incompletely and irregularly recorded seismic data always have severely negative effects on the subsequent data processing workflow, such as surface related multiple elimination, velocity analysis, The associate editor coordinating the review of this manuscript and approving it for publication was Shuangqing Wei . full waveform inversion, wave equation migration, timefrequency analysis, amplitude versus offset inversion and seismic interpretation [21], [24]- [26].
Seismic data reconstruction, therefore, has drawn a large amount of attention from both academia and industry in the past several decades. It can be divided into two very basic categories: model-driven methods and data-driven methods. Although model driven methods are powerful and could be highly accurate, a-priori knowledge about the model needs to be provided, while it is not simple at all to build such an accurate model. In addition to a-priori knowledge about the model, computational cost is another big issue for model-driven reconstruction methods. On the contrary, data-driven reconstruction methods are significantly popular among seismic exploration community because of their higher efficiency and less dependency on model. Usually, VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ simple mathematical assumptions are utilized in data-driven methods for interpolating seismic data, which turn out to be more solid and effective. Since the last decade, a majority of researchers have made lots of efforts to investigate all different types of data-driven reconstruction methods. In recent years, dictionary learning and machine learning are applied to the reconstruction of 5D simple data [28]- [33]. Data driven tight frame (DDTF) is a kind of dictionary-learning method, which can simultaneously denoise and interpolate 5D seismic data [28]. In DDTF, a sparsity-promoting algorithm is used to build the dictionary which can represent the observed data and estimate the complete data. Reference [30] combined the DDTF with a classic machine learning method named support vector regression (SVR) to optimize the learning, which obtain better performance than the Gauss SVR method.
With the continuous improvement of intelligent methods, learning-based 5D data reconstruction will be a hot research topic in the future.
In the seismic data processing literature, it is also noticed that the challenge of random noise always comes along with data reconstruction. In general, random noise suppression is of significant importance for both increasing the signal-tonoise ratio (SNR) [35] and improving seismic data reconstruction performance. Random noise is capable of severely affecting the final performances of reconstruction methods, and thus, much attention has been drawn to random noise attenuation algorithms [35]- [38]. While most random noise attenuation approaches utilize the signal predictability, the irregularity of seismic data can directly influence this predictability and result in unsatisfactory denoising performance. Because of these mutual influences mentioned above between random noise suppression and data reconstruction, simultaneous reconstruction and denoising approaches are widely studied in [42]- [45], [46] and [47]. Among all these approaches, low-rank approximation methods are most commonly studied.
Meanwhile in recent years, simultaneous 5D seismic data reconstruction and denoising has become very popular due to the fact that it is capable of taking all physical dimensions of seismic data into consideration [49]- [51], and thus, taking advantage of more data constraints and correlations to improve the reconstruction performance even at an extremely high data missing ratio. Especially when compared with 2D and 3D data cases, 5D data interpolation and denoising can always provide much better results. There exist two different categories for 5D interpolation: Fourier based approaches and low-rank approximation based approaches. Nowadays, 5D Fourier based interpolation methods have already become a standard for industry due to its stability and efficiency. Nevertheless, it is reported that the reconstruction performance of 5D Fourier based methods is less accurate than the low-rank approximation based approaches especially when dealing with curved events. On the other hand, low-rank approximation methods for simultaneous 5D seismic data reconstruction and denoising are widely studied during the last decade.
Low-rank approximation approaches still have two subcategories: tensor based methods and block Hankel/ Toeplitz matrix based methods. The first category regards multi-dimensional data as multi-linear arrays and applies dimensionality reduction on multi-linear arrays, in which some folding and unfolding operations are utilized. These methods are usually referred to as tensor completion. The other category transforms seismic data into multi-level block Hankel/Toeplitz matrix and a rank reduction algorithm can be used to recover the data, which can also be named multichannel singular spectrum analysis (MSSA) or Cadzow filtering. Although simultaneous 5D seismic data reconstruction and denoising via low-rank approximation becomes more attractive recently, there are still many challenges within this framework. Major issues in low-rank approximation approaches for 5D data interpolation are [53]: (1) high computational cost during SVD process, (2) the inevitable residual noise [19] and (3) the rank inconsistency problem [71].
Firstly, high computational cost during SVD process is the biggest problem we encounter especially for 5D dataset [54], [55]. Since our target multi-level block Hankel/Toeplitz matrix or tensor array tends to be extremely large, the truncated SVD process consumes huge amounts of time. Many researchers have made lots of efforts to solve this problem. Reference [42] proposed to replace truncated SVD with randomized SVD to accelerate the MSSA based rank reduction method. [43] introduced a fast rank reduction approach along four spacial dimensions based on Lanczos bidiagonalization. Instead of accelerating low-rank approximation for multi-level Hankel/Toeplitz matrix, [47] presented a tensor completion based fast algorithm via parallel matrix factorization to speed up the 5D interpolation process. Recently, [52] developed a parallel square matrix factorization method for efficient 5D data reconstruction without SVD calculation, which however could be difficult to implement.
Secondly, the recently discovered inevitable residual noise is another big issue preventing low-rank approximation approaches from having better performance [56]- [59]. Usually, conventional low-rank approximation methods can obtain reasonably good results while there are still amount of visible residual noise left, which is more severe for land seismic data. The reason is simply due to the incomplete decomposition between signals and noise. Reference [61] first noticed the inevitable residual noise during the truncated SVD process in MSSA for random noise attenuation and it is proved that the residual noise comes from the signalplus-noise subspace, which is decomposed by the truncated SVD. Therefore, damped MSSA is proposed to mitigate this noise leakage problem by introducing a damping factor to better decompose the signal and noise [61]. Note that the later proposed double least-squares projection has more or less the similar concept with damped MSSA [62]. Reference [63] proposed a multi-step damped MSSA algorithm to further improve the SNR during reconstruction. Straightforwardly, [19] extended damped MSSA to simultaneous 5D seismic data reconstruction and denoising in the condition with extremely noisy land seismic data. A much cleaner and more robust 5D interpolation result can be obtained by imposing more accurate signal and noise space decomposition. In addition, [64] presented a novel hybrid rank sparsity constraint on seismic data to improve the interpolation performance with less residual noise. This hybrid constraint can take advantage of both constraints to upgrade the final results. A similar hybrid model was proposed in [65] based on a different framework called the three-operator proximal splitting scheme.
Thirdly, an easily neglected but highly important issue is the rank inconsistency problem during low-rank approximation [27], [66], [67]. Seismic events are usually complicated, therefore, it is difficult to select a single appropriate rank value for the whole dataset. If the rank is too large, little noise will be removed and if the rank is too small, even useful signals are lost. An effective strategy to mitigate this issue is to process the data within local windows, which better satisfies the basic linear assumption of low-rank approximation methods. However, it still fails due to highly nonstationary property of seismic data in both time and space, for example, when dealing with crossing events where the rank should be higher than other local windows. Reference [66] proposed empirical low-rank approximation, in which the multi-dip data are decomposed to multiple single-dip data so that the optimal rank (e.g., one) can be applied to each local window. More recently, [71] developed an adaptive rank thresholding method to make the damped rank-reduction method suitable for local processing.
In this paper, we address the first issue, i.e., the problem of low computational efficiency, of the low-rank approximation methods. We discuss an SVD-free low-rank approximation method for efficient and effective reconstruction and denoising of 5D seismic data that contains four spatial dimensions. Our SVD-free rank constraint model is based on an alternating minimization strategy, which updates one variable each time while fixing the other two. For each update, we only need to solve a linear least-squares problem with much more efficient QR factorization. We compare the performance of the traditional SVD-based and the proposed SVD-free lowrank approximation methods in detail. Extensive results show that the proposed SVD-free low-rank approximation method can obtain similar performance as the traditional SVD-based method, but with a much higher computing speed. The computational efficiency becomes more noticeable when the data size becomes larger.

II. LOW-RANK APPROXIMATION FRAMEWORK FOR 5D SEISMIC DATA A. THE INVERSE PROBLEM
The 5D reconstruction problem aims to solve the following inverse problem [42], [43], [68]: where D obs denotes the observed data, D ideal denotes the complete data, and S denotes the sampling operator, which can be more specifically displayed as follows: represents the subset of observed indices. S equals 1 at the observation point while 0 at the missing traces. • denotes element-wise product. Both D obs and D ideal are considered in the frequency domain.
The inverse problem can be solved via the rank-reduction method, based on a weighted iterative algorithm [42], [43]: where D 0 = D obs . The operator M represents the low-rank approximation operator that is based on truncated singular value decomposition (TSVD) for rank reduction. a n is an iteration-dependent scalar that linearly decreases from a 1 = 1 to a n max = 0. The algorithm stops when either a maximum number of iterations n max has been reached or tol denotes a very small tolerance value. In fact, equation 3 represents the general low-rank approximation for simultaneous seismic data reconstruction and denoising. Next, we will introduce the rank constraint model that is connected with the rank-reduction operator.

B. THE RANK CONSTRAINT MODEL
Presuming the original data matrix X is of low-rank structure, due to the theory of matrix completion, the randomly missing data can be reconstructed by a rank-minimizing model: where rank(HX) is the number of non-zero singular values of HX. H denotes some mapping operations, for instance, the block Hankelization/Toeplitization operator in equation 6. Y is the observed seismic data. The constrained minimization problem in equation 4 can be converted to the unconstrained minimization problem in equation 5 by introducing the balancing parameter: where · F denotes the Frobenius norm. S denotes trace sampling. µ denotes balancing coefficient, which can balance the weight between rank constraint and data misfit. In practice, this balancing parameter is explicitly connected to the estimated rank value in the rank constraint. The transformation of the 4D data matrix into a level-four block Hankel/Toeplitz matrix can be represented in operator notation as follows: where H denotes the level-four block Hankelization/ Toeplitization operator. VOLUME 8, 2020 Both missing traces and additive noise increase the rank of the matrix M, which originally possesses low-rank structure. We assume that M has full rank, rank(M) = J and the desired low-rank matrix M K has deficient rank, rank(M K ) = K < J . Therefore, the desired low-rank approximation can be obtained by a rank reduction method using truncated singular value decomposition (TSVD): where we use R K as the rank reduction operator via TSVD. U K resulting from TSVD denotes the first K singular vectors of matrix M. The symbol [·] H represents the conjugate transpose of a matrix. Due to the TSVD step, the low-rank approximation method discussed in this paper can also be understood as a PCA method [62]. After rank reduction, the filtered data can be recovered with random noise attenuated and missing traces reconstructed via properly averaging along the anti-diagonals of the target low-rank approximation matrix M K : where A denotes the averaging operator.
Here is a complete and detailed description of low-rank approximation framework for 5D seismic data. We adopt the following algorithm for reconstructing 5D noisy seismic data with missing traces: The iteration terminates after all F frequencies are finished. For each frequency slice, the iteration terminates after N iterations or upon reaching convergence to the specified tolerance tol.

III. FAST AND ROBUST LOW-RANK APPROXIMATION APPROACH
Both conventional rank constraint and its solution are valid and robust. However, the computational cost is a major issue due to SVD in the low-rank approximation framework. SVD consumes most computational time during the whole process and it is especially severe for large scale problem, such as 5D seismic data reconstruction. In order to avoid SVD-related computations, we introduce a new SVD-free rank constraint model following Wen et al. (2012) [69].

A. SVD-FREE RANK CONSTRAINT MODEL
Still presuming the original data matrix X is of low-rank structure, our target level-four block Hankel/Toeplitz matrix after mapping HX possess the same property. It is common sense that m by n matrix HX with rank up to K can be represented by a matrix product HX = PQ where the size of P is m by K and the size of Q is K by n. A new SVD-free rank constraint model can be expressed as: where R is introduced mainly for a computational reason. P and Q are the two low-rank decomposition matrices using the SVD-free algorithm, such that the product of the estimated P and Q is the low-rank approximated result. The advantage of this new SVD-free rank constraint model is the higher computational efficiency without costly SVD-related process. More specifically speaking, normal QR factorization instead of full or partial SVD is the principle components for computation. Similar to other algorithms, the solution to our SVD-free rank constraint model should be naturally based on an alternating minimization strategy. That is to say, updating one variable each time while fixing the other two. For each update, we only need to solve a linear least-squares problem with much less expensive QR factorization. In addition, a more complex nonlinear successive over-relaxation scheme with dynamically adjusted relaxation weight is utilized [69]. The detailed updates in each iteration are displayed as follows: where (.) † denotes the Moore-Penrose pseudo-inverse of a target matrix, which is implemented by a normal QR algorithm. ω is the relaxation weight parameter and S H represents sampling on the Hankel/Toeplitz matrix instead of the original data matrix. P * , Q * and R * are the updated versions of P, Q and R.
To improve the efficiency of the introduced algorithm, an update on the relaxation weight parameter is used based on the calculation of residual ratio γ [69]: where if γ > 1, we simply reset ω to 1 and if γ < 1, we should keep it unchanged. However, if γ < 1 but is close to 1, we need to increase ω to min(ω + τ, ω max ). Here, τ is a tiny increment and ω max is the maximum value that ω can be.

B. FAST 5D RECONSTRUCTION
Therefore, according to equations 10, 11 and 12, the new SVD-free rank constraint model can be effectively solved.
Note that a proof of the convergence of SVD-free lowrank approximation algorithm can be found in [69]. Due to the absence of SVD-related calculation, the new model is generally able to have a large speedup compared to the conventional model. Besides, it is extremely easy to apply the new model to the low-rank approximation framework. The whole solution to new model can be regarded as a low-rank approximation operator, and then equation 8 can be revised to:D where R L denotes the rank reduction operator based on the SVD-free rank constraint model and L is its corresponding low-rank approximation operator. The final fast and robust low-rank approximation framework for 5D seismic data is introduced as follows: it is clear that the only difference between equations 3 and 17 is the SVD-free rank constraint based low-rank approximation operator. It is the first time that the fast SVD-free low-rank decomposition method [69] is used in the 5D seismic reconstruction problem. It is known that the matrix completion is a common mathematical model. The method developed in Wen et al. (2012) [69] is a general solution to accelerate a matrix completion problem. It was developed to accelerate a variety of real-world applications, like the one in this paper. Here, we leverage this fast algorithm to accelerate the stateof-the-art 5D seismic reconstruction framework.

IV. SYNTHETIC EXAMPLES
In this section, we use several synthetic examples to demonstrate that the proposed SVD-free low-rank approximation method can significantly accelerate the reconstruction of 5D seismic data without loss of accuracy. Although Hankel and Toeplitz matrices [72], [73] seem to be similar as the low-rank matrix, but in some cases one can lead us to better results in comparison to the other one. According to our experience, the Toeplitz matrix results in a slightly better reconstruction performance. So, here, we prefer to use the Toeplitz matrix as the low-rank matrix. The first example is a synthetic example with linear events. We first construct a clean synthetic data, then we add some random noise to the data and also randomly remove some traces from the noisy data to simulate the incomplete seismic data. The size of this synthetic data is 100 × 10 × 10 × 10 × 10. A comparison of the clean data, noisy data, and the incomplete data is presented in Figure 1, from which we can see that the added random noise is so strong that most useful seismic signals are not visible. Since 80% of the total traces are removed from the original noisy data, the observed data is considered to be highly incomplete. We compare the reconstruction performance between the SVD-based low-rank approximation method and the SVD-free version in Figure 2. The top row in Figure 2 corresponds to the results from the SVD-based low-rank approximation method and the bottom row in Figure 2 corresponds to the results from the SVD-free version. In Figure 2, we show the reconstructed data in common offset gathers in the left column, the removed noise from the noisy data in the middle column, and the reconstruction error in the right column. The results show that the SVD-based and SVD-free low-rank approximation methods both obtain very successful reconstructions. The gaps in the observed incomplete data have been filled with seismic signals for both methods. The removed noise cubes of both methods do not contain obvious signal patterns, indicating that there are no significant damages to the useful signals during the reconstruction. The reconstruction error cubes also demonstrate that the error is negligible for both methods. However, for this example, the proposed SVD-free method is much faster than the SVD-based method. As tested on a Mac Pro Laptop equipped with an Intel Core i7 CPU clocked at 2.5 GHz and 16 GB of RAM, the SVD-based method takes 603 s to finish the calculation while the proposed one only takes 153 s. The computational speedup is almost four times. A comparison of the computational efficiency of different methods for different data sizes is shown in Table 1. The denoising performance in the common midpoint gathers is presented in Figures 3 and 4. Figure 3 plots the clean synthetic data, noisy data, and decimated data with 80% traces randomly removed. Figure 4 shows the reconstructed data, removed noise, and reconstruction error cubes using the SVD-based and SVD-free low-rank approximation methods. A single-slice comparison in the common offset gather is presented in Figure 5. From the single-slice comparison, we can see that both methods obtain very similar results (Figures 5(d) and 5(e)), which are very close to the clean data ( Figure 5(a)). To numerically compare the reconstruction VOLUME 8, 2020 TABLE 1. Comparison of computational cost of both the SVD-based and SVD-free low-rank approximation for different data sizes based on an Intel Core i7 CPU. performance, we use the signal-to-noise ratio (SNR) defined as follows [61], [70]: where s denotes the vectorized clean data andŝ denotes the vectorized reconstructed data.
To further compare the amplitude difference among different datasets, we extract a single trace from each data cubes of the clean, noisy, observed, and reconstructed data using two methods. The trace to be compared is the middle trace in each dataset. The single-trace comparison is plotted in Figure 6. The clean data is plotted as the black line. The noisy data is plotted as the red line. The observed data is not plotted because in this position it is a blank trace. The blue line plots the trace using the SVD-based low-rank approximation method. The green line plots the trace using the SVD-free method. It is apparent that both SVD-based and SVD-free methods obtain very successful recovery of the exact solution, i.e., the clean trace. Although there are tiny differences  between the blue and the green lines, they are overall very similar.
The calculated SNR of the incomplete data is −8.08 dB. The SNRs of the SVD-based and SVD-free low-rank approximation methods are 9.53 dB and 9.79 dB, respectively. We find that the proposed SVD-free low-rank approximation can even obtain a slight improvement in terms of SNR.  Besides, we also use the local similarity metric [35] to compare the reconstruction performance of different methods. Figure 7 shows the comparison of local similarity maps of different datasets. Figure 7(a) plots the local similarity between clean data and the noisy data, where we can see that the high local similarity is consistent with the distribution of the useful signals. Although, it is difficult to observe signals from the very noisy data, it is possible to use local similarity to highlight the distribution of the useful signals in the presence of extremely strong random noise. Figure 7(b) shows the local similarity map between the observed incomplete data and the clean data. The right side of the section has distinctly zero local similarity value, due to the missing traces.  7(c) and 7(d) show the local similarity maps corresponding to the SVD-based and SVD-free low-rank approximation methods. The very close local similarity maps of two methods indicate that the reconstruction performance is nearly the same.
Next, we use a synthetic example with hyperbolic events to test the performance of the proposed method in the case of curved events. The clean data, noisy data, and incomplete data with 80% missing traces are shown in Figure 8. The size of this synthetic data is 101×32×32×5×5. The comparison of reconstruction results using the SVD-based and SVD-free methods is presented in Figure 9. The left column in Figure 9 plots the reconstructed data. The middle column in Figure 9 plots the removed noise. The right column in Figure 9 plots the reconstruction error. From Figure 9, it is clear that both methods obtain very similar performance. However, the computational costs for the SVD-based and SVD-free methods are 7226 s, and 1543s, respectively. Figure 10 shows the comparison of local similarity cubes. Figure 10(a) plots the local similarity between the noisy data and the clean data. Figure 10(b) plots the local similarity between the observed incomplete data and the clean data. Figure 10(c) plots the local similarity between the reconstructed data and the clean data using the SVD-based low-rank approximation method. Figure 10(d) plots the local similarity between the VOLUME 8, 2020  reconstructed data and the clean data using the SVD-free lowrank approximation method. From the similarity cubes, we see that the reconstructions from both methods help recover most of the useful signals and obtain high local similarity measures. The performance of the SVD-free method is further confirmed to be close to the SVD-based method in terms of local similarity measure. The SNRs of the reconstructed data using the SVD-based method and the SVD-free method are 1.75 dB and 1.92 dB, respectively. Note that the SNR of the incomplete data of this example is extremely low, i.e., −13.73 dB. This example also demonstrates the effectiveness of the SVD-free low-rank approximation method in the case of ultra-low SNR.

V. REAL DATA EXAMPLE
Finally, we apply the proposed SVD-free low-rank approximation method to a 5D field data. The data have been binned onto a regular grid and a common offset gather of the field data is shown in Figure 11. In Figure 11, the colored stripes are the recorded seismic traces. The white blanks denote the missing traces, which means that we do not observe seismic data in these positions. The size of this field data is 301×10× 10 × 20 × 20. Because of the difficulty in displaying a 5D data set, we only show a common mid-point gather here. The 3D common midpoint gather is re-arranged into a 2D matrix for a better view. The transparent colored windows denote zooming areas for an amplified comparison. In this example, roughly 80 percent traces are missing from the regular grids. Because of the high ratio of missing traces, the observed seismic traces do not show any spatial coherency. It is difficult to see the waveforms from the raw data. The results from the two aforementioned methods, i.e., SVD-based and SVD-free methods, are shown in Figures 11(b) and 11(c). After 5D reconstruction, the white blanks in the raw data have been filled with seismic traces. The waveforms become well aligned along the spatial direction. Compared with the raw data, all methods seem to obtain a dramatic improvement on FIGURE 11. Comparison of reconstruction performance in common offset gathers for the 5D field data example. (a) Observed data with roughly 80% traces missing. (b) Reconstructed data using the SVD-based low-rank approximation method. (c) Reconstructed data using the SVD-free low-rank approximation method. The frame box areas are zoomed and highlighted in FIGURE 12. the data quality. It is salient that both low-rank approximation methods obtain much smoother and cleaner results. When zooming the data in the transparent red rectangles in both Figures 11(b) and 11(c), the comparison among different datasets becomes much clearer. From Figure 12, we observe that both low-rank approximation methods obtain very similar results. In this example, the SVD-free method takes 15378 s while the SVD-based method takes 92340 s.
The speedup using the proposed SVD-free low-rank approximation is nearly six times.

VI. CONCLUSION
The traditional singular value decomposition (SVD) based low-rank approximation method suffers from the bottleneck of computational cost due to many SVD calculations. In order to relieve the computation overburden, we develop an SVD-free low-rank approximation method in this paper. The SVD-free matrix decomposition is based on an alternating minimization strategy. During the iterative update, we only need a much faster QR factorization for solving a linear least-squares problem. The SVD-free low-rank approximation can be easily inserted into the singular spectrum analysis (SSA) framework. We test the effectiveness of the SVD-free method on synthetic examples with linear and curving events, and a real 5D seismic data. The results show that the proposed SVD-free method does not degrade the reconstruction performance compared with the traditional method. However, the proposed SVD-free method can greatly improve the computational efficiency. He has authored or coauthored more than 100 internationally renowned journal articles and more than 80 international conference papers. His research interests include seismic signal analysis, seismic modeling and inversion, simultaneous-source data deblending and imaging, global adjoint tomography, machine learning, and high-performance computing.