3.2.1. Beamforming

Wu et al. [22] studied one straight cable of the PRMRA and showed the good effect of the adapted eigenvalue-based filter. Since the PRMRA covers an area of 50 km2, local environments for different sensors and the coherence of the noise sources around different sensors may be different. It is necessary to study the application of the adapted eigenvalue-based filter to the entire array further. The study is also valuable for the shear-wave tomography in the following section.

To investigate the directionality of the wavefield, conventional beamforming is applied to the data recorded by different gathers (choose 30 continuously adjacent sensors in one cable as a gather, as is marked in red in Figure 1, for example). If the wavefield is diffuse and equipartitioned, the beamforming energy should be evenly distributed in all directions. An empirical estimation of the average phase speed 900 m/s is used in the beamforming. Although the estimation is rough, it does not make a big difference and is good enough to show the distribution of incoming energy. In addition, an iterative study can be done after the dispersion curves of the surface wave are extracted if it is necessary.

Figure 4 shows the beamforming result of one gather which is marked in red in Figure 1. Figure 4a,c show the beamforming energy distribution of the original wavefield. The energy is higher in the directions of −20◦ to 20◦ in Figure 4a than the background level. The energy becomes extraordinarily high from 30 to 40 min at 0◦ and can be recognized as Event B considering the high consistency in time and frequency. It indicates that Event B comes from the normal direction of the array. In Figure 4c, a source dominates the beam power with an almost constant direction of about 40◦. It can be recognized as Event A considering the high consistency in time. The travel time study of Event A shows the incident angle of the airgun shots is approximately 45◦. The misfit is caused by the usage of averaged phase speed in the beamforming. Figure 4a,c indicate that the original wavefield is affected by several strong, directional sources at both 2 Hz and 4 Hz. Event C, however, is not visible in the beamforming results because they have relatively lower energy than other events.

**Figure 4.** Beamforming results. (**<sup>a</sup>**,**b**) show the beamforming energy distribution of 2 Hz; (**<sup>c</sup>**,**d**) show the beamforming energy distribution of 4 Hz; (**<sup>a</sup>**,**<sup>c</sup>**) show the beamforming energy distribution of the original wavefield; (**b**,**d**) show the beamforming energy distribution of the filtered wavefield using the adapted eigenvalue-based filter.

Figure 4b,d illustrate the beamforming energy distribution of the filtered wavefield using the adapted eigenvalue-based filter. Figure 4b,d show relatively uniform distributed energy compared with Figure 4a,c. It indicates that the adapted eigenvalue-based filter can suppress most of the high energy and partially recover the expected diffuse and equipartitioned wavefield.

Keeping all the parameters (especially the weight) of the adapted eigenvalue-based filter the same, similar results have been obtained for other gathers. It indicates that the influence of non-microseismic events on the receiver array over 50 km<sup>2</sup> is more or less the same and that the adapted eigenvalue-based filter and the parameters used keep effective at the scale of 50 km2.

#### 3.2.2. The NCC Retrieval

Selecting 30 continuously adjacent sensors in one cable as a gather and 25 sensors as the overlap between the adjacent gathers, 573 gathers are obtained from the entire PRMRA. We have retrieved the NCC for all gathers using Equation (4).

Figure 5 displays the NCC results for two gathers marked in red (Figure 5a,b) and green (Figure 5c,d) in Figure 1. All the NCCs are bandpass filtered between 0.2 and 4.5 Hz. Figure 5a,c displays the NCC retrieved from unfiltered SCM. The waveforms are apparently contaminated by directional, energetic events. The NCCs in Figure 5a uses the same gather with beamforming in Figure 4. The narrow peaks (close to 0 s) traveling with a velocity of about 1450 m/s are nondispersive and can be recognized as Event A. However, the effect of Events B and C is not very tyical in the NCC, making it hard to recognize. The existence of these directional, energetic events shields the expected surface waves severely. Figure 5b,d displays the NCC retrieved from the filtered SCM. Most of the influence of the directional, energetic events has been removed from the NCC and the expected surface waves are significantly enhanced. Figure 5b,d also show a clear dispersion feature, which is a typical characteristic for surface waves.

