Shallow Shear-Wave Velocity Structure beneath the West Lake Area in Hangzhou, China, from Ambient-Noise Tomography

Urban geophysical exploration plays an important role in the sustainable development of and the mitigation of geological hazards in metropolitan areas. However, it is not suitable to implement active seismic methods in densely populated urban areas. The rapidly developing ambient-noise tomography (ANT) method is a promising technique for imaging the near-surface seismic velocity structure. We selected the West Lake area of the city of Hangzhou as a case study to probe the shallow subsurface shear-wave velocity (Vs) structure using ANT. We conducted seismic interferometry on the ambient-noise data recorded by 28 seismograph stations during a time period of 17 days. Fundamental-mode Rayleigh-wave group- and phase-velocity dispersion data were measured from cross-correlation functions and then inverted for a 3D Vs model of the uppermost 1 km that covers an area of about 7 km × 8 km. The tomographic results reveal two prominent anomalies, with high velocities in the southwest and low velocities in the northeast. The fast anomaly corresponds to the presence of limestone and sandstone, whereas the slow anomaly is due to the relatively low-velocity rhyolite and volcanic tuff in the area. The boundary between the two anomalies lies to the NE of an NW–SE trending fault, indicating that the fault dips toward the NE. In addition, the pronounced low-velocity anomalies appear under the Baoshi mountain, likely due to the thick rhyolite and volcanic tuff beneath the extinct volcano. Our results correlate well with regional geological features and suggest that ANT could be a promising technique for facilitating the exploration of urban underground space.


