Depth migration of crustal-scale seismic reflection profiles: A case study in the Bohai Bay basin

The application of the pre-stack depth migration (PSDM) method using deep seismic reflection data is still a challenge. In this paper, a deep seismic reflection dataset collected from the Bohai Bay Basin is used to test the applicability of the pre-stack depth migration method. After improving the signal-to-noise ratio of the prestack gathers and building an accurate velocity model capturing lateral and vertical variations, a high-quality pre-stack depth migration profile is obtained. This profile shows detailed crustal structures and an undulating Moho. The depth difference between the east and west sides of the Moho imaged in the profile is about 6 km. The Moho is disrupted by several discontinuities. Compared with conventionally stacked time-migration profiles, the pre-stack depth migration profile provides a higher-resolution subsurface image with accurate depth information. Combining geological and geophysical data, the pre-stack depth migration profile reveals geologic structures responsible for crust thinning during the destruction of the North China Craton. The profile also provides more accurate depth domain information for further study of crust-mantle interaction in this area.


Introduction
Crustal-scale seismic reflection imaging is effective in providing high-resolution information about the continental lithosphere, which has been proven successful all over the world (Clowes et al., 1968;Wang et al., 2010). The evolution of the lithosphere can be established by determining subsurface structures obtained by deep seismic reflection profiles. Since the 1970s, seismic reflection projects for continental studies have been carried out across the globe. These projects help to understand crustal deformation and continental crustal evolution. However, the depth information in the current interpretation of deep seismic reflection profiles is usually obtained by converting two-way traveltime using the average crustal velocity (Potter et al., 1986;Brown, 1991;Cook, 1995;Cook et al., 1999), or by time migration results (Calvert and Doublier, 2018) to the depth. This inevitably leads to geometrical distortion of the actual geological structures. Meanwhile, the application of depth migration for producing deep seismic reflection profiles is still a challenge. Specially, the depth migration profiles often suffer from severe arc reflection when using a reasonable crustal velocity to implement depth migration. This has a serious negative impact on the extraction of effective information from the depth migration profile. Warner (1987) even asserts that depth migration is not applicable to deep seismic reflection profiles. Apparently, obtaining a deep seismic reflection profile with reliable depth information is not an easy task.
Post-stack depth migration has been applied to the study of oceanic lithosphere structures, which have simple subsurface structures and smaller lateral velocity variation compared to that for the continental crust (e.g., Kopp et al., 2009). However, the continental lithosphere is difficult to image because of the complex subsurface structures such as undulating topographic surface, more complex fault-system geometries, and heterogeneous seismic wave velocities. Post-depth migration has also been attempted to produce continental-lithosphere seismic reflection profiles. The COCORP project used a relatively reasonable crustal velocity structure: 3-5 km/s in the sedimentary layer, 6 km/s in the basement, and 7 km/s in the lower crust, and conducted a post-stack depth migration test with distorted results (Lynn et al., 1983). In order to overcome the distortion problem and obtain a suitable profile, Warner (Warner, 1987) points out that the optimal velocity for migration is about 50% of the reflection or refraction velocity, which is obviously not in line with the actual underground conditions. In the early part of this century, post-stack depth migration of deep seismic reflection profiles was successfully applied to works done in Alberta of Canada, where a gradient velocity model (surface velocity value from 5.5 km/s to 7.1 km/s in the lower crust) was used to Frontiers in Earth Science frontiersin.org 02 establish the varying crustal thickness (Christensen and Mooney, 1995;Bouzidi et al., 2002).
In this century, pre-stack depth migration (PSDM) has become the main technology to image the high-resolution subsurface structure. Obviously, it is an important step to test and promote the application of the PSDM method in the continental deep seismic reflection profile to adapt to the development of technology. It can reveal the depth domain structure of the continental crust. In this paper, high-quality deep seismic reflection data (Li et al., 2022) from the Bohai Bay Basin in eastern China is selected for research experiments. The Kirchhoff pre-stack depth migration method is used to produce the depth migration profile. By comparing the results of PSDM and pre-stack time migration (PSTM), we discuss the feasibility of applying PSDM in deep seismic reflection and its significance for deep crustal structure interpretation.

