Fast Least-Squares Reverse Time Migration of OBN Down-Going Multiples

Nowadays the ocean bottom node (OBN) acquisition is widely used in oil and gas resource exploration and seismic monitoring. Conventional imaging algorithms of OBN data mainly focus on the processing of up-going primaries and down-going first-order multiples. Up-going multiples and higher-order down-going multiples are generally regarded as noise and should be eliminated or ignored in conventional migration methods. However, multiples carry abundant information about subsurface structures where primaries cannot achieve. To take full advantage of multiples, we propose a migration method using OBN down-going all-order multiples. And then the least-squares optimization algorithm is used to suppress crosstalks. Finally, a phase-encoding-based migration algorithm is developed to cut down the computational cost by blending several common receiver gathers together using random time delays and polarity reversals. Numerical experiments on the complex Marmousi model illustrate that the developed approach can enlarge the imaging area evidently, reduce the computational cost effectively, and enhance the image quality by suppressing crosstalks and improving the resolution.


INTRODUCTION
Many advanced imaging algorithms have been developed over the past few decades, among which reverse time migration (RTM) can provide quality images of subsurface complex structures aided by the two-way wave equation. Conventional RTM only accounts for primaries and regards multiples as noise. RTM using the raw data containing primaries and multiples will produce serious crosstalks in non-structural subsurface areas, which results in the reduction of the signal-to-noise (S/N) ratio of final imaging results and the incorrect geological interpretation. Therefore, many researchers concentrated on the suppression of multiples (Verschuur et al., 1992;Dragoset et al., 2010;Liu et al., 2010). However, multiples also carry lots of useful information about strata. Compared with primaries, multiples can provide extra illumination in areas where primaries cannot achieve for their smaller reflection angles and longer travel paths (Berkhout and Verschuur, 2006;Liu et al., 2011;Lu et al., 2015). To obtain higher-precision imaging results of subsurface structures, geophysical researchers have developed a variety of imaging approaches using multiples (Jiang et al., 2007;Muijs et al., 2007;Liu et al., 2011Liu et al., , 2015Verschuur and Berkhout, 2011;Lu et al., 2015;Liu X. et al., 2016;Zhang et al., 2019). Among these techniques, reverse time migration of multiples (RTMM) Liu et al. (2011) combining the advantages of multiples and RTM can correctly migrate the multiple reflections and image the subsalt structures. RTMM is currently one of the easiest and most widely used algorithms to migrate multiples (Zuberi and Alkhalifah, 2013;Liu Y. et al., 2016;Li et al., 2017;Zhang et al., 2020).
Due to the existence of interfaces with large impedance differences such as the air-water surface, seabed, and salt boundaries, marine seismic data often contain plenty of multiples with strong energy. The imaging results with satisfactory accuracy cannot be obtained if multiples are not effectively processed. Traditionally, marine seismic data are acquired using the towed streamer, ocean bottom cable (OBC), or OBN. Equipped with a geophone and a hydrophone, OBN is a seismometer that is placed at the seabed. For the flexible deployment and the abundant wavefield information, the OBN acquisition is widely used in oil and gas resource exploration, seismic monitoring, detection of deep structures, etc. (Katzman et al., 1994;Dash et al., 2009;Mathias et al., 2009). In the OBN acquisition, the record is composed of up-going and down-going components, as shown in Figure 1. The up-going record includes the up-going primaries and multiples, while the down-going record contains the direct wave and down-going multiples. And the two components can be separated using the hydrophone recording (P) and the vertical component of the geophone (Z) aided by PZ summation technique (Barr and Sander, 1989;Schalkwijk et al., 2003).
Traditional imaging algorithms using up-going primaries have limits in imaging the seabed and the illumination range. The migration using the down-going first-order multiples is developed to solve the problems mentioned above. And it is also called mirror migration (Grion et al., 2007;Dash et al., 2009;Wong et al., 2015). In mirror migration, the virtual OBN position is obtained by turning the real OBN up with the sea surface as a mirror, followed by implementing conventional migration procedures. Regardless of the conventional migration or the mirror migration, the up-going multiples and down-going high-order multiples are regarded as noise and should be suppressed before or during migration. Compared with primaries and first-order down-going multiples, higher-order multiples possess longer travel paths and sufficient smaller reflection angles, which can guarantee that imaging using multiples can supply quality results with a wider illumination range (Liu et al., 2011;Lu et al., 2015;Wong et al., 2015). To make full use of the structural information contained in multiples and further increase the illumination range, we proposed RTM using down-going multiples for OBN acquisition, in which the lowerorder multiples are recognized as the secondary sources of the higher-order multiples. However, the biggest challenge of RTMM is the suppression of crosstalks (Liu et al., 2011). Further, the least-squares algorithm is introduced to suppress the crosstalk noise (Nemeth et al., 1999;Dai et al., 2012;Liu et al., 2017;He et al., 2019).
The calculation cost of the migration process is linear with the total number of shot gathers. Therefore, combining several shot gathers into a supergather can effectively reduce the calculation unit and improve the calculation efficiency. However, when migrating one supergather, the interactions between the unrelated sources wavefield and the receivers wavefield will cause serious crosstalk noise, which decreases the S/N ratio of the imaging results tremendously (Romero et al., 2000;Dai et al., 2012). For example, by combining two shot gathers into one supergather and then implementing the migration, the crosscorrelation between the source wavefield activated by the first source and the receiver wavefield generated by the unrelated second backward-penetrating receiver data will produce strong crosstalks. To mitigate the above-mentioned crosstalk noise, Romero et al. (2000) first applied the phase encoding strategies in the migration procedures and showed a variety of encoding strategies. Among the commonly used encoding strategies, the random phase encoding technology Schuster et al. (2011), combining random time delays and random polarity reversals, has a better crosstalk suppression effect because it can shift and discrete crosstalk noise randomly. Krebs et al. (2009) introduced the random phase encoding scheme into the full waveform inversion and obtained results of the same quality as conventional methods, while the amount of calculation is reduced by an order of magnitude. The multisource least-squares reverse time migration (LSRTM) proposed by Dai et al. (2012) is even more efficient than RTM, although there are some random artifacts in the imaging results. It is worth noting that to avoid convergence problems caused by the incomplete matching between the predicted data and observed data, the application of random phase encoding is generally limited to the fixed-spread acquisition geometries (Krebs et al., 2009;Dai et al., 2012;. In the RTM of OBN down-going multiples, the lower-order multiples are recognized as the sources of the higher-order multiples, and the lower-and higher-order multiples share the same receivers. Therefore, the random phase encoding scheme can be introduced into the migration of multiples for OBN acquisition conveniently. To save the computational cost, we developed a phase-encodingbased LSRTM of OBN down-going multiples.
In this paper, we first specify the principle of RTM using OBN down-going multiples and then introduce least-squares technology to suppress crosstalks, followed by the explicit description of the phase-encoding-based LSRTM of OBN downgoing multiples. Numerical experiments on the Marmousi model are used to verify the effectiveness and efficiency of the proposed approach, and the results demonstrate that the proposed approach can effectively suppress crosstalks, improve the imaging resolution, and save the computational cost. Finally, the discussion and conclusion are provided.

Common Receiver Domain RTM of OBN Down-Going Multiples
The seismic record in OBN acquisition consists of up-going and down-going components, which can be separated using PZ summation: where U and D denote the up-going and down-going records, respectively. P and Z represent the pressure recorded on the hydrophone and particle velocity recorded on the vertical geophone, respectively. ρ and c are the density and the acoustic velocity of water, respectively. As shown in Figure 2, the basic principle of PZ summation is that the geophone and the hydrophone data record opposite polarity for the down-going propagating wavefield. If the wavefield is propagating upward, they will record the same polarity.
In OBN acquisition, receivers are placed at the seabed sparsely, while sources are activated on the sea surface densely. For the sake of computational efficiency, the seismic reciprocal principle Knopoff and Gangi (1959) is implemented and migrations for OBN data are usually performed in the common receiver gather (CRG) domain. The data is converted from the common shot gather (CSG) into the CRG by exchanging the source and receiver locations during imaging. After conversion, the sources and nodes are regarded as the virtual receivers and the virtual sources, respectively. The conventional migration algorithms for OBN data often account for up-going primaries or down-going first-order multiples, as shown in Figures 3A,B. For the imaging of up-going primaries, the virtual sources are arranged at the seabed, which results in serious aliasing noise around the real node locations. And the subsurface coverage area of primaries is limited to the vicinity of the node positions. The mirror migration, utilizing down-going first-order multiples, can image the seabed well and enlarge the imaging range by turning the OBNs to the mirror OBNs taking the sea surface as the axis of symmetry. Both the abovementioned migration approaches regard the up-going multiples and down-going high-order multiples as the data noise and suppress them before migration. The structural information of high-order multiples has not been fully exploited. To take full advantage of higher-order multiples, we propose RTM using down-going multiples for the OBN acquisition. Prior to the migration procedure, the down-going record is separated using PZ summation, and then the down-going multiples are acquired by muting the direct wave from the down-going record. The RTM of down-going multiples can be accomplished through the following three steps: 1 Producing the source wavefield by forward-propagating the separated down-going record including the direct wave and multiples: where x r is the horizontal location of virtual receivers, v 0 (x, z) is the migration velocity model, and P s (x, z, t) is the generated source wavefield. d(x r , t) and m(x r , t) are the direct wave and the down-going multiples in the CRG domain, respectively.
2 Backward-extrapolating the down-going multiples to obtain the receiver wavefield M r (x, z, t): 3 Implementing the zero-lag crosscorrelation imaging condition to the source and receiver wavefield: where I(x, z) represents the migrated image and t max is the total recording time. Compared to the imaging point of mirror Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730476 migration (the red solid point in Figure 3B), the imaging point of the proposed RTM using down-going multiples in Figure 3C is located at the left of that in Figure 3B (the shadow zone of mirror migration), which means wider illumination range. And theoretically with the order of multiples gets larger, the illumination can be further expanded. Unlike conventional RTM, RTMM correlates the complex forward-and backward-penetrated records, thus the interactions among unrelated multiples will cause undesired crosstalk artifacts which critically degrade the S/N ratio of the final imaging results. Figure 4 illustrates the representative crosstalks in RTMM for OBN acquisition. In Figure 4, the red solid lines, the black solid lines, and the red dotted lines represent the wavepaths of the forward-propagating data, the backwardpropagating data, and the crosscorrelation, respectively. As shown in Figure 4, the forward-propagating direct wave recorded at S3 can be regarded as the secondary source of the backward-propagating second-order multiples recorded at S1, and the crosscorrelation between the two multiples will produce crosstalk point A1 denoted by the magenta point. In Figure 4B, the crosscorrelation between the forward-propagating first-order multiples recorded at S2 and the backward-propagating firstorder multiples recorded at S1 will cause the crosstalk point A2.

LSRTM of OBN Down-Going Multiples
Compared with RTM, the LSRTM has been extensively researched to enhance the imaging quality by suppressing RTM imaging artifacts, improving the imaging resolution, and balancing the energies of imaging events. LSRTM also can be extended to the migration of multiples to suppress the severe crosstalks of RTMM. Unlike the conventional LSRTM using primaries, LSRTM of multiples (LSRTMM) Zhang and Schuster (2014), Wong et al. (2015), Liu X. et al. (2016) regards multiples data as the secondary source and tries to find a reflectivity model that minimizes the misfit function between predicted and observed multiple reflections. The Born modeling operator (see Appendix A) that corresponds to the multiples can be represented as FIGURE 4 | Illustrations of representative crosstalks in RTMM. As shown in panel A, the crosscorrelation between the forward-propagating direct wave recorded at S3 and the backward-propagating second-order multiples recorded at S1 will generate the crosstalk point A1 denoted by the magenta point. In panel B, the first-order multiple recorded at S2 can be regarded as the secondary source of the backward-propagating first-order multiple recorded at S3, and the crosscorrelation of the two multiples will cause crosstalk point A2.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730476 Eq. (5) contains three steps: 1) The down-going record of OBN data d(x r , t) + m(x r , t) are treated as the virtual sources and injected into the subsurface background model v 0 , which will produce the incident wavefield M.
2) The interaction between the incident wavefield M and the subsurface reflectivity R will generate the scattered wavefield M r . 3) Through saving the boundary wavefield of the scattered wavefield M r , we can obtain the predicted down-going multiples m pre . Eq. (5) can be depicted compactly using a matrix-vector notation: where L i and m pre,i represent the linear Born modeling operator and the predicted OBN down-going multiples for the ith CRG, respectively. The RTM operator can be treated as the adjoint of the forward-modeling operator, thus the final migrated image R mig can be obtained by: where the superscript T is the conjugate transpose operation. L T i and m i are the migration operator and the observed down-going multiples for the ith CRG, respectively. Based on the Born modeling operator defined in Eq. 5 and its adjoint RTM operator in Eq. (7), we can minimize the difference between the modeled/predicted multiples and the observed multiples via the misfit function f (R), which can be iteratively solved aided by any gradient-based inversion approach.