Introduction
Hangzhou is the capital city of Zhejiang province, southeastern China, with a population of over 9.8 million. The area around the West Lake in Hangzhou, a tourist attraction well-known nationwide, has suffered from tens of geological hazards, such as landslides and rockfall, which are in connection with the mountainous environment and complex geostructures such as the West Lake synclinorium. In addition, widespread unconsolidated sediments in the city may increase earthquake risk by amplifying the ground motion from potential adjacent earthquakes. Hence, there is an urgent demand to conduct urban geophysical surveys in this area, which are necessary for the related exploitation of underground resources and the mitigation of geological hazards. However, active (controlled-source) seismic methods are not suitable due to environmental issues in the urban area. To this end, we employed an environmentally friendly imaging technique, that is, ambient-noise tomography (ANT), to image the shallow shear-wave velocity (Vs) structure beneath the West Lake area.
ANT is based on the principle of generating empirical Green's functions (EGFs) by cross-correlating ambient noise recorded at different receivers [1][2][3][4], which are in turn used for tomography. Previous theoretical work [1][2][3][4] showed that under the assumption of a homogenous source distribution, the EGF between stations A and B (G AB ) can be expressed as where C AB (t) represents the correlation functions. Campillo and Paul [5] first demonstrated that EGFs between two receivers can be retrieved by cross-correlation of the seismic coda waves. Then, Shapiro and Campillo [6] successfully extracted Rayleigh waves by correlating ambient seismic noise. Shapiro et al. [7] performed the first ambient-noise Rayleigh-wave tomography in South California. A large number of subsequent applications of ANT have been focused on continental-or regional-scale crustal structure [8][9][10][11][12][13].
The aforementioned studies pave the way for the application of ANT in urban geophysical exploration. We selected the West Lake area as a case study to image the shallow Vs structure and to verify the applicability of ambient-noise tomography in urban underground space exploration. Our study region (Figure 1) is to the northwest of the Jiangshan-Shaoxing fault zone, the suture between the Yangtze Craton and the Cathaysia Terrane in the South China Block. The Cathaysia Terrane was subducted under the Yangtze Craton during the Neoproterozoic [18,19]. Subsequently, the Mesozoic subduction of the paleo-Pacific plate beneath the Eurasian plate led to the formation of several NE-SW trending faults, folds, and volcanos, which form the current tectonic frameworks in Hangzhou [20]. The Baoshi mountain (BSM in Figure 2) is an extinct volcano and once erupted tuff that covered the northeastern WLS during the early Cretaceous. During the neotectonics, the northeastern part of the West Lake synclinorium (WLS) experienced significant subsidence while the southwestern WLS was uplifted [21,22]. From the geological map of the West Lake area (Figure 2), it can be seen that the study area is mainly composed of Paleozoic sandstone and limestone with Cretaceous rhyolite and volcanic tuff in the northeastern part.
Remote Sens. 2021, 13, x FOR PEER REVIEW 2 of 12 used for tomography. Previous theoretical work [1][2][3][4] showed that under the assumption of a homogenous source distribution, the EGF between stations A and B ( ) can be expressed as where represents the correlation functions. Campillo and Paul [5] first demonstrated that EGFs between two receivers can be retrieved by cross-correlation of the seismic coda waves. Then, Shapiro and Campillo [6] successfully extracted Rayleigh waves by correlating ambient seismic noise. Shapiro et al. [7] performed the first ambient-noise Rayleigh-wave tomography in South California. A large number of subsequent applications of ANT have been focused on continental-or regional-scale crustal structure [8][9][10][11][12][13].
The aforementioned studies pave the way for the application of ANT in urban geophysical exploration. We selected the West Lake area as a case study to image the shallow Vs structure and to verify the applicability of ambient-noise tomography in urban underground space exploration. Our study region ( Figure 1) is to the northwest of the Jiangshan-Shaoxing fault zone, the suture between the Yangtze Craton and the Cathaysia Terrane in the South China Block. The Cathaysia Terrane was subducted under the Yangtze Craton during the Neoproterozoic [18,19]. Subsequently, the Mesozoic subduction of the paleo-Pacific plate beneath the Eurasian plate led to the formation of several NE-SW trending faults, folds, and volcanos, which form the current tectonic frameworks in Hangzhou [20]. The Baoshi mountain (BSM in Figure 2) is an extinct volcano and once erupted tuff that covered the northeastern WLS during the early Cretaceous. During the neotectonics, the northeastern part of the West Lake synclinorium (WLS) experienced significant subsidence while the southwestern WLS was uplifted [21,22]. From the geological map of the West Lake area (Figure 2), it can be seen that the study area is mainly composed of Paleozoic sandstone and limestone with Cretaceous rhyolite and volcanic tuff in the northeastern part.

