Depth Migration Based on Two-Way Wave Equation to Image OBS Multiples: A Case Study in the South Shetland Margin (Antarctica)

With the development of marine seismic exploration, the ocean bottom seismometer (OBS) as a new seismic acquisition technology has been widely concerned. Although multiple waves are frequently viewed as noises, they may carry a wealth of subsurface information and produce a broader illumination than primary waves. To perform multiple wave imaging, we propose to utilize a two-way wave equation depth wavefield extrapolation method which is rarely used in this field. A simple dipping model is imaged by using primary and multiple waves, which proves the superiority of multiple waves in imaging over the primary waves and lays a foundation for practical application. Moreover, the comparison of multiple imaging results by reverse time migration and those by our proposed method demonstrates that our proposed method requires less storage space. In this study, we apply this migration method to actual OBS data collected in the South Shetland margin (Antarctica), where gas hydrates have been well documented. Firstly, the wavefield separation method is adopted to process the OBS data, so as to produce reliable primary and multiples waves; secondly, the ray-tracing method is used to derive the velocity field; and finally, the depth wavefield extrapolation method based on the two-way wave equation is applied to image primary and multiple waves. Migration results show that multiple waves provide a broader illumination and a clearer sediment structure than primary waves, especially for the highly shallow reflections.


Introduction
In recent years, a multicomponent ocean bottom seismometer (OBS) technique has been widely used in gas hydrate exploration. OBS data have the following advantages over conventional near-surface recording: (1) OBS is deployed on the seafloor and can therefore record both P-waves and converted PS-waves; (2) OBS provides the geometry of wide azimuth, which is of great significance for obtaining good imaging under some complex overburden; (3) OBS data have a higher signal to noise ratio. Even so, OBS data processing still faces many difficulties. For example, the primary waves are usually used to image OBS data. However, due to the ray path geometry, OBSs generally produce a narrow illumi-nation by traditional migration, particularly for those shallow reflections with a depth below the seafloor less than the OBS interval. If OBSs are located at a great spacing, this problem will be aggravated and a large gap between two OBSs could occur on the imaging. It can be improved by dense surface acquisition, but the cost of acquisition is very high. For this reason, migration with OBS multiple waves was proposed by Dash et al. [1]. As multiple waves penetrating into the subsurface several times have longer traveling paths and contain abundant reflection information, they can produce a broader illumination of subsurface than primary waves. Water-layer multiples have been used by several authors in the migration of OBS data (i.e., [2][3][4][5][6]). Recently, Kazuya et al. [7] proposed to use seismic interferometry to address the OBS data for wide-angle imaging. In order to improve the imaging quality, Tong et al. [8] used the OBS data to image the deep water area and achieved a better result compared with conventional methods. However, most researches have used the conventional migration methods such as the Kirchhoff one-way wave equation migration to deal with OBS data. Few researches have introduced the newest migration into the application of OBS data, which is the topic to be discussed in this paper.
Seismic migration is one of the vital data processing in the academic and industrial circles. Researchers have carried out a branch of works on seismic migration. At present, migration based on wave equation has been widely used and can be split into three categories, among which, the most deeply researched one is the one-way wave equation method, which offers a kinematic solution for the acoustic wave equation, such as the split-step Fourier method [9][10][11][12][13]. Because the conventional migration methods of one-way wave equation are frequently based on the limited Taylor expansion, their shortcoming is the limited imaging angles. To overcome this restriction, the optimization methods are utilized to recalculate coefficients of partial fraction, thus achieving better imaging performance [14][15][16]. But this improvement has its limits, and these methods seem to make it difficult to make the imaging angle reach 90 degrees. As a very important migration method, reverse time migration (RTM) is a hot issue in imaging researches, because it can handle various waves in arbitrarily velocity models including turning waves. More importantly, it has no dip limitations [17,18], so it can image complicated models that are unable to be imaged by one-way wave equation. Although RTM has its advantages, it also has weak points. The most challenging problem it faced is storage. In order to compute the imaging amplitudes, RTM needs to restore all of the wavefields, which will require a huge memory cost. Researchers are seeking different techniques to solve this problem, such as proposing the random boundary method to avoid the storage of wavefields [19]. Another big problem of the RTM method is the lowfrequency artifacts existing in the imaging sections. To address this issue, Yoon and Marfurt [20] proposed a Laplacian filter. Although RTM has some weaknesses and it still has a long way to go to solve them completely, it is still a concern in the field of imaging.
A two-way wave equation-based depth migration method as the latest development is a second-order partial differential equation in respect of coordinates, so two boundary conditions are needed to make it solvable in the domain of depth [21,22]. It is very hard to meet the requirements of boundary conditions, so researchers have separated the two-way wave equation mathematically and developed one-way wave equation migration methods. In order to handle the boundary conditions, You et al. [23] presented a dual-sensor seismic acquisition system and performed the wavefield depth migration on the base of the two-way wave equation. Because the two-way wave equation depth migration method is on the base of fullwave equation, theoretically it also has no limited imaging angles. However, compared with the RTM method, it needs less storage memory and produces less low-frequency artifacts. Therefore, it is a promising migration method.
In this study, we propose to image OBS multiple waves with the two-way wave equation depth migration method, which is the first application to OBS data. Here, we present an application of this method to OBS data collected in the South Shetland margin (Antarctica) (Figure 1). The South