**Figure 5.** The NCC results. Note that the NCCs are bandpass filtered between 0.2 and 4.5 Hz. (**<sup>a</sup>**,**b**) displays the NCC of the gather marked in red in Figure 1; (**<sup>c</sup>**,**d**) displays the NCC of the gather marked in green in Figure 1; (**<sup>a</sup>**,**<sup>c</sup>**) displays the NCC retrieved from unfiltered SCM while (**b**,**d**) displays the NCC retrieved from the filtered SCM.

The values of basic parameters used in the AEF are displayed in Figures 6 and 7. The value of *N* shown in Figure 6 is a function of frequency when gathers with the same shape are used and an averaged slowness of the medium 1.1 s/km is chosen. The weight *w* = 0.2 is used for this study. Note that *K* may be different for different gathers. Figure 7 displays the values of *K* of the gather marked in red in Figure 1. It shows increased values for larger frequencies because energetic events dominate higher frequencies over almost the entire time, as is indicated in Figure 4c. Between 30–40 min at around 1.5–3 Hz, Figure 7 shows larger values of *K*. It corresponds well with Figure 4a, where Event B happens at 30–40 min. For frequencies lower than 1.5 Hz, the values of *K* are very small, indicating a better equipartitioned wavefield for lower frequencies.

**Figure 6.** The value of *N* as a function of frequency.

**Figure 7.** The value of *K* as a function of frequency and time. The chosen gather is marked in red in Figure 1.

## **4. Shear-Wave Tomography**

The dispersive behavior of the retrieved surface waves carry useful information about the subsurface. The dispersion property of the retrieved surface waves is highly sensitive to shear-wave velocities at different depths and is used to perform 3D shear-wave tomography of the area.

The slowness-frequency transform (SFT) [31] is applied to the retrieved NCC of each gather. The result of the SFT can be expressed as a Gabor diagram, which shows signal energy as a function of frequency and phase velocity. Then, the phase–velocity dispersion curves can be extracted from the resulting Gabor diagram of the SFT. Examples of the retrieved dispersion curves can be found in [22].

#### *4.1. The ASSA Inversion*

A nonlinear optimization method, adaptive simplex simulated annealing algorithm (ASSA) [32], is applied to the fundamental mode of extracted dispersion curves of all gathers for estimating shear-wave velocity as a function of depth and range in the sea bottom.

The ASSA combines the simulated annealing and downhill simplex method. It operates on a simplex of *M* models and randomly perturbs the parameters after a downhill simplex step. The obtained trial model is subsequently accepted or rejected according to the Metropolis criteria

$$P(\Delta E) = \exp(-\Delta E/T),\tag{5}$$

where *T* denotes the temperature and *E* denotes the mismatch. The temperature *T* is decreased every *Np* accepted perturbations using

$$T\_{\bar{j}} = \beta^{\bar{j}} T\_{0\prime} \tag{6}$$

where *T*0 denotes the initial temperature and *β* denotes a temperature reduction factor satisfying 0 < *β* < 1. The procedure continues until the models in the simplex satisfy the convergence condition

$$\frac{E\_{\text{high}} - E\_{\text{low}}}{\left(E\_{\text{high}} + E\_{\text{low}}\right)/2} < \varepsilon\_{\prime} \tag{7}$$

where *<sup>E</sup>*high and *E*low represent the highest and lowest mismatch model in the simplex, respectively. Assuming that the measured dispersion data errors (including measurement error and theory error) are independent, Gaussian-distributed random processes, the likelihood function can be written as

$$L(\mathbf{m}) = \frac{1}{(2\pi)^{N/2} \prod\_{i=1}^{N} \sigma\_i} \exp\left\{-\sum\_{i=1}^{N} \frac{\left[d\_i - d\_i(\mathbf{m})\right]^2}{2\sigma\_i^2}\right\},\tag{8}$$

where **d** denotes the dispersion data vector with *N* entries, **d**(**m**) denotes the modeled data vector and *σ* is the standard deviation of **d**. The data misfit

$$E(\mathbf{m}) = \sum\_{i=1}^{N} \frac{\left[d\_i - d\_i(\mathbf{m})\right]^2}{2\sigma\_i^2} \tag{9}$$