Data and Method
We deployed 28 broadband seismograph stations (equipped with Trillium Horizon 120 seismometer and Centaur digitizer) in the West Lake area to acquire ambient-noise data with a sampling rate of 200 Hz and an average interstation spacing of~0.5-1 km ( Figure 2). Data acquisition lasted for 17 days, from 29 July to 14 August 2018. The study region covers an area of 7 km × 8 km. Figure 3a,b show a one-day vertical-component waveform recorded by station WL.B04 ( Figure 2) and its spectrogram, respectively. The amplitude of the spectrogram is higher during the daytime than during the night, which is due to the influence of anthropogenic activities.  Only vertical-component data were processed to obtain Rayleigh-wave EGFs using procedures similar to those documented in Reference [23]. First, we truncated the data into hour-long segments and decimated the recordings at a sampling rate of 50 Hz. We removed the mean and trend from the data. Then, the data were bandpass-filtered between 0.2 and 5 s periods. Both temporal and spectral normalizations were applied before crosscorrelation to reduce the influence of earthquakes and non-stationary noise energy. Since running-absolute-mean normalization of the recordings can lead to faster convergence of noise cross-correlation functions (NCFs) [24], we implemented running-absolute-mean normalization instead of 1-bit normalization for the purpose of temporal normalization. Finally, cross-correlation was performed on hour-long recordings for each station pair with lag times ranging from −30 s to 30 s. These hourly results were stacked to obtain the final NCFs corresponding to the whole recording duration for each station pair. In total, we obtained 378 NCFs.
Conventional linear stacking and time-frequency domain phase-weighted stacking (tf-PWS) techniques were compared to increase the signal-to-noise ratio (SNR) of NCFs. Tf-PWS calculates the phase coherence coefficients and downweights the inconsistent signals in the linear stacked NCFs [25,26]. Figure 4a,b show the NCFs for one station pair recovered by two stacking schemes, respectively, illustrating the considerable increase in the SNR by the tf-PWS. More NCFs with higher SNRs can be obtained from the tf-PWS method than the linear stacking method (Figure 4c,d). Hence, NCFs constructed from the tf-PWS method were used for later processing. As seen in Figure 4c,d, the coherent Rayleigh waves elucidate an apparent velocity of approximately 2 km/s. removed the mean and trend from the data. Then, the data were bandpass-filtered between 0.2 and 5 s periods. Both temporal and spectral normalizations were applied before cross-correlation to reduce the influence of earthquakes and non-stationary noise energy.
Since running-absolute-mean normalization of the recordings can lead to faster convergence of noise cross-correlation functions (NCFs) [24], we implemented running-absolutemean normalization instead of 1-bit normalization for the purpose of temporal normalization. Finally, cross-correlation was performed on hour-long recordings for each station pair with lag times ranging from −30 s to 30 s. These hourly results were stacked to obtain the final NCFs corresponding to the whole recording duration for each station pair. In total, we obtained 378 NCFs. Conventional linear stacking and time-frequency domain phase-weighted stacking (tf-PWS) techniques were compared to increase the signal-to-noise ratio (SNR) of NCFs. Tf-PWS calculates the phase coherence coefficients and downweights the inconsistent signals in the linear stacked NCFs [25,26]. Figure 4a,b show the NCFs for one station pair recovered by two stacking schemes, respectively, illustrating the considerable increase in the SNR by the tf-PWS. More NCFs with higher SNRs can be obtained from the tf-PWS method than the linear stacking method (Figure 4c,d). Hence, NCFs constructed from the tf-PWS method were used for later processing. As seen in Figure 4c,d, the coherent Rayleigh waves elucidate an apparent velocity of approximately 2 km/s. Since the uneven distribution of noise energy would distort the NCFs and affect the reliability of dispersion measurements [27], we calculated the SNR for both acausal and causal components in each NCF as a proxy to identify the azimuth dependence of noise energy. SNR is defined by the ratio of the absolute maximum amplitude within the signal window and the root mean square value in the noise window (−30 to −15 s and 15 to 30 s). The analysis of noise directionality is similar to the technique adopted by Yang and Since the uneven distribution of noise energy would distort the NCFs and affect the reliability of dispersion measurements [27], we calculated the SNR for both acausal and causal components in each NCF as a proxy to identify the azimuth dependence of noise energy. SNR is defined by the ratio of the absolute maximum amplitude within the signal window and the root mean square value in the noise window (−30 to −15 s and 15 to 30 s). The analysis of noise directionality is similar to the technique adopted by Yang and Ritzwoller [28], in which the strength of directional noise energy from two opposite azimuths was estimated by the SNR values. The azimuthal distribution of SNRs for both acausal and causal components of all NCFs is shown in Figure 5. We normalized the SNRs Ritzwoller [28], in which the strength of directional noise energy from two opposite azimuths was estimated by the SNR values. The azimuthal distribution of SNRs for both acausal and causal components of all NCFs is shown in Figure 5. We normalized the SNRs by the maximum SNR value and found that approximately equal power of noise energy arrived from variant directions, indicating that the noise source is nearly homogeneous.