Geofluids
Shetland margin lies in the northeastern tip of the Pacific margin of the Antarctic Peninsula. It marks the convergence zone boundary where the Antarctic and the former Phoenix plates are subducting beneath the South Shetland microcontinental block. The continental margin shows a complicated tectonic setting, where a trench-accretionary prism-fore-arc basin sequence can be identified [24,25]. The subduction process is now believed to take place due to sinking and rollback of the oceanic plate, coupled with the extension of the Bransfield Strait marginal basin [25][26][27][28]. In this area, seismic data collected in the Italian Antarctic cruises have shown the presence of a possible gas hydrate reservoir along this margin [29][30][31][32]. Like tight sandstone, gas hydrate is an unconventional energy source with great economic potential [33,34].
Recent studies have pointed out that long-term ocean warming may have an influence on the stability of gas hydrate reservoir [35]. If gas hydrates decompose and release methane, it could lead to climate change. As a result, it is of great significance to increase our information on gas hydrates, including a good image of the gas hydrate system in this area.

OBS Data Acquisition and Processing.
In this paper, the seismic data were collected along the South Shetland margin in the austral summer 2003/2004. The data were concentrated on the continental slope with a strong bottom simulating reflection (BSR) recognized in the second survey. The two major objectives of this survey were to confirm the presence of a possible hydrate reservoir and to rebuild the tectonic background of this area. There were two OBSs located along the multichannel seismic line BSRstar8, as shown in Figure 1.
Two GI guns were used as the seismic source with an overall volume of 3.5 L. They were dragged to a depth of 8 m and fired at an interval of 50 m. The OBS comprised a threecomponent geophone and a hydrophone. Vertical and horizontal components are usually used to study compressional and shear waves, respectively. The data were recorded up to 20 s every 2 ms.
During the OBS investigation, its real position on the seafloor is uncertain, since it might displace quite a few hundreds of meters away from the stationing position where it was designed to be on the sea surface, considering oceanic currents and water depths. Recognizing the accurate seafloor position of an OBS is a precondition to a precise velocity field and a satisfying image section because velocity analysis is quite susceptible to the error of OBS position. Thus, the first procedure we have to operate is OBS relocation. In most cases, we utilize direct waves to estimate the OBS location on the seafloor. The relocation of OBS is an inverse issue, the purpose of which is to locate the OBS position as well as to reduce the residual between calculated and observed travel times of direct waves as much as possible. The accurate seafloor positions of OBSs were calculated by the direct arrivals from the different shots, presuming a consistent water column velocity, while the depth of seafloor was obtained from bathymetric data. The relocation result suggests that OBS6 hovered approximately 750 m away, at right angles to the shooting line, at a depth of 1790 m, which might be associated with powerful water currents [37]. To test the quality of the relocated OBS position, the linear normal travel time correction for the water column was conducted on OBS data. If the position was accurate, direct waves should appear flat after the correction as a result of the elimination of the water column with offset ( Figure 2).
We applied an amplitude correction for spherical divergence as well as a bandpass filtering (10-100 Hz) to the vertical and hydrophone component, so as to improve the data quality and strengthen the recognition of phases. Moreover, to reduce the ringing noises, we performed a predictive deconvolution with an operator length of 160 ms and a step of 25 ms to the vertical component data. Nevertheless, the deconvolution was not quite effective, still accompanied by noticeable ringing noises. In this case, we chose the hydrophone component instead of the vertical component to recognize phases. The vertical component, instead of the hydrophone component, showed an obvious ringing phenomenon which was the result of the effect of coupling between the seafloor and the instrument. Furthermore, the