Geological setting
The seismic reflection profile is located in the central part of the Jizhong Depression in the Bohai Bay Basin. The basin contains Mesozoic-Cenozoic strata on top of the North China basement rocks . The depressions can be further divided into the Baoding Sag, the Rongcheng Uplift, the Niutuozhen Uplift, and the Southern Niubei Slope ( Figure 1B). The depression is located on the eastern continental block of the North China Craton (NCC) (Zhao et al., 2005). The crystalline basement of theNCC was formed in the Archean-Paleoproterozoic (Kusky and Li, 2003;Zhai, 2010). In the Late Mesozoic, the closure of the Mongol-Okhotsk Ocean and the subduction of the Paleo-Pacific plate caused a significant change in the structure of the NCC, from a nearly north-south contraction to a NWW-SEE extension (Zhu et al., 2011;Zhu et al., 2012). From late Mesozoic to Cenozoic, due to extensive magmatic activity and largescale tectonic extension, the ancient, cold and thick cratonic lithosphere changed into the thin and hot lithosphere (Chen et al., 2010). At the same time, the Bohai Bay Basin began to develop . Under the control of the Taihang Mountain Piedmont Faults and a series of secondary faults, intense extensional activity occurred and lasted until the end of the Paleogene. During the Neogene-Quaternary, the back-arc extension in North China weakened, the Jizhong Depression entered the depression evolution stage, dominated by overall subsidence, and the entire Bohai Bay Basin entered the post-rift thermal subsidence stage (Ma et al., 1983;Qi et al., 2003).
The seismic reflection line passes through the Xushui Depression, Rongcheng Uplift, Niubei Slope and Niutuozhen Uplift, and ends in the Baxian Depression from west to east. The surface of the study area is covered by a Quaternary mudstone layer with a variable thickness. In order to obtain lithospheric structures at a high resolution with a high signalto-noise ratio, our seismic surveying was conducted using a shot spacing of 40 m, a receiver spacing of 20 m, and 1000 receiving channels. A total of 1850 high-quality single records were obtained by means of large, small, and medium charges combined with explosive source excitation, with a full coverage length of 60 km ( Figure 1). More details on the data acquisition can be found at Table 1 or Li et al. (2022).

Data and methods
The testing process of the PSDM method includes data preprocessing, velocity-model building, and the application of the PSDM method. Figure 2 illustrates the processing workflow: the left part is the processing flow of pre-processing to generate effective CMP gathers, whereas the right part is the processing employing the PSDM method.

Pre-processing for pre-stack depth migration
Pre-processing for PSDM includes tomographic static correction, pre-stack multi-domain denoise, amplitude compensation and surface consistency deconvolution. Tomographic static correction can resolve the distortion problem caused by the low-velocity zone and the undulating surface topography. Abnormal noises and surface waves are attenuated by multi-channel identification, single-channel suppression, and adaptive attenuation methods. The records before and after noise attenuation are shown in Figure 3. Surface-consistent amplitude compensation and trace-to-trace amplitude normalization were used to reduce lateral amplitude variation. Surface-consistent deconvolution is used to improve the wavelet consistency of the wavelet and broaden the frequency band. The detailed parameters of data pre-processing is listed in Table 2. Frontiers in Earth Science frontiersin.org 03

FIGURE 2
The deep seismic reflection data processing flow. The left part represents the pre-processing scheme for PSDM. The right part indicates the PSDM processing.

FIGURE 3
The comparison before and after noise attenuation. (A) The record before noise attenuation with time function gain. (B) The record after noise attenuation with the same time function gain.

Frontiers in Earth Science
frontiersin.org 04