LSRTM of OBN Down-Going Multiples With Random Phase Encoding
However, in LSRTM, each iteration contains migration, demigration and gradient calculation steps, the calculation cost of which is several times of that of RTM. And the computational cost of LSRTM for OBN down-going multiples may be larger for more iterations are needed to achieve better crosstalk artifacts suppression effect. Therefore, we introduce the random phase encoding scheme into LSRTM of OBN down-going multiples to improve the computational efficiency. By encoding the forward modeling operators (in Eq. (6)), OBN down-going multiples of different CRGs can be simultaneously modeled as [α 1 + α 2 + / + α n ]m pre,i [α 1 L 1 + α 2 L 2 + / + α n L n ]R where α n is the phase encoding function. Eq. (9) can be written in a compact form of m pre LR where L and m pre represent the encoded forward modeling operator and predicted data for one supergather, respectively. The image of a supergather of observed OBN down-going multiples m can be produced by where R mig represents the phase-encoded image of OBN down-going multiples. Eq. (12) is the expanded form of where for any integer value i, we assume that α p i α i 1. In Eq. (12), the valid term, which is identical with Eq. (7), will produce the correct imaging result and crosstalks. The noise term only generates crosstalk noise, which is caused by the interaction between the encoded CRGs. For phase-encodingbased RTM of OBN down-going multiples, the crosstalk noise generated by the noise term is much stronger than the noise generated by the valid term . In the numerical experiment section, we will illustrate the different crosstalks caused by the valid and noise terms of Eq. (12). Through encoding ns supergathers from original CRGs and based on Eqs (10) and (11), the objective function of phaseencoding-based LSRTM of OBN down-going multiples is expressed as The implementation details of the proposed migration approach are illustrated in Figure 5. One of the key steps to the proposed approach is the choice of the phase encoding functions. In this paper, the phase encoding function is designed as where τ n and S n denote the random time delay and random polarity reversal for a blended CRG, respectively. The random time delays τ n follows the normal distribution with a standard deviation of σ and a mean of zero, expressed as: The probability that the polarity reversal function S n takes a value of 1 or -1 for each blended CRG is 1/2. According to the property of normal distribution, the polarity reversal function S n can be defined by the following: where sgn represents the sign function. If τ n ≥ 0 S n 1, otherwise S n −1.