Wavefield Separation.
Before migration, the primaries and multiples of OBS data need to be separated. The most commonly used method is to separate the wavefields into upgoing as well as downgoing wavefields by merging the hydrophone with the vertical components (dual-sensor or PZ summation; [38]). In practice, the hydrophone and the vertical geophone usually show different instrument responses and instrument sensitivities, as well as different coupling characteristics with the seafloor; thus, the two component data have differences in amplitudes and phases. Therefore, it is necessary to perform a match on the two components prior to PZ summation [39]. After calibrating the vertical component to the hydrophone component, the upgoing and downgoing wavefields can be separated as where U is upgoing wavefield, D is downgoing field, P is hydrophone component, Z cal is vertical component after calibration, ρ is water density, and c is acoustic velocity. OBS is deployed at the seafloor, the boundary between acoustic and elastic mediums, so the wavefields just above and below the seafloor needed to be separated [39][40][41]. After the initial signal is generated by a source and travels via the water column, it is partly reflected when meeting the watersediment interface and partly transferred into the seafloor. For instance, the travel time of the direct arrival is almost equal to the time it takes for the wave to be reflected at the seafloor and arrives at the receiver just above the seafloor.
Above the seafloor, the direct wave consists of upgoing and downgoing wavefields, while it shows only the downgoing wavefield beneath the seafloor, as well as the multiples at receiver sides. Analogously, when a reflection event propagates upward and meets the water-sediment interface, it is partly transferred into the water column and partly reflected into the seafloor. Accordingly, the reflected signal consists of upgoing and downgoing components just beneath the seafloor, but it has only an upgoing component just above the seafloor. In this study, we performed the wavefield separation on the base of the method as described in Ref. [1].

OBS Multiple Wave Imaging Based on Two-Way Wave Equation Migration
Method. In a 2D case, the acoustic wave equation with initial boundary conditions can be written as follows [22]: where uðx, z, tÞ is the wavefield at time t, vðx, zÞ is the medium velocity, Δz is the grid space and the continuing step, dðx, 0, tÞ and dðx, Δz, tÞ are the recorded seismic data at z = 0 and z = Δz depth levels, and u z is the derivative wavefield. It is well to know that two initial boundary conditions are required because the two-way wave equation is a second-order partial derivative differential equation.
In conventional marine exploration, we usually collect wavefields at the surfaces. Therefore, only one boundary condition can be provided. In order to handle this issue, because the velocity of sea body is constant, the common approach is to use a migration method based on one-way wave equation to calculate the derivative wavefields, which can be expressed as where k z = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðω 2 /c 2 Þ + ð∂ 2 /∂x 2 Þ p and dz is the extrapolated step. C is the velocity of sea body.ũðx, z = 0, ωÞ andũ z ðx, z = 0, ωÞ are the results of Fourier transform of uðx, z = 0, tÞ and u z ðx, z = 0, tÞ on time, respectively.
It seems that directly solving equation (2) in the depth domain is very hard. We attempt to transfer the time variable of equation (1) into the frequency domain, which can be expressed as where ω is angular frequency and uðx, z, ωÞ is frequencydomain wavefield.   4 Geofluids Combining equation (4) and two boundary conditions, the two-way wave equation used for depth extrapolation scheme can be expressed as wheredðx, 0, ωÞ anddðx, Δz, ωÞ are the results of Fourier transform of dðx, 0, tÞ and dðx, Δz, tÞ on time, respectively [22]. All in all, we can use the wavefields recorded at the surface to compute the derivative wavefields and then utilize equation (5) to perform wavefield depth extrapolation recursively. Obviously, equation (5) is a first-order partial derivative differential equation, and a classical Runge-Kutta method is used to approximate it.
To perform wavefield depth extrapolation, how to calculate L operator is a key core. In order to handle this issue, we need to discretize the Helmholtz operator L as a matrix form at a certain frequency point for a medium with arbitrary changes in velocity: where , Analyzing equation (6), we can find that L operator is a symmetric matrix. So matrix decomposition can be implemented, which can be written as Seafloor OBS Source Mirror OBS Sea surface Figure 4: Schematic diagram of the mirror principle. For migration, multiple waves can be regarded as primary waves supposing that the OBS station is deployed on one surface above the seafloor, whose distance to the seafloor is twice the water column.