Velocity building
The velocity model directly affects the quality of the PSDM method. The migration process requires a rational velocity model corresponding to the underground location. The standard industrial migration process employs the interval velocity converted from the stacking velocity with normal move-out correction. It is effective for the shallow-crust data within reasonable parameters. In this paper, the initial depth domain layer velocity model is mainly divided into shallow part (0-7 km) and deep part (7-50 km). The shallow velocity model is a fusion of velocity model obtained by firstarrival wave tomography and interval velocity model obtained using constrained velocity inversion. The velocity of the first arrival wave tomography inversion is only about 2 km, for reference only. In the depth where the first arrival wave cannot be retrieved, the velocity in the depth domain is transformed by the velocity inversion method of the constrained layer. And the deep velocity model is from wide-angle reflection and refraction, which can break through the limitations on reflected deep velocity analysis and thus optimize the imaging effect (Li, 2013). Several deep seismic sounding profiles in the study area intersect the deep seismic reflection profiles ( Figure 1B) so the deep part of initial velocity model is the velocity from the wide-angle reflection and refraction (Duan et al., 2016). The interval velocity model is shown in Figure 4A.
After the first application of the PSDM method using the initial velocity, the velocity error was estimated according to the depth deviation of the CRP gather. The Residual Moveout (RMO) and grid   Figure 4B; Figure 4C shows that in addition to the vertical velocity variation, the final velocity model can show a sharp lateral velocity difference (more obvious around 8 km). The blue solid line with a sharp increase at a depth of 8 km represents the velocity of CDP 1500 which is modified by tomographic inversion. The velocity deeper than 8 km in the western profile (west of CDP 2700) is significantly higher than the eastern profile. The velocity value decreases from the location at CDP 2000.

Pre-stack depth migration
The Kirchhoff pre-stack depth migration is an efficient and practical PSDM method with a high computational efficiency, high accuracy and good computational stability. It has been widely used in the petroleum industry (Yang et al., 1996). The equation for producing imaging results can be expressed as follows: where R(r) is the reflection coefficient; r, r s and r g represent the position of the imaging point, source, and receiver, respectively; t(r, r s ) and t(r, r g ) represent the travel time from r to r s and r g , respectively; A(r, r s ) and A(r, r g ) stand for the amplitude of r to r s and r g , respectively; z is the depth, and p is the pressure wave field. This method first calculates the travel time of seismic waves and then realizes the migration. The calculation of travel time is based on the dynamic ray tracing algorithm on wavefront construction. And the algorithm selects and saves the propagation time of the event with the largest amplitude and offset to that point depth. Theoretically, the Kirchhoff pre-stack depth migration can achieve true amplitude imaging.

Data reliability analysis
In order to analyze the authenticity of the migration results and verify the reliability of the lateral variation of velocity, this study conducted reliability analysis using the first-arrival tomographic inversion and analysis of shot records. The first arrival wave of the original unprocessed shot records carries the direct wave, refracted wave, or wave combination that reaches the surface first. The velocity structure can be inverted by using the first arrival wave ray travel time and path. Figure 5 is the velocity model obtained from the inversion of the first arrival, which shows a sharp lateral velocity change at a shallow depth of 1.8 km. We can see that at a depth of 2 km or less, the layer velocity value reaches 7 km/s. From the first arrival waves (direct arrival waves, refraction waves) recorded by shots (Figure 6), it can be found that the crustal velocity structure from west to east has huge heterogeneity and systematic changes. The velocity of the first arrival wave of a single shot varies drastically, and there are even as high as 7 km/s or more. The shot records and first-arrival inversion velocity both show that the velocity increases sharply from the middle crust on the west side, and the vertical velocity increases slowly in the east. This results in a more stable moho reflection time in the time profile between the thicker western part and the thinner eastern part. These comparisons demonstrate the reliability and stability of prestack depth migration imaging structures. Although the absolute values of layer velocities corresponding to different analysis methods are slightly different, both analyses show that lateral velocity variation is always visible.

Subsurface structure revealed by the PSDM profile
The final PSDM profile is shown in Figure 7A. In the profile, the sedimentary layers are present in the top 10 km, and the basement morphology is continuous. The middle and lower crust is dominated by strong reflections and weak refections. There are longitudinal alternating strong and weak reflected domains. The reflection of the crust-mantle transition zone has a high signal-to-noise ratio. The

Structure of the upper crust
The tectonic pattern of two depressions and two uplifts in the upper crust is consistent with the regional geological data. Based on the drilling data (He et al., 2018;Li et al., 2022), it is inferred that the Quaternary and Neogene sedimentary strata are within 2 km deep, displaying multiple groups of near-horizontal reflection strata with a thickness of 200-2500 m. In the depression, there are three groups of half-graben structures in the basement controlled by the hidden faults F1, F2, and F3 as the western boundary. Where F1 is the Baoding-Shijiazhuang fault, F2 is the Rongcheng fault and F3 is the Niudong fault. There are also several secondary faults of these three basement-involved faults that can be clearly observed in the locally enlarged profile shown in Figure 8. Figure 9 shows the comparison between the PSDM profile and PSTM profile in the Baxian fault-bounded depression. Both profiles reveal a buried hill internal which is located at the depth of 4 km in the PSDM profile and at the two-way time of 3 s in the PSTM profile. The PSDM profile reveals the buried hill internals structure more clearly with better reflection continuity and a higher signal-to-noise ratio.