Numerical Examples
In this section, we apply the proposed phase-encoding-based LSRTM of OBN down-going multiples on the Marmousi model. A high-order staggered-grid finite-difference algorithm with PML absorbing boundaries is utilized to solve the first-order velocitystress acoustic wave equations and produce the OBN data. The results of conventional RTM using primaries and mirror migration are also used to verify the effectiveness of the proposed approach. The Behavior of Conventional LSRTM of OBN Down-Going Multiples Figure 6 shows the Marmousi velocity model. The model is discretized into 351 grids in the vertical direction and 1,001 grids in the horizontal direction both with a 10 m grid interval. The seismic source is selected as the Ricker wavelet with a peak frequency of 20 Hz.
A set of 1,001 CSGs with a source interval of 10 m are simulated to certify the performance of the novel approach. The total recording time of each CSG is 4s with a sampling interval of 4 m s. First, we make use of 20 OBNs, distributed evenly from 2.98 to 6.78 km in the horizontal direction with a 200-m node interval, to test the imaging performances of RTM and LSRTM of down-going multiples. The pressure components of a CSG and a converted CRG are illustrated in Figure 7A,B respectively. The up-going and down-going records after PZ summation are displayed in Figure 8A,B, respectively. The downgoing multiples after muting the direct wave from the down-going record (Figure 8) are shown in Figure 8C.
In the CRG domain, since the seismometers are placed at the seabed and they act as the virtual sources, there is strong aliasing noise around the seabed in the image of the conventional RTM using up-going primaries as shown by the arrow in Figure 9A. Moreover, the illumination range of Figure 9A is limited to the vicinity of the node positions due to the insufficient subsurface coverage of upgoing primaries. To settle the above-mentioned problems, the RTM using the down-going first-order multiples (mirror migration) is developed. As shown in Figure 9B, the imaging quality of the seabed has been greatly improved pointed out by the arrow, and the illumination range has also been expanded compared with the imaging result of up-going primaries in Figure 9A. To make full advantage of higher-order multiples, the proposed RTM of OBN down-going multiples algorithm is implemented, and the imaging result is displayed in Figure 9C. Compared with the results in Figures 9A,B, the imaging result of all-order multiples in Figure 9C shows a much wider illumination range. However, the interferences among unrelated multiples cause crosstalk artifacts which degrade  the imaging fidelity as denoted by the arrows in Figure 9C. The kind of crosstalks is inherent in conventional RTMM results (Eq. (7) or the valid term of Eq. (12)) (Liu et al., 2011). Compared with RTM of down-going multiples, the developed LSRTM of down-going multiples can suppress the crosstalks effectively as pointed out by the arrows in Figures 9C,D, and produce a higher-resolution image ( Figure 9D).

