Seismic imaging of incomplete data and simultaneous-source data using least-squares reverse time migration with shaping regularization

Simultaneous-source acquisition improves the efficiency of the seismic data acquisition process. However, direct imaging of simultaneous-source data may introduce crosstalk artifacts in the final image. Likewise, direct imaging of incomplete data avoids the step of data reconstruction, but it can suffer from migration artifacts. We have proposed to incorporate shaping regularization into least-squares reverse time migration (LSRTM) and use it for suppressing interference noise caused by simultaneous-source data or migration artifacts caused by incomplete data. To implement LSRTM, we have applied lowrank one-step reverse time migration and its adjoint iteratively in the conjugate-gradient algorithm to minimize the data misfit. A shaping operator imposing structure constraints on the estimated model was applied at each iteration. We constructed the shaping operator as a structure-enhancing filtering to attenuate migration artifacts and crosstalk noise while preserving structural information. We have carried out numerical tests on synthetic models in which the proposed method exhibited a fast convergence rate and was effective in attenuating migration artifacts and crosstalk noise.


INTRODUCTION
Seismic migration can be regarded as a linear inverse problem, and its solution can be obtained by iteratively seeking a model that best fits the observed data.This approach is often called leastsquares migration (LSM) (Ronen and Liner, 2000).LSM's original implementation was with Kirchhoff migration (Nemeth et al., 1999;Duquet et al., 2000); then, it was developed for one-way waveequation migration (Kuehl and Sacchi, 2003;Clapp et al., 2005;Wang et al., 2005), and more recently, it has been implemented with reverse time migration (RTM) (Dai et al., 2012;Liu et al., 2013;Zhang et al., 2013;Hou and Symes, 2015;Sun et al., 2015b;Zhao et al., 2015).It has been demonstrated that LSM helps to attenuate migration artifacts and provides a sharper, better frequency distribution of estimated reflectivity (Nemeth et al., 1999;Kuehl and Sacchi, 2003;Tang, 2006;Dutta et al., 2014).LSM can also be applied to simultaneous-source data (Beasley, 2008;Hampson et al., 2008) or blended data (Berkhout, 2008) in which it not only enhances the resolution of the image, but also suppresses the crosstalk without any preseparation of the blended data (Dai and Schuster, 2009;Tang and Biondi, 2009).Although LSM using gradient-type iterative methods is able to attenuate migration artifacts and crosstalk, its expensive computational cost necessitates a regularization approach that can help us to get a noise-free imaging result more efficiently.
Angle-domain common-image gathers generated with the extended imaging condition (Sava and Fomel, 2003;Biondi and Symes, 2004;Sava and Fomel, 2006) are perhaps the most straightforward dimension for regularizing LSM because the reflection angle dimension is generally continuous.Prucha and Biondi (2002) and Kuehl and Sacchi (2003) apply smoothing regularization along the ray parameter axis in LSM to suppress migration artifacts caused by irregular sampling.Salomons et al. (2014) add a penalization term in LSM for the difference between adjacent angle traces in common-image gathers to suppress artifacts caused by incomplete illumination.Stanton and Sacchi (2015) regularize LSM of elastic data by dip filtering in the angle domain to reduce the effect of source/receiver sampling, noise, and crosstalk artifacts.Although common-image gathers provide a more straightforward domain to impose regularization for LSM, this requires more computation to construct common-image gathers and significantly larger memory to store the image volume.Fomel (2007) proposes to incorporate a shaping operator enforcing model compliant to the specified shape at each iteration into a conjugate-gradient algorithm and shows that shaping regularization provides better control on the estimated model and often leads to a faster iterative convergence (Fomel, 2009).We propose to incorporate shaping regularization into least-squares reverse time migration (LSRTM).We choose structure-enhancing filtering (Liu et al., 2010;Swindeman and Fomel, 2015) as a shaping regularization operator for effectively removing noise while preserving structural information.The extra computational cost caused by this model constraint is negligible when compared with the cost of LSM.For our linear modeling/migration operator pair, we choose lowrank one-step modeling and its adjoint (Sun et al., 2014).
The paper is organized as follows: First, we present the theory of LSRTM with shaping regularization (LSRTM-SR) in applica- tion to incomplete data and blended data.Then, we describe two synthetic examples of suppressing migration artifacts introduced by incomplete data and interference caused by blended data.These experiments demonstrate that the proposed method is highly effective in removing migration artifacts and crosstalk noise.