Rayleigh-Wave Dispersion Curve Measurement
We stacked the causal and acausal parts of the NCFs linearly to form the symmetric component of NCFs. Then, the dispersion curve was extracted from every symmetric NCF using the multiple filter analysis package in Computer Program in Seismology (CPS) [29,30]. We checked the quality of every symmetric NCF carefully and chose dispersion curves manually. Two selection criteria were applied to avoid biased dispersion measurements. First, interstation distance should be at least 1.5 times the wavelength to satisfy the surface-wave far-field approximation condition. Second, the SNR of NCFs should be greater than 10 to retain the most reliable measurements. The extracted group-and phasevelocity dispersion curves for the fundamental-mode Rayleigh wave in the period range from 0.35 s to 2 s are shown in Figure 6a

Rayleigh-Wave Dispersion Curve Measurement
We stacked the causal and acausal parts of the NCFs linearly to form the symmetric component of NCFs. Then, the dispersion curve was extracted from every symmetric NCF using the multiple filter analysis package in Computer Program in Seismology (CPS) [29,30]. We checked the quality of every symmetric NCF carefully and chose dispersion curves manually. Two selection criteria were applied to avoid biased dispersion measurements. First, interstation distance should be at least 1.5 times the wavelength to satisfy the surfacewave far-field approximation condition. Second, the SNR of NCFs should be greater than 10 to retain the most reliable measurements. The extracted group-and phase-velocity dispersion curves for the fundamental-mode Rayleigh wave in the period range from 0.35 s to 2 s are shown in Figure 6a Ritzwoller [28], in which the strength of directional noise energy from two opposite azimuths was estimated by the SNR values. The azimuthal distribution of SNRs for both acausal and causal components of all NCFs is shown in Figure 5. We normalized the SNRs by the maximum SNR value and found that approximately equal power of noise energy arrived from variant directions, indicating that the noise source is nearly homogeneous.

Rayleigh-Wave Dispersion Curve Measurement
We stacked the causal and acausal parts of the NCFs linearly to form the symmetric component of NCFs. Then, the dispersion curve was extracted from every symmetric NCF using the multiple filter analysis package in Computer Program in Seismology (CPS) [29,30]. We checked the quality of every symmetric NCF carefully and chose dispersion curves manually. Two selection criteria were applied to avoid biased dispersion measurements. First, interstation distance should be at least 1.5 times the wavelength to satisfy the surface-wave far-field approximation condition. Second, the SNR of NCFs should be greater than 10 to retain the most reliable measurements. The extracted group-and phasevelocity dispersion curves for the fundamental-mode Rayleigh wave in the period range from 0.35 s to 2 s are shown in Figure 6a