The Behavior of the Proposed Phase-Encoding-Based LSRTM of OBN Down-Going Multiples
Subsequently, the node interval is reduced from 200 to 100 m, and a total of 40 nodes are evenly placed at the seabed from 2.98 to 6.88 km.
After the simulation and PZ summation, 40 CRGs are encoded into 1,  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730476 8 2, 4, 8, 20 supergathers to verify the effectiveness and efficiency of the proposed phased-encoding-based algorithm. The random encoding function which combines the random time delays and polarity reversals is utilized to blend different CRGs. The random time delay sequences distribute normally with a 0.5 standard deviation and zero mean. The blended supergathers of 1, 2, 5,10, 20, and 40 CRGs of down-going multiples are illustrated in Figure 10, respectively. Images generated by RTM of up-going primaries, down-going firstorder multiples (mirror migration), and down-going multiples are shown in Figures 11A-C, respectively. Figure 11D represents the result of LSRTM of down-going multiples. Compared with the images with the 200 m node interval in Figure 9, the corresponding images with the 100 m node interval provide more sufficient subsurface coverage, more continuous strata, and a higher S/N ratio.
In the following, we use five sets of supergathers that regard the number of encoding CRGs as the variable to present the effectivity and efficiency of the proposed phase-encoding-based LSRTM of down-going multiples. Figures 12-16 display the RTM and LSRTM (10 iterations) images using the supergathers in Figures  10B-F, respectively. Compared with the RTM, the LSRTM can effectively suppress the noise caused by the interferences among encoded CRGs (the noise term of Eq. 12) and significantly improve the imaging resolution as shown by the arrows in panels A and B of Figures 12-16. From the five pairs of migration results, as the number of encoding CRGs decreases, the imaging results are getting closer to those of conventional migrations ( Figures  11C,D). The trend is caused by the decrease of the energy proportion of the noise term in Eq. (12). However, as the number of encoding CRGs decreases, the computational cost will increase. In the last two experiments, the original 40 shot gathers are encoded into 8 and 20 supergathers, respectively. And the final results (Figures 15, 16) show better imaging qualities and a higher S/N ratio compared with those results in Figures 12-14. Therefore, we choose them to conduct the calculation comparison with conventional LSRTM of down-going multiples. Scaling all calculations to the forward modeling, each iteration of the LSRTM requires 6 forward modelings including two for the RTM, two for the de-migration, and two for the gradient calculation. The conventional LSRTM of down-going multiples with 40 supergathers and 10 iterations (see Figure 11D) requires a computational cost of 40 × 6 × 10 2,400 forward modelings.
The computational costs for LSRTM images in Figure 15B (8 supergathers) and Figure 16B (20 supergathers) are 8 × 6 × 10 480 and 20 × 6 × 10 1,200 forward modelings, respectively. The comparison of computational cost demonstrates that the proposed phased-encoded LSRTM method can provide the results which can meet the requirement of interpretation and improve the computational efficiency manyfold.  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730476 10 DISCUSSION Up to now, there are two widely used crosstalk suppression methods in the migration of multiples. The first one is the LSRTMM which introduces least-squares optimal algorithms to suppress crosstalks and obtains higher-precision images of multiples migration. The second one is the RTM of controlledorder multiples (RTM-CM) in which the different orders of multiples are extracted before migration procedures and the separated nth-order and (n+1)th-order multiples respectively   Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 730476 11 act as the source and receiver data. RTM-CM can significantly reduce the crosstalks due to the simplifies of the forward-and backwardpropagating data. Liu Y. et al. (2016) separated different orders of multiples using surface-related multiple elimination (SRME) and then proposed LSRTM-CM. Compared with LSRTMM, LSRTM-CM obtains images with a higher S/N ratio. When extracting multiples of different orders by conventional SRME, the locations of sources and receivers must keep coincident. However, in OBN acquisition, shots are activated densely on the sea surface and seismometers are sparsely placed at the seabed, which makes it difficult to isolate different-order up-or down-going multiples. Therefore, in this study, LSRTMM is utilized to suppress crosstalks rather than RTM-CM. The least-squares algorithms can not only suppress the inherent crosstalks of multiple imaging but also reduce the noises caused by interactions among the encoding CRGs.
In this paper, we select the static phase encoding strategy, that is, during the entire least squares implementation the encoding function remains constant, which will result in the undesired suppression of the crosstalk noise generated by the encoded CRGs. In the following research, we will design an optimal dynamic phase encoding scheme to achieve a better noise suppression effect. It should be illustrated that we assumed that the geophones and the hydrophones are placed above the seabed. Consequently, the acquired synthetic data only contain P-wave components. And the pressure component P and the particle velocity component Z directly can be obtained directly by using the first-order velocity-stress acoustic wave equation. S-wave information is the key to the data processing of OBN acquisition, and the availability and efficiency of the proposed method on S-wave remain to be verified. In the following researches, we will concentrate on the exploitation of S-wave information.

CONCLUSION
Conventional imaging approaches of OBN data account for the up-going primaries and the down-going first-order multiples. However, OBN data contain the up-going multiples and down-going higher-order multiples, and the information carried by those multiples could also be used to enhance the imaging quality. In this paper, we developed a CRG-domain LSRTM of OBN down-going multiples with a random phase encoding scheme in which the full downgoing data and down-going multiples are recognized as secondary sources and observed data. The encoding function is composed of random time delay sequences and random polarity reversals which are designed as a normal distribution. Synthetic experiments demonstrate that the proposed approach can efficiently produce images with a wide illumination range and high resolution. Our approach has the potential to serve as a powerful tool for imaging OBN data of complex subsurface structures.

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.

AUTHOR CONTRIBUTIONS
YZ and YL contributed as the major authors of the manuscript. JY did a part of writing and coding works. XL provided some constructive ideas and polished the manuscript. All authors contributed to the article and improved the submitted version.