Structure of crust-mantle transition zone
The crust-mantle transition zone is defined by the lower-crust bottom reflectors and a strong reflective Moho interface. The characteristics in PSDM profile are as follows.

1) Generally continuous seismic reflections of crust-mantle
transition are well developed in the whole profile except the section from CDP 4200 to CDP 5000. To the west of CDP 4200, the crust-mantle transition zone is composed of dense, strong reflections with a downward convex shape. These reflection events are of a certain thickness and has a longitudinal sustained depth of about 3-4 km in the profile. The Moho depth of CDP 2500 and CDP 4200-5000 has changed. To the east of CDP 5000, the crust-mantle reflection zone is spread horizontally, and the thickness is about 1.5-2.5 km. Compared with the western area, the reflection amplitude of the eastern crust-mantle transition zone is weaker and the thickness is thinner. 2) As a whole, the crust is thicker in the west and thinner in the east, with the section between CDP 4200 and CDP 5000 as the transition zone.  Frontiers in Earth Science frontiersin.org 2000-3000 and CDP 4200-5000, which are low reflectivity areas. These transparent reflections of the lower crust correspond to deep Moho discontinuities. It is worth noting that east of CDP 4900, the lower crust has an anticlinal reflection (marked in the blue dashed box).

Factors of PSDM
The predominant challenge to applying the PSDM method is to get a reasonable depth-velocity structure. Compared to shallow-crustal seismic reflection profiles used in oil exploration, which can be validated using drill-hole data, lithospheric scale reflection profiles focusing on structures in the lower crust and the upper mantle and having lower signal-to-noise ratios do not have this advantage. The non-single-valued error of velocity-depth increases linearly with depth in shallow layers and squarely with depth in deep layers (Ross, 1994;Bloor et al., 1997). This means that the error between velocity and depth is significant when depth exceeds the maximum offset. To improve the accuracy of seismic reflection velocity is to increase the offset, and the offset of wide-angle reflection and refraction can reach tens to hundreds of kilometers. Therefore, it is of great significance to use wide-angle reflection and refraction to constrain deep seismic reflection profiles for depth migration processing. The necessary condition for accurate reflection wave inversion is that there is adequate effective reflection in the lower crust so that the reliability and rationality of the velocity and depth values can be judged according to the RMO semblance of the CRP gathers. Figure 10 shows the final velocity model and the shear wave velocity structure of the crust and upper mantle through receiving function waveform inversion. The comparison shows that the velocity model in this paper not only has the same variation trend, but also is more refined.
Another factor for the reasonable deep velocity model obtained by this test inversion is that there is an obvious wave impedance Frontiers in Earth Science frontiersin.org difference in the middle and lower crust, which is consistent with the requirements of the tomographic inversion theory. It also shows that velocity modeling of the lithosphere is a process of combining the understanding of local geology with adequate seismological methods.
Migration aperture is one of the important parameters affecting the PSDM of deep seismic reflection. In this study area, if the aperture is less than 15 km, it will limit the summation of the diffraction hyperbolas, making it impossible to image the lower crustal steep interface. If the aperture is greater than 19 km, a "smiley-face" section which is antiformal migration artifacts appears in the boundary of the transparent area and high reflectivity area below the lower crust. PSDM is a process of stacking all signals. The aperture can be appropriately reduced to avoid wasting calculation time and reducing the signal-to-noise ratio of imaging when the structural fluctuation is not large. However, the aperture should be larger than the length of the receiving length containing effective reflection at least. Or the best is the maximum horizontal displacement of the phase axis in the migration with the steepest angle of the profile (Yilmaz, 2001). In practical application,  Figure 7).

FIGURE 9
Local enlarged profile comparison of the Baxian Sag of the PSDM profile (A) and the PSTM profile (B). (The enlarged range is marked in Figure 7).