THEORY
Linear seismic modeling can be expressed in a compact linearalgebra notation as Lm ¼ d; (1) a) b) where L is the modeling operator, m is the reflectivity model, and d is the observed seismic reflection data.
Conventional migration corresponds to approximating the inverse of L with its adjoint (Claerbout, 1992): (2) where L T represents the migration operator.To implement L and L T , we use lowrank one-step modeling and RTM (Sun et al., 2014(Sun et al., , 2015a)).
The objective function of regularized LSRTM can be expressed as (3) where D is a regularization operator.Alternatively, Fomel (2007) defines the least-squares solution with shaping regularization as follows: where S is a model shaping operator, λ is a scaling parameter, and I is the identity operator.
The goal of regularization is to impose constraints on the estimated model.Shaping regularization implements it by explicit mapping of the estimated model to the space of the admissible model (Fomel, 2007;Xue and Zhu, 2015).In our approach to image simultaneous-source data, the estimated model contains crosstalk.The admissible model is the one that has reduced crosstalk.The two models can be bridged by a shaping regularization operator, which we construct as a structure-enhancing filter that attenuates noise along the dominant dip direction.The connection between equations 4 and 5 can be formally represented as When the shaping operator is symmetric and representable in the form S ¼ HH T , equation 5 can be rewritten as Following Fomel (2007), we perform inversion in equation 7 iteratively using the conjugate-gradient algorithm (Hestenes and Stiefel, 1952).

LSRTM with shaping regularization for incomplete data
Due to many physical and economic constraints, seismic data may have limited aperture, coarse wavefield sampling, and gaps (Nemeth et al., 1999).Direct imaging of incomplete data avoids the steps of interpolation, but it may suffer from intense artifacts.
In the case of irregularly sampled data, the forward operator in equation 1 is a combination of forward modeling operator L and sampling operator M. The new form of equation 1 is In our synthetic experiments, we define M as a simple masking operator.

LSRTM with shaping regularization for simultaneous-source data
Simultaneous-source acquisition allows more than one source to be activated at the same time regardless of the interference between different sources.This approach has attracted much attention from industry because it can help to reduce the acquisition cost and improve the quality of acquired seismic data (Berkhout, 2008).
For a constant-receiver survey, the simultaneous-source data can be expressed as where d is the blended data; M is the number of simultaneous sources; P M j¼1 T j denotes the blending operator (Berkhout, 2008), also known as the encoding function (Romero et al., 2000) or dithering (Chen et al., 2014); and d j represents an individual unblended shot gather, which can be modeled with the modeling operator L j associated with the jth shot: Thus, in the case of simultaneous-source data, the forward operator L in equation 1 becomes a combination of modeling operators L j and the blending operator P M j¼1 T j .

Review of structure-enhancing filtering
Structure-enhancing filtering is a nonstationary smoothing operator that follows local structural dips.Liu et al. (2010) realize it by applying an averaging filter along an extended dimension to remove the structurally nonconforming noise.The extended dimension is acquired by plane-wave prediction, which predicts a trace from its neighbors according to the local slopes of seismic events (Fomel, 2010).
Figure 1 illustrates the shaping process of structure-enhancing filtering.Figure 1a shows a simple velocity model.Figure 1b is a local dip estimate obtained using plane-wave destruction (Fomel, 2002).Figure 1c shows 16 points picked from velocity interfaces.We use the local dips shown in Figure 1b to smooth the points with structure-enhancing filtering and to generate the corresponding impulse responses.Figure 1d-1f shows three groups of impulse responses with increasing shaping strength controlled by the length of prediction (ns in Figure 1).As the length increases, the impulse responses become longer, and finally they blend together and reconstruct the structure of the velocity model.The directions of the impulse responses follow exactly the structural directions at the corresponding locations in Figure 1a, which justifies the effectiveness of structure-enhancing filtering as a shaping regularization operator.
In application to 3D imaging, dips of inline and crossline directions can be estimated, and 3D shaping regularization becomes a chain of shaping operations in inline and crossline directions (Fomel, 2007).

