Seismic data reconstruction based on low dimensional manifold model

Nan-Ying Lan，F(xiàn)an-Chang Zhang，Xing-Yao Yin

School of Geosciences，China University of Petroleum (East China)，Qingdao，Shandong 266580，China

Keywords:Seismic data reconstruction Low dimensional manifold model Regularization Low-rank approximation

ABSTRACT Seismic data reconstruction is an essential and yet fundamental step in seismic data processing workflow，which is of profound significance to improve migration imaging quality，multiple suppression effect，and seismic inversion accuracy.Regularization methods play a central role in solving the underdetermined inverse problem of seismic data reconstruction.In this paper，a novel regularization approach is proposed，the low dimensional manifold model (LDMM)，for reconstructing the missing seismic data.Our work relies on the fact that seismic patches always occupy a low dimensional manifold.Specifically，we exploit the dimension of the seismic patches manifold as a regularization term in the reconstruction problem，and reconstruct the missing seismic data by enforcing low dimensionality on this manifold.The crucial procedure of the proposed method is to solve the dimension of the patches manifold.Toward this，we adopt an efficient dimensionality calculation method based on low-rank approximation，which provides a reliable safeguard to enforce the constraints in the reconstruction process.Numerical experiments performed on synthetic and field seismic data demonstrate that，compared with the curvelet-based sparsity-promoting L1-norm minimization method and the multichannel singular spectrum analysis method，the proposed method obtains state-of-the-art reconstruction results.

1.Introduction

Seismic data is an important instrument for studying the subsurface geological structure and hydrocarbon distribution，and it contains considerable geophysical information.However，due to the constraints of obstacles or other forbidden areas such as reservoirs，levees，villages，mines，etc.，the collected seismic data are often under-sampled or aliased，which leads to the loss of geophysical information and thus is not conducive to the subsequent research.To restore complete geophysical information，reconstructing the missing seismic data becomes an essential and critical task (Stolt，2002;Trad，2009).

Methods for seismic data reconstruction can be divided into two categories:wave equation methods and signal processing methods.The wave equation methods use the physics of wave propagation to recover missing observations (Ronen，1987;Trad，2003;Malcolm et al.，2005;Zhao et al.，2021).These methods need to understand the velocity distribution in the subsurface and perform full wavefield calculations，which limit their practical applicability.The signal processing methods can be further subdivided into:prediction error filtering methods based on linear event hypothesis，transform domain methods based on sparsity promotion，and lowrank methods.Predictive error filtering methods typically utilize low frequency data components in a regular space grid to estimate the predictive filter required to interpolate high frequency data components.Spitz(1991)first developed a frequency(f-x)domain interpolation method based on a unit-step prediction filter.However，using the unit-step prediction filter to reconstruct the missing data is time-consuming and complex because it requires solving two systems of linear equations.To solve this problem，Porsani(1999) provided a method based on a half-step prediction filter，which requires solving only one system and thus greatly improves the efficiency and tractability of interpolation.Subsequently，Gülünay (2003) extended the idea of prediction filtering to the frequency wavenumber (f-k) domain and gave an efficient f-k interpolation method.Another improved f-x interpolation method has been proposed by Naghizadeh and Sacchi(2009)，which adopts an exponentially weighted recursive least-squares method to estimate the prediction filters，further improving the reconstruction efficiency.