is used as the objective function in the inversion. Therefore, the ASSA minimizes the misfit and yields an maximum-likelihood solution.

The chosen values of the mentioned parameters in this paper are listed in Table 1. Note that *M* is chosen based on the dimensional of the search space, *E*0 is the initial misfit of the initialized model and *σi* is assumed to be the same for all entries of *σ*.

**Table 1.** Values of the parameters used in the ASSA inversion.


The forward problem is solved using the DISPER80 subroutines [33], which can calculate phase–velocity dispersion curves given a geoacoustic model.

#### *4.2. Seabed Parameterizations*

A 6-layer geoacoustic model (shown in Figure 8) is used in the inversion. The 6th layer is a semi-infinite space. Considering that the dispersion property of the surface wave is not sensitive to the density *ρ* and the compressional velocity *vP* in the frequency range used (0.2 to 4.5 Hz), we use the following empirical equation [34,35]

$$wp = 1.16 \times v\_{\\$} + 1.36,\tag{10}$$

$$
\rho = 1.74 v\_{\text{P}}^{0.25},
\tag{11}
$$

to evaluate the values of them during the inversion.


**Figure 8.** Geoacoustic model used in the inversion.

The shear-wave velocity *vS* and the thickness of the layer *h* are variable in the iteration of the inversion process. Prior bounds are applied to the geoacoustic parameters for excluding physically unreasonable values. Table 2 shows the chosen bounds for different parameters.


**Table 2.** Search bound of seabed parameters.

## *4.3. Tomography Result*

The ASSA provides the maximum-likelihood solution of the data. We have performed 100 independent ASSA inversions on each gather to form a mean shear-wave velocity profile. A 3D shear-wave structure (Figures 9a and 10) can be reconstructed based on the solution.

**Figure 9.** (**a**) a vertical slice of the resolved mean shear-wave subsurface structure under the 10th cable counted from left to right of the PRMRA in Figure 1. The top white layer represents the water layer. The dashed lines indicate the horizontal slices used in Figure 10; (**b**) a vertical slice of the shear-wave subsurface structure from active-source reflection tomography [36]; (**c**) a vertical slice of the shear-wave subsurface structure from 1D Monte Carlo ambient noise tomography [37]. The circles in (**<sup>a</sup>**,**b**) mark the high-velocity zones.

Figure 9a displays a vertical slice of the shear-wave subsurface structure under the 10th cable counted from left to right of the PRMRA in Figure 1. Four subsurface layers (excluding the water layer) are well resolved from a holistic perspective, although a 6-layer model is used in the inversion. The top white layer (0–125 m) denotes the water layer. The first layer (125–140 m) is a low-velocity layer (lower than 300 m/s) which is very thin. This layer locates near the ocean bottom and is thought to be very soft. The second layer (140–240 m) is relatively thicker and harder, and contains a velocity of about 500 to 600 m/s. As the depth goes deeper, the layers become much thicker and contain larger velocities. The third layer (240–500 m) contains a velocity of about 800 to 900 m/s. The interfaces for the top three layers below the water layer are clear and the velocity values distribute relatively homogeneous over range. The bottom layer (>500 m) is a high-velocity layer. The velocity values of it mostly span from 900 to 1300 m/s (can be even larger in several areas), indicating a more complicated structure

compared with the top three layers. Note that the resolution in range of Figure 9 is 200 m, which is dependent on the chosen overlap between the adjacent gathers.