Frontiers in Earth Science
frontiersin.org the selection of aperture size needs to be combined with the depth of the target layer and the degree of inclination. The aperture should be increased to cover the extensive range of high and steep inclines when the structure fluctuates greatly. At the same time, a too-large aperture will reduce the resolution of imaging results (Ma, 2002), so different apertures should be selected according to the complexity of the structure to obtain the best imaging.

Artifacts of moho horizontal features
The Moho in the stack profile and PSTM profile is nearly horizontal. It is placed at approximately 11s in two-way time. However, we can see that the Moho in the west is about 6 km deeper than in the east. The boundary is positioned from CDP 4200 to CDP 4500, corresponding to an obvious weak reflection of the lower crust. The updated velocity model shows that the velocity of the lower crust is higher in the west and lower in the east. As can be seen from the original records, the apparent velocity of the refracted wave at the west is about 7 km/s at about 3s in twoway time. Although the apparent velocity cannot be directly regarded as the propagation velocity, it can qualitatively reflect the true velocity. The shallow velocity structure obtained from the first-arrival wave inversion verifies the existence of a highvelocity layer in this area. This is the main factor leading to the difference in Moho morphology between the time-domain section and depth-domain section.

Structural implications of the PSDM profile
The comparison of the PSDM profile and the stack profile of the same survey line is shown in Figures 11A,B. The sedimentary layer

FIGURE 10
The velocity profile in the study area. (A) Shear wave velocity obtained by receiver function waveform inversion of the corresponding area (Zheng et al., 2009) (B) Final velocity model used by the PSDM method.

FIGURE 11
The comparison of the PSDM profile and stack profile. (A) PSDM profile (B) Stack profile of the same survey line (Li et al., 2022).

Frontiers in Earth Science
frontiersin.org in the time domain profile accounts for a large proportion due to the slow propagation speed of the sedimentary layer. And the velocity of the middle and lower crust is higher, so the middle and lower crust occupy a higher proportion in the PSDM profile. From the perspective of the overall lithosphere structure, the Moho offset of the PSDM profile is more obvious than that of the stack profile, and the crust on the west side of the survey line is thicker than that on the east side. Crustal thickness is an important basis for judging the evolution history of the crust, and thinning and thickening generally indicate crustal extension or compression respectively (Meissner and Mooney, 1998). Similar features are also obtained from the dense array receiver function study (Zheng et al., 2009), further confirming the reliability of the pre-stack depth migration profile. The Moho discontinuities and offset are more clearly revealed by the high-resolution PSDM profile which corresponds to weak reflection areas in the lower crust and may represent channels for new mantle material to flow into the crust after the destruction of the thick North China lithosphere. It also shows that the formation time of the Moho predates the mantle upwelling event that did not create a new Moho. On the west side of the survey line, the reflection is short and horizontal in the middle and lower crust. On the east side, it is an anticlinal deformation pattern. The obvious differences in the structure on both sides may prove that the western lower crust once belonged to hard and ancient block, which was closer to the characteristics of the central orogen of the North China Craton (Zhao et al., 2005). It is probably the deep boundary of the central orogeny and eastern block of NCC in the lithosphere scale.

Conclusion
The deep seismic reflection profile produced by the PSDM method accurately depicts the crustal structures in the Bohai Bay Basin; sediments thicken from west to east controlled by three major faults. Multiple groups of secondary faults and fault breakpoints are shown on the PSDM profile. The velocity varies laterally, higher in the west and lower in the east. The crustal thickness defines the basic structural feature of being thick in the west and thin in the east. The Moho displays several discontinuities along the profile. On both sides of the CDP 4200-CDP 5000 section in the profile, the Moho offset is up to 6 km. These discontinuities correspond to weakly reflective domains in the lower crust and may represent channels for new mantle material to flow into the crust during or after the destruction of the earlier, thicker North China lithosphere.
There are differences between the PSDM results and timedomain imaging results. The PSDM method can provide more reliable imaging results than time-domain sections. The PSDM method can not only resolve the problem of fault imaging but also make the contact relationship of shallow faults more reliable. At the same time, it can better locate lower crustal complex structures and provide the depth position information of reflection interfaces. The Moho discontinuities are clear, depicting the thinning of the crust due to the interaction between the crust and the mantle. For the imaging results of the lower crust and lithospheric mantle, other geological and geophysical data in this area confirm the rationality of the pre-stack depth migration results in this area.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.