The transform domain methods achieve the reconstruction of under-sampled data by applying the sparsity-promoted minimization to the coefficients in a transform domain.Nowadays，many mathematical transformations have been introduced in seismic data reconstruction，such as Fourier transform(Zwartjes and Gisolf，2006;Tang and Yang，2010;Naghizadeh and Innanen，2011;Zhang et al.，2013)，Radon transform (Tang and Mao，2014;Wang et al.，2017;Ibrahim et al.，2018;Tang et al.，2020)，Curvelet transform(Hennenfent et al.，2010;Zhang and Chen，2013;Bai et al.，2014;Zhang et al.，2017)，Shearlet transform (Liu et al.，2018;Yang et al.，2020)，Seislet transform (Liu et al.2013，2017;Gan et al.，2015)，Dreamlet transform (Wang et al.，2015)，and dictionary-learningbased transforms (Yu et al.，2015;Siahsar et al.，2017a;Ma and Yu，2017;Lan et al.2020，2021).Sparsity-promoted minimization belongs to a classical sparse optimization problem that requires an appropriate regularization tool to obtain a unique and stable solution.Regularization tools used to promote sparsity include L0 norm (Chen et al.，2013);L1 norm (Wang et al.，2011;Yin et al.，2015);Cauchy norm (Sacchi and Ulrych，1996);or mixed L1_L2 norm(Li et al.，2012).Besides，the non-convex Lp(0 ＜p ＜1)norm is also applied to seismic data reconstruction.Zhong et al.(2015)applied the L1/2 norm as the constraint specification to reconstruct 3D irregularly sampled data.Recently，Zhang et al.(2019)proposed a smoothed L1/2 norm and applied it to restore the under-sampled seismic data，which yielded satisfactory reconstruction results.

Low-rank regularization is another most popular reconstruction approach.Trickett et al.(2010) gave a multidimensional interpolation algorithm built on the matrix-rank reduction of constantfrequency slices.Oropeza and Sacchi (2011) rearranged the seismic data into block Hankel matrix in the frequency domain，and then reduced the rank of the block Hankel matrix by multichannel singular spectrum analysis (MSSA) to gain the reconstructed seismic data.Further，Kreimer and Sacchi (2012) provided a fivedimensional seismic data reconstruction method by combining the Hankel operation and high-order singular value decomposition.Jia et al.(2016) presented the orthogonal matrix pursuit Hankel reconstruction (OMPHR)，which used the rank-one orthogonal matching pursuit method instead of the singular value decomposition to reduce the rank and greatly improved the computational efficiency.Alternately，Chen et al.(2016)applied the damped multichannel singular spectrum analysis (DMSSA) operator to reconstruct highly incomplete 5D seismic data.Siahsar et al.(2017b)proposed an optimal singular value contraction method for seismic data recovery，which generates the optimal low-rank estimator according to the random matrix theory，and utilizes the attenuation factor to constrain the singular value to improve the robustness of the rank estimation.More recently，a fast and memory-efficient implementation of the MSSA method is developed by Cheng et al.(2019)for seismic data reconstruction.It takes advantage of random projections and the structure of Hankel matrices to avoid the construction of large Hankel trajectory matrices and thus significantly reduces the computational costs.

Recently，a novel regularization method based on patches manifold had been put forward for general image processing(Osher et al.，2017)，where it has achieved state-of-the-art results.This method held that patches manifold of general images is close to low dimensional，and this low dimensional manifold is an effective geometric prior in the image processing field.The superior performance of this low dimensional manifold model (LDMM) was also verified in other different fields (Lai and Li，2018;Zhu et al.，2018;Abdullah et al.，2019).Latterly，this model has been used in geophysics.Yu et al.(2017) firstly proved that the manifold of seismic data is close to low dimensional，and then applied the LDMM to attenuate seismic noise.In this paper，motivated by the inherent low-dimensional properties of seismic data，we exploit the LDMM regularization to reconstruct the missing seismic data from under-sampled observed data.The article is organized as follows.In the first section，we review the mathematical model of seismic data reconstruction.The subsequent section presents the details of the proposed LDMM-based seismic reconstruction model.A fast and efficient algorithm for solving the proposed reconstruction model is given in the third section.The fourth section performs numerical experiments and compares various methods.The discussion and conclusion are given in the fifth and sixth sections，respectively.

2.Methodology

2.1.Review of seismic data reconstruction model

The mathematical model of multi-channel seismic data reconstruction can be described as:

where y denotes the recorded wave-field data，which often is under-sampled;f denotes the complete wave-field data，Ψ denotes the sampling operator，which is a matrix of diagonal structure，consisted of zero and identity operator:

where each I corresponds to a trace sampled from the complete data，and each O corresponds to a missing trace.Since the sampling operator is highly singular，Equation (1) is highly underdetermined.To solve the under-determined problem shown in Equation (1)，the popular strategy is to convert it into an optimization problem with a regularization term

where R（·）denotes the regularization operator that contains prior information on the modeldenotes the square of L2 norm，λ is the regularization parameter weighing the fidelity termand the model constraint term R（f）.

The most extensively used model constraint methods in the reconstruction field include:sparsity constraints and low-rank constraints.The regularization optimization problem based on sparsity constraints can be gave as(Hennenfent et al.，2010;Wang et al.，2011;Bai et al.，2014):

where D is a sparse transform，such as Fourier transform，Curvelet transform，etc.X is the coefficients of the signal in the transform domain.‖·‖1denotes the L1 norm.Another regularization optimization problem based on low-rank constraints can be described as follows:

where T（·）represents the mapping operations，such as Hankel/Toeplitz operation(Oropeza and Sacchi，2011;Jia et al.，2016).‖·‖*denotes the nuclear norm.

Essentially，both sparsity constraints and low-rank constraints are based on certain assumptions(sparsity or low-rank)of seismic data in the transformation space (sparse space or rank space) to reconstruct the missing data.These model constraint methods have demonstrated their unique properties in seismic data reconstruction.However，when seismic data do not meet these assumptions，they may be difficult to obtain satisfactory reconstruction results.Unlike the aforementioned model-constrained approaches，this paper exploits the inherent low-dimensional properties of seismic data to achieve high-quality reconstruction.More specifically，relying on the fact that seismic patches occupy a low-dimensional manifold，we propose a new model constraint approach，the lowdimensional manifold model，for reconstructing the missing seismic data.The specific details are provided in the next section.

2.2.Seismic data reconstruction based on LDMM

For given seismic data f，definition the matrix G（f）is a set of all patches of size r×r extracted from the seismic image

where P［·］denotes the patch operator，x is center of the patch，and Θ is an index set that describes the location of the patches.The patches set G（f）can be regarded as a point cloud in，containing a large number of patch points.The basic idea of LDMM is that this point cloud is usually close to the smooth manifold embedded in，which is called the patches manifold of seismic data and denotes as M .Fig.1 shows the relationship between patches，patches set，and patches manifold.Due to the low-dimensional properties of seismic patches manifold，the dimension of patches manifold can be exploited a regularization term to reconstruct missing seismic data.More specifically，the reconstruction model based on LDMM regularization can be described as:

where dim（M ）denotes the dimension of patches manifold.

The optimization problem of Equation (7) has strong nonconvex and nonlinear，and its effective solution strategy is iterative method.The crucial step of iterative method is to calculate the dimension of manifold on the patch point cloud utilizing the lowrank approximation method，and the details are introduced in the next section.

2.3.Optimization by alternating iteration method

In most cases，patches manifold is not a single manifold with a specific dimension，but a set of several manifolds with different sub-dimensions.Thus，the patches manifold dimension can be rewritten into the form of integral:

where dim（M（Gx（f）））represents the dimension of the submanifold which the patch Gx（f）belongs.The key point to solving the optimization problem in Equation (8) is how to calculate the dimension of each sub-manifold.In the original LDMM，the dimension of the sub-manifold can be calculated as the Dirichlet energy of the coordinate function on the manifold，and the corresponding Laplace-Beltrami equation can be solved by the point integration method (Osher et al.，2017;Shi and Sun，2017).However，the original LDMM has the drawback of poor computational efficiency (Shi et al.，2018).