Geofluids
where matrix Λ includes the eigenvalues of L, Q stands for its eigenvectors, and superscript T denotes transposition.
However, when performing the wavefield depth extrapolation, one unavoidable question is the evanescent waves. Sandberg and Beylkin [22] analyzed that the evanescent waves are produced because of negative values of the L operator. Based on this understanding, we introduce a matrix decomposition method to remove the negative values of the L operator so as to suppress the evanescent waves; more details can be found in Ref. [42].
For imaging multiple waves by using the OBS data, we prefer to utilize the mirror principle; the schematic diagram of the mirror principle is drawn in Figure 4. OBS is usually located at the seafloor, and its position can be settled by the direct wave. Assuming the sea surface is flat, therefore we can find a mirrored OBS which has the same distance to surface as the actual OBS. Viewing Figure 4, we can find that the mirrored OBS shares the same traveling time with the actual OBS for multiple waves. In summary, we can use the mirrored OBS data to image multiple waves.

Results and Discussion
In the numerical experiment, we designed a dipping velocity model with three layers ( Figure 5); the velocity of which was 1500 m/s, 2000 m/s, and 2500 m/s, respectively, to test the imaging performance by means of primary and multiple waves. The sources were put within the range of 0~2000 m at the sea level, and the OBS was put on the seafloor at the middle of model. The wavefield was calculated by a fullwave equation wavefield simulation method, as shown in Figure 6. Viewing the wavefield propagation, we can clearly observe the primary waves, because the longer traveling paths of multiple waves, the weaker the amplitudes of multiple waves than that of primary waves. Then, the primary waves were used to perform depth migration and the mirror principle was utilized to perform the imaging of multiple waves. Figure 7 shows the imaging results of primary and multiple waves. It can be clearly seen that the primary waves cannot image some parts of the interface such as the boundary of the model and the seafloor, while the multiple waves allow imaging the boundary of the model successfully and provide a boarder imaging illumination, especially for the seafloor reflector. This is because multiple waves penetrate into the subsurface several times and have longer traveling 6 Geofluids paths as well as contain abundant reflection information ( Figure 8). This simple experiment testifies the advantages of multiple waves in imaging over primary waves. As discussed before, as a conventional migration method, reverse time migration has been widely used in the subsurface imaging. In order to testify the imaging performances of the two-way wave equation depth migration method, the reverse time migration was also applied to image OBS multiples. By comparing the imaging results produced by our proposed migration and reverse time migration methods, we can see that their imaging results are very similar (Figure 7). But the RTM method needs about 6.4 Gb memory to store wavefields while our proposed migration method only needs about 46 Mb memory. That is an outstanding advantage of the depth migration method but can achieve a similar imaging result.
For the real OBS data, we used the mirror principle to image the multiple waves based on the two-way wave equation depth migration method. But for a fair comparison, the imaging result of primary waves was also included. The velocity field used to calculate travel times was determined by the travel time inversion of refractions and reflections from the two OBS data in this area [36]. Considering the quality of the OBS data, only one OBS data along the multichannel seismic line was used for imaging. The width of illumination is very limited by using the primary waves. However, when we observe the propagation of multiple waves, we can see that the multiple waves have a longer traveling path than the primary waves. Theoretically, multiple waves bring more information about submarines.
We zoom in the imaging zone under OBS for a detailed comparison, as shown in Figure 9. Observing the imaging results of primary and multiple waves, it is clear to see that the imaging quality of multiple waves is improved dramatically in contrast to the imaging result of primary waves. The most apparent place is the imaging result of the seafloor and the shallow reflections just below the seafloor, which are not imaged by primary waves but well imaged by multiple waves. Moreover, multiple waves produce a much better image and lateral illumination is much enhanced. Due to the long traveling time of multiple waves, it can illumine a wider region than primary waves. The shallow structures under the seafloor are imaged clearly. This practical application of OBS data also proves the advantages of multiple waves in imaging over the primary waves.

Conclusions
Multiples are common in marine seismic data, and they are often discarded as noises. However, they might carry 7 Geofluids abundant information on subsurface and provide additional information on the subsurface parts not illuminated by primaries. In this study, the gas hydrate system along the South Shetland margin (Antarctica) was imaged by using multiple waves of OBS data and the imaging result was compared with common imaging using primaries. Because of its efficiency, a two-way wave equation depth wavefield extrapolation method was used to obtain the migration imaging, which is the first application on the OBS data. The imaging results show that the migration imaging with multiples provides a wider illumination and a clearer structure of the subsurface than traditional imaging with primary waves, particularly for the shallow near-seafloor reflected events, which can provide a reference for the study and exploration of gas hydrates. Moreover, this case study proves the power of the two-way wave equation depth wavefield extrapolation method in imaging OBS multiples.

Data Availability
The data are available from the first author on request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.