Estimating the Shear-Wave Velocities of Shallow Sediments in the Yellow Sea Using Ocean-Bottom-Seismometer Multicomponent Scholte-Wave Data

Scholte-wave dispersion analysis is effective at imaging the relatively low shear-wave velocity of shallow marine sediments in marginal seas. The combination of a four-component ocean-bottom-seismometer (OBS) and a towed air-gun source can economically and effectively acquire the marine dispersive seismic data. Extracting higher-order dispersive Scholte wave modes is the most critical problem in the dispersion analysis method. The extremely low shear-wave velocity and severe attenuation in the top hundreds of meters of marginal sea sediment provide an uneven dispersive energy distribution for the four components of the Scholte wave data. The fundamental mode dispersive energy dominates in the vertical component and higher-order modes dominate in the horizontal component. We developed the method of the four-component OBS Scholte velocity-spectra stacking, which can effectively, rapidly, and robustly extract higher-order modes. We imaged the shear-wave velocity structure of complicated shallow marine sediment in the North Yellow Sea using an active OBS seismic profile with a large-volume air-gun array. The fourth higher-order Scholte wave mode can be imaged with the four-component velocity-spectra stacking method with a lower frequency range of 1.0–7.0 Hz. Only the second-order mode can be recognized from the dispersion energy image of the single vertical component. The joint inversion of multimode dispersion curves can provide more accuracy and deeper constraints for the inverted model; thus, the constraint depth with five modes increases by a factor of 1.9 compared with single fundamental mode inversion. The inverted profile suggests a low shear-wave velocity of 123–670 m/s and strong lateral variations within 350 m. The main regional geological structures are shown by the inverted shear-wave velocity structure.


INTRODUCTION
Imaging suboceanic shear-wave velocity (V s ) structures is fundamental in the geophysical investigations of marginal seas (Ewing et al., 1992;Klein et al., 2005;Kugler et al., 2007;Wang et al., 2016;Wang et al., 2020). The V s for shallow marine sediment is treated as an important parameter in many marine engineering construction activities, such as offshore platforms, wind parks, and pipelines. Moreover, the V s structure of shallow marine sediment provides a reliable reference model for offshore multicomponent seismic exploration (Bohlen et al., 2004;Kugler et al., 2007), ocean lithosphere shear-wave tomography (Zhao et al., 2008;Zhao et al., 2010), underwater acoustic wave attenuation (Rauch, 1980;Hughes et al., 1990;Bradshaw, 2015), and very-lowfrequency geoacoustic studies (Du et al., 2020).
The dispersive characteristics of interfacial and surface waves are widely used to estimate the V s structure of shallow marine sediment (Bohlen et al., 2004). Compared with the marine converted-body-wave survey approach, marine interfacial waves show greater promise due to the low signal-to-noise ratio of ocean-bottom body shear waves from the weak P-to-S wave conversion coefficient (Kugler et al., 2007). The dispersion curves for interfacial waves are sensitive to shear-wave velocities, which is the precondition and fundamental advantage of using interfacial waves to invert S-wave structures.
The vertical or radial component of an interfacial wave that travels on the water-sediment interface in the marine environment is known as a Scholte wave and is analogous to Rayleigh waves on land (Scholte, 1947). Its transverse component is known as a Love wave (Love, 1911), both in this case and on land. Scholte waves can be detected through a combination of a single deployed underwater seismic station and multiple air-gun shots on the sea surface, which is known as the "stationaryreceiver" method . An ocean-bottom seismometer (OBS) (also known as an ocean-bottom hydrophone or ocean-bottom node) is generally used as the ocean-bottom seismic station as it is cheaper and technically simpler for marine fieldwork. However, acquiring marine Lovewave data is more complex because of the need for an oceanbottom P-SH type of source (e.g., shear-wave vibrator; Socco et al., 2011).
The single vertical component of OBS active source Scholtewave data was used to investigate the V s structure of shallow marine sediments (Shtivelman, 2003;Shtivelman, 2004;Bohlen et al., 2004;Klein et al., 2005;Kugler et al., 2005;Kugler et al., 2007;Nguyen et al., 2009;Socco et al., 2011;Vanneste et al., 2011;Dong et al., 2012;Wang et al., 2016). Wang et al. (2020) noted that the three-component Scholte-wave dispersive analysis can extract additional higher-order modes compared with using a single vertical component. Compared with the fundamentalmode inversion, joint higher-order modes can effectively increase the inversion depth and resolution (Xia et al., 2003;Luo et al., 2007). Wang et al. (2020) demonstrated that the joint inversion of the four-mode dispersion curves can decrease the maximum inverted error by a factor of 16 compared with one fundamental mode inversion. However, the possibility of and method to image higher-order dispersive modes from the fourcomponent (4C with three seismic components and one hydrophone component) OBS Scholte-wave data remain unknown and are worth exploring.
This study considers the possibility and potential of using the 4C OBS Scholte-wave data to image higher-order dispersive modes. Moreover, the feasibility of imaging the complex shallow marine sediment of marginal seas using the Scholtewave dispersion inversion method is verified. We begin by presenting the acquisition and geometry of the OBS activesource Scholte-wave field data from the Yellow Sea and then describe the detailed 4C Scholte-wave analysis procedure from the raw OBS data to the final V s structure. We next show the dispersion energy images obtained from the 4C Scholte-wave analysis method and finally illustrate the inverted onedimensional (1D) and pseudo-two-dimensional (pseudo-2D) V s structures.