To improve work efficiency，we calculate the dimension of each sub-manifold by using low-rank approximation on its manifold.Since the dimension of the sub-manifold which each patch belongs is same as the dimension of its tangent space，it can be approximated as the rank of the covariance matrix generated by the similar set of each patch Gx（f）on M .On the discrete case，the patch set G（f）belongs to the sampling of the manifold M，so the low-rank property of the covariance matrix formed by the similar set of Gx（f）on M can be linearly approximated by the matrix that is generated for the K-nearest-neighbors of corresponding patch in G（f）.Therefore，if we define the matrix FM（Gx（f））as the Knearest-neighbors of patch Gx（f）in the patch set G（f），then there are

Fig.1.The relationship diagram of patches，patches set，and patches manifold.

where rank（·）denotes the rank operator.Note that the matrix FM（Gx（f））is obtained by concatenating the K-nearest-neighbors patches as columns.Given the above definition，Equation(8)can be rewritten as:

The alternating iterative method is a common technique for solving Equation (10)，which needs to divide this optimization problem into two sub-problems:first，fixed the manifold Mk，update the seismic datafk+1;second，the seismic data fk+1is fixed，update the manifold Mk+1.Specifically，these two sub-problems can be stated as:

where k is the number of outer iterations.

For the first sub-problem，by introducing the auxiliary variableit can be expressed as follows:

Using the augmented Lagrangian method，the constrained optimization problem in Equation (12) can be reformulated the following unconstrained form:

where μ represents penalty factor，Q is the scaled Lagrange multipliers，and ‖·‖F(xiàn)denotes the Frobenius norm.Similar to the split-Bregman iteration (Goldstein and Osher，2009;Osher et al.，2017;Lai and Li，2018;Liu et al.，2018)，the solution of optimization Equation (13) can be characterized by the following iterative format:

where j is the number of inner iterations.Equation(14)is a classical nuclear norm optimization problem，which can be solved by the singular value thresholding(Cai et al.，2010).Specifically，its closedform solution can be described as:

where D1/μ（·）presents the soft-thresholding operator，which defines as the following equation:

The optimization problem in Equation (15) is a least squares problem whose solution can be formulated as:

For the second sub-problem，the manifold Mk+1can be directly updated through Equation (11b).Thus，by combining equations(11b)，(16) and (17) and (21)，the LDMM-based reconstruction model shown in Equation(7)can be solved efficiently.Algorithm 1 summarizes the overall algorithm workflow of the proposed LDMM-based seismic data reconstruction method.

Algorithm 1.Seismic data reconstruction based on LDMM

3.Examples

In this section，we perform numerical experiments to test the proposed LDMM method on synthetic and field seismic data.Note that the number of inner and outer iterations of the proposed method is set to 3 and 10 in all the examples，respectively.For comparison，we also present the reconstruction results with the Curvelet-based sparsity-promoting L1-norm minimization method(Wang et al.，2011)and the multichannel singular spectrum analysis method (Oropeza and Sacchi，2011).To quantitatively evaluate the reconstruction performance，this paper introduces two measurement parameters，which defined as follows:

where f*denotes the reconstructed data，f denotes the complete data.Here，the higher the SNR value and the lower the ERR value means that the result has better reconstruction quality.

3.1.Synthetic example

In order to analyze the feasibility and stability of the proposed method，we forward a synthetic record (see Fig.2a) using Marmousi model and finite-difference approach with the following parameter settings:the temporal sampling interval is 2 ms，the time samples are 1001，the trace interval is 25 m，and the trace number is 500.By irregularly removing 50% traces from the complete data，we simulate the under-sampled data and show it in Fig.2b.Fig.3a illustrates the reconstruction results using the LDMM algorithm with patch size 16 × 16 and patch overlap 15，which is substantially close to the exact result shown in Fig.2a.We compared the proposed LDMM method with the Curvelet and MSSA methods.For Curvelet method，the maximum number of iterations is set to 300，and the Lagrangian parameter is set to 0.15.In MSSA method，the processing frequency range is 0-125 Hz，the size of the local window is set to 100，the window overlap is set to 70 in each dimension，the rank is determined by the adaptive strategy (Wang et al.，2020)，and the rank reduction process is implemented by the rank-one orthogonal matching pursuit algorithm (Jia et al.，2016).We reveal the reconstruction results of the Curvelet method and the MSSA method，as shown in Fig.3c and e.From Fig.3c and e，it can be observed that Curvelet method and MSSA method all effectively restore the removing traces，but these are slightly inferior to LDMM in the continuity of seismic events.In addition，we also draw the difference between the reconstructed results of three methods and the original data，as shown in Fig.3b，d，and 3f，respectively.As can be seen from Fig.3b，d，and 3f，the LDMM method has a smaller residual signal than the other two methods.This clearly confirms that the proposed method obtains a more accurate reconstruction result than other two methods.