Construction of 3D Vs Model
A two-step inversion scheme was implemented to obtain the 3D Vs model. We first inverted the group-and phase-velocity dispersion measurements for 2D group-and phase-Remote Sens. 2021, 13, 2845 6 of 11 velocity maps, respectively, on a 0.003 • × 0.003 • grid at different periods using a non-linear tomographic technique [31]. For the forward problem, a fast-marching method [32,33] was employed to calculate the travel-time field by solving the grid-based Eikonal equation through a wave-front tracking method, and then the travel-time gradient was calculated to determine the ray paths between station pairs. The final model was solved by an iterative subspace inversion scheme [34]. The optimum damping and smoothing weights were obtained from the standard trade-off curves for travel-time residuals and model roughnesses. We also performed over-smoothed tomographic inversions and rejected dispersion measurements with residuals greater than 2.5 standard deviations from the mean residual. The resulting group-velocity maps were plotted in the period range of 0.46 s to 1.1 s and phase-velocity maps from 0.40 s to 0.95 s (Figure 7). We further evaluated the spatial resolution of our model through synthetic checkerboard tests. The synthetic initial model consists of high-and low-velocity checkers with the maximum velocity perturbation of 10% relative to the average velocity. The checkerboard size was three times the grid size, and randomly distributed noise was added to the synthetic data. Here, the randomly distributed noise level was set to be 1% of the theoretical travel time. The recovered phase-velocity checkerboards at three specific periods are shown in Figure 8, suggesting that the tomographic results are capable of resolving structures with a lateral resolution of about 1 km × 1 km. board tests. The synthetic initial model consists of high-and low-velocity checkers with the maximum velocity perturbation of 10% relative to the average velocity. The checkerboard size was three times the grid size, and randomly distributed noise was added to the synthetic data. Here, the randomly distributed noise level was set to be 1% of the theoretical travel time. The recovered phase-velocity checkerboards at three specific periods are shown in Figure 8, suggesting that the tomographic results are capable of resolving structures with a lateral resolution of about 1 km × 1 km. Rayleigh-wave group-and phase-velocity dispersion curves at every grid point were combined to jointly invert for the isotropic 1D Vs profile using a linearized procedure [30]. The average velocity for the uppermost 40 m from the borehole data (yellow rectangles in Figure 8a) and the deep crustal Vs model from Bao et al. [10] were integrated to generate the initial model with a layer thickness of 100 m. Surface waves at different periods are sensitive to structures at different depth ranges. The peak of the sensitivity kernel for the Rayleigh wave is around one-third of its wavelength. For a given period, Rayleigh-wave group velocity is sensitive to shallower structures than phase velocity. The inverted 1D Vs model at one grid point is used to calculate the sensitivity kernel using CPS package [30]. From the sensitivity kernels shown in Figure 9, our dataset can effectively constrain the Vs structure down to a depth of ~1 km. Rayleigh-wave group-and phase-velocity dispersion curves at every grid point were combined to jointly invert for the isotropic 1D Vs profile using a linearized procedure [30]. The average velocity for the uppermost 40 m from the borehole data (yellow rectangles in Figure 8a) and the deep crustal Vs model from Bao et al. [10] were integrated to generate the initial model with a layer thickness of 100 m. Surface waves at different periods are sensitive to structures at different depth ranges. The peak of the sensitivity kernel for the Rayleigh wave is around one-third of its wavelength. For a given period, Rayleigh-wave group velocity is sensitive to shallower structures than phase velocity. The inverted 1D Vs model at one grid point is used to calculate the sensitivity kernel using CPS package [30]. From the sensitivity kernels shown in Figure 9, our dataset can effectively constrain the Vs structure down to a depth of~1 km.
We assembled all the 1D Vs models at each grid point to construct the final 3D Vs model. Several representative horizontal slices were plotted in Figure 10, showing an obvious bimodal distribution of Vs in the area as observed in Rayleigh-wave group-and phase-velocity maps (Figure 7): high-velocity anomalies in the southwest and low-velocity anomalies in the northeast. The fast anomalies mainly occur in the southwestern WLS, whereas the slow anomalies in the northeast generally correspond to the BSM and the West Lake region. The boundary separating the two anomalies lies to the NE of a NW-SE trending fault, suggesting that the fault dips towards NE. Furthermore, four velocity profiles in Figure 11 traversing the main geological structures ( Figure 2) were plotted to display the lateral velocity variations. It is noticeable again that the boundary between the low-velocity structure in the northeast and the high-velocity structure in the southwest dips toward NE (AA' and BB' profiles). As shown in profiles AA' and DD', relatively thick low-velocity layers appear beneath the BSM. Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 12 (a) (b) Figure 9. The depth sensitivity kernels for group (a) and phase (b) velocities at different periods.
We assembled all the 1D Vs models at each grid point to construct the final 3D Vs model. Several representative horizontal slices were plotted in Figure 10, showing an obvious bimodal distribution of Vs in the area as observed in Rayleigh-wave group-and phase-velocity maps (Figure 7): high-velocity anomalies in the southwest and low-velocity anomalies in the northeast. The fast anomalies mainly occur in the southwestern WLS, whereas the slow anomalies in the northeast generally correspond to the BSM and the West Lake region. The boundary separating the two anomalies lies to the NE of a NW-SE trending fault, suggesting that the fault dips towards NE. Furthermore, four velocity profiles in Figure 11 traversing the main geological structures ( Figure 2) were plotted to display the lateral velocity variations. It is noticeable again that the boundary between the low-velocity structure in the northeast and the high-velocity structure in the southwest dips toward NE (AA' and BB' profiles). As shown in profiles AA' and DD', relatively thick low-velocity layers appear beneath the BSM.