Seismic Data Acquisition
As shown in Figure 1, a north-northwest-oriented OBS activesource survey profile with a total length of 342.2 km was acquired in the Yellow Sea of China in 2013. The geological location of the survey area is the Yangtze-Gyeonggi plate (Li et al., 2012). The survey profile crosses the continental shelf and northern basin of the South Yellow Sea and ends at the north rim of its southern basin. A CAS-Micro-4C-OBS was used as the marine seismic station, which was developed by the Institute of Geology and Geophysics of the Chinese Academy of Sciences. The CAS-Micro- Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 812744 4C-OBS assemble had a 13-inch glass sphere (Vitrovex) with a maximum deployment depth of 6,000 m. It can record data for a period of 3 months with a maximum recovery period of up to 6 months. Three geophones with a frequency bandwidth of 0.1-150 Hz were installed orthogonally in the glass sphere. A low frequency (1 Hz-15 kHz) hydrophone was set up in the OBS and placed outside the glass sphere. The dynamic range of the OBS was greater than or equal to 120 dB. A 24-bit A/D converter was used with a sampling frequency freely set from 1 to 250 Hz. A 32 GB FAT32 flashcard was used to store the data. In total, 16 OBSs with an average deployment distance of~12 km were deployed along the profile and were designated as C02-C32 with an interval of 2. However, as C02 and C04 were relatively close to the coastline where there was significant additional disturbing noise, they are excluded from this study due to low data quality; thus, only 14 OBSs were used for the Scholte-wave analysis. The air-gun shots each had a volume of 6,060 in 3 , an average interval of 136 m, and were approximately 10 m below the sea surface. A differential global positioning system (DGPS) was used to control the air-gun shot timing and vessel navigation.
Combined with the shot information, the common receiver gather (CRG) was cut from the raw OBS continuous binaryformat data. As an example, the 4C CRGs of the OBS from station number C22 are shown in Figure 2. After obtaining the traces by subtracting their means, a 0.5-6.3 Hz bandpass filter was applied.
The pre-processed seismogram shows a body wave (green arrows in Figure 2) with a strong reflection and refraction and a visible Scholte wave (black arrows in Figure 2) for both the seismic components (Figures 2A-C) and hydrophone component ( Figure 2D). The raw CRGs also show that the Scholte-wave energy is distributed unevenly between the seismic and hydrophone components. The Scholte-wave energy from the seismic components is stronger than that of the hydrophone component. This may be because the amplitude of the Scholtewave decays exponentially with distance from the interface, and the geophones of the OBS instrument are directly coupled with the seafloor while the hydrophone floats in the water at 0.8 m from the seabed.
The three seismic CRGs show three primary Scholte-wave features. First, both the seismic components and hydrophone component record the Scholte wave. Second, the seismic components with different directions exhibit various energy amplitudes. A dominant Scholte wave is visible in the vertical BHZ component ( Figure 2C), whereas a relatively weak Scholte wave is visible in the random-azimuth horizontal BH1 component ( Figure 2A) and nearly disappears in the other horizontal BH2 component ( Figure 2B). Third, the Scholte wave modes recorded from the different seismic components are distinct, with the BH1 component recording a higher-order mode ( Figure 2A) and the BHZ component recording the dominant fundamental mode ( Figure 2D).