Fig.2.The synthetic data for testing.(a) Original synthetic record.(b) Decimated data with irregularly missing 50% traces.

Fig.3.Reconstructed results of three methods on the decimated data with irregularly missing 50% of traces.(a) Reconstruction by LDMM，SNR=20.2231 dB.(b) Reconstruction error of LDMM.(c) Reconstruction by Curvelet，SNR=18.5131 dB.(d)Reconstruction error of Curvelet.(e) Reconstruction by MSSA，SNR=16.2443 dB.(f)Reconstruction error of MSSA.

To show the details of the reconstruction in different methods，we present a comparison of the 217th trace in Fig.4.In Fig.4，the blue line is the original seismic trace.The red，magenta and green line is the recovery result of the LDMM，Curvelet，and MSSA，respectively.As can be seen from Fig.4，the blue and red lines are still very close to each other while the magenta line and green line deviate from the blue line a lot as shown in arrows，meaning that the proposed LDMM method can better preserve the amplitude.We summarize the evaluation metrics and calculation time of the above three reconstruction methods in Table 1.From the evaluation metrics in Table 1，we can note that the SNR value of reconstructed data of the LDMM is 20.2231 dB，which is higher than other methods.At the same time，the LDMM method yields a lower ERR value than the other methods.It is apparent that the LDMM method achieves a higher reconstruction quality than the other methods.From the time cost comparisons in Table 1，although the efficient iterative algorithm is utilized to solve the LDMM-based reconstruction model，its efficiency is still slightly inferior to the other two methods.

Table 1The evaluation metrics and time cost of three methods in synthetic example.

To test the sensitivities of different methods to the undersampling rate，we consider several sampling rates vary from 30%to 80%.For a fair comparison，we repeat the test for each sampling ratio with 10 irregularly sampling operators，and then take the average value of the 10 indicators (SNR and ERR) as the metric at this sampling ratio.Fig.5a and b respectively provide a comparison of the SNR values and ERR values for all of the above sampling ratios.From Fig.5，the SNR values of the LDMM method are always the largest and the ERR values are always the smallest among these methods.Meanwhile，it also can be found that the reconstruction accuracy of the three approaches is not much different at high sampling rates，while the proposed LDMM method has a significant improvement at low sampling rates.This demonstrates that the proposed LDMM method can be potentially advantageous in reconstructing seismic data with low sampling rates.To illustrate this point more visually，F(xiàn)ig.6 shows an example of the above synthetic data reconstruction where 30% of traces are sampling.Fig.6a and b plot the complete and under-sampled data，respectively.Fig.6c-e presents the reconstructed data of the LDMM，Curvelet，and MSSA methods，and Fig.6f-g draw their reconstruction errors，respectively.As shown in Fig.6，the Curvelet and MSSA methods recover the missing traces to a certain extent，while introducing a large number of amplitude artifacts in reconstructed profiles.In contrast，the LDMM method can produce much better results because it recovers the effective signal well and greatly avoids the generation of additional artifacts.Using Curvelet method and MSSA method in this case，we could achieve SNR=11.2358 dB(ERR=0.2743)，and SNR=9.1107 dB(ERR=0.3503)，respectively，while LDMM achieved SNR=13.4563 dB(ERR=0.2124).Obviously，the LDMM method does superb work of reconstructing the seismic data with a low sampling rate.