Discussion
In this study, using the 17-day continuous ambient-noise data, we extracted high-SNR Rayleigh waves within the period range of 0.35-2 s, which were inverted to constrain the subsurface Vs structure to~1 km depth. Higher-and lower-frequency surface waves with high SNRs cannot be recovered, likely due to the limitation of inter-station spacing and high attenuation of the shallow medium. Generally speaking, we can obtain higher.
Frequency surface waves with shorter station-pair distance and vice versa. For example, using data from the extremely dense seismic array in Long Beach, California, with an average inter-station spacing of~100 m, Lin et al. [14] extracted Rayleigh waves with frequencies of up to 4 Hz, which can improve the resolution of the near-surface structure.
From the sensitivity kernels of the extracted Rayleigh waves (Figure 9), we can see that our data are highly sensitive to the subsurface structure below 100 m depth. Thus, we suggest that the shallow structures located above 50 m depth are not the likely cause of the imaged velocity variations. Hence, the predominant velocity anomalies are mainly associated with the variations in the underlying bedrock. From the geological map (Figure 2), the western WLS comprises two synclines and is mainly composed of Paleozoic sandstone and limestone, while the subsided northeastern WLS is covered by the Cretaceous rhyolite and tuff. The limestone and sandstone have higher shear velocity (about 2.0-3.3 km/s) than the rhyolite and tuff (about 1.5-2.4 km/s). Thus, our tomographic models correlate well with the known geological features. The high-velocity anomaly in the southwest corresponds to the sandstone and limestone, while the slow anomaly in the northeast is related to the low-velocity rhyolite and volcanic tuff. It can be seen that the low-velocity area shrinks at greater depths and appears only beneath the BSM at depths of 850-950 m, likely due to the thick rhyolite and tuff beneath the extinct volcano.
Based on the correlation between our tomographic model and the geologic features, we further demonstrate here the applicability of ambient-noise tomography in delineating the shallow Vs structure in urban areas for facilitating the exploration of urban underground structures. In addition, we also reveal the geometry of a regional fault. The structural variations across other faults have also been illuminated by previous studies, such as the Newport-Inglewood fault in Southern California [14] and the Shushan fault in the city of Hefei [15], emphasizing the promising future of ANT in studying the fault zone structure and for seismic hazard and risk assessment, since the ground motion of an earthquake is significantly influenced by the shallow Vs structure. With the potential of extracting body waves [35], more detailed structures of the shallow crust could also be resolved by ambient-noise interferometry.

Conclusions
In this study, we applied the environmentally friendly ANT method to investigate the shallow Vs structure surrounding the West Lake. The SNRs of NCFs were increased by the tf-PWS method. We measured fundamental-mode Rayleigh-wave dispersion curves within 0.35-2 s period band. A two-step inversion scheme was implemented to construct a 3D Vs model with a horizontal resolution of about 1 km × 1 km, which correlates well with regional geological features. Our tomographic model illustrates two prominent anomalies with high velocities in the southwest and low velocities in the northeast. The former anomaly corresponds to the presence of limestone and sandstone, whereas the latter is due to the relatively low-velocity rhyolite and volcanic tuff in that region. In addition, the boundary between the two anomalies lies to the NE of an NW-SE trending fault, suggesting that the fault dips toward NE. The Baoshi mountain is underlain by a thick, pronounced low-velocity structure, likely due to the thick rhyolite and volcanic tuff beneath the extinct volcano. Our results further suggest that ambient-noise tomography is a promising technique in facilitating the exploration of urban underground space and seismic hazard and risk assessment.