4C Scholte-Wave Analysis Procedure
The OBS 4C Scholte-wave dispersion energy imaging (DEI) and inversion method for each OBS station involved the following four main steps. FIGURE 2 | Pre-processed common receiver gathers (CRGs) recorded at OBS station C22. Panels (A, B) show two horizontal seismic components with deployed random azimuths. Panels (C, D) show the CRGs of the vertical and hydrophone components, respectively. The red and blue dashed lines denote apparent velocities of 100 and 160 m/s, respectively. The pre-processing involved trace normalization using the mean with bandpass filtering at 0.5-6.3 Hz, and finally Wiener (least-squares) predictive error filtering. The green and black arrows denote the body and Scholte waves, respectively.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 812744 Step 1: Preparation of 4C CRGs. Transfer the raw continuous binary OBS data to the 4C CRGs with the SEGY data format, which is followed by pre-processing the raw CRGs with trace normalization, bandpass filtering, Wiener predictive error filtering, and cutting the traces with an absolute offset of 400-5,000 m.
Step 2: 4C DEI. Calculate the velocity spectrum (VS) of each component using a phase-shift algorithm.
Step 3: 4C velocity-spectra stacking. Obtain the final stacked VS by stacking the VSs from the three seismic components with the hydrophone component.
Step 4: Extraction and inversion of the dispersion curve. Pick the multimode dispersion curves from the stacked VS and finally invert the structure for the shallow marine sediment.

Preparation of 4C CRGs
As an example, we show the 4C Scholte-wave dispersion analysis method for the OBS station C22. We begin by pre-processing the 4C CRGs following step 1, as shown in Figure 2. The entire CRG is split into two branches: positive and negative offsets. For each CRG branch, the minimum and maximum offset windows are chosen to cut the traces, which guarantees a high signal-to-noise ratio for use in the next step. Here, the near-and far-field effects of the CRG are the primary references used to select the offset window (from 400 to 5,000 m).