Figure 9b displays a similar slice of the shear-wave subsurface structure from active-source reflection tomography at the same field from [36]. The active-source has higher energy and the resolution of the reflection tomography is better. Thus, we consider it as a good approximation to the true model of the seabed. The top thin layer (125–140 m) of Figure 9a can not be recognized in Figure 9b. It might be caused by the regularization method which smooths the layers. However, it can be recognized easily in the P-wave velocity structure in [36]. The existence of the top thin layer in Figure 9a is reasonable because the top layer of shallow ocean is usually very soft with low shear-wave velocity. Figure 9b combines the first and second layer ( 125–240) of Figure 9a and shows a velocity of around 500 m/s at this layer. When the depth goes deeper in Figure 9b (240–500 m, corresponding to the third layer of Figure 9a), the velocity increases from 500 to 900 m/s. The third layer of Figure 9a can be considered as a transition layer at this place, compared with Figure 9b. In the depth range of 500 to 1700 m, although Figure 9a generally shows higher velocities than Figure 9b, they do have some similarities in the structure. At the depth of 600 m, Figure 9b shows two high-velocity zones A and B. We can find a similar zone A in Figure 9a, while zone B is not visible here possibly because of limited resolution of inversion using ambient noise. If we look further at the horizontal slice of the subsurface structure at the depth of 600 m in Figure 10c, both the high-velocity zones A and B can be recognized. At the depth of 1000 to 1600 m, Figure 9b shows another large high-velocity zone C. We can find a similar zone C in Figure 9a, although it is not as large as the zone in Figure 9b. The horizontal slice of Figure 9a at the depth of 1100 m demonstrates the distribution of this zone in Figure 10d, which is larger than it shows in Figure 9a. The high similarity in these zones of Figure 9a,b indicates that these high velocity anomalies are reliable and can be recognized as real features of the seabed. Below the depth of 1700 m, Figure 9b indicates a layer with very high velocities, which can not be recognized in Figure 9a. One possible reason is that surface wave can not penetrate to this depth, so that the inversion of Figure 9a is not sensitive to the shear-wave velocity at this depth.

Figure 9c displays a similar slice of the shear-wave subsurface structure from 1D Monte Carlo ambient noise tomography using 12 h of continuous ambient noise records at the same field from [37]. Although Zhang et al. [37] also did 2D and 3D tomography with the same data and obtained smoother results, Figure 9c is more comparable with the result of the present paper other than 2D and 3D results because they both use 1D parametrization. Figure 9c illustrates continuously increasing velocities from the depth of 125 to 500 m. The general structure corresponds well with Figure 9a,b. When the depth goes deeper than 500 m, Figure 9c shows complicated structure. These complicated features are thought to be artifacts due to the inversion strategy. Although lateral spatial constraints in 2D and 3D tomography make the structure smoother [37], the zones recognized in Figure 9a,b are not visible. One possible reason is that the authors used a narrower frequency band (0.35–1.5 Hz) compared with what we have used, and the information of the dispersion curves is limited. If more modes are included in the inversion, the frequency band needs to be wider. On the other side, people have to deal with the problem of the interferences included when a wider frequency band is used. At a depth deeper than 500 m, Figure 9c seems to ge<sup>t</sup> lower mean velocities than Figure 9a. The difference can be explained by several factors. One factor is the measurement errors of the dispersion data. The depth sensitivity (described in the following section) indicates that only lower frequency components can penetrate to a deeper depth so that the measurement errors at low frequencies may result in large bias to the velocities of deeper layers. The other factor is the large uncertainty at deeper layers caused by the inversion method, which can be quantified by statistical methods. The uncertainty for depth deeper than 500 m of Figure 9c is about 250 m/s [37] and the uncertainty for depth deeper than 500 m of Figure 9a is about 100 m/s (described in the following section).

**Figure 10.** (**<sup>a</sup>**–**d**) show horizontal slices of the mean shear-wave subsurface structure at the depth of 150, 300, 600 and 1100 m, respectively. The black line in (**a**) shows the location of Figure 9. The letters in (**<sup>c</sup>**,**d**) correspond to the letters in Figure 9.

Figure 10 displays four horizontal slices of the shear-wave subsurface structure at different depths. The distribution is roughly 'flat' for Figure 10a,b except for some small areas with a little bit high shear-wave velocities. Figure 10c,d show the horizontal distribution of the shear-wave velocities corresponding to the straight dashed lines marked in Figure 9a. In Figure 10c, zones A and B are a little discrete. This is probably because the results slightly overfit the local data. Adding lateral constraints to the parametrization of the inversion may help to improve it. The high velocity patches in Figure 10c,d correspond to the high velocity patches in Figure 9b well. It indicates that there may be similar geological features at this field although the ground-truth is not known.