3.2.Field data example

To verify the applicability of the proposed method，we choose a field pre-stack gather as the experimental data，as shown in Fig.7a.Fig.7b presents the under-sampled data，which regularly removes 50%traces of Fig.7a.The reconstructed data obtained by the LDMM，Curvelet，and MSSA methods are depicted in Fig.8a，c，8e，and their corresponding reconstruction errors are shown in Fig.8b，d，8f，respectively.In this example，the maximum number of iterations of Curvelet method is 300，and the Lagrange parameter is 0.062.For the MSSA method，we linearly increase the frequency band from 0 to 125 Hz at an interval of 0.05 Hz and implement the reconstruction on each frequency slice.Meanwhile，the local window strategy is still applied，and its specific parameters are set as follows:the window size is 60，the window overlap is 35，and the rank is also determined by the adaptive strategy.For the proposed LDMM method，the parameters are set up in the same way as the synthetic example.From Fig.8c and d，we can observe that the Curvelet method successfully recovers missing data，but it also produces certain amplitude artifacts.It is found from Fig.8e and f that the MSSA method performs slightly worse on the wave-field boundaries，and there are serious energy residues in these areas.From Fig.8a and b，it can be observed that，the continuity of seismic events is well enhanced and the amplitude variation is effectively maintained after reconstruction by the present method，while the reconstruction result is closer to the original field data，demonstrating that the proposed LDMM method still has higher accuracy in field data reconstruction.

Fig.4.Comparison of the 217th seismic trace in Fig.3.In which the blue line is the original seismic trace.The red，magenta and green line is the recovery result of the LDMM，Curvelet，and MSSA，respectively.

Fig.5.The SNR (a) and ERR (b) diagrams of the three approaches with different sampling rate.

Fig.6.Reconstructed results of three methods on the under-sampled data with missing 70% of traces.(a) Complete data.(b) Under-sampled data.(c) Reconstruction by LDMM，SNR=13.4563 dB.(d) Reconstruction by Curvelet，SNR=11.2358 dB.(e)Reconstruction by MSSA，SNR=9.1107 dB.(f) Reconstruction error of LDMM.(g) Reconstruction error of Curvelet.(h) Reconstruction error of MSSA.

Fig.7.The first field data for testing.(a) The complete field data.(b) The under-sampled data with regularly removing 50% traces.

In addition，we also compared the F-K spectrum of the second field data before and after reconstruction.Fig.9a and b are the F-K representations of the data in Fig.7a and b，respectively.From Fig.9b，we can see that severe spatial aliasing that is attributable to the regular removal of seismic traces.Fig.10 illustrates the F-K spectrum of the data in Fig.8，in which Fig.10a，c，10e are the F-K spectrum of the reconstructed results of three methods，and Fig.10b，d，10f are the F-K spectrum of the reconstructed errors.From F-K spectrum of the reconstruction results，we observe that the most of aliased energy presented in the F-K spectrum has been effectively eliminated by three methods.However，from the F-K spectrum of the reconstruction errors，we can find that the more residual energy is distributed around the effective signal in both the Curvelet method and the MSSA method (as shown the white arrows)，and the proposed LDMM method has fewer residual energy around the effective signal.This result further proves that the proposed LDMM method has distinct strengths in aspect of preserving the reconstructed amplitude.Table 2 presents the reconstruction evaluation metrics and time cost of the above three methods in this example.From Table 2，it is clear that the proposed LDMM method outperforms the other two methods in terms of accuracy，but is slightly inferior to them in terms of efficiency.

Table 2The evaluation metrics and time cost of three methods in first field data example.

The second field example is shown in Fig.11，which intends to further test the effectiveness of the proposed LDMM method for reconstructing the mixed missing data.Fig.11a plots the original complete seismic gather with a temporal sampling interval of 2 ms and a trace interval of 10 m.We regularly remove 50%of the traces and irregularly remove 10%of the traces from Fig.11a to obtain the mixed missing data，as shown in Fig.11b.We implement the proposed LDMM method to reconstruct this decimated data.In this case，the patch size is 20 × 20，and the patch overlap is 19.The reconstruction result and error are shown in Fig.11c and d，for which we can see that the proposed LDMM method successfully restore the removing seismic traces，and the signal leakage is quite small.To more vividly show the reconstruction effect，we have drawn the F-K spectrum of Fig.11.Fig.12a displays the F-K spectrum of the complete data.Fig.12b displays the F-K spectrum of the decimated data and clearly shows that severe noise and spectrum aliasing are produced in the F-K domain due to the mixed missing.Fig.12c displays the F-K spectrum of the reconstruction data performed by the LDMM method.As Fig.12c exhibits，after reconstruction by the proposed LDMM method，the noise and aliasing in F-K spectrum are well suppressed，and the reconstructed spectrum is essentially consistent with the spectrum of the original data.In this example，the value of SNR and ERR is equal to 15.1765 dB and 0.1742，respectively.All the evidence demonstrates that the proposed LDMM method has good adaptability for reconstructing the mixed missing data.

The last example is employed to validate the practicality of this method in actual processing.Fig.13a shows an acquired pre-stack gather，which contains 150 traces with a spatial sampling interval of 10 m.Fig.13c shows an enlargement of the rectangular area in Fig.13a.As seen in Fig.13a and c，there is an obvious sawtooth phenomenon in the acquired gather，which is caused by too large spatial sampling interval.Using the presented LDMM method for reconstruction，we obtained reconstructed gather with 300 traces and 5 m spatial sampling interval，as shown in Fig.13b.It can be seen from Fig.13b that，the sawtooth phenomenon is well alleviated after reconstruction by the LDMM method，which can also be clearly seen in Fig.13d.Fig.14 plots the F-K spectrum of the gather before and after reconstruction using the proposed LDMM method.Fig.14a pictures the F-K spectrum of the seismic gather before reconstruction，in which the spatial aliasing caused by big trace interval is visible，as shown at the arrow.Fig.14b is the F-K spectrum of the seismic gather reconstructed by the LDMM method presented in this paper.After reconstruction，trace interval is halved so that the spatial aliasing phenomenon is effectively suppressed.This result further demonstrates the effectiveness and practicality of the proposed method for seismic data reconstruction.

Fig.8.Reconstructed results of three methods on the first field data.(a) Reconstruction by LDMM，SNR=16.7316 dB.(b) Reconstruction error of LDMM.(c) Reconstruction by Curvelet，SNR=14.5530 dB.(d) Reconstruction error of Curvelet.(e)Reconstruction by MSSA，SNR=13.7563 dB.(f) Reconstruction error of MSSA.

4.Discussion

The low dimensional manifold model is a patch-based regularization method，thus the characteristics (patch size and overlap degree) of the patch will directly determine the reconstruction accuracy and efficiency of this method.Herein，we will examine the impact of different patch features on the reconstruction results.We first analyze the sensitivity of the LDMM method to the patch size parameter using Fig.7b as the test data.Fig.15 shows the reconstruction accuracy and efficiency of LDMM with patch sizes from 6 × 6 to 24 × 24 (the patch overlap is set to maximum).It can be observed from Fig.15a that，as the patch size increases，the SNR of the reconstruction result increases first and then decreases slightly，and the ERR decreases rapidly and then increases slightly.This indicates that when the patch size is too small，the complete structural information contained in the patch will be very scarce，which is unfavorable for the recovery of missing data.Conversely，when the patch size is too large，over-fitting will make the result too smooth and lose more detail features，resulting in a slight decrease in reconstruction quality.In Fig.15b，we can see that the consuming time exhibits an exponential increase with the patch size increases.Comprehensive considering of the relationship between reconstruction quality and computational time，we hold that the patch size of 16 × 16 is the most appropriate in this example.

Here we discuss how overlap degree affects the accuracy and computational efficiency of LDMM.Similarly，taking Fig.7b as an example，the relationship between patch overlap and the reconstruction result is studied where the patch size is equal to 16×16.Fig.16a presents the variation of the number of patches with the patch overlap degree and clearly shows that the higher patch overlap can produce a greater number of patches.Fig.16b portrays the change curve between the patch overlap and the evaluation metrics of reconstruction result，and Fig.16c demonstrates how the computational time changes with patch overlap.As Figs.16b and c shown，the reconstruction accuracy and computational time both increase as the patch overlap increases.The reason is perhaps that the highly overlapped patches can generate more points (patches)to approximate the manifold，thus improving the reconstruction quality.However，using more points (patches) to approximate the manifold also brings a huge amount of computation.In seismic data reconstruction work，its goal is to obtain optimal reconstruction result in order to recover the lost geophysical information as much as possible.Therefore，we recommend that the overlap degree of patch should be set to the maximum when using the LDMM method for seismic data reconstruction.

Fig.9.The F-K spectrum of Fig.7.(a) The F-K spectrum of Fig.7a.(b) The F-K spectrum of Fig.7b.

Fig.10.The F-K spectrum of Fig.8.(a)The F-K spectrum of Fig.8a.(b)The F-K spectrum of Fig.8b.(c)The F-K spectrum of Fig.8c.(d)The F-K spectrum of Fig.8d.(e)The F-K spectrum of Fig.8e.(f) The F-K spectrum of Fig.8f.

Fig.11.The second field example for testing.(a)The original field gather.(b)The mixed missing data with removing 60%of traces.(c)The reconstructed result using the proposed LDMM method.(d) The reconstructed residual.

Fig.12.The F-K spectrum of Fig.11.(a) The F-K spectrum of Fig.11a.(b) The F-K spectrum of Fig.11b.(c) The F-K spectrum of Fig.11c.

Fig.13.The last field example for testing.(a)Original field gather with 10 m trace interval.(b)Reconstructed field gather using proposed LDMM method with 5 m trace interval.(c)Enlargement of the rectangular area in Fig.13a.(d) Enlargement of the rectangular area in Fig.13b.

Fig.14.The comparisons of the F-K spectrum.(a) The F-K spectrum of Fig.13a.(b) The F-K spectrum of Fig.13b.

Fig.15.The change curves of evaluation metrics (a) and the computational time (b) with patch size r × r.

5.Conclusions

In this paper，we have presented a novel framework for reconstructing the under-sampled seismic data based on low dimensional manifold model.In a manner，the proposed method can be regarded as a generalization of conventional regularized reconstruction method.Specifically，the proposed method replaces the transformation or mapping operations of the conventional regularized reconstruction method with a patches manifold of the seismic data，and utilizes the dimension of the patches manifold as the regularizer to reconstruct the under-sampled seismic data.The crucial procedure of the proposed method is to solve the dimension of the patches manifold.Toward this，we adopt an efficient dimensionality calculation method based on low-rank approximation，which provides a reliable safeguard to enforce low dimensionality in the reconstruction process.We use one synthetic data set and three field gathers to test the performance of the proposed method.Numerical experiments show that the proposed method has good flexibility and adaptability.Meanwhile，the proposed method also achieves state-of-the-art reconstruction results，compared to the Curvelet-based sparsity-promoting L1-norm minimization method and the multichannel singular spectrum analysis method.In future work，we will focus on the expansion of this method for reconstructing high-dimensional seismic data，as well as investigate the parallelization to further speed up its execution efficiency.

Fig.16.The change curves of the number of patches (a)，the evaluation metrics (b) and the computational time (c) with different patch overlap.

Acknowledgments

We are grateful to all the reviewers for their helpful comments on this paper.This work is supported by National Natural Science Foundation of China (Grant No.41874146 and No.42030103) and Postgraduate Innovation Project of China University of Petroleum(East China) (No.YCX2021012).