NUMERICAL EXAMPLES
Our first example is based on the velocity model shown in Figure 1a.We use forward modeling to generate a synthetic data set under conventional and simultaneous-source acquisitions.For simultaneous-source acquisition, two shots were excited at nearly the same time.Figure 2a shows the incomplete data obtained after we randomly removed 80% traces.Figure 2b indicates the simultaneous-source data.
In Figure 3, we compare the effects of different methods on migrating incomplete data.Figure 3a  constraints is shown in Figure 3b.It contains less artifacts, and the effect of poor illumination has been partially compensated.Figure 3c shows the fifth iteration result of LSRTM-SR.After five iterations, most of the noise has been removed; however, the amplitude of the events is still unbalanced.The 20th iteration result is shown in Figure 3d.We can find that the image becomes cleaner, and the energy of the weak reflectors has been compensated.
The structure-enhancing filtering used in LSRTM-SR requires local dip information.In our experiments, the dip is iteratively estimated at each iteration from the result of the previous iteration.Figure 4a and 4b shows the dips used at the first and last iterations of LSRTM-SR.The initial dip has some large values caused by migration artifacts on the left and right top corners.These erroneous values gradually disappeared over iterations, which is verified by the fact that the dip for the last iteration has the same pattern as the dip directly estimated from the exact velocity model, shown in Figure 1b. Figure 4c shows the direct application of the shaping operator to the adjoint RTM.The remaining artifacts are easily observed.Figure 4d shows the convergence curves of normalized data errors against the iteration number for LSRTM and LSRTM-SR.The two curves appear similar and exhibit a fast convergence rate.
Figure 5 indicates a series of results for the blended data case.Blended data cause strong interference between layers shown in Figure 5a.LSRTM without shaping regularization can remove part of the interference and compensate for the amplitude of reflectivity (Figure 5b).A crosstalk-free profile is obtained by LSRTM-SR after 20 iterations (Figure 5d).Its convergence curve (Figure 6d) is similar to that of LSRTM and reveals that the data error has decreased by more than 90% after 20 iterations.Figure 6a   In comparison, the image obtained by applying shaping operator to the output of RTM still has strong remaining crosstalk artifacts and cannot get any benefits of LSRTM (Figure 6c).Next, we use a more realistic Sigsbee 2A model (Paffenholz et al., 2002) to test the proposed method.Figure 7 shows the blended data set, which has 16 super shot gathers, with each super shot gather generated by firing eight sources at the same time.If directly using RTM to migrate the data, the result contains significant crosstalk, as shown in Figure 8a. Figure 9a shows the LSRTM result after 10 iterations, which has strong remaining crosstalk.Then, we switch to using LSRTM-SR.The tenth iteration result (Figure 9b) shows that the crosstalk has been greatly eliminated, while the structural information has been preserved.Figure 8b shows the dip estimated from the RTM image and used for the first iteration of LSRTM-SR.
In this experiment, we enforce different shaping strengths through the inversion: a stronger shaping (larger radius) at the first two iterations and a weaker shaping (smaller radius) at the later iterations.We do that to remove the strong crosstalk first, and then we use the subsequent inversion process to reduce the data misfit and to recover small-scale features, such as faults and diffraction points, that may get smeared by the shaping operator at early iterations.Figure 10 shows how this idea works.Comparing Figure 10a and 10b, we find that the crosstalk has been totally suppressed at the second iteration; however, shaping somewhat smoothed out the diffraction points.As shown in Figure 10c, this issue has been partially resolved by the subsequent inversion process.The diffraction points became focused at around the 10th iteration.
Figure 11 shows the incomplete data with 60% of its traces missing.Just as in the case of simultaneous-source data, the same conclusions can be reached for the incomplete data case.A comparison between Figure 12a and 12b demonstrates that LSRTM-SR can effectively suppress artifacts caused by the incomplete data.
Figure 13a shows the convergence history of the normalized data misfit against the iteration number for LSRTM and LSRTM-SR in the simultaneous-source data case.The two curves are very close, which means that the two methods have a similar data misfit decrease rate.However, we can deduce that LSRTM-SR has a faster model misfit convergence from the comparison shown in Figure 9. Figure 13b shows the convergence history curves for the incomplete data case.We can find that the data misfit of LSRTM-SR decreased by approximately 90% after 10 iterations.

CONCLUSIONS
We have proposed to incorporate shaping regularization imposing structural constraints on the estimated reflectivity in the process of LSM.Our numerical tests on synthetic models confirmed the applicability of the proposed method and showed that least-squares RTM with shaping regularization, while having a similar convergence rate to LSRTM, can effectively suppress migration artifacts caused by incomplete data and interference effects caused by simultaneous-source data.

Figure 1 .
Figure 1.Shaping by structure-enhancing filtering: (a) an example velocity, (b) local dip map, (c) 16 points selected as input, (d) impulse responses of structure-enhancing filtering for the 16 points (ns ¼ 20), (e) impulse responses for a larger prediction radius (ns ¼ 40), and (f) impulse responses for an even larger prediction radius (ns ¼ 90).

Figure 2 .Figure 3 .
Figure 2. Synthetic data for model in Figure 1a: (a) incomplete data with 80% traces missing and (b) blended shot gathers.

Figure 4 .
Figure 4. Dip estimation and convergence curve for an incomplete data case: (a) initial dip estimated from the RTM image (Figure 3a), which is used at the first iteration of LSRTM-SR, (b) dip used at the last iteration of LSRTM-SR, (c) image obtained by directly applying a shaping operator to the RTM image with the initial dip in Figure 4a, and (d) data-misfit convergence curves of LSRTM (dashed line) and LSRTM-SR (solid line).

Figure 6 .
Figure 6.Dip estimation and convergence curve for a simultaneous-source case: (a) initial dip estimated from the RTM image (Figure 5a), which is used at the first iteration of LSRTM-SR, (b) dip used at the last iteration of LSRTM-SR, (c) image obtained by directly applying the shaping operator to the RTM image with the initial dip in Figure 6a, and (d) convergence curves of LSRTM (dashed line) and LSRTM-SR (solid line).

Figure 7 .
Figure 7. Synthetic simultaneous-source data for the Sigsbee model.