4C DEI
The pre-processed CRG branch is then transferred from the t-x domain to the f-v domain using the phase-shift algorithm (Park et al., 1998) following step 2. For one trace of the CRG data r i (i = 1, 2, . . . ,N) in the t-x domain, its Fourier transform R i in the frequency domain is obtained by: Where ω is the angular frequency, and A i and P i denote the amplitude and phase terms, respectively. After removing the A i via normalization: We define a trial phase-velocity (c ω ) range with small increments, and the velocity spectrum (VS) in the frequencyphase-velocity domain is calculated by: Where x i denotes the offset from the source to the receiver. As shown in Figure 3, the 4C VSs are imaged using the phaseshift algorithm for the positive CRG branch of the OBS station C22. First, the 4C dispersion energy images show that both the seismic and hydrophone components can be used to image the dispersive VSs. Second, both the fundamental mode and higherorder modes are obtained simultaneously. Third, the mode energy is distributed unequally between the different components. The fundamental mode is dominant in the VS of the vertical component ( Figure 3C) but missing in those for both horizontal components ( Figures 3A,B). Correspondingly, the first higher-order mode is dominant in the two horizontal components. Finally, the recognizable highest dispersive Scholte-wave mode differs between the four components. Thus, the fourth higher-order mode can be distinguished from the BH1 component, the weak third higher-order mode is present in the BHZ component, and the first higher-order mode can be distinguished from the hydrophone component.
There is a significant difference between the VSs of the two horizontal components; only the first higher-order mode can be recognized from the VS of the BH2 component, whereas the VS of the BH1 component clearly shows the first four higher-order modes. Figure 4 shows hodograms of the BH1-BH2 ( Figure 4A) and BH2-BHZ ( Figure 4C) components plotted using a time window of 1.7-2.2 s for a seismic trace with an offset of 1,499 m for the OBS station C22. In the BH1-BH2 hodogram, the particle vibration trajectory is characterized as a straight line along the BH1 axis, while the BH2-BHZ particle motion trajectory is a vertical straight line. These two main features of the hodograms indicate that the BH1 component corresponds almost directly to the shot-line azimuth, whereas the BH2 component is perpendicular to the shot line. For this OBS station, BH1 can be treated as the inline component, with BH2 treated as the corresponding crossline component. This explains reasonably well why there is only one relatively weak dispersive energy mode (M1 in Figure 3B) in the VS of BH2.

4C VS Stacking
We stack the 4C VSs together after completing the above DEI as described in step 3. The 4C stacked VSs (VS stack ) are obtained by Where VS BH1 , VS BH2 , VS BHZ , and VS HYD denote the imaged VSs of the BH1, BH2, BHZ, and HYD components, respectively. Here, the precondition for the 4C VSs being directly linearly stacked is that they contain the same dispersion energy as the Scholte-wave source. As the source-to-receiver seismic wavefield propagation path can be treated as a pseudo-2D profile, only the P-SV type of Scholte waves exist along this inline profile. No P-SH-type Love waves exist as the P-wave-type acoustic waves are excited from the air-gun sources. This can also be demonstrated inversely from the concordance comparison of the 4C dispersion curves. If the VSs contain Love waves, then the dispersion curves picked from the vertical component differ from those from the horizontal component. Conversely, consistent dispersion curves from the two components suggest that only the same type of Scholte waves exist and that there are no Love waves. Therefore, we regard the two horizontal components as projections of the inline component and a direct linear stacking of their spectra without azimuthal correction.
As an example, the dispersion curves picked from the 4C VSs as shown in Figure 5 are compared in detail. As shown in Figure 5A, the dispersion curves show great consistency between the different components, from the fundamental to the fourth higher-order dispersive modes. The mean of the cross-correlation coefficients between the 4C first-higher-Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 812744 order-mode dispersion curves approaches 99.88%. The correlation coefficient of the first-higher-order-mode dispersion curves produced by the BH1 and BHZ components is the highest and reaches 99.98%. The BH1 and BH2 components have the lowest correlation coefficient of 99.78%, while the other correlation coefficients between BH1 and HYD, BH2 and BHZ, BH2 and HYD,and BHZ and HYD are 99.90,99.93,99.80,and 99.88%, respectively. These higher correlation coefficients for the first higher-order mode between different components inversely prove that the dispersive energies come from the same wave source, which is the P-SV type of Scholte wave. The stacked VS (SVS) is obtained from the normalized 4C VS stacking results. As shown in Figure 5B, the final full-mode SVS can be imaged using the proposed 4C VS stacking method, including the fundamental mode (M0 in Figure 5B) and the first to fourth higher-order modes. Compared with the traditional single vertical component VS, in which only the fundamental and first higher-order modes can be recognized, the fundamental to fourth higher-order modes are observed in the 4C SVS. The missed fundamental mode in the BH1 component appears in the SVS, and the missed second to fourth higher-order modes in the BHZ component are observed from the full-mode SVS.
We directly stack the 4C velocity spectra linearly using the same weights for each component. This approach has the significant advantages of being fast, stable, and effective. To verify whether the weights of the four components affect the extraction of higher-order modes for the Scholte-wave 4C VS stacking method, we chose two sets of random weights as a comparison. We denote the same weight stacking velocity spectrum as VS, which is calculated from Eq. 4 and is shown in Figure 6C. Two sets of different random weights with the 4C stacking velocity spectrum are obtained by: and shown in Figures 6A,B, respectively. All three SVSs show no visible differences and the dispersive characteristics are the same. That is, a total of five modes can be distinguished and the first two lower modes have much stronger energies with the last three higher modes being weaker. The absolute difference between every two SVSs is calculated. Thus, the differences between the VS and VS 0 , VS and VS 1 , and VS 0 and VS 1 are shown in Figures 6D-F, respectively. The average absolute difference is 0.028 for all three SVSs, which indicates that the deviation between the VSs as calculated with the same and different weights is approximately equal to or less than 3%. The average absolute difference is 0.035 and 0.036 between the VS and VS 0 and the VS and VS 1 , respectively. The absolute difference between the VSs calculated from two different weights is 0.116, which is the minimum for all three differences. This further verifies that different weights have little impact on the 4C dispersion spectrum stacking. The 4C dispersive VS stacking method for Scholte waves can effectively image the full-mode dispersion energy. Compared with the single-vertical-component DEI method, this procedure either increases the number of recognizable higherorder modes or widens the selectable frequency range of the dispersion curves for higher-order modes. For the singlehorizontal-component DEI processing, our 4C stacking method entirely avoids mistaken mode misidentification. Thus, the first higher-order mode is treated as the missed fundamental mode in the horizontal component. This also inversely demonstrates the necessity of the 4C Scholte-wave analysis and the failings of the single-component Scholte-wave imaging method.

Extraction and Inversion of Dispersion Curves
As shown in Figure 7, all dispersion curves were extracted manually from the 4C SVS of all 14 OBS stations. The manual selection error can be controlled within 5 m/s, which was considered in the subsequent inversion step. The numbers of extracted Scholte-wave modes differ between the OBS stations. Most stations can extract the third higher-order mode, whereas only three stations can extract the fourth higher-order mode. The highest number of extracted modes varies irregularly between the OBSs. A reasonable explanation for this phenomenon of unequal mode numbers is that it is caused by the lateral geological variations in the shallow sediment along the OBS survey line, or by differences in the degree of OBS coupling. Overall, all OBS stations can extract the available multimode dispersion curves for the subsequent inversion procedure.
The implicit relationship between the earth structure and the dispersion curves is described as (Xia et al., 2003).
Where n is the total number of points on the dispersion curve, (f i , c i ) is the frequency-phase velocity pair of the previous extracted dispersion curves, and v s , v p , ρ, and h are the geophysical model parameters (vector of shear-wave velocity, compressional wave velocity, density, and layer thickness, respectively). We use the multimode-dispersion-curve joint inversion method (Xia et al., 2003;Luo et al., 2007) to image the geological structures of shallow sediment. The final convergent optimal solution was obtained using the damped iterative leastsquares algorithm and was implemented in the Computer  Programs in Seismology (CPS) (Herrmann, 2013). During the iterations, only v s is considered as a free parameter, v p and h are fixed, and ρ is related to v s via the empirical relation ρ (g/cm 3 ) 0.8 log(v s ) + 0.23, with units of m/s (Herrmann, 2013). This inversion strategy has been proven reliable for the inversion of Scholte-wave dispersion curves (Wang et al., 2016;Wang et al., 2020). Wang et al. (2016)  To determine the deepest constraint depth during inversion, we calculate the sensitivity kernels for the five modes of the Scholte-wave dispersion curves based on the initial geophysical model. As shown in Figure 8, v s is the most sensitive to the dispersion curves, v p is the least sensitive, and ρ is in between. The statistical results show that the absolute mean sensitivities of v s , ρ, and v p are 3.62, 0.01, and 0.91, respectively. The mean sensitivity of v p is only 0.38% that of v s , and the corresponding ratio for ρ is higher at 25.05%. The sensitivity degrees for the different modes also show differences as the fundamental mode is more sensitive than the higher-order modes. The fundamental dispersive mode sensitivity of v s is 4.16, whereas the corresponding sensitivity of the higher-order modes is approximately 3.5. However, the higher-order modes can constrain the deeper layers more than the fundamental mode. Therefore, combining the higher-order modes in the inversion allows "seeing" deeper than only using the fundamental mode.
The sensitivity kernel results illustrate this constraint in Figure 8. The higher-order modes have greater sensitivities to deeper layers than the fundamental mode for the same frequency.  For example, when selecting 10% of the mean sensitivity (0.362) as the threshold, the greatest depths for the fundamental to the fourth higher-order mode at 6.0 Hz are 107,157,195,230,272, and 315 m, respectively. These simple statistical results demonstrate that for the same frequency, the constraint depth of the fourth higher-order mode is 2.94 times that of the fundamental mode. Correspondingly, the higher-order Scholtewave modes allow "seeing" deeper with an approximately linear trend compared with the fundamental mode. Finally, we choose 10% of the mean sensitivity (0.362) as the threshold to determine the deepest confidence depth for the inverted models, which is 347 m. We first obtain all maximum depth values from the sensitivity kernel of each mode for each frequency based on the threshold (Figure 9), and the final deepest inversion depth is the mean of these maximum depths.
We implement the five-mode joint dispersion-curve inversion to obtain the best-fit marine sedimentary model. All modes of the phase-velocity dispersion curves converge well after 20 iterations with a damping factor of 0.1 (Figure 10). The average signal power fit between the inverted and selected dispersion curves for all stations reaches 99.98%. The average standard fitting error for the multimode phase velocities is 4.13 m/s, where the corresponding average fitting residual is 2.46 m/s. The final best-fit v s model is very closer to the initial model of this station, while the corresponding density model differs considerably from the initial model. To illustrate how much the Scholte-wave phase velocity changes due to the density alone, another group of predicted dispersion curves (green solid line in Figure 10C) were calculated based on the inverted density and the initial v s profile. The result shows that the dispersion curves are almost the same between the initial model (blue dashed line in Figure 10C) and the inverted density and the initial v s model. Furthermore, this demonstrates that the density has an ignorable influence on Scholte-wave dispersion curves inversion. In particular, the shallowest layers of the inverted model have much lower shearwave velocities of 109.9 m/s to 137.9 m/s for the top two sedimentary layers.

RESULTS
In total, 26 dispersion-curve groups were successfully extracted from 13 OBS stations, with one station (C18) having data errors. Of the 26 dispersion-curve groups, three groups extracted a FIGURE 8 | Sensitivity kernels of the five dispersive Scholte-wave modes based on the initial marine sedimentary model. The red, blue, and black lines denote the sensitivities of v s , ρ, and v p , respectively. Modes 0-4 denote the Scholte-wave dispersive fundamental to the fourth higher-order modes, respectively.
FIGURE 9 | The maximum constrained depth obtained from the given threshold of the sensitivity kernel for each mode, where m0-m4 denote the fundamental to fourth-higher-order mode.
Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 10 | Article 812744 maximum of five modes, five groups extracted a maximum of four modes, 13 groups extracted a maximum of three modes, four groups extracted a maximum of two modes, and one group extracted a maximum of one mode. The phase velocities vary from 150 to 500 m/s for all dispersion-curve groups (Figure 7). The lower velocities belong to classical soft sedimentary geological conditions. Along the survey profile, the Scholtewave phase velocities show lateral differences, which indicates strong lateral heterogeneity. All 26 of the 1D shear-wave velocity structures were inverted from the multimode dispersion curves. As shown in Figure 11, the pseudo-2D shear-wave velocity profile was constructed from all the 1D velocity structures. The adjustable tension continuous curvature spline algorithm (Smith and Wessel, 1990) with a tension factor of 0.25 was adopted to interpolate the 1D V s profiles. Relatively low shear-wave velocities of 123-257 m/s were obtained in the top 100-m-deep sedimentary layers, which increased with an average highest rate of 3.47 m/s. The bottom sedimentary layer with a depth of 280-350 m had higher S-wave velocities of 500-670 m/s, with a mid-level velocity gradient of 1.25 m/s. The middle sedimentary layer showed mid-scale shear-wave velocities that varied from 257 to 500 m/ s with the lowest gradient of 0.77 m/s. A large lateral-velocity variation appeared in the constructed 2D S-wave velocity profile. There is a strong high-velocity zone in the middle and deep areas of the northern end (left side of Figure 11) of the survey line, which corresponds to the geological unit of the Qianliyan uplift. The apparent velocity difference between stations C08 and C10 may be indicated as an active fault, which is part of the Qianliyan fault zone (Wu et al., 2020). This corresponds to regional deep large faults with scales that cut through multiple sedimentary layers over a long development time. The V s profile between C10 and C22 shows the classical Quaternary sedimentary layers. At the right-hand section of the V s profile, there is a strong "piggyback pattern" thrust nappe structure induced from the reverse active fault, which indicates this type of extrusion effect is continuous, as shown from previous reflection seismic profiles (Wu et al., 2020). Some faults are converted from a normal type in the lower part to a reverse fault in the upper part, which may indicate that the stress environment changes from tension to compression and has a typical normal inverted structure (Wu et al., 2020).

CONCLUSION
The proposed OBS multicomponent Scholte-wave dispersion analysis method can effectively image the complex, extremely low, and strong lateral velocity variation V s structure in the shallow sediment of marginal seas. A complete, fast, and effective 2D shear-wave imaging approach using 4C OBS active-source Scholte-wave data was developed and demonstrated. The lower-frequency Scholte waves can be recorded from the three-component geophones and onecomponent hydrophone of the OBS. The four VSs imaged from the 4C Scholte-wave CRGs can be directly stacked linearly based on the conditions that all components record the same source of P-SV-type Scholte waves and no Love waves exist. Compared with the traditional single-verticalcomponent Scholte-wave DEI method, the stacking procedure can effectively image higher weak dispersive modes. Thus, the highest recognizable dispersive mode was the fourth from the SVS, whereas it was only the second higher-order mode from the single vertical component. Furthermore, the stacking procedure makes the dominant fundamental mode of the vertical component compensate perfectly for the missing fundamental mode in the horizontal components. Similarly, the dominant higher-order mode in the two horizontal components equivalently compensates for the weak energy modes in the vertical and hydrophone components. This kind of mutual benefit and disadvantage is complementarily for horizontal, vertical, and hydrophone components and illustrates the 4C Scholte-wave imaging mechanism and key points. In total, 26 dispersion-curve groups were extracted from the negative and positive branches of 13 OBS CRGs in the North Yellow Sea. Lower phase velocities of 150-500 m/s were obtained from all dispersion-curve groups. The damped iterative leastsquares algorithm effectively inverted the multimode Scholtewave dispersion curves (>99.98% goodness-of-fit) based on the water-sedimentary layered homogenous earth model. A maximum inversion-constraint model depth of 350 m was accurately defined by the sensitivity kernel threshold (10% of maximum). The pseudo-2D shear-wave velocity profile was constructed from the 26 1D inverted models. The key geological features of the shallow sediment appeared in the shear-wave velocity profile, including the uplift, active fault, and thrust nappe structures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YW, TH, and QY all contributed as major authors of the manuscript. YW wrote the text, processed the data, and created the figures. TH collected the marine OBS active data and performed the geological interpretation. QY developed the OBS instruments used to collect the seismic